东南大学计算方法实验报告(共16页).doc

《东南大学计算方法实验报告(共16页).doc》由会员分享，可在线阅读，更多相关《东南大学计算方法实验报告(共16页).doc（16页珍藏版）》请在淘文阁 - 分享文档赚钱的网站上搜索。

1、精选优质文档-倾情为你奉上 计算方法与实习实验报告学 院：电气工程学院指导老师：李翠平班 级：姓 名：黄芃菲学 号：实习题一实验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

2、(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函数语句与形参说明程序源代码

3、如下：形参与函数类型参数意义int n节点的个数double xn（double *x）存放n个节点的值double yn（double *y）存放n个节点相对应的函数值double p指定插值点的值double fun()函数返回一个双精度实型函数值，即插值点p处的近似函数值#include#includeusing namespace std;#define N 100double fun(double *x,double *y, int n,double p);void main()int i,n;coutn;double xN, yN,p;coutplease input xiangl

4、iang x= endl;for(i=0;ixi;coutplease input xiangliang y= endl;for(i=0;iyi;coutplease input LagelangrichazhiJieDian p= p;coutThe Answer= fun(x,y,n,p)endl;system(pause) ;double fun(double x,double y, int n,double p)double z=0,s=1.0;int k=0,i=0;double LN;while(kn) if(k=0) for(i=1;in;i+)s=s*(p-xi)/(x0-xi

5、); L0=s*y0; k=k+1; else s=1.0; for(i=0;i=k-1;i+)s=s*(p-xi)/(xk-xi); for(i=k+1;in;i+) s=s*(p-xi)/(xk-xi); Lk=s*yk;k+;for(i=0;in;i+)z=z+Li;return z;五实验分析n=2时，为一次插值，即线性插值n=3时，为二次插值，即抛物线插值n=1，此时只有一个节点，插值点的值就是该节点的函数值n1时，结果都是返回0的；这里做了n=0和n=-7两种情况3n100时，也都有相应的答案常用的是线性插值和抛物线插值，显然，抛物线精度相对高些n次插值多项式Ln（x）通常是次数为

6、n的多项式，特殊情况可能次数小于n.例如：通过三点的二次插值多项式L2(x),如果三点共线，则y=L2(x)就是一条直线，而不是抛物线，这时L2(x)是一次式。拟合曲线光顺性差实验2 牛顿插值法一、方法原理及基本思路在拉格朗日插值方法中，若增加一个节点数据，其插值的多项式需重新计算。现构造一个插值多项式Nn(x)，只需对Nn-1(x)作简单修正(如增加某项)即可得到，这样计算方便。利用牛顿插值公式，当增加一个节点时，只需在后面多计算一项，而前面的计算仍有用；另一方面Nn(x)的各项系数恰好又是各阶差商，而各阶差商可用差商公式来计算。由线性代数知，对任何一个不高n次的多项式P(x)=b0b1xb

7、2x2bnxn (幂基) 也可将其写成P(x)=a0a1(xx0)a2(xx0) (xx1)an(xx0) (xxn-1)其中ai为系数，xi为给定节点，可由求出ai 一般情况下，牛顿插值多项式Nn(x)可写成：Nn(x)= a0a1(xx0)a2(xx0) (xx1)an(xx0) (xxn-1)只需求出系数ai，即可得到插值多项式。二、计算方法及过程1.先后输入节点个数n和节点的横纵坐标，插值点的横坐标，最后输入精度e2. 用do-while循环语句得到跳出循环时k的值3.将k值与n-1进行比较，若在达到精度时kn-1,则输出计算结果；若此时k=n-1,则计算失败！函数语句与形参说明形参与

8、函数类型参数意义int n节点的个数float xMAX存放n个节点的值（从小到大）Float yMAX;存放n个节点相对应的函数值float x0,y0指定插值点的横纵坐标float e精度程序源代码如下：#include#includeusing namespace std;#define MAX 100void main() float xMAX,yMAX; float x0,y0,e,N1; float N0=0; int i,k,n; coutn; cout请输入插值节点！endl; for(i=0;ixi; cout请输入对应的函数值！endl; for(i=0;iyi; cout

9、x0; coute; y0=1; N1=y0; k=0; do k+; N0=N1; y0=y0*(x0-xk-1); for(i=0;i e & k(n-1); if (k=(n-1) cout计算失败！; if (k(n-1) cout输出结果yx0=N1endl;system(pause);三运行结果测试：题一：已知f(x)=sh(x)的函数表如下：计算f(0.23)的近似值xi00.200.300.50Sh(xi)00.201340.304520.52110实习题二1. 用牛顿法求下列方程的根：3）实验代码：#include #include #define N 100#define

10、eps 1e-6#define eta 1e-8float Newton(float(*f)(float),float(*f1)(float),float x0)float x1,d;int k=0;dox1=x0-(*f)(x0)/(*f1)(x0);if(k+N|fabs(*f1)(x1)eps)printf(n Newton 迭代发散);break;d=fabs(x1)eps&fabs(*f)(x1)eta);return x1;float f(float x)return log10(x)+x-2;float fl(float x)return 1.0/(x*log(10)+1;voi

11、d main()float x0,y0;printf(请输入迭代初值 x0n);scanf(%f,&x0);printf(x(0)=%fn,x0);y0=Newton(f,fl,x0);printf(方程的根为: %fn,y0);运行窗口：实习题三4编写用追赶法解三对角线性方程组的程序，并解下列方程组：2）l 实验代码#include#includeusing namespace std;void main()float a=1;float b=-4;float c=1;float d11=0,-27,-15,-15,-15,-15,-15,-15,-15,-15,-15;float l11;

12、float bb11;float y11;float x11;bb1=b;y1=d1;int i;for(i=2;i0;i-)xi=(yi-c*xi+1)/bbi;for(i=1;i11;i+)coutxi:xiendl;l 运行窗口#include#include#include#include#include#define N 11main()float aN=0,0,1,1,1,1,1,1,1,1,1;float bN=0,-4,-4,-4,-4,-4,-4,-4,-4,-4,-4;float cN=0,1,1,1,1,1,1,1,1,1,0;float dN=0,-27,-15,-15

13、,-15,-15,-15,-15,-15,-15,-15;float xN=0,0,0,0,0,0,0,0,0,0,0;float rN=0,0,0,0,0,0,0,0,0,0,0;float yN=0,0,0,0,0,0,0,0,0,0,0;float q;/clrscr();int k;r1=c1/b1;y1=d1/b1;for(k=2;k=1;k-)xk=yk-rk*xk+1;for(k=1;kN;k+)printf(x%d=%fn,k,xk);getch();return 0 ;l 运行窗口5 雅克比迭代法l 实验代码#include#include#define eps 1e-6#d

14、efine max 100void Jacobi(float *a,int n,float x) int i,j,k=0; float epsilon,s; float *y=new floatn; for(i=0;in;i+)xi=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)*xj; yi=(*(a+i*(n+1)+n)-s)/(*(a+i*(n+1)+i); epsilon+=fabs(yi-xi); for(i=0;in;i+)xi=yi;

15、if(epsilon=max) printf(迭代发散);return; delete y;void main() int i; float a45=10,-1,2,0,-11,0,8,-1,3,-11,2,-1,10,0,6,-1,3,-1,11,25; float x4; Jacobi(a0,4,x); for(i=0;i4;i+)printf(x%d=%fn,i,xi);l 运行窗口高斯-塞德尔迭代法l 实验代码#include#include#define N 500void main()int i;float x4;float c45=10,-1,2,0,-11,0,8,-1,3,-

16、11,2,-1,10,0,6,-1,3,-1,11,25;void GaussSeidel(float *,int,float);GaussSeidel(c0,4,x);for(i=0;i4;i+)printf(x%d=%fn,i,xi);void GaussSeidel(float *a,int n,float x)int i,j,k=1;float d,dx,eps;for(i=0;in;i+)xi=0.0;while(1)eps=0;for(i=0;in;i+)d=0;for(j=0;jn;j+)if(j=i)continue;d+=*(a+i*(n+1)+j)*xj;dx=(*(a+i

17、*(n+1)+n)-d)/(*(a+i*(n+1)+i);eps+=fabs(dx-xi);xi=dx;if(epsN)printf(迭代发散n);return;k+;l 运行窗口实习题四2. 按下列数据Xi0.300.420.500.580.660.72Yi1.044031.084621.118031.156031.198171.23223作5次插值，并求X1=0.46，X2=0.55，X3=0.60时的函数近似值。l 实验代码#include #define N 5void Difference(float x,float y,int n)float *f=new floatn+1;int

18、 k,i;for(k=1;k=n;k+)f0=yk;for(i=0;i=0;i-)a=a*(varx-xi)+yi;for(i=N-1;i=0;i-)b=b*(vary-xi)+yi;for(i=N-1;i=0;i-)c=c*(varz-xi)+yi;printf(Nn(%f)=%fn,varx,a);printf(Nn(%f)=%fn,vary,b);printf(Nn(%f)=%fn,varz,c);l 运行窗口实习题51.1） 试分别用抛物线y=a+bx+cxx和指数曲线y=ae(bx)拟合下列数据：Xi11.522.533.5Yi33.479.50122.65159.05189.152

19、14.15Xi44.555.566.5Yi238.65252.50267.55280.50296.65301.40XI77.58YI310.40318.15325.15实验运行代码：/最小二乘法拟合曲线#include#includeusing namespace std;const int n=15;/数据对个数const int m=3;/线性无关函数个数class ColPivotfriend class Approx;private:double cmm+1;double xm;public:ColPivot()for(int i=0;im;i+)for(int j=0;jm+1;j+

20、) cij=0;xi=0;void Cal();void ColPivot:Cal()int i,j,t,k;double p;for(i=0;im-1;i+)k=i;for(j=i+1;jfabs(cki) k=j;if(k!=j)for(j=i;jm+1;j+)p=cij;cij=ckj;ckj=p;for(j=i+1;jm;j+)p=cji/cii;for(t=i;t=0;i-)p=cim;for(j=m-1;ji;j-)p-=cij*xj;if(cii=0) xi=0;else xi=p/cii;class Approxprivate:double xn,yn,a3,b2;public

21、:void GetValue(double,double);void Para();void Expo();void Show();void Approx:GetValue(double an,double bn)for(int i=0;in;i+)xi=ai;yi=bi;void Approx:Para()int i,j,t;ColPivot col;for(i=0;i3;i+)for(j=0;j3;j+)for(t=0;tn;t+) col.cij+=pow(xt,i)*pow(xt,j);for(t=0;tn;t+) col.ci3+=yt*pow(xt,i);col.Cal();for

22、(i=0;i3;i+) ai=col.xi;void Approx:Expo()int i,j,t;ColPivot col;for(i=0;i2;i+)for(j=0;j2;j+)for(t=0;tn;t+) col.cij+=pow(xt,i)*pow(xt,j);for(t=0;tn;t+) col.ci3+=log(yt)*pow(xt,i);col.Cal();b0=exp(col.x0);b1=col.x1;void Approx:Show()cout用抛物线拟合：endl;couta=a0endl;coutb=a1endl;coutc=a2endl;cout用指数曲线拟合：endl;couta=b0endl;coutb=b1endl;int main() Approx ap;double x15;for(int i=0;i15;i+) xi=1+i*0.5;double y15=33.4,79.5,122.65,159.05,189.15,214.15,238.65,252.5,267.55,280.5,296.65,301.4,310.4,318.15,325.15;ap.GetValue(x,y);ap.Para();ap.Expo();ap.Show();return 0;专心-专注-专业