main.cpp

多项式插值的震荡现象与三次样条插值实验

#include<iostream>
#include<math.h>
#include<string>
using namespace std;
double f(double x);
double L(int n,double a,double b,double x);
void S(int n,double X[],double Y[],double y0,double yn,int x_num,double x[]);
int main()
{
	cout<<"说明：由于输出表达式形式效果不好也比较麻烦，所以对于L10(x)、"<<"\n"
		<<"L20(x)、S10(x)、S20(x)这里只输出公式中所需的所有常数项。"<<endl;
	cout<<"\n精确值f(4.8)="<<f(4.8)<<endl;
	cout<<"\nL10的情况："<<endl;
	double L10=L(10,-5,5,4.8);
	cout<<"L(4.8)="<<L10<<endl;
	cout<<"\nL20的情况："<<endl;
	double L20=L(20,-5,5,4.8);
	cout<<"L(4.8)="<<L20<<endl;
	double x10[11]={-5,-4,-3,-2,-1,0,1,2,3,4,5};
	double y10[11];
	double x_1[1]={4.8};
	for(int i=0;i<11;i++)
		y10[i]=f(x10[i]);
	S(11,x10,y10,2/125,-2/125,1,x_1);
	double x20[21]={-5,-4.5,-4,-3.5,-3,-2.5,-2,-1.5,-1,-0.5,0,0.5,1,1.5,2,2.5,3,3.5,4,4.5,5};
	double y20[21];
	for(int i=0;i<21;i++)
		y20[i]=f(x20[i]);
	S(21,x20,y20,2/125,-2/125,1,x_1);
	double x18[19]={0.520,3.1,8.0,17.95,28.65,39.62,50.65,78,104.6,156.6,208.6,260.7,312.50,364.4,416.3,468,494,507,520};
	double y18[19]={5.288,9.4,13.84,20.20,24.90,28.44,31.10,35,36.9,36.6,34.6,31.0,26.34,20.9,14.8,7.8,3.7,1.5,0.2};
	double x_2[5]={2,30,130,350,515};
	S(19,x18,y18,1.86548,-0.046115,5,x_2);
	system("pause");
	return 0;
}

double f(double x)
{
	return 1/(1+x*x);
}
double L(int n,double a,double b,double x)
{
	double x_[30],y_[30],l_[30],Ln=0.0;
	for(int i=0;i<=n;i++)
	{
		x_[i]=a+i*(b-a)/n;
		y_[i]=f(x_[i]);
		cout<<"X"<<i<<"="<<x_[i]<<"\tY"<<i<<"="<<y_[i]<<endl;
	}
	for(int i=0;i<=n;i++)
	{
		l_[i]=1.0;//初始化为1.0，否则会溢出
		for(int j=0;j<=n;j++)
			if(i!=j)
				l_[i]*=(x-x_[j])/(x_[i]-x_[j]);
	}
	for(int i=0;i<=n;i++)
		Ln+=y_[i]*l_[i];
	return Ln;
}
void S(int n,double X[],double Y[],double y0,double yn,int x_num,double x[])
{	
	double h[30],u[30],l[30],d[30];
	for(int i=0;i<n;i++)
		h[i]=X[i+1]-X[i];
	for(int i=1;i<n;i++)
	{	
		u[i]=h[i-1]/(h[i-1]+h[i]);
		l[i]=h[i]/(h[i-1]+h[i]);
		d[i]=((f(X[i])-f(X[i+1]))/(X[i]-X[i+1])-(f(X[i])-f(X[i-1]))/(X[i]-X[i-1]))/(h[i]+h[i-1]);
	}	
	l[0]=1;
	d[0]=6*((f(X[0])-f(X[1]))/(X[0]-X[1])-y0)/h[0];
	u[n]=1;
	d[n]=6*(yn-(f(X[n-1])-f(X[n]))/(X[n-1]-X[n]))/h[n-1];
	double M[30],B[30],y[30];
	//追赶法解出M
	B[0]=l[0]/2;
	for(int i=1;i<n-1;i++)
		B[i]=l[i]/(2-u[i]*B[i-1]);
	y[0]=d[0]/2;
	for(int i=1;i<n;i++)
		y[i]=(d[i]-u[i]*y[i-1])/(2-u[i]*B[i-1]);
	M[n-1]=y[n-1];
	for(int i=n-2;i>=0;i--)
		M[i]=y[i]-B[i]*M[i+1];
	//输出M
	cout<<"\nS"<<n-1<<"的情况："<<endl;
	for(int i=0;i<n;i++)
		cout<<"M"<<i<<"="<<M[i]<<endl;
	for(int i=0;i<x_num;i++)
		for(int j=0;j<n;j++)
			//判断所要求的值在哪段区间内
			if(x[i]>=X[j]&&x[i]<=X[j+1])
			{
				double s;
				s=M[j]*(X[j+1]-x[i])*(X[j+1]-x[i])*(X[j+1]-x[i])/(6*h[j])+M[j+1]*(x[i]-X[j])*(x[i]-X[j])*(x[i]-X[j])
					/(6*h[j])+(Y[j]-M[j]*h[j]*h[j]/6)*(X[j+1]-x[i])/h[j]+(Y[j+1]-M[j+1]*h[j]*h[j]/6)*(x[i]-X[j])/h[j];
				cout<<"S("<<x[i]<<")="<<s<<endl;
				if(x_num!=1)
				{
					double s1;
					s1=-M[j]*(X[j+1]-x[i])*(X[j+1]-x[i])/(2*h[j])+M[j+1]*(x[i]-X[j])*(x[i]-X[j])
					/(2*h[j])+(Y[j+1]-Y[j])/h[j]-(M[j+1]-M[j])/(6*h[j]);
					cout<<"S'("<<x[i]<<")="<<s1<<endl;
					double s2;
					s2=M[j]*(X[j+1]-x[i])/h[j]+M[j+1]*(x[i]-X[j])/h[j];
					cout<<"S''("<<x[i]<<")="<<s2<<endl;
				}
			}
}