lagrange.cpp

对函数构造Lagrange插值多项式

#include<iostream>
#include<math.h>
#include<iomanip>
using namespace std;

const int N = 40;
double PI = 3.1415926;
double x1[N+1],x2[N+1];
double f1[N+1],f2[N+1];
double c[101],f[101];

double lagrange(double t,int n,double x[],double y[])
{
	double fx = 0.0;
	for(int i=0;i<=n;i++)
	{
		int j;
		double tmp = 1.0;
		for(j=0;j<i;j++)
			tmp = tmp*(t-x[j])/(x[i]-x[j]);
		for(j=i+1;j<=n;j++)
			tmp = tmp*(t-x[j])/(x[i]-x[j]);
		fx = fx + tmp * y[i];
	}
	return fx;
}

void Newton_1(int n,double x[],double y[])
{
	int i;
	for(i=1;i<=n;i++)
	{
		for(int j=n;j>=i;j--)
			y[j] = (y[j]-y[j-1])/(x[j]-x[j-i]);
	}
}

double Newton_2(double t,int n,double x[],double y[])
{
	double fx = y[n];
	for(int i=n;i>=1;i--)
	{
		fx = y[i-1]  + (t-x[i-1])*fx;
	}
	return fx;
}

int create(int n)
{
	int i;
	for(i=0;i<=n;i++)
	{
		x1[i] = 5-(10.0/n)*i;
		x2[i] = -5*cos((2*i+1)*PI/(2*n+2));
		f1[i] = 1.0/(1+x1[i]*x1[i]);
		f2[i] = 1.0/(1+x2[i]*x2[i]);
	}
	for(i=0;i<=100;i++)
	{
		c[i] = i/10.0 - 5;
		f[i] = 1.0/(1+c[i]*c[i]);
	}
	return 1;
}

int main()
{
	int i;
	double p1,p2,max1,max2;
	max1 = max2 = 0;

	double p3,p4,max3,max4;
	max3 = max4 = 0;

	create(N);
	for(i=0;i<=100;i++)
	{
		p1 = lagrange(c[i],N,x1,f1);
		p2 = lagrange(c[i],N,x2,f2);
	
		if(f[i]-p1>max1)
			max1 = f[i]-p1;
		else if(f[i]-p1<-1*max1)
			max1 = p1 - f[i];
		if(f[i]-p2>max2)
			max2 = f[i]-p2;
		else if(f[i]-p2<-1*max2)
			max2 = p2-f[i];
	}
	
	Newton_1(N,x1,f1);
	Newton_1(N,x2,f2);
	
	for(i=0;i<=100;i++)
	{
		p3 = Newton_2(c[i],N,x1,f1);
		p4 = Newton_2(c[i],N,x2,f2);

		if(f[i]-p3>max3)
			max3 = f[i]-p3;
		else if(f[i]-p3<-1*max3)
			max3 = p3 - f[i];
		if(f[i]-p4>max4)
			max4 = f[i]-p4;
		else if(f[i]-p4<-1*max4)
			max4 = p4-f[i];
	}

	cout.precision(18);
	cout.setf(ios::fixed);
	cout<<"N = "<<N<<endl;
	cout<<"\nLagrange"<<endl;
	cout<<"Max Error of grid(1) : "<<max1<<endl;
	cout<<"Max Error of grid(2) : "<<max2<<endl;
	cout<<"\nNewton"<<endl;
	cout<<"Max Error of grid(1) : "<<max3<<endl;
	cout<<"Max Error of grid(2) : "<<max4<<endl;

	return 0;
}