东南大学计算方法实验报告(共16页).doc
精选优质文档-倾情为你奉上 计算方法与实习实验报告学 院：电气工程学院指导老师：李翠平班 级：姓 名：黄芃菲学 号：实习题一实验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.输入节点的个数n 2.输入各个节点的横纵坐标 3.输入插值点 4.调用函数，返回z函数语句与形参说明程序源代码如下：形参与函数类型参数意义int n节点的个数double xn（double *x）存放n个节点的值double yn（double *y）存放n个节点相对应的函数值double p指定插值点的值double fun()函数返回一个双精度实型函数值，即插值点p处的近似函数值 #include<iostream>#include<math.h>using namespace std;#define N 100 double fun(double *x,double *y, int n,double p);void main()int i,n;cout<<"输入节点的个数n：" cin>>n;double xN, yN,p;cout<<"please input xiangliang x= "<<endl;for(i=0;i<n;i+)cin>>xi;cout<<"please input xiangliang y= "<<endl;for(i=0;i<n;i+)cin>>yi;cout<<"please input LagelangrichazhiJieDian p= "<<endl; cin>>p;cout<<"The 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(k<n) if(k=0) for(i=1;i<n;i+)s=s*(p-xi)/(x0-xi); 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;i<n;i+) s=s*(p-xi)/(xk-xi); Lk=s*yk; k+;for(i=0;i<n;i+)z=z+Li;return z;五实验分析n=2时，为一次插值，即线性插值n=3时，为二次插值，即抛物线插值n=1，此时只有一个节点，插值点的值就是该节点的函数值n<1时，结果都是返回0的；这里做了n=0和n=-7两种情况3<n<100时，也都有相应的答案常用的是线性插值和抛物线插值，显然，抛物线精度相对高些n次插值多项式Ln（x）通常是次数为n的多项式，特殊情况可能次数小于n.例如：通过三点的二次插值多项式L2(x),如果三点共线，则y=L2(x)就是一条直线，而不是抛物线，这时L2(x)是一次式。拟合曲线光顺性差 实验2 牛顿插值法 一、方法原理及基本思路在拉格朗日插值方法中，若增加一个节点数据，其插值的多项式需重新计算。现构造一个插值多项式Nn(x)，只需对Nn-1(x)作简单修正(如增加某项)即可得到，这样计算方便。利用牛顿插值公式，当增加一个节点时，只需在后面多计算一项，而前面的计算仍有用；另一方面Nn(x)的各项系数恰好又是各阶差商，而各阶差商可用差商公式来计算。由线性代数知，对任何一个不高n次的多项式P(x)=b0b1xb2x2bnxn (幂基) 也可将其写成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进行比较，若在达到精度时k<n-1,则输出计算结果；若此时k=n-1,则计算失败！函数语句与形参说明形参与函数类型参数意义int n节点的个数float xMAX存放n个节点的值（从小到大）Float yMAX;存放n个节点相对应的函数值float x0,y0指定插值点的横纵坐标float e精度 程序源代码如下：#include<iostream>#include<math.h>using namespace std;#define MAX 100 void main() float xMAX,yMAX; float x0,y0,e,N1; float N0=0; int i,k,n; cout<<"输入节点的个数,n=" cin>>n; cout<<"请输入插值节点！"<<endl; for(i=0;i<n;i+) cin>>xi; cout<<"请输入对应的函数值！"<<endl; for(i=0;i<n;i+) cin>>yi; cout<<"请输入插值点,x0=" cin>>x0; cout<<"请输入精度,e=" cin>>e; y0=1; N1=y0; k=0; do k+; N0=N1; y0=y0*(x0-xk-1); for(i=0;i<k;i+) yk=(yk-yi)/(xk-xi); N1=N0+y0*yk; while (fabs(N1-N0) > e && k<(n-1); if (k=(n-1) cout<<"计算失败！" if (k<(n-1) cout<<"输出结果yx0="<<N1<<endl;system("pause"); 三运行结果测试：题一：已知f(x)=sh(x)的函数表如下：计算f(0.23)的近似值xi00.200.300.50Sh(xi)00.201340.304520.52110 实习题二1. 用牛顿法求下列方程的根：3）实验代码：#include <stdio.h>#include <math.h>#define N 100#define 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)<1?x1-x0:(x1-x0)/x1;x0=x1;printf("x(%d)=%ft",k,x0);while(fabs(d)>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;void 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<iostream>#include<cmath>using 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;float bb11;float y11;float x11;bb1=b;y1=d1;int i;for(i=2;i<11;i+)li=a/bbi-1;bbi=b-li*c;yi=di-li*yi;x10=y10/bb10;for(i=9;1>0;i-)xi=(yi-c*xi+1)/bbi;for(i=1;i<11;i+)cout<<'x'<<i<<':'<<xi<<endl;l 运行窗口#include<stdio.h>#include<string.h>#include<math.h>#include<conio.h>#include<stdlib.h>#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,-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<N-1;k+)q=bk-rk-1*ak;rk=ck/q;yk=(dk-yk-1*ak)/q;yN-1=(dN-1-yN-2*aN-1)/(bN-1-rN-2*aN-1);xN-1=yN-1;for(k=N-2;k>=1;k-)xk=yk-rk*xk+1;for(k=1;k<N;k+)printf("x%d=%fn",k,xk);getch();return 0 ;l 运行窗口5 雅克比迭代法l 实验代码#include<stdio.h>#include<math.h>#define eps 1e-6#define max 100void Jacobi(float *a,int n,float x) int i,j,k=0; float epsilon,s; float *y=new floatn; for(i=0;i<n;i+)xi=0; while(1) epsilon=0; k+; for(i=0;i<n;i+) s=0; for(j=0;j<n;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;i<n;i+)xi=yi; if(epsilon<eps) printf("迭代次数为:%dn",k);return; if(k>=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;i<4;i+)printf("x%d=%fn",i,xi);l 运行窗口高斯-塞德尔迭代法l 实验代码#include<stdio.h>#include<math.h>#define N 500void main()int i;float x4;float c45=10,-1,2,0,-11,0,8,-1,3,-11,2,-1,10,0,6,-1,3,-1,11,25;void GaussSeidel(float *,int,float);GaussSeidel(c0,4,x);for(i=0;i<4;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;i<n;i+)xi=0.0;while(1)eps=0;for(i=0;i<n;i+)d=0;for(j=0;j<n;j+)if(j=i)continue;d+=*(a+i*(n+1)+j)*xj;dx=(*(a+i*(n+1)+n)-d)/(*(a+i*(n+1)+i);eps+=fabs(dx-xi);xi=dx;if(eps<1e-6)printf("迭代次数为:%dn",k);return;if(k>N)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 <stdio.h>#define N 5void Difference(float x,float y,int n)float *f=new floatn+1;int k,i;for(k=1;k<=n;k+)f0=yk;for(i=0;i<k;i+)fi+1=(fi-yi)/(xk-xi);yk=fk;delete f;return;void main()int i;float a,b,c,varx=0.46,vary=0.55,varz=0.60;float xN+1=0.30,0.42,0.50,0.58,0.66,0.72;float yN+1=1.04403,1.08462,1.11803,1.15603,1.19817,1.23223;Difference(x,y,N);a=yN;b=yN;c=yN;for(i=N-1;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.15214.15Xi44.555.566.5Yi238.65252.50267.55280.50296.65301.40XI77.58YI310.40318.15325.15实验运行代码：/最小二乘法拟合曲线#include<iostream>#include<cmath>using 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;i<m;i+)for(int j=0;j<m+1;j+) cij=0;xi=0;void Cal();void ColPivot:Cal()int i,j,t,k;double p;for(i=0;i<m-1;i+)k=i;for(j=i+1;j<m;j+)if(fabs(cji)>fabs(cki) k=j;if(k!=j)for(j=i;j<m+1;j+)p=cij;cij=ckj;ckj=p;for(j=i+1;j<m;j+)p=cji/cii;for(t=i;t<m+1;t+)cjt-=p*cit;for(i=m-1;i>=0;i-)p=cim;for(j=m-1;j>i;j-)p-=cij*xj;if(cii=0) xi=0;else xi=p/cii;class Approxprivate:double xn,yn,a3,b2;public:void GetValue(double,double);void Para();void Expo();void Show();void Approx:GetValue(double an,double bn)for(int i=0;i<n;i+)xi=ai;yi=bi;void Approx:Para()int i,j,t;ColPivot col;for(i=0;i<3;i+)for(j=0;j<3;j+)for(t=0;t<n;t+) col.cij+=pow(xt,i)*pow(xt,j);for(t=0;t<n;t+) col.ci3+=yt*pow(xt,i);col.Cal();for(i=0;i<3;i+) ai=col.xi;void Approx:Expo()int i,j,t;ColPivot col;for(i=0;i<2;i+)for(j=0;j<2;j+)for(t=0;t<n;t+) col.cij+=pow(xt,i)*pow(xt,j);for(t=0;t<n;t+) col.ci3+=log(yt)*pow(xt,i);col.Cal();b0=exp(col.x0);b1=col.x1;void Approx:Show()cout<<"用抛物线拟合："<<endl;cout<<"a="<<a0<<endl;cout<<"b="<<a1<<endl;cout<<"c="<<a2<<endl;cout<<"用指数曲线拟合："<<endl;cout<<"a="<<b0<<endl;cout<<"b="<<b1<<endl;int main() Approx ap;double x15;for(int i=0;i<15;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;专心-专注-专业
