中北大学 数值分析1

中北大学理学院实验报告实验类别：数值分析专业：信息与计算科学班 级：13080241学 号：1308024120姓名：杨燕中北大学理学院实验一函数插值方法实验内容】给定一元函数y二f(x)的n+1个节点值y二f(x)(j二0,1,n)，数据如下:jjxj0.40.550.650.800.951.05yj0.410750.578150.696750.901.001.25382求五次Lagrange多项式L5(x)或分段三次插值多项式或Newton插值多项式，并计f(0.596),f(0.99)的值。(提示：结果为f(0.596)~0.625732,f(0.99)~1.05423)实验方法与步骤】利用Lagrange插值公式x利用Lagrange插值公式x-x

i

x-xI=0k,i丰kyk，用C语言编写出插值多项式程序如下：#include<stdio.h>#defineN5floatx[]={0.4,0.55,0.65,0.80,0.95,1.05};floaty[]={0.41075,0.57815,0.69675,0.90,1.00,1.25382};floatp(floatxx){inti,k;floatpp=0,m1,m2;for(i=0;i<=N;i++){m1=1;m2=1;for(k=0;k<=N;k++)if(k!=i){m1*=xx-x[k];m2*=x[i]-x[k];}pp+=y[i]*m1/m2;}returnpp;}main(){printf("f(0.596)=%lf\n",p(0.596));printf("f(0.99)=%lf\n",p(0.99));}【实验结果】【思考】1、给出的程序求f(1.06)行不行，精度高不高？2、五次Lagrange多项式与Newton插值多项式是同一个多项式吗？五次Lagrange多项式与Newton插值多项式是同一个多项式。3、为什么高次插值不能令人满意？一般来说，节点个数越多，插值函数和被插值函数就有越多的地方相等。但是随着插值节点个数的增加，两个插值节点之间插值函数并不一定能够很好地逼近被插值函数。再次，从舍入误差看，高次插值由于计算量大，可能会产生更严重的误差积累，所以，稳定性得不到保证。此时就会出现龙格现象。中北大学理学院实验报告实验类别：数值分析专业：信息与计算科学班 级：13080241学 号：1308024120姓名：杨燕中北大学理学院实验二函数逼近与曲线拟合实验内容】1、编写出Legendre、Chebyshev多项式的程序；2、从随机的数据中找出其规律性，给出其近似表达式，在生产实践和科学实验中大量存在，通常利用数据的最小二乘法求得拟合曲线。例如在某冶炼过程中根据统计数据的含碳量与时间关系，试求含碳量y与时间t的拟合曲线。（t分）0510152025303540455055yC10一4)01.72.162.863.443.874.154.374.514.584.024.64i111111111111实验方法与步骤】1、用C编写语言出Legendre多项式的程序如下:#include<stdio.h>doublep(intn,doublex){if(n==0)return1;elseif(n==1)returnx;elsereturn((2*n-1)*x*p(n-1,x)-(n-1)*p(n-2,x))/n;}intmain(){intn;doublex;doubley;printf("inputn,x:\n");scanf("%d%lf",&n,&x);y=p(n,x);printf("%lf\n",y);return0;}2、给出表格数据的近似解析表达式为；申（t）=a+bt，选取基函数申o（x）二1,申1（x）=x，则得到(p,p)=12，(甲,*)=(*,*)=丈x,(甲,甲)=艺x2，(cp,f)=艺f，G,f)=£fx‘001001i11i0i1iii=0 i=0 i=0 i=0于是得方程组険P0)a+%P1)b=%f),用C语言编写曲线拟合的程序如I(p,p丿a+(p,p)b=(p,f)10111下：#include<stdio.h>#defineMax_N25main(){inti,n;doublex[Max_N],y[Max_N];doublea11,a12,a21,a22,d1,d2;doublea,b;printf("\nInputnvalue:");do{scanf("%d",&n);if(n>Max_N)printf("\npleasere_inputnvalue:");}while(n>Max_N||n<=0);printf("inputx[i],i=0,...%d:\n",n-1);for(i=0;i<n;i++)scanf("%lf",&x[i]);printf("inputy[i],i=0,...%d:\n",n-1);for(i=0;i<n;i++)scanf("%lf",&y[i]);for(i=0;i<n;i++){a21+=x[i];a22+=x[i]*x[i];d1+=y[i];d2+=x[i]*y[i];}a12=a21;a11=n;a=(d1*a22-d2*a12)/(a11*a22-a12*a21);b=(d1*a21-d2*a11)/(a21*a12-a22*a11);printf("slove:P(x)=%f+%fx\n",a,b);getchar();return0;}3、为作比较，用MATLAB进行曲线拟合，编写的程序如下:x=[0510152025303540455055];y=[01.272.162.863.443.874.154.374.514.584.024.64];p=polyfit(x,y,1);disp([num2str(p(1)),'*x+',num2str(p(2))]);xx=linspace(0,60,60);yy=polyval(p,xx);plot(x,y,'rx',xx,yy)【实验结果】1、C语言编写的Legendre多项式的程序结果如下inputn,x：43321.000000Pressanykeytocontinue2、C语言拟合的曲线结果如下：Inputnualue-12inputxti1,1=0,...11：0ii1315翎2S303540455055inputyti.],1=0,...11：Q1.272.162.863.443.874.154.374.514.584.924.64Sloue：P<x>=l.321539+0.672762x3、MATLAB拟合的曲线结果和误差分析结果如下:拟合的方程为：y=0.072762x+1.3215【思考】1、 最佳一致逼近与最佳平方逼近的区别是什么？最佳一致逼近中的范数取的是无穷大范数，公式为||f(x)-P*(x) =minmaxf(x)-P(x)。gP^Hna<x<b最佳平方逼近中的范数取的是二范数，公式为||f(x)—P*(x)|2=minJb[f(x)—P(x)1dx。2PwHan2、 曲线拟合与最佳平方逼近的区别是什么？曲线拟合中的f(x)是[a,b]上的一个列表函数，公式为：||f-P*||2=min£[f(x)-P(x)1。2唄i=o i i最佳平方逼近中的f(x)是一个未知函数，公式为：f(x-P*x(2=) Jmfixn-Px(I2dx。()2PeHan3、能用什么方法确定最小二乘法的拟合函数？拟合的方法除了最小二乘法外，还有拉格朗日插值法、牛顿插值法、区间二分法、雅克比迭代法和牛顿科特斯数值积分法等，都可以确定拟合函数。中北大学理学院实验报告实验类别：数值分析专业：信息与计算科学班 级：13080241学 号：1308024120姓名：杨燕中北大学理学院实验三数值积分与数值微分实验内容】选用复合梯形公式，复合Simpson公式，Romberg算法、高斯算法计算(1)I二宀小-sin2xdx(2)I=J1 ―dx(3)I=11旦匕卫dXo o4+x2 o1+x2实验方法与步骤】1、用C语言编写Romberg算法数值积分的程序如下：#include<stdio.h>#include<math.h>#defineepsilono.ooo1/*求基函数*/floatf(floatx){return(sqrt(4-sin(x)*sin(x)));} /*梯形公式*/floatRomberg(floata,floatb){intk=1;floatS,x,T1,T2,S1,S2,C1,C2,R1,R2,h=b-a;T1=h/2*(f(a)+f(b));while(1){S=o;x=a+h/2;do{S+=f(x);x+=h;}while(x<b);T2=(T1+h*S)/2.o;if(fabs(T2-T1)<epsilon)return(T2);S2=T2+(T2-T1)/3.0;if(k==1){T1=T2;S1=S2;h/=2;k+=1;continue;}if(fabs(S2-S1)<epsilon)return(S2);C2=S2+(S2-S1)/15.0;if(k==2){C1=C2;T1=T2;S1=S2;h/=2;k+=1;continue;}if(fabs(C2-C1)<epsilon)return(C2);R2=C2+(C2-C1)/63.0;if(k==3){R1=R2;C1=C2;T1=T2;S1=S2;h/=2;k+=1;continue;}if(fabs(R2-R1)<epsilon)return(R2);R1=R2;C1=C2;T1=T2;S1=S2;h/=2;k+=1;}}main(){inti;floata,b,S;printf("\nInputbeginandend:");scanf("%f%f",&a,&b);S=Romberg(a,b);printf("Solveis:%f",S);scanf("%",S);}2、再用复合Simpson公式计算同一个积分，并比较其结果，且用C语言编写的程序如下：#include<stdio.h>#include<math.h>floatf(floatx){return(sqrt(4-sin(x)*sin(x)));}floatSimpson(floata,floatb){intk=1,n=20;floats,x1,x2,h=(b-a)/n,c=a+h/2;s=h/6*(f(a)+f(b)+4*f(c));for(k=1;k<=n-1;k++){x1=a+k*h;x2=x1+h/2;s=s+(2*h/3)*f(x2)+(h/3)*f(x1);

