计算方法与实习实验报告学院：电气工程学院指导老师：***班级：160093******学号：********实习题一实验1 拉格朗日插值法一、方法原理n次拉格朗日插值多项式为：L n(x)=y0l0(x)+y1l1(x)+y2l2(x)+…+y n l n(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)二、主要思路使用线性方程组求系数构造插值公式相对复杂，可改用构造方法来插值。

对节点x i(i=0,1,…,n)中任一点x k(0<=k<=n)作一n 次多项式l k(x k)，使它在该点上取值为1，而在其余点x i(i=0,1,…,k-1,k+1,…,n)上为0，则插值多项式为L n(x)=y0l0(x)+y1l1(x)+y2l2(x)+…+y n l n(x) 上式表明：n 个点x i(i=0,1,…,k-1,k+1,…,n)都是l k(x)的零点。

可求得l k三．计算方法及过程：1.输入节点的个数n2.输入各个节点的横纵坐标3.输入插值点4.调用函数，返回z函数语句与形参说明程序源代码如下：#include<iostream>#include<math.h>using namespace std;#define N 100double fun(double *x,double *y, int n,double p);void main(){int i,n;cout<<"输入节点的个数n：";cin>>n;double x[N], y[N],p;cout<<"please input xiangliang x= "<<endl;for(i=0;i<n;i++)cin>>x[i];cout<<"please input xiangliang y= "<<endl;for(i=0;i<n;i++)cin>>y[i];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 L[N];while(k<n){ if(k==0){ for(i=1;i<n;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;i<n;i++) s=s*((p-x[i])/(x[k]-x[i]));L[k]=s*y[k];k++;}}for(i=0;i<n;i++)z=z+L[i];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)=b0＋b1x＋b2x2＋…＋bnxn (幂基) ①也可将其写成P(x)=a0＋a1(x－x0)＋a2(x－x0) (x－x1)＋…＋a n(x－x0) …(x－x n-1)其中ai为系数，xi为给定节点，可由①求出ai 一般情况下，牛顿插值多项式Nn(x)可写成：Nn(x)= a0＋a1(x－x0)＋a2(x－x0) (x－x1)＋…＋a n(x－x0) …(x－x n-1))只需求出系数ai，即可得到插值多项式。

