//本题是上机题的第五大题
#include"iostream.h"
#include"stdlib.h"
#include"math.h"
void main()
{
void fx(int N,double *p);
void dxs(int N,double *p1);
void scytiao(int N,double *p,double *p1);
void calcu(int N,double *p,double *p1,double *p2,double *p3);//计算各拟合函数
void Esss(double *p,double *p1);
double a5[6][7],a10[11][12],a20[21][22];//计算所得xi 与yi的值，以及差商表
double M5[6],M10[11],M20[21];//三次样条差值，利用M样条插值，开辟各M存储空间
double v[20];//保存第三小题的自变量值
double r5[3][20],r20[3][20]; //开辟各拟合函数值的空间，第一行为已知函数值，第二行为牛顿拟合函数值，第三行为三次样条拟合函数值
double es5[2],es20[2];//误差，数组第一个数为真值与牛顿的最大误差，第二项为真值与三次插值的最大误差
int a[3]={5,10,20},i,j,k;
for(i=1;i<=10;i++) v[i-1]=(double)2*i/100-1;
for(i=90;i<=99;i++) v[i-80]=(double)2*i/100-1;
for(i=0;i<=2;i++)
{
j=a[i];
switch(j)
{
case 5:
{
fx(5,a5[0]);//N为5的第一问
dxs(5,a5[0]);//N为5的差商表
scytiao(5,a5[0],M5);//样条参数计算函数
calcu(5,a5[0],M5,r5[0],v);//计算各值
Esss(r5[0],es5);
}
continue;
case 10:
{
fx(10,a10[0]);
dxs(10,a10[0]);
scytiao(10,a10[0],M10);
}
continue;
case 20:
{
fx(20,a20[0]);
dxs(20,a20[0]);
scytiao(20,a20[0],M20);
calcu(20,a20[0],M20,r20[0],v);//计算各值
Esss(r20[0],es20);
}
continue;
}
}
}
void fx(int N,double *p)
{
double x,y;
int i,j,k;
for(i=0;i<=N;i++)
{
x=(-1)+(double)(2*i)/N;
*(p+(N+2)*i+0)=x;
y=1/(1+25*x*x);
*(p+(N+2)*i+1)=y;
}
cout<<"######################################################################################################"<<endl;
cout<<"当N="<<N<<"时，自变量xi与对应函数值yi为："<<endl;
for(i=0;i<=N;i++)
{
cout<<"x"<<i<<"="<<*(p+(N+2)*i+0)<<endl;
}
for(i=0;i<=N;i++)
{
cout<<"y"<<i<<"="<<*(p+(N+2)*i+1)<<endl;
}
}
void dxs(int N,double *p1)
{
double x,y;
int i,j,k;
for(i=2;i<=N+1;i++)
{
for(j=i-1;j<=N;j++)
{
x=(*(p1+(N+2)*j+(i-1))-*(p1+(N+2)*(j-1)+(i-1)))/(*(p1+(N+2)*j+0)-*(p1+(N+2)*(j-i+1)+0));
*(p1+(N+2)*j+i)=x;
}
}
/* cout<<"这是N="<<N<<"的差商表："<<endl;
for(i=0;i<=N;i++)
{
for(j=0;j<=N+1;j++)
{
cout<<*(p1+(N+2)*i+j)<<" ";
}
cout<<endl;
}*/
cout<<"现在写出他的多项式系数:"<<endl;
for(i=1;i<=N+1;i++)
{
j=i-1;
cout<<*(p1+(N+2)*j+i)<<" ";
}
}
void scytiao(int N,double *p,double *p1)
{
int i,j,k;
double x,y,*pt[30],*pt1[30],*pt2[30];
for(i=0;i<=N-2;i++)
pt[i]=new double[N];//弯矩矩阵形式，即追赶法形式，最后一列存储di的值共N-1行N列即为增广矩阵
for(i=0;i<=N-2;i++)
pt1[i]=new double[N-1];//追赶法的L对角下三角单位矩阵
for(i=0;i<=N-2;i++)
pt2[i]=new double[N];//追赶法的R单位上三角单位矩阵
x=*(p+(N+2)*0+0);
y=x*x;
y=(4*(1+25*y)*(1+25*y)-200*y*(1+25*y))/((1+25*y)*(1+25*y)*(1+25*y)*(1+25*y));//计算M0的二阶导数
*(p1+0)=y;//M0
x=*(p+(N+2)*N+0);
y=(4*(1+25*x*x)*(1+25*x*x)-200*x*x*(1+25*x*x))/((1+25*x*x)*(1+25*x*x)*(1+25*x*x)*(1+25*x*x));//计算M0的二阶导数
*(p1+N)=y;//计算Mn的二阶导数值
x=(*(p+(N+2)*1+0)-*(p+(N+2)*0+0))/(*(p+(N+2)*2+0)-*(p+(N+2)*0+0));//计算u1的值
y=6*(*(p+(N+2)*2+3));//计算d1的值
y=y-x*(*(p1+0));//计算新的的d1的值
pt[0][N-1]=y;//保存新的d1的值
pt[0][0]=2.000;//
pt[0][1]=(*(p+(N+2)*2+0)-*(p+(N+2)*1+0))/(*(p+(N+2)*2+0)-*(p+(N+2)*0+0));//计算兰布拉1的值并保存
x=(*(p+(N+2)*N+0)-*(p+(N+2)*(N-1)+0))/(*(p+(N+2)*N+0)-*(p+(N+2)*(N-2)+0));//计算兰布拉N-1的值
y=6*(*(p+(N+2)*N+3));//计算dn-1的值
y=y-x*(*(p1+N));//计算新的的dn-1的值
pt[N-2][N-1]=y;//保存新的dn-1值
pt[N-2][N-3]=(*(p+(N+2)*(N-1)+0)-*(p+(N+2)*(N-2)+0))/(*(p+(N+2)*N+0)-*(p+(N+2)*(N-2)+0));//计算un-1的值保存
pt[N-2][N-2]=2.00;//
for(i=1;i<=N-3;i++)
{
pt[i][i]=2.00;//
x=(*(p+(N+2)*(i+1)+0)-*(p+(N+2)*i+0))/(*(p+(N+2)*(i+2)+0)-*(p+(N+2)*i+0));//计算ui+1的值保存
pt[i][i-1]=x;
pt[i][i+1]=1-x;//计算兰布拉i+1的值并保存
pt[i][N-1]=(*(p+(N+2)*(i+2)+3))*6;//计算di+1的值，
}// 这样追赶法弯矩矩阵就计算完了
pt1[0][0]=pt[0][0];
pt2[0][N-1]=pt[0][N-1]/pt1[0][0];
for(i=1;i<=(N-2);i++)
{
pt2[i-1][i]=pt[i-1][i]/pt1[i-1][i-1];//计算r的值
pt1[i][i]=pt[i][i]-pt[i][i-1]*pt2[i-1][i];//计算l的值
pt2[i][N-1]=(pt[i][N-1]-pt[i][i-1]*pt2[i-1][N-1])/pt1[i][i];//计算y的值
} //计算完成追赶矩阵的各个分解矩阵，下面要计算Mi的各个数值
*(p1+(N-1))=pt2[N-2][N-1];
for(i=N-2;i>=1;i--)
*(p1+i)=pt2[i-1][N-1]-pt2[i-1][i]*(*(p1+(i+1)));
cout<<endl;
cout<<"显示各个Mi的值："<<endl;
for(i=0;i<=N;i++)
cout<<*(p1+i)<<endl;
cout<<endl;
/*cout<<"输出上三角阵"<<endl;
for(i=0;i<=N-2;i++)
{ for(j=0;j<=N-1;j++)
{cout<<pt2[i][j]<<" ";}
cout<<endl;
}
cout<<"输出下三角阵"<<endl;
for(i=0;i<=N-2;i++)
{ for(j=0;j<=N-2;j++)
{cout<<pt1[i][j]<<" ";}
cout<<endl;
}*/
for(i=0;i<=N-2;i++)
{
delete []pt[i];
delete []pt1[i];
delete []pt2[i];
}
} //三次样条插值的各M值计算完毕
void calcu(int N,double *p,double *p1,double *p2,double *p3)
{
double sum,Nx,hi,h1,h2;
int i,j,k;
cout<<"输出第三那小题的各自变量x的值"<<endl;
for(i=0;i<=19;i++)
cout<<*(p3+i)<<" ";
cout<<endl;
for(i=0;i<=19;i++)
{
sum=*(p+(N+2)*0+1);
Nx=1;
for(j=0;j<=N-1;j++)
{
Nx=Nx*(*(p3+i)-*(p+(N+2)*j+0));//*(*(p+(N+2)*(j+2)+(j+1)));//(*(p+(N+2)*(j+1)+j));//由前一项推出后一项各项再相加就可得出
sum=sum+Nx*(*(p+(N+2)*(j+1)+(j+2)));
// cout<<(*(p3+i)-*(p+(N+2)*j+0))<<" ";
}
*(p2+20*0+i)=1/(1+25*(*(p3+i))*(*(p3+i)));//计算已知函数的真值
*(p2+20*1+i)=sum;//得出牛顿的值
// cout<<sum<<" ";
// break;
}//下面将做三才样条插值法
for(i=0;i<=19;i++)
{
for(j=1;j<=N;j++)
{
if(!((*(p3+i)>=(*(p+(N+2)*(j-1)+0)))&&(*(p3+i)<=(*(p+(N+2)*j+0)))))
{ // cout<<j<<" "<<endl;
continue;
}
else
{
hi=*(p+(N+2)*j+0)-*(p+(N+2)*(j-1)+0);
h1=*(p+(N+2)*j+0)-*(p3+i);
h2=*(p3+i)-*(p+(N+2)*(j-1)+0);
sum=(h1*h1*h1*(*(p1+(j-1)))+h2*h2*h2*(*(p1+j)))/(6*hi)+(*(p+(N+2)*(j-1)+1)-hi*hi*(*(p1+(j-1)))/6)*h1/hi+(*(p+(N+2)*j+1)-hi*hi*(*(p1+j))/6)*h2/hi;
// cout<<*(p1+(j-1))<<" "<<endl;
*(p2+20*2+i)=sum;
break;
}
}
}
cout<<endl;
cout<<"输出真值、牛顿值的数据和三项插值的函数值"<<endl;
for(i=0;i<=2;i++)
{
for(j=0;j<=19;j++)
cout<<*(p2+20*i+j)<<" ";
cout<<endl;
}
}
void Esss(double *p,double *p1)
{
double max,x;
int i,j,k;
max=0;
for(i=0;i<=19;i++)
{
x=fabs(*(p+20*0+i)-*(p+20*1+i));
// cout<<x<<endl;
if(x>=max) max=x;
}
*(p1+0)=max;
max=0;
for(i=0;i<=19;i++)
{
x=fabs(*(p+20*0+i)-*(p+20*2+i));
if(x>=max) max=x;
}
*(p1+1)=max;
cout<<endl;
cout<<"输出误差："<<endl;
cout<<"真值与牛顿的最大误差为："<<*(p1+0)<<endl;
cout<<"真值与三次样条插值的最大误差为："<<*(p1+1)<<endl;
}
