东南大学计算方法实验报告

计算方法与实习实验报告学院：电气工程学院指导老师：李翠平班级：160093姓名：黄芃菲学号：16009330实习题一实验1拉格朗日插值法一、方法原理n次拉格朗日插值多项式为：Ln(x)=y0l0(x)+y1l1(x)+y2l2(x)+…+ynln(x)n=1时，称为线性插值，L1(x)=y0(x-x1)/(x0-x1)+y1(x-x0)/(x1-x0)=y0+(y1-x0)(x-x0)/(x1-x0)n=2时，称为二次插值或抛物线插值，精度相对高些L2(x)=y0(x-x1)(x-x2)/(x0-x1)/(x0-x2)+y1(x-x0)(x-x2)/(x1-x0)/(x1-x2)+y2(x-x0)(x-x1)/(x2-x0)/(x2-x1)二、主要思路使用线性方程组求系数构造插值公式相对复杂，可改用构造方法来插值。对节点xi(i=0,1,…,n)中任一点xk(0=k=n)作一n次多项式lk(xk)，使它在该点上取值为1，而在其余点xi(i=0,1,…,k-1,k+1,…,n)上为0，则插值多项式为Ln(x)=y0l0(x)+y1l1(x)+y2l2(x)+…+ynln(x)上式表明：n个点xi(i=0,1,…,k-1,k+1,…,n)都是lk(x)的零点。可求得lk三．计算方法及过程：1.输入节点的个数n2.输入各个节点的横纵坐标3.输入插值点4.调用函数，返回z函数语句与形参说明程序源代码如下：#includeiostream#includemath.husingnamespacestd;#defineN100doublefun(double*x,double*y,intn,doublep);voidmain(){inti,n;cout输入节点的个数n：;cinn;doublex[N],y[N],p;coutpleaseinputxiangliangx=endl;for(i=0;in;i++)cinx[i];coutpleaseinputxiangliangy=endl;for(i=0;in;i++)ciny[i];coutpleaseinputLagelangrichazhiJieDianp=endl;cinp;形参与函数类型参数意义intn节点的个数doublex[n]（double*x）存放n个节点的值doubley[n]（double*y）存放n个节点相对应的函数值doublep指定插值点的值doublefun()函数返回一个双精度实型函数值，即插值点p处的近似函数值coutTheAnswer=fun(x,y,n,p)endl;system(pause);}doublefun(doublex[],doubley[],intn,doublep){doublez=0,s=1.0;intk=0,i=0;doubleL[N];while(kn){if(k==0){for(i=1;in;i++)s=s*(p-x[i])/(x[0]-x[i]);L[0]=s*y[0];k=k+1;}else{s=1.0;for(i=0;i=k-1;i++)s=s*((p-x[i])/(x[k]-x[i]));for(i=k+1;in;i++)s=s*((p-x[i])/(x[k]-x[i]));L[k]=s*y[k];k++;}}for(i=0;in;i++)z=z+L[i];returnz;}五．实验分析n=2时，为一次插值，即线性插值n=3时，为二次插值，即抛物线插值n=1，此时只有一个节点，插值点的值就是该节点的函数值n1时，结果都是返回0的；这里做了n=0和n=-7两种情况3n100时，也都有相应的答案常用的是线性插值和抛物线插值，显然，抛物线精度相对高些n次插值多项式Ln（x）通常是次数为n的多项式，特殊情况可能次数小于n.例如：通过三点的二次插值多项式L2(x),如果三点共线，则y=L2(x)就是一条直线，而不是抛物线，这时L2(x)是一次式。拟合曲线光顺性差实验2牛顿插值法一、方法原理及基本思路在拉格朗日插值方法中，若增加一个节点数据，其插值的多项式需重新计算。现构造一个插值多项式Nn(x)，只需对Nn-1(x)作简单修正(如增加某项)即可得到，这样计算方便。利用牛顿插值公式，当增加一个节点时，只需在后面多计算一项，而前面的计算仍有用；另一方面Nn(x)的各项系数恰好又是各阶差商，而各阶差商可用差商公式来计算。由线性代数知，对任何一个不高n次的多项式P(x)=b0＋b1x＋b2x2＋…＋bnxn(幂基)①也可将其写成P(x)=a0＋a1(x－x0)＋a2(x－x0)(x－x1)＋…＋an(x－x0)…(x－xn-1)其中ai为系数，xi为给定节点，可由①求出ai一般情况下，牛顿插值多项式Nn(x)可写成：Nn(x)=a0＋a1(x－x0)＋a2(x－x0)(x－x1)＋…＋an(x－x0)…(x－xn-1))只需求出系数ai，即可得到插值多项式。二、计算方法及过程1.先后输入节点个数n和节点的横纵坐标，插值点的横坐标，最后输入精度e2.用do-while循环语句得到跳出循环时k的值3.将k值与n-1进行比较，若在达到精度时kn-1,则输出计算结果；若此时k=n-1,则计算失败！函数语句与形参说明程序源代码如下：#includeiostream#includemath.husingnamespacestd;#defineMAX100voidmain(){floatx[MAX],y[MAX];floatx0,y0,e,N1;floatN0=0;inti,k,n;cout输入节点的个数,n=;cinn;cout请输入插值节点！endl;for(i=0;in;i++)cinx[i];cout请输入对应的函数值！endl;for(i=0;in;i++)ciny[i];cout请输入插值点,x0=;cinx0;cout请输入精度,e=;cine;y0=1;N1=y[0];k=0;do{k++;N0=N1;y0=y0*(x0-x[k-1]);for(i=0;ik;i++)y[k]=(y[k]-y[i])/(x[k]-x[i]);N1=N0+y0*y[k];}while(fabs(N1-N0)e&&k(n-1));形参与函数类型参数意义intn节点的个数floatx[MAX]存放n个节点的值（从小到大）Floaty[MAX];存放n个节点相对应的函数值floatx0,y0指定插值点的横纵坐标floate精度if(k==(n-1))cout计算失败！;if(k(n-1))cout输出结果y[x0]=N1endl;system(pause);}三．运行结果测试：题一：已知f(x)=sh(x)的函数表如下：计算f(0.23)的近似值xi00.200.300.50Sh(xi)00.201340.304520.52110实习题二1.用牛顿法求下列方程的根：3）lg20xx实验代码：#includestdio.h#includemath.h#defineN100#defineeps1e-6#defineeta1e-8floatNewton(float(*f)(float),float(*f1)(float),floatx0){floatx1,d;intk=0;do{x1=x0-(*f)(x0)/(*f1)(x0);if(k++N||fabs((*f1)(x1))eps){printf(\nNewton迭代发散);break;}d=fabs(x1)1?x1-x0:(x1-x0)/x1;x0=x1;printf(x(%d)=%f\t,k,x0);}while(fabs(d)eps&&fabs((*f)(x1))eta);returnx1;}floatf(floatx){returnlog10(x)+x-2;}floatfl(floatx){return1.0/(x*log(10))+1;}voidmain(){floatx0,y0;printf(请输入迭代初值x0\n);scanf(%f,&x0);printf(x(0)=%f\n,x0);y0=Newton(f,fl,x0);printf(方程的根为:%f\n,y0);}运行窗口：实习题三4．编写用追赶法解三对角线性方程组的程序，并解下列方程组：2）Axb实验代码#includeiostream#includecmathusingnamespacestd;voidmain(){floata=1;floatb=-4;floatc=1;floatd[11]={0,-27,-15,-15,-15,-15,-15,-15,-15,-15,-15};floatl[11];floatbb[11];floaty[11];floatx[11];bb[1]=b;y[1]=d[1];inti;for(i=2;i11;i++){l[i]=a/bb[i-1];bb[i]=b-l[i]*c;y[i]=d[i]-l[i]*y[i];}x[10]=y[10]/bb[10];for(i=9;10;i--){x[i]=(y[i]-c*x[i+1])/bb[i];}for(i=1;i11;i++)cout'x'i':'x[i]endl;}运行窗口#includestdio.h#includestring.h#includemath.h#includeconio.h#includestdlib.h#defineN11main(){floata[N]={0,0,1,1,1,1,1,1,1,1,1};floatb[N]={0,-4,-4,-4,-4,-4,-4,-4,-4,-4,-4};floatc[N]={0,1,1,1,1,1,1,1,1,1,0};floatd[N]={0,-27,-15,-15,-15,-15,-15,-15,-15,-15,-15};floatx[N]={0,0,0,0,0,0,0,0,0,0,0};floatr[N]={0,0,0,0,0,0,0,0,0,0,0};floaty[N]={0,0,0,0,0,0,0,0,0,0,0};floatq;//clrscr();intk;r[1]=c[1]/b[1];y[1]=d[1]/b[1];for(k=2;kN-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;kN;k++)printf(x[%d]=%f\n,k,x[k]);getch();return0;}运行窗口5．12310211xxx2348311xxx1232106xxx123431125xxxx雅克比迭代法实验代码#includestdio.h#includemath.h#defineeps1e-6#definemax100voidJacobi(float*a,intn,floatx[]){inti,j,k=0;floatepsilon,s;float*y=newfloat[n];for(i=0;in;i++)x[i]=0;while(1){epsilon=0;k++;for(i=0;in;i++){s=0;for(j=0;jn;j++){if(j==i)continue;s+=*(a+i*(n+1)+j)*x[j];}y[i]=(*(a+i*(n+1)+n)-s)/(*(a+i*(n+