计算方法 实验一 插值.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;i2for（j2=0;j2a2[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;i2b2[i2]=d[i2+1];
/*courant*/
for（i2=0;i2for（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;j2x2[i2]-=（u2[i2][j2]*x2[j2]）;
}
for（i2=0;i2i=0;
while（x>=points[i].x&&iif（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;
}
升级会员 