lagrage.cpp

lagrange插值算法

#include<iostream.h>
#include<iomanip.h>
#include<fstream.h>
#include<iomanip.h>
#include<math.h>
#include<malloc.h>
#define PI 3.14159265358979
double fx( double x )
{
	double tmp=x*x;
	tmp+=1.0;
	return( 1.0/tmp );
}
double Lagrange( int n, double x[], double y[], double val )
{
	double lx=0.0;
	double tmp;
	int i;
	int j;
	double mid1,mid2,mid3;
	for( i=0; i<=n; i++ )
	{
		tmp=1.0;
		for( j=0; j<i; j++ )
		{
			mid1=val-x[j];
			mid2=x[i]-x[j];
			mid3=mid1/mid2;
			tmp=tmp*mid3;
		}
		for( j=i+1; j<=n; j++ )
		{
			mid1=val-x[j];
			mid2=x[i]-x[j];
			mid3=mid1/mid2;
			tmp=tmp*mid3;
		}
		tmp*=y[i];
		lx+=tmp;
	}
	return( lx );
}
double err( int n, double x[], double y[], double val )
{
	double tmp=fx( val );
	double mid=Lagrange( n, x, y, val );
	tmp-=mid;
	if( tmp>0 )
		return tmp;
	else
		return -tmp;
}
void main( void )
{
	double val;
	double max;
	double error;
	int i,j;
	double mid;
	ofstream file( "result.txt" );
	cout<<"第1组：插值节点为 xi=-5+10*i/N"<<endl;
	file<<"第1组：插值节点为 xi=-5+10*i/N"<<endl;
	for( j=5; j<=40; j*=2 )
	//j=20;
	{
		max=0;
		double *x=( double * )malloc( ( j+1 )*sizeof( double ) );
		double *y=( double * )malloc( ( j+1 )*sizeof( double ) );
		for( i=0; i<=j; i++ )
		{
			mid=10.0*i;
			mid/=j;
			x[i]=-5.0+mid;
			y[i]=fx( x[i] );
		}
		for( i=0; i<=100; i++ )
		{
			val=i/10.0-5.0;
            error=err( j, x, y, val );
			if( max<error )
				max=error;
		}
		cout<<"N="<<j<<"\t Err="<<setprecision( 13 )<<max<<endl;
		file<<"N="<<j<<"\t Err="<<setprecision( 13 )<<max<<endl;
		free( x );
		free( y );
	}
	cout<<"第2组：插值节点为 xi=-5cos((2i+1)/(2N+2)PI)"<<endl;
	file<<"第2组：插值节点为 xi=-5cos((2i+1)/(2N+2)PI)"<<endl;
	for( j=5; j<=40; j*=2 )
	//j=20;
	{
		max=0;
		double *x=( double * )malloc( ( j+1 )*sizeof( double ) );
		double *y=( double * )malloc( ( j+1 )*sizeof( double ) );
		for( i=0; i<=j; i++ )
		{
			mid=( 2.0*i+1.0 )/( 2.0*j+2.0 );
			mid*=PI;
			mid=cos( mid );
			x[i]=-5.0*mid;
			y[i]=fx( x[i] );
		}
		for( i=0; i<=100; i++ )
		{
			val=i/10.0-5.0;
            error=err( j, x, y, val );
			if( max<error )
				max=error;
		}
		cout<<"N="<<j<<"\t Err="<<setprecision( 13 )<<max<<endl;
		file<<"N="<<j<<"\t Err="<<setprecision( 13 )<<max<<endl;
		free( x );
		free( y );
	}
}