计算方法 实验一 插值.docx

1、计算方法 实验一 插值实验一 插值专业班级： 姓名： 学号： 1）实验目的 熟悉C语言编程； 学习拉格朗日插值,多项式插值,和自然样条插值方法及程序设计算法2）实验题目下面给出美国从1920年到1970年的人口表：1. 用表中数据构造一个5次拉格朗日插值多项式，并用此估计1910，1965和2002年的人口。在1910年的实际人口约为91772000，请判断插值计算得到的1965年和2002年的人口数据准确性是多少？2. 用牛顿插值估计：（1）1965年的人口数；（2）2002年的人口数。3. 用自然样条函数估计在1910，1965和2002年的人口数。4. 请比较以上三种方法所求值的效果。那

2、一种方法最优？3）实验原理与理论基础在本次实验当中主要运用了以下原理及公式： 拉格朗日插值 插值基函数: 牛顿插值牛顿插值公式：Pn (x)=a0+a1(x-x0)+ an(x-x0)(x-x1)(x-xn-1) 自然样条函数 4）实验内容联系所学的知识，按照实验题目要求进行编程以及对数值的分析，能够更加深刻的理解拉格朗日插值,多项式插值,和自然样条插值方法及程序设计算法。5）实验结果 拉格朗日插值 牛顿插值 自然样条函数6）实验结果分析与小结1. 使用拉格朗日插值进行计算时，首先计算插值基函数，再构造拉格朗日插值多项式对数据进行计算。进行准确性分析，采用事后误差估计。1910年的实际人口约为

3、91772000，判断插值计算得到的1965年和2002年的人口数据准确性。 共有7组数据则 取前6组得 L（1965）178340.5 L（2002） -3964442.345取后6组得L（1965）=193082.5则R（x） R（1965） 55/（-60）*（193081.5 - 178340.5）-13512.5x=2002时，所得结果明显偏离事实甚远。可见拉格朗日插值的求值效果不理想。2. N(1965) 193096.277344 与拉格朗日插直球的结果相近 N(2002)15823.279616 仍然偏离事实3. S(1910)88219.0 S(1965)191873.430

4、622 与前两个插值求得的结果相近S(2002)157868.777837 比前两个插值求得的结果靠近与事实在自然样条插值计算的编程中，我用到了Courant法求解线性方程组的方法，主要是用来解决Mn的 求值问题。4. 比较以上实验结果，并加以分析可得，自然样条插值相对其他更好。但是，脱离了给出的x的区间，在两端延伸的区间内进行数值计算所得到的数值依然不够精确，不够理想。只有在已知Xi，Xi+1的区间内计算，才能够得到比较满意的数值。程序详见附录。附录: 编程语言采用c语言，运行环境为win-tc。1. 拉格朗日插值#include #define MAX_N 20typedef struct

5、 tagPOINT double x; double y;POINT;int main() int n; int i,j; POINT pointsMAX_N+1; double diffMAX_N+1; double x,tmp,newton=0; printf(nInput n value:); scanf(%d,&n); if(nMAX_N) printf(The input n is larger then MAX_N,please redefine the MAX_N.n); return 1; if(n=0) printf(Please input a number between

6、 1 and %dn,MAX_N); return 1; printf(Now input the (x_i,y_i),i=0,.,%d:n,n); for(i=0;i=n;i+) scanf(%lf%lf,&pointsi.x,&pointsi.y); printf(Now input the x value:); scanf(%lf,&x); for(i=0;i=n;i+)diffi=pointsi.y; for(i=0;ii;j-) diffj=(diffj-diffj-1)/(pointsj.x-pointsj-1-i.x); tmp=1; newton=diff0; for(i=0;

7、in;i+) tmp=tmp*(x-pointsi.x); newton=newton+tmp*diffi+1; printf(newton(%f)=%fn,x,newton); getch(); return 0;2. 牛顿插值#include #define MAX_N 20typedef struct tagPOINT double x; double y;POINT;int main() int n; int i,j; POINT pointsMAX_N+1; double diffMAX_N+1; double x,tmp=0,lagrange=0,tx,ty; printf(nIn

