计算方法 实验一 插值.docx
- 文档编号:9414412
- 上传时间:2023-02-04
- 格式:DOCX
- 页数:16
- 大小:156.22KB
计算方法 实验一 插值.docx
《计算方法 实验一 插值.docx》由会员分享,可在线阅读,更多相关《计算方法 实验一 插值.docx(16页珍藏版)》请在冰豆网上搜索。
计算方法实验一插值
实验一插值
专业班级:
姓名:
学号:
1)实验目的
•熟悉C语言编程;
•学习拉格朗日插值,多项式插值,和自然样条插值方法及程序设计算法
2)实验题目
下面给出美国从1920年到1970年的人口表:
1.用表中数据构造一个5次拉格朗日插值多项式,并用此估计1910,1965和2002年的人口。
在1910年的实际人口约为91772000,请判断插值计算得到的1965年和2002年的人口数据准确性是多少?
2.用牛顿插值估计:
(1)1965年的人口数;
(2)2002年的人口数。
3.用自然样条函数估计在1910,1965和2002年的人口数。
4.请比较以上三种方法所求值的效果。
那一种方法最优?
3)实验原理与理论基础
在本次实验当中主要运用了以下原理及公式:
Ø拉格朗日插值
插值基函数:
:
Ø牛顿插值
牛顿插值公式:
Pn(x)=a0+a1(x-x0)+…+an(x-x0)(x-x1)…(x-xn-1)
Ø自然样条函数
Ø
4)实验内容
联系所学的知识,按照实验题目要求进行编程以及对数值的分析,能够更加深刻的理解拉格朗日插值,多项式插值,和自然样条插值方法及程序设计算法。
5)实验结果
Ø拉格朗日插值
Ø牛顿插值
Ø自然样条函数
6)实验结果分析与小结
1.使用拉格朗日插值进行计算时,首先计算插值基函数,再构造拉格朗日插值多项式对数据进行计算。
进行准确性分析,采用事后误差估计。
1910年的实际人口约为91772000,判断插值计算得到的1965年和2002年的人口数据准确性。
~~
共有7组数据则取前6组得L(1965)≈178340.5L(2002)≈-3964442.345
取后6组得L(1965)=193082.5
则R(x)≈
R(1965)≈[55/(-60)]*(193081.5-178340.5)≈-13512.5
x=2002时,所得结果明显偏离事实甚远。
可见拉格朗日插值的求值效果不理想。
2.N(1965)≈193096.277344与拉格朗日插直球的结果相近
N(2002)≈15823.279616仍然偏离事实
3.S(1910)≈88219.0
S(1965)≈191873.430622与前两个插值求得的结果相近
S(2002)≈157868.777837比前两个插值求得的结果靠近与事实
在自然样条插值计算的编程中,我用到了Courant法求解线性方程组的方法,主要是用来解决Mn的求值问题。
4.比较以上实验结果,并加以分析可得,自然样条插值相对其他更好。
但是,脱离了给出的x的区间,在两端延伸的区间内进行数值计算所得到的数值依然不够精确,不够理想。
只有在已知[Xi,Xi+1]的区间内计算,才能够得到比较满意的数值。
程序详见附录。
附录:
编程语言采用c语言,运行环境为win-tc。
1.拉格朗日插值
#include
#defineMAX_N20
typedefstructtagPOINT
{
doublex;
doubley;
}POINT;
intmain()
{
intn;
inti,j;
POINTpoints[MAX_N+1];
doublediff[MAX_N+1];
doublex,tmp,newton=0;
printf("\nInputnvalue:
");
scanf("%d",&n);
if(n>MAX_N)
{
printf("TheinputnislargerthenMAX_N,pleaseredefinetheMAX_N.\n");
return1;
}
if(n<=0)
{
printf("Pleaseinputanumberbetween1and%d\n",MAX_N);
return1;
}
printf("Nowinputthe(x_i,y_i),i=0,...,%d:
\n",n);
for(i=0;i<=n;i++)
scanf("%lf%lf",&points[i].x,&points[i].y);
printf("Nowinputthexvalue:
");
scanf("%lf",&x);
for(i=0;i<=n;i++)diff[i]=points[i].y;
for(i=0;i { for(j=n;j>i;j--) { diff[j]=(diff[j]-diff[j-1])/(points[j].x-points[j-1-i].x); } } tmp=1; newton=diff[0]; for(i=0;i { tmp=tmp*(x-points[i].x); newton=newton+tmp*diff[i+1]; } printf("newton(%f)=%f\n",x,newton); getch(); return0; } 2.牛顿插值 #include #defineMAX_N20 typedefstructtagPOINT { doublex; doubley; }POINT; intmain() { intn; inti,j; POINTpoints[MAX_N+1]; doublediff[MAX_N+1]; doublex,tmp=0,lagrange=0,tx,ty; printf("\nInputnvalue: "); scanf("%d",&n); if(n>MAX_N) { printf("TheinputnislargerthenMAX_N,pleaseredefinetheMAX_N.\n");getch(); return1; } if(n<=0) { printf("Pleaseinputanumberbetween1and%d\n",MAX_N);getch(); return1; } printf("Nowinputthe(x_i,y_i),i=0,...,%d: \n",n); for(i=0;i<=n;i++) scanf("%lf%lf",&points[i].x,&points[i].y); printf("Nowinputthexvalue: "); scanf("%lf",&x); for(i=0;i<=n;i++) { diff[i]=0; tx=1; ty=1; for(j=0;j<=n;j++) { if(i! =j) { tx=tx*(x-points[j].x); ty=ty*(points[i].x-points[j].x); } } diff[i]=tx/ty; } for(i=0;i<=n;i++) { tmp=points[i].y*diff[i];printf("%f",tmp); lagrange+=tmp; } printf("lagrange(%f)=%f\n",x,lagrange); getch(); return0; } 3.自然样条函数 #include #include #defineMAX_N20 typedefstructtagPOINT { doublex; doubley; }POINT; intmain() { intn; inti,j; POINTpoints[MAX_N+1]; doublex,tmp,sa=0,sb=0,s,t[3]; doubleh[MAX_N],v[MAX_N],u[MAX_N],M[MAX_N+1],d[MAX_N]; /*courant中所用变量*/ intn2; inti2,j2,k2; intmi2; doublemx2,tmp2; staticdoublea2[MAX_N][MAX_N],b2[MAX_N],x2[MAX_N],y2[MAX_N]; staticdoublel2[MAX_N][MAX_N],u2[MAX_N][MAX_N]; printf("\nInputnvalue: "); scanf("%d",&n); if(n>MAX_N) { printf("TheinputnislargerthenMAX_N,pleaseredefinetheMAX_N.\n"); return1; } if(n<=0) { printf("Pleaseinputanumberbetween1and%d\n",MAX_N); return1; } printf("Nowinputthe(x_i,y_i),i=0,...,%d: \n",n); for(i=0;i<=n;i++) scanf("%lf%lf",&points[i].x,&points[i].y); printf("Nowinputthexvalue: "); scanf("%lf",&x); for(i=0;i { h[i]=points[i+1].x-points[i].x; } for(i=1;i { v[i]=h[i]/(h[i]+h[i-1]); u[i]=1-v[i]; } for(i=1;i { tmp=0; tmp+=points[i-1].y/((points[i-1].x-points[i].x)*(points[i-1].x-points[i+1].x)); tmp+=points[i].y/((points[i].x-points[i-1].x)*(points[i].x-points[i+1].x)); tmp+=points[i+1].y/((points[i+1].x-points[i].x)*(points[i+1].x-points[i-1].x)); d[i]=tmp*6; } /*线性方程组求解解放在M[i]中,M[0]=M[n]=0*/ n2=n-1; M[0]=0; M[n]=0; for(i2=0;i2 for(j2=0;j2 a2[i2][j2]=0; for(i2=0;i2 { if(i2-1>=0)a2[i2][i2-1]=u[i2+1]; a2[i2][i2]=2; if(i2+1 } for(i2=0;i2 b2[i2]=d[i2+1]; /*courant*/ for(i2=0;i2 for(k2=0;k2 { for(i2=k2;i2 { l2[i2][k2]=a2[i2][k2]; for(j2=0;j2<=k2-1;j2++) l2[i2][k2]-=(l2[i2][j2]*u2[j2][k2]); } for(j2=k2+1;j2 { u2[k2][j2]=a2[k2][j2]; for(i2=0;i2<=k2-1;i2++) u2[k2][j2]-=(l2[k2][i2]*u2[i2][j2]); u2[k2][j2]/=l2[k2][k2]; } } for(i2=0;i2 { y2[i2]=b2[i2]; for(j2=0;j2<=i2-1;j2++) y2[i2]-=(l2[i2][j2]*y2[j2]); y2[i2]/=l2[i2][i2]; } for(i2=n2-1;i2>=0;i2--) { x2[i2]=y2[i2]; for(j2=i2+1;j2 x2[i2]-=(u2[i2][j2]*x2[j2]); } for(i2=0;i2 i=0; while(x>=points[i].x&&i if(i>0)i--; sa=points[i+1].x-x; sb=x-points[i].x; s=M[i]*sa*sa*sa/(h[i]*6)+M[i+1]*sb*sb*sb/(h[i]*6)+ (points[i].y-M[i]*h[i]*h[i]/6)*sa/h[i]+ (points[i+1].y-M[i+1]*h[i]*h[i]/6)*sb/h[i]; printf("S(%f)=%f\n",x,s); return0; }
- 配套讲稿:
如PPT文件的首页显示word图标,表示该PPT已包含配套word讲稿。双击word图标可打开word文档。
- 特殊限制:
部分文档作品中含有的国旗、国徽等图片,仅作为作品整体效果示例展示,禁止商用。设计者仅对作品中独创性部分享有著作权。
- 关 键 词:
- 计算方法 实验一 插值 实验