insertvalue.cpp

本人的计算方法课程作业

//insertvalue.cpp--牛顿和三次样条多项式插值
#include<iostream>
#include<cmath>
using namespace std;
const int N=20;
float func(float x);
void print_array(float a[][N+1]);
float Newton(float x[N+1],float bb[][N+1],float xx);
float EVASPLINE(float x[N+1],float y[N+1],float m[N+1],float xx);
void main()
{
	static float x[N+1],y[N+1],aa[N+1][N+1],bb[N+1][N+1];
    float a[N+1],b[N+1],c[N+1],d[N+1],m[N+1],h[N+1];
	int i,j,k;
    for(i=0;i<=N;i++)       //求函数f(x)在各节点处的值
	{
        x[i]=i*(2.0/N)-1;
	    y[i]=func(x[i]);
	}
    cout<<"函数f(x)在各节点处的值为:"<<endl;
	for(i=0;i<=N;i++)       //输出函数f(x)在各节点处的值
	   cout<<"i="<<i<<"  xi="<<x[i]<<"  f("<<x[i]<<")="<<y[i]<<endl;
    for(i=0;i<=N;i++)
		aa[i][0]=y[i];
	for(j=1;j<=N;j++)
	{
		for(i=j;i<=N;i++)
			aa[i][j]=(aa[i][j-1]-aa[i-1][j-1])/(x[i]-x[i-j]);  //求牛顿插值的各阶差商		 
	}                                                          //结果存放在二维数组aa[][]中
    print_array(aa);        

//计算三次样条插值的参数Mi，存放于数组m[]中,其它各参数根据第三边界条件确定
	for(i=0;i<=N;i++)
		m[i]=y[i];
	for(k=1;k<=2;k++)
	{
		for(i=N;i>=k;i--)
			m[i]=(m[i]-m[i-1])/(x[i]-x[i-k]);
	}
	h[1]=x[1]-x[0];
	for(i=1;i<=N-1;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;
		d[i]=6*m[i+1];
	}
	d[0]=(-12)*h[1]*(m[3]-m[2])/(x[3]-x[0]);
	d[N]=12*h[N]*(m[N]-m[N-1])/(x[N]-x[N-3]);
	c[0]=-2;
	b[0]=2;
	a[N]=-2;
	b[N]=2;

//调用TSS算法计算Mi，系数存放于数组a,b,c中，
//d中是右端向量，计算结果存放于数组m[]中
	m[0]=d[0];
	for(k=2;k<=N+1;k++)
	{
		a[k-1]=a[k-1]/b[k-2];
		b[k-1]=b[k-1]-a[k-1]*c[k-2];
		m[k-1]=d[k-1]-a[k-1]*m[k-2];
	}
	m[N]/=b[N];
	for(k=N;k>=1;k--)
		m[k-1]=(m[k-1]-c[k-1]*m[k])/b[k-1];
	cout<<"--------------------------------------------------------"<<endl;
	cout<<"输出三次样条插值的Mi参数:"<<endl;
	for(i=0;i<=N;i++)
		cout<<"M"<<i<<"="<<m[i]<<endl;
//在得到牛顿和三次样条插值多项式之后，输出运算结果
	cout<<"利用各插值函数，计算给定节点处的相应函数值如下:"<<endl;
	float xi;
	for(i=1;i<=10;i++)
	{
		xi=i*1.0/50-1;
		cout<<"k="<<i<<"  xi="<<xi<<"  yk="<<func(xi)
			<<"  nk="<<Newton(x,aa,xi)<<"  sk="<<EVASPLINE(x,y,m,xi)<<endl;
	}
	cout<<"--------------------------------------------------------"<<endl;
	for(i=90;i<=99;i++)
	{
		xi=i*1.0/50-1;
		cout<<"k="<<i<<"  xi="<<xi<<"  yk="<<func(xi)
			<<"  nk="<<Newton(x,aa,xi)<<"  sk="<<EVASPLINE(x,y,m,xi)<<endl;
	}
	float En=0,Es=0;
    for(i=1;i<=10;i++)        //求牛顿和三次样条插值的最大误差
	{
		xi=i*1.0/50-1;
		if(En<fabs(Newton(x,aa,xi)-func(xi)))
			En=fabs(Newton(x,aa,xi)-func(xi));
		if(Es<fabs(EVASPLINE(x,y,m,xi)-func(xi)))
			Es=fabs(EVASPLINE(x,y,m,xi)-func(xi));
	}
	cout<<"牛顿插值和三次样条插值的最大误差分别为:"<<endl;
	
	cout<<"E(Nn)="<<En<<endl;
	cout<<"E(Sn)="<<Es<<endl;
}
//func()给出原函数的定义
float func(float x)
{
	float f;
	if(x>=-1&&x<=1) 
		f=1/(25*x*x+1);
	return f;
}
//print_array()用于输出牛顿插值的差商表
void print_array(float aa[N+1][N+1])
{
	int i,j;
	cout<<"--------------------------------------"<<endl;
	cout<<"输出牛顿插值差商表:"<<endl;
	for(i=0;i<=N;i++)
	{
	    for(j=0;j<=i;j++)
		{
			cout.precision(6);
		    cout<<aa[i][j]<<"  ";
		}
		cout<<endl;
	}
}
//Newton()用牛顿插值公式计算给定插值数据点的值
float Newton(float x[N+1],float bb[][N+1],float xx)
{
	int i,j;
	float n=bb[0][0],t;
    for(i=1;i<=N;i++)
	{   
		t=bb[i][i];
        for(j=1;j<=i;j++)
            t=t*(xx-x[j-1]);
	    n+=t;
	}
    return n;
}
//FINDK()找出给定插值数据点所在的区间(x[k-1],x[k]),结果返回下标值k
int FINDK(float xx,float x[])
{
	int k=1,i;
	for(i=1;i<=N;i++)
	{
		if (xx<=x[i])
		{ k=i; goto end;}
	    else
		{ k=i+1;}
	}
end:    ;
	  return k;
}
//EVASPLINE()用三次样条插值函数计算各给定数据点对应的函数值
float EVASPLINE(float x[N+1],float y[N+1],float m[N+1],float xx)
{
	int k;
	float s,h;
	k=FINDK(xx,x);
	h=x[k]-x[k-1];
	if(k>=1)
		s=(m[k-1]*pow((x[k]-xx),3)/6
		    +m[k]*pow((xx-x[k-1]),3)/6
		    +(y[k-1]-m[k-1]*h*h/6)*(x[k]-xx)
			+(y[k]-m[k]*h*h/6)*(xx-x[k-1]))/h;
	
	return s;
}