二、计算方法及过程1.先后输入节点个数n和节点的横纵坐标，插值点的横坐标，最后输入精度e2. 用do-while循环语句得到跳出循环时k的值3.将k值与n-1进行比较，若在达到精度时k<n-1,则输出计算结果；若此时k=n-1,则计算失败！函数语句与形参说明程序源代码如下：#include<iostream>#include<math.h>using namespace std;#define MAX 100void main(){ float x[MAX],y[MAX];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>>x[i];cout<<"请输入对应的函数值！"<<endl;for(i=0;i<n;i++) cin>>y[i];cout<<"请输入插值点,x0=";cin>>x0;cout<<"请输入精度,e=";cin>>e;y0=1; N1=y[0]; k=0;do { k++;N0=N1;y0=y0*(x0-x[k-1]);for(i=0;i<k;i++) y[k]=(y[k]-y[i])/(x[k]-x[i]);N1=N0+y0*y[k];} while (fabs(N1-N0) > e && k<(n-1));if (k==(n-1)) cout<<"计算失败！";if (k<(n-1)) cout<<"输出结果y[x0]="<<N1<<endl; system("pause");}三．运行结果测试：实习题二1.用牛顿法求下列方程的根：3）lg20 x x+-=实验代码：#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;do{x1=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)=%f\t",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("请输入迭代初值x0\n");scanf("%f",&x0);printf("x(0)=%f\n",x0);y0=Newton(f,fl,x0);printf("方程的根为: %f\n",y0);}运行窗口：实习题三4．编写用追赶法解三对角线性方程组的程序，并解下列方程组：2）Ax b实验代码#include<iostream>#include<cmath>using namespace std;void main(){float a=1;float b=-4;float c=1;float d[11]={0,-27,-15,-15,-15,-15,-15,-15,-15,-15,-15};float l[11];float bb[11];float y[11];float x[11];bb[1]=b;y[1]=d[1];int i;for(i=2;i<11;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;1>0;i--){x[i]=(y[i]-c*x[i+1])/bb[i];}for(i=1;i<11;i++)cout<<'x'<<i<<':'<<x[i]<<endl;}运行窗口#include<stdio.h>#include<string.h>#include<math.h>#include<conio.h>#include<stdlib.h>#define N 11main(){float a[N]={0,0,1,1,1,1,1,1,1,1,1};float b[N]={0,-4,-4,-4,-4,-4,-4,-4,-4,-4,-4};float c[N]={0,1,1,1,1,1,1,1,1,1,0};float d[N]={0,-27,-15,-15,-15,-15,-15,-15,-15,-15,-15};float x[N]={0,0,0,0,0,0,0,0,0,0,0};float r[N]={0,0,0,0,0,0,0,0,0,0,0};float y[N]={0,0,0,0,0,0,0,0,0,0,0};float q;//clrscr();int k;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]=%f\n",k,x[k]);getch();return 0 ; }● 运行窗口5．12310211x x x -+=-2348311x x x -+=-1232106x x x -+=123431125x x x x -+-+=雅克比迭代法● 实验代码#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 float[n];for(i=0;i<n;i++)x[i]=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)*x[j];}y[i]=(*(a+i*(n+1)+n)-s)/(*(a+i*(n+1)+i));epsilon+=fabs(y[i]-x[i]);}for(i=0;i<n;i++)x[i]=y[i];if(epsilon<eps){printf("迭代次数为:%d\n",k);return;}if(k>=max){printf("迭代发散");return;}}delete y;}void main(){int i;float a[4][5]={10,-1,2,0,-11,0,8,-1,3,-11,2,-1,10,0,6,-1,3,-1,11,25};float x[4];Jacobi(a[0],4,x);for(i=0;i<4;i++)printf("x[%d]=%f\n",i,x[i]);}运行窗口高斯-塞德尔迭代法实验代码#include<stdio.h>#include<math.h>#define N 500void main(){int i;float x[4];float c[4][5]={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(c[0],4,x);for(i=0;i<4;i++)printf("x[%d]=%f\n",i,x[i]);}void GaussSeidel(float *a,int n,float x[]){int i,j,k=1;float d,dx,eps;for(i=0;i<n;i++)x[i]=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)*x[j];}dx=(*(a+i*(n+1)+n)-d)/(*(a+i*(n+1)+i));eps+=fabs(dx-x[i]);x[i]=dx;}if(eps<1e-6){printf("迭代次数为:%d\n",k);return;}if(k>N){printf("迭代发散\n");return;}k++;}}●运行窗口实习题四Xi 0.30 0.42 0.50 0.58 0.66 0.72Yi 1.04403 1.08462 1.11803 1.15603 1.19817 1.23223●实验代码#include <stdio.h>#define N 5void Difference(float x[],float y[],int n){float *f=new float[n+1];int k,i;for(k=1;k<=n;k++){f[0]=y[k];for(i=0;i<k;i++)f[i+1]=(f[i]-y[i])/(x[k]-x[i]);y[k]=f[k];}delete f;return;}void main(){int i;float a,b,c,varx=0.46,vary=0.55,varz=0.60;float x[N+1]={0.30,0.42,0.50,0.58,0.66,0.72};float y[N+1]={1.04403,1.08462,1.11803,1.15603,1.19817,1.23223};Difference(x,y,N);a=y[N];b=y[N];c=y[N];for(i=N-1;i>=0;i--)a=a*(varx-x[i])+y[i];for(i=N-1;i>=0;i--)b=b*(vary-x[i])+y[i];for(i=N-1;i>=0;i--)c=c*(varz-x[i])+y[i];printf("Nn(%f)=%f\n",varx,a);printf("Nn(%f)=%f\n",vary,b);printf("Nn(%f)=%f\n",varz,c);}运行窗口实习题51.Xi 1 1.5 2 2.5 3 3.5Yi 33.4 79.50 122.65 159.05 189.15 214.15 Xi 4 4.5 5 5.5 6 6.5Yi 238.65 252.50 267.55 280.50 296.65 301.40实验运行代码：//最小二乘法拟合曲线#include<iostream>#include<cmath>using namespace std;const int n=15;//数据对个数const int m=3;//线性无关函数个数class ColPivot{friend class Approx;private:double c[m][m+1];double x[m];public:ColPivot(){for(int i=0;i<m;i++){for(int j=0;j<m+1;j++) c[i][j]=0;x[i]=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(c[j][i])>fabs(c[k][i])) k=j;}if(k!=j){for(j=i;j<m+1;j++){p=c[i][j];c[i][j]=c[k][j];c[k][j]=p;}}for(j=i+1;j<m;j++){p=c[j][i]/c[i][i];for(t=i;t<m+1;t++){c[j][t]-=p*c[i][t];}}}for(i=m-1;i>=0;i--){p=c[i][m];for(j=m-1;j>i;j--){p-=c[i][j]*x[j];}if(c[i][i]==0) x[i]=0;else x[i]=p/c[i][i];}}class Approx{private:double x[n],y[n],a[3],b[2];public:void GetValue(double[],double[]);void Para();void Expo();void Show();};void Approx::GetValue(double a[n],double b[n]) {for(int i=0;i<n;i++){x[i]=a[i];y[i]=b[i];}}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.c[i][j]+=pow(x[t],i)*pow(x[t],j);}for(t=0;t<n;t++) col.c[i][3]+=y[t]*pow(x[t],i);}col.Cal();for(i=0;i<3;i++) a[i]=col.x[i];}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.c[i][j]+=pow(x[t],i)*pow(x[t],j);}for(t=0;t<n;t++) col.c[i][3]+=log(y[t])*pow(x[t],i);}col.Cal();b[0]=exp(col.x[0]);b[1]=col.x[1];}void Approx::Show(){cout<<"用抛物线拟合："<<endl;cout<<"a="<<a[0]<<endl;cout<<"b="<<a[1]<<endl;cout<<"c="<<a[2]<<endl;cout<<"用指数曲线拟合："<<endl;cout<<"a="<<b[0]<<endl;cout<<"b="<<b[1]<<endl;}int main(){Approx ap;double x[15];for(int i=0;i<15;i++) x[i]=1+i*0.5;doubley[15]={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;}。