数值分析报告Lagrange差值和牛顿插值

实验一一、实验名称Lagrange插值多项式和牛顿插值多项式二、实验目的与要求：实验目的：掌握Lagrange插值多项式和牛顿插值多项式的算法。实验要求：1.给出Lagrange插值和牛顿插值算法思路，2.用C语言实现算法，运行环境为MicrosoftVisualC++，3.计算误差（这里只要求给出（-5,5）内101个点的误差）。三、实验内容：1.对Lagrange插值多项式算法作编程练习和上机运算，2.对牛顿插值多项式算法作编程练习和上机运算，3.比较两种方法。算法思路：1.Lagrange算法是把多项式p写成如下形式：0011()()()+()nnpxylxylxylx……，其中0()njijijjixxlxxx称为Lagrange基函数。计算Lagrange基函数的方法：fx=0.0;for(i=0;i=n;i++){tmp=1.0;for(j=0;ji;j++)tmp=tmp*(x-x[j])/(x[i]-x[j]);for(j=i+1;jn;j++)tmp=tmp*(x-x[j])/(x[i]-x[j]);fx=fx+tmp*y[i];}return(fx);2.牛顿算法是把多项式p写成如下形式：010101-1()()+()()()nnpxccxxcxxxxxx……其中01,,nxxx…,是插值点，011,,nccc…,是待定系数。可以通过插值点01,,nxxx…,和插值点处的函数值01,,nyy…,y算出待定系数011,,nccc…,，方法如下：(1)01010k-1()()+()()kkkkkkpxccxxcxxxx……(2)101020k-2()()+()()kkkkkkpxccxxcxxxx……将(1)-(2)并利用()kkkypx，得1011()()()()()kkkkkkkkkpxpxcxxxxxx…1011()()()()kkkkkkkypxxxxxxx…算法：c[0]=y1[0];for(k=1;k=N;k++){d=x1[k]-x1[k-1];u=c[k-1];for(i=k-2;i=0;i--){u=u*(x1[k]-x1[i])+c[i];d=d*(x1[k]-x1[i]);}c[k]=(y1[k]-u)/d;}3.这里，我们先编写两个函数px1和px2，分别用来计算插值节点取105ixiN和Chebyshevpoint215cos()22iixN，i=0,1,…，N时的Lagrange多项式函数或牛顿多项式函数，然后在main函数中，主要目的是计算误差|f(x)-p(x)|在（-5,5）内101个点处的最大值，其中当然要引用前面编写的两个函数px1和px2，最后将N=5,10,20,40时两种情况下的最大误差输出，并分析结果。四、实验题目：对函数21(),[5,5]1fxxx，构造Lagrange插值多项式和牛顿插值多项式，插值节点取为105ixiN，i=0,1,…，N和215cos()22iixN，i=0,1,…，N（Chebyshevpoint），并计算如下误差0100max{|()()|},5,10iiiiifypyy对N=5,10,20,40比较以上两组节点的结果。编写程序（程序见后面附录），输出结果如下：1.Lagrange插值法结果：可以看出：（1）对于该函数而言，N值越大，即插值点数越多，相对误差（大概地）越小，但是要求得严格的误差，应该取|f(x)-p(x)|在区间（-5,5）上的最大值。这里认为该区间上101个点足够密集，而且f(x)-p(x)在（-5,5）上是一致连续的，所以将其看作所求误差是合理的。（2）MaxErrorofgrid（1）总是小于MaxErrorofgrid（2），而且MaxErrorofgrid（1）迅速收敛到零，然而MaxErrorofgrid（2）则很缓慢地向零靠拢，这说明该函数f（x）用（-5,5）上平均分布的N个点得到的Lagrange函数比用（-5,5）上N个Chebyshevpoint更靠近原函数f(x)。2.牛顿插值法得到的结果：可以看出：（1）单独地观察这组数据，会得到和Lagrange插值法一样的结论，对于该函数而言，N值越大，即插值点数越多，相对误差（大概地）越小，MaxErrorofgrid（1）总是小于MaxErrorofgrid（2），且MaxErrorofgrid（1）更迅速地收敛到零。（2）和Lagrange插值法数据比较，我们会发现当取相同的插值点时，两种方法得到的误差结果是一样的。追究其原因知，这是多项式插值定理决定两组数据必定是相同的，多项式插值定理是这样阐述的，若01,,nxxx…,是不同的实数，则对任意数值01,,nyy…,y，存在唯一的次数至多是n次的多项式P（x），使得(),(0)niipxyin，所以两种方法算得的Lagrange多项式和牛顿多项式实质上是同一个函数P（x），故0100max{|()()|},510iiiiifypyy当然是相同的了。五、附录：实验编程，运行环境为MicrosoftVisualC++A．Lagrange插值多项式法：#includestdio.h#includemath.hdoublepx1(floatw,intN)/*构造变量为w的Lagrange函数，xi=5-10i/N为插值点*/{inti,j;doublez1=0.0,x1[50],y1[50];for(i=0;i=N;i++){x1[i]=5.0-10.0*i/N;y1[i]=1/(1+x1[i]*x1[i]);}for(i=0;i=N;i++){doubletmp=1.0;for(j=0;ji;j++)tmp=tmp*(w-x1[j])/(x1[i]-x1[j]);for(j=i+1;j=N;j++)tmp=tmp*(w-x1[j])/(x1[i]-x1[j]);z1=z1+tmp*y1[i];}return(z1);}doublepx2(floatw,intN)/*构造变量为w的Lagrange函数，Chebyshevpoint*/{inti,j;doublez2=0.0,x2[50],y2[50];for(i=0;i=N;i++){x2[i]=-5.0*cos((2*i+1)*3.1415926/(2*N+2));y2[i]=1/(1+x2[i]*x2[i]);}for(i=0;i=N;i++){doubletmp=1.0;for(j=0;ji;j++)tmp=tmp*(w-x2[j])/(x2[i]-x2[j]);for(j=i+1;j=N;j++)tmp=tmp*(w-x2[j])/(x2[i]-x2[j]);z2=z2+tmp*y2[i];}return(z2);}voidmain()/*计算-5至5中101个点上的最大误差fmax*/{intk,N;doublem1[101],m2[101],z[101],fmax1=0.0,fmax2=0.0;/*z[101]为101个点，m[101]为相应差值*/printf(pleseinputthenumberN:);scanf(%d,&N);printf(N=%d\n,N);for(k=0;k=100;k++){z[k]=k/10-5;m1[k]=fabs(1/(1+z[k]*z[k])-px1(z[k],N));m2[k]=fabs(1/(1+z[k]*z[k])-px2(z[k],N));}fmax1=m1[0];for(k=1;k=100;k++){if(m1[k]m1[k-1])fmax1=m1[k];}printf(MaxErrorofgrid(1)：%6.12f\n,fmax1);fmax2=m2[0];for(k=1;k=100;k++){if(m2[k]m2[k-1])fmax2=m2[k];}printf(MaxErrorofgrid(2)：%6.12f\n,fmax2);}B．牛顿插值多项式法：#includestdio.h#includemath.hdoublepx1(floatw,intN)/*构造变量为w的牛顿插值多项式函数，xi=5-10i/N为插值点*/{inti,k;doubled,u,z1,tmp,x1[50],y1[50],c[50];for(i=0;i=N;i++){x1[i]=5.0-10.0*i/N;y1[i]=1/(1+x1[i]*x1[i]);}c[0]=y1[0];for(k=1;k=N;k++){d=x1[k]-x1[k-1];u=c[k-1];for(i=k-2;i=0;i--){u=u*(x1[k]-x1[i])+c[i];d=d*(x1[k]-x1[i]);}c[k]=(y1[k]-u)/d;}z1=c[0];tmp=1.0;for(k=1;k=N;k++){tmp=tmp*(w-x1[k-1]);z1=z1+c[k]*tmp;}return(z1);}doublepx2(floatw,intN)/*构造变量为w的牛顿函数，Chebyshevpoint*/{inti,k;doubled,u,z2,tmp,x2[50],y2[50],c[50];for(i=0;i=N;i++){x2[i]=-5.0*cos((2*i+1)*3.1415926/(2*N+2));y2[i]=1/(1+x2[i]*x2[i]);}c[0]=y2[0];for(k=1;k=N;k++){d=x2[k]-x2[k-1];u=c[k-1];for(i=k-2;i=0;i--){u=u*(x2[k]-x2[i])+c[i];d=d*(x2[k]-x2[i]);}c[k]=(y2[k]-u)/d;}z2=c[0];tmp=1.0;for(k=1;k=N;k++){tmp=tmp*(w-x2[k-1]);z2=z2+c[k]*tmp;}return(z2);}voidmain()/*计算-5至5中101个点上的最大误差fmax*/{intk,N;doublem1[101],m2[101],z[101],fmax1=0.0,fmax2=0.0;/*z[101]为101个点，m[101]为相应差值*/printf(pleseinputthenumberN:);scanf(%d,&N);printf(牛顿插值法：\nN=%d\n,N);for(k=0;k=100;k++){z[k]=k/10-5;m1[k]=fabs(1/(1+z[k]*z[k])-px1(z[k],N));m2[k]=fabs(1/(1+z[k]*z[k])-px2(z[k],N));}fmax1=m1[0];for(k=1;k=100;k++){if(m1[k]m1[k-1])fmax1=m1[k];}printf(MaxErrorofgrid(1)：%6.12f\n,fmax1);fmax2=m2[0];for(k=1;k=100;k++){if(m2[k]m2[k-1])fmax2=m2[k];}printf(MaxErrorofgrid(2)：%6.12f\n,fmax2);}