8、put n value:); scanf(%d,&n); if(nMAX_N) printf(The input n is larger then MAX_N,please redefine the MAX_N.n);getch(); return 1; if(n=0) printf(Please input a number between 1 and %dn,MAX_N); getch(); return 1; printf(Now input the (x_i,y_i),i=0,.,%d:n,n); for(i=0;i=n;i+) scanf(%lf%lf,&pointsi.x,&poi

9、ntsi.y); printf(Now input the x value:); scanf(%lf,&x); for(i=0;i=n;i+) diffi=0; tx=1; ty=1; for(j=0;j=n;j+) if(i!=j) tx=tx*(x-pointsj.x); ty=ty*(pointsi.x-pointsj.x); diffi=tx/ty; for(i=0;i=n;i+) tmp=pointsi.y*diffi;printf(%f,tmp); lagrange+=tmp; printf(lagrange(%f)=%fn,x,lagrange); getch(); return

10、 0;3. 自然样条函数#include #include #define MAX_N 20typedef struct tagPOINT double x; double y;POINT;int main() int n; int i,j; POINT pointsMAX_N+1; double x,tmp,sa=0,sb=0,s,t3; double hMAX_N,vMAX_N,uMAX_N,MMAX_N+1,dMAX_N; /*courant中所用变量*/ int n2; int i2,j2,k2; int mi2; double mx2,tmp2; static double a2MA

11、X_NMAX_N,b2MAX_N,x2MAX_N,y2MAX_N; static double l2MAX_NMAX_N,u2MAX_NMAX_N; printf(nInput n value:); scanf(%d,&n); if(nMAX_N) printf(The input n is larger then MAX_N,please redefine the MAX_N.n); return 1; if(n=0) printf(Please input a number between 1 and %dn,MAX_N); return 1; printf(Now input the (

12、x_i,y_i),i=0,.,%d:n,n); for(i=0;i=n;i+) scanf(%lf%lf,&pointsi.x,&pointsi.y); printf(Now input the x value:); scanf(%lf,&x); for(i=0;in;i+) hi=pointsi+1.x-pointsi.x; for(i=1;in;i+) vi=hi/(hi+hi-1); ui=1-vi; for(i=1;in;i+) tmp=0; tmp+=pointsi-1.y /(pointsi-1.x-pointsi.x)*(pointsi-1.x-pointsi+1.x); tmp

13、+=pointsi.y /(pointsi.x-pointsi-1.x)*(pointsi.x-pointsi+1.x); tmp+=pointsi+1.y /(pointsi+1.x-pointsi.x)*(pointsi+1.x-pointsi-1.x); di=tmp*6; /*线性方程组求解 解放在Mi中 , M0=Mn=0 */ n2=n-1; M0=0; Mn=0; for(i2=0;i2n2;i2+) for(j2=0;j2n2;j2+) a2i2j2=0; for(i2=0;i2=0) a2i2i2-1=ui2+1; a2i2i2=2; if(i2+1n2) a2i2i2+1=

14、vi2+1; for(i2=0;i2n2;i2+) b2i2=di2+1; /*courant*/for(i2=0;i2n2;i2+)u2i2i2=1; for(k2=0;k2n2;k2+) for(i2=k2;i2n2;i2+) l2i2k2=a2i2k2; for(j2=0;j2=k2-1;j2+) l2i2k2-=(l2i2j2*u2j2k2); for(j2=k2+1;j2n2;j2+) u2k2j2=a2k2j2; for(i2=0;i2=k2-1;i2+) u2k2j2-=(l2k2i2*u2i2j2); u2k2j2/=l2k2k2; for(i2=0;i2n2;i2+) y2i2=b2i2; for(j2=0;j2=0;i2-) x2i2=y2i2; for(j2=i2+1;j2n2;j2+) x2i2-=(u2i2j2*x2j2); for(i2=0;i2=pointsi.x&i0)i-; sa=pointsi+1.x-x; sb=x-pointsi.x ; s=Mi*sa*sa*sa/(hi*6)+Mi+1*sb*sb*sb/(hi*6)+ (pointsi.y-Mi*hi*hi/6)*sa/hi+ (pointsi+1.y-Mi+1*hi*hi/6)*sb/hi ; printf(S(%f)=%fn,x,s); return 0;