return(s);main(){inti;floata,b,s;printf(“请输入a,b:");scanf("%f%f",&a,&b);s=Simpson(a,b);printf("解为:%f",s);scanf("%",s);}3、分别取不同步长h二(b-a)/n，试比较计算结果(如n二10,20等);4、给定精度要求e，试用变步长算法，确定最佳步长。实验结果】1、所给三个定积分的结果如下：Inputbeginandend：00.25Solueis=0.498711Inputbeginandend：01Bolueis：0.2721972、对第一个定积分I=I‘4J4-sin2xdx用Romberg算法和复合Simpson公式计0算的结果分别如下：=00-25解为：0・498711Inputbeginandend=00-25解为：0・4987113、步长为h二(b-a)/10时，所给三个定积分的结果如下:Inputbeginandend：00.25Solueis：0.498708Inputbeginandend：01Solveis：0.2721784、步长为h二(b-a)/20时，所给三个定积分的结果如下:Inputbeginandend：Inputbeginandend：00.25Solueis=0.498520Inputbeginandend：01olueis：0.390654Inputbeginandend：01Solueis=0.272082【思考】1、 复合求积公式的优点是什么？复合求积公式通过把区间分成若干个子区间，再在每个子区间上用低阶求积公式来提高求积的精度。2、 Romberg算法与与牛顿-柯特斯求积算法的联系是什么？Romberg算法与与牛顿-柯特斯求积算法的节点都是等距选取均匀分布的，且n>8的牛顿-柯特斯公式因为误差太大而不能使用。3、 高斯算法的求积节点如何确定？以这些节点为零的多项式①(x)二(x-x)(x-x)…(x-x)与任何次数不超过TOC\o"1-5"\h\zn+1 0 1 nn的多项式p(x)带权p(x)正交，即：Ibp(x)e(x)p(x)dx=0。a n+14、 牛顿-柯特斯求积与高斯算法的节点分布有什么不同？牛顿-柯特斯求积的节点是等距选取均匀分布的，而高斯算法的节点是满足以这些节点为零的多项式3(x)二(x-x)(x-x)•••(x-x)与任何次数不超过n的多n+1 0 1 n项式p(x)带权p(x)正交，即是不均匀分布的。中北大学理学院实验报告实验类别：数值分析专业：信息与计算科学班 级：13080241学 号：1308024120姓名：杨燕中北大学理学院

实验四线性方程组的直接解法实验内容】2x一x=312用追赶法）一x+3x一2x=1用追赶法）1 2 3一2x+4x一2x=0234一2x+5x=一534实验方法与步骤】1、对上述方程组用追赶法求解，用C语言编程如下：#include<stdio.h>#include<math.h>#include<string.h>#defineN5main(){floata[N]={0,0,-1,-2,-2};floatb[N]={0,2,3,4,5};floatc[N]={0,-1,-2,-2,0};floatd[N]={0,3,1,0,-5};floatx[N]={0,0,0,0};floatr[N]={0,0,0,0};floaty[N]={0,0,0,0};floatq;intk;r[1]=c[1]/b[1];y[1]=d[1]/b[1];for(k=2;k<=N-1;k++){q=b[k]-r[k-1]*a[k];r[k]=c[k]/q;y[k]=(d[k]-y[k-1]*a[k])/q;}y[N-1]=(d[N-1]-y[N-2]*a[N-1])/(b[N-1]-r[N-2]*a[N-1]);x[N-1]=y[N-1];for(k=N-2;k>=1;k--)x[k]=y[k]-r[k]*x[k+1];for(k=1;k<N;k++)printf("x[%d]=%lf\n",