t3_1.cpp

求解三次样条函数思路为：由 连续 预设 , 为一次多项式：故 积分2次 积分常数（2个） 由插值条件 确定 积分常数 得 （含预设的 ） 利用 连续： 确定 的 个方程 + 2边界条件 确定

//  cpp  t3_1源代码  上机练习题3_1"
#include <iostream.h>
#include <iomanip.h>
#include <math.h>
const int N=20;
static double x[N+1],fx[N+1],yN[N+1],xc[N+1],yc[N+1],m[N+1];
void Calculate(int n);      //函数声明
void Splinem(int n);             //函数声明
int main()
{
	cout<<"上机练习题3_1"<<endl;
	cout<<"当n=5时的计算结果为:////////////////////////////////////"<<endl;
	Calculate(5);
	cout<<"当n=10时的计算结果为://////////////////////////////////"<<endl;
	Calculate(10);
	cout<<"当n=20时的计算结果为://////////////////////////////////"<<endl;
	Calculate(20);
	return 1;
}
void Calculate(int n)
{
	int i,j,k;
	double nN[N+1],s[N+1];
	double x1,x2,h;
	double en,es;
	cout<<"a:"<<endl;            /////////////////////////a
	for(i=0;i<=n;i++)
	{
		x[i]=-1+2.0*i/n;
		fx[i]=1/(1+25*x[i]*x[i]);
		yN[i]=fx[i];
		cout<<"f("<<x[i]<<")="<<fx[i]<<";";
		if((i+1)%5==0)cout<<endl;
	}
	cout<<endl;
	cout<<"b:"<<endl;         //////////////////////////////////b	 
	for(k=1;k<=n;k++)        //////////////Newton插值
		for(i=n;i>=k;i--)
			yN[i]=(yN[i]-yN[i-1])/(x[i]-x[i-k]);	
	cout<<"Newton插值多项式为:"<<endl;
	for(i=0;i<=n;i++)
	{
		cout<<"yN["<<i<<"]="<<yN[i]<<";";
		if((i+1)%5==0)cout<<endl;
	}
	cout<<endl;
	Splinem(n);   ////样条插值
	if(n==5||n==20)                 /////////////////c
	{
		cout<<"c:"<<endl;
		for(k=1;k<=20;k++)     //////计算f(xk)
			if(k>=11)
			{
				xc[k]=-1+0.02*(k+79);
				yc[k]=1/(1+25*xc[k]*xc[k]);


			}
			else
			{
				xc[k]=-1+0.02*k;
				yc[k]=1/(1+25*xc[k]*xc[k]);
			}
		for(k=1;k<=20;k++)         //////////计算N(xk)
		{
			nN[k]=yN[n];
				for(i=n-1;i>=0;i--)
					nN[k]=nN[k]*(xc[k]-x[i])+yN[i];			
		}		
		for(i=1;i<=20;i++)      ////////计算s(xk)
		{
			k=1;
			for(j=1;j<n;j++)
				if(xc[i]<=x[j])
				{
					k=j;
					break;
				}
				else
					k=j+1;		
			h=x[k]-x[k-1];
			x1=x[k]-xc[i];
			x2=xc[i]-x[k-1];
			s[i]=m[k-1]*(x1*x1*x1)/6+m[k]*(x2*x2)/6;
			s[i]+=(fx[k-1]-m[k-1]*h*h/6)*x1;
			s[i]+=(fx[k]-m[k]*h*h/6)*x2;
			s[i]/=h;
		}

		cout<<"输出y(xk):"<<endl;
		for(i=1;i<=20;i++)                   ////输出y(xk)
		{                           
			cout<<"yc["<<setw(5)<<xc[i]<<"]="<<setprecision(5)<<yc[i]<<";";			
			if(i%4==0)cout<<endl;
		}
		cout<<"输出N(xk):"<<endl;
		for(i=1;i<=20;i++)                     ////输出N(xk)
		{
			cout<<"N["<<setw(5)<<xc[i]<<"]="<<setw(9)<<nN[i]<<";";		
			if(i%4==0)cout<<endl;
		}
		cout<<"输出s(xk):"<<endl;
		for(i=1;i<=20;i++)                      	////输出s(xk)
		{
			cout<<"s["<<setw(5)<<xc[i]<<"]="<<setw(9)<<s[i]<<";";	
			if(i%4==0)cout<<endl;
		}
		cout<<"d:"<<endl; 
		en=fabs(yc[1]-nN[1]);
		es=fabs(yc[1]-s[1]);
		for(k=2;k<=20;k++)
		{
			if(en<fabs(yc[k]-nN[k]))
				en=fabs(yc[k]-nN[k]);
			if(es<fabs(yc[k]-s[k]))
				es=fabs(yc[k]-s[k]);
		}
		cout<<"E(N)="<<en<<endl;
		cout<<"E(s)="<<es<<endl;
	}
}
void Splinem(int n)     /////三次样条插值
{
	int i,k;
	double h[N+1];
	double a[N+1],b[N+1],c[N],u[N+1],l[N+1],yy[N+1];
	for(i=0;i<=n;i++)
		m[i]=fx[i]; 
	for(k=1;k<=2;k++)
		for(i=n;i>=k;i--)
			m[i]=(m[i]-m[i-1])/(x[i]-x[i-1]);
	h[1]=x[1]-x[0];
	for(i=1;i<n;i++)
	{
		h[i+1]=x[i+1]-x[i];
		c[i]=h[i+1]/(h[i]+h[i+1]);
		a[i]=1-c[i];
		b[i]=2;
		m[i]=6*m[i+1];
	}
		
	m[0]=0;m[n]=0;c[0]=0;
	b[0]=2;b[n]=2;a[n]=0;
	///////////追赶法求解
	u[0]=b[0];
	yy[0]=m[0];
	for(k=1;k<=n;k++)
	{
		l[k]=a[k]/u[k-1];
		u[k]=b[k]-l[k]*c[k-1];
		yy[k]=m[k]-l[k]*yy[k-1];	
	}
	m[n]=yy[n]/u[n];
	for(k=n-1;k>=0;k--)
	{
		m[k]=(yy[k]-c[k]*m[k+1])/u[k];	
	}
	cout<<"样条插值多项式s(x)为:"<<endl;
	for(i=0;i<=n;i++)
	{
		cout<<"M["<<i<<"]="<<m[i]<<";";	
		if((i+1)%5==0)cout<<endl;
	}
	cout<<endl;
}