3insert.c

三次样条插值算法和拉格朗日插值算法的实现。

#include"math.h"
#include "stdio.h"
main()
{
	double sx(double,double*,double*);
	double lx(double);
	int i,j=0;
	int n;
	double p;
	double a[4][4];/*追赶法矩阵*/
	double x[6]={0.03846,0.05882,0.1,0.2,0.5,1.0};/*插值点数值*/
	double fx[5];/*求值*/
	double g[4];
	double u,v;
	double m[6]={0.01479,0,0,0,0,0};	
	double y[4];
	double b[4];
	//double x[4];
	for(i=0;i<4;i++)
		for(j=0;j<4;j++)
		{
		    	a[i][j]=0;
		}
	for(i=0;i<5;i++)
	{
		fx[i]=x[i+1]-x[i];
	}
	j=0;
	for(i=0;i<4;i++)/*构建追赶法中的矩阵a[][]*/
	{
		v=0.5;u=0.5;
		g[i]=3*(v*fx[i]+u*fx[i+1]);
		a[j][j]=2;
		if(j!=3)
			a[j][j+1]=u;
		if(j!=0)
			a[j][j-1]=v;
		j++;
	}
	g[0]=g[0]-0.5*m[0];
	g[3]=g[3]-0.5*m[5];
	/*生成了追赶法中的系数矩阵a[][]以及g[]*/
	/*以下求解追赶法方程的解*/
		
	/*求解追赶法中的矩阵得到m[]*/
	
	b[0]=a[0][1]/a[0][0];
	for(i=1;i<4;i++)
		b[i]=a[i][i+1]/(a[i][i]-a[i][i-1]*b[i-1]);
	y[0]=g[0]/a[0][0];
	for(i=1;i<4;i++)
		y[i]=(g[i]-a[i][i-1]*y[i-1])/(a[i][i]-a[i][i-1]*b[i-1]);
	m[4]=y[3];
	for(i=3;i>0;i--)
		m[i]=y[i-1]-b[i-1]*m[i+1];
	for(j=0;j<6;j++)
	{
		printf("m[%d]=%f,",j,m[j]);
		printf("\n");
	}
	printf("计算插值的点n=");
	scanf("%d",&n);
	//n=20;
	for(i=0;i<n;i++)
	{
		p=-5+(5*1.0/n)*i;
		printf("s[%f]=%f,  ",p,sx(p,x,m));
		printf("  f[%f]=%f,  ",p,1/(1+p*p));
		printf("  l[%f]=%f\n",p,lx(p));
	}
}
double sx(double x,double *y,double *m)/*三湾矩求值函数,对插值函数求值*/
{
/* a,b为小区间插值点,x为需要求的某点x坐标,sx为要求的某点的值*/

	int a,b;
	double sx;
	double y1,y2,m1,m2;
	if(x<0)
	{
		a=(int)x-1;
		b=a+1;
		y1=*(y+5+a);
		y2=*(y+5+b);
		m1=*(m+5+a);
		m2=*(m+5+b);
		
	}
	else
	{
		a=(int)x;
		b=a+1;
	}
	sx=(x-b)*(x-b)*(1+2*(x-a))*y1+(x-a)*(x-a)*(1+2*(b-x))*y2+(x-b)*(x-b)*(x-a)*m1+(x-a)*(x-a)*(x-b)*m2;
	return sx;
}
double lx(double x)/*拉格朗日插值求值函数*/
{
	int i,j;
	double l=1;
	double lx=0;
	for(i=0;i<6;i++)
	{
		for(j=0;j<6;j++)
		{
			if(j!=i)
			l=l*(x+j)/(-i+j);
		}
		lx+=l/(1+i*i);
		l=1;
	}		
	return lx;		
}