线性插值.cpp

线性数值积分c++代码已通过上机验证适合于数值计算学习者编程欢迎下载

#include<stdio.h>
#include<math.h>
#define PI 3.141592653
double p(int i,int N)
{  double pi;
   pi=5-10*i/N;
   return pi;
}

double q(int i,int N)
{	double qi;
	qi=-5*cos((2*i+1)*PI/(2*N+2));
	return qi;
}

double qiuji(int N,int i,double y,double(*f)(int))
{	double x[N+1];
	int j,k;
	double temp;
	for(j=0;j<=N;j++)
		x[j]=(*f)(j);
	temp=1.0;
	for(k=0;k<i;k++)
		temp*=(y-x[k])/(x[i]-x[k]);
	for(k=i+1;k<=N;k++)
		temp*=(y-x[k])/(x[i]-x[j]);
	return temp;
}

double lagrange(int N ,double y,double(*f)(int))
{	int i;double m,n,s;
	s=0.0;
	for(i=0;i<=N;i++)
	{m=qiuji(N,i,y,f);
	n=1/(1+y*y);
	s+=m*n;
	}
    return s;
}

double max(double x[])
{	double temp;
	temp=x[0];
	int i;
	for(i=1;i<=101;i++)
	{
	if(x[i]>temp)
		temp=x[i+1];
	}
	return temp;
}

double shulie(int n,double(*f)(int))
{	double temp[n+1],grid[n+1];
	int i; double x=0.0,y=0.0;
	double *p;
	for(i=0;i<=100;i++)
	{y=i/10-5;
	x=lagrange(n,y,f);
	temp[i]=x;
	grid[i]=abs(temp-1/(1+y*y));
	}
	p=grid;
	temp=max(p);
	return temp;
}

void main()
{
	int x[4]={5,10,20,40};
	int j,n;
	for(j=0;j<=3;j++)
	{n=x[j];
	double a1=0.0,a2=0.0;
	a1=shulie(n,p);
	a2=shulie(n,q);
	printf("N:%d\n",n);
	printf("Max Error of grid1:%16.12f\n",a1);
	printf("Max Error of grid2:%16.12f\n",a2);
	}
}