计算数值实验报告(太原理工大学)

v1.0可编辑可修改v1.0可编辑可修改本科实验报告课程名称：计算数值方法 实验地点：综合楼五层506室专业班级：计科1002学号：J4学生姓名：xxx 指导教师：王峥2012年6月20

太原理工大学学生实验报告学院名称计算机科学与技术专业班级计科1002学号14学生姓名xxx实验日期成绩课程名称计算数值方法实验题目实验一方程求根一、课题名称方程求根：熟悉使用、迭代法、牛顿法、割线法等方法对给定的方程进行根的求解。选择上述方法中的两种方法求方程： 一分法f(x)=x+4x-10=0在[1,2]内的一个头根，且要求满足精度|x*-xn|<X10-5迭代法：用迭代公式x=f(x)进行迭代计算，直到满足|x-xn|<Xl0—5为止。二分法：设f(x)在[a,b]上连续，且f(al)*f(x1)<0,记(a2,b2)=(x1,bl)带入计算式进行计算直到|x*-xn|<X10-5为止。二、目的和意义(1)了解非线性方程求根的常见方法，如二分法、迭代法、牛顿法、割线法。(2)加深对方程求根方法的认识，掌握算法。会进行误差分析，并能对不同方法进行比较。三、计算公式(1)迭代法).首先对给定的计算公式进行变形使其能够迭代或者找出相应迭代速度较快的式子。).带入求好的式子到循环中去比如：xk1 (xk) (k0,1,2,L)(2)二分法：f(x)在区间(x,y)上连续.先找到a、b属于区间(x,y),使f(a),f(b)异号，说明在区间(a,b)内一定后零点，然后求f[(a+b)/2],.如果f[(a+b)/2]=0 ，该点就是零点，如果f[(a+b)/2]<0,则在区间((a+b)/2,b)内肩零点,反之在(a,(a+b)/2)内后零点带入1)中继续。四、主要仪器设备Vc++C-freeCodeBlocks五、结构程序设计迭代法：#include<>#include<>main(){inti;doublexn[15],y,x1,x2,m;printf("请输入x1,x2的值：\n");scanf("%lf%lf",&x1,&x2);printf(" 请输入精度要求：\n");scanf("%lf",&m);printf("nxn'n");i=0;do{xn[0]=(x1+x2)/2;xn[i+1]= sqrt(10/(4+xn[i]));__ _150.31059.1431X159.174215X1123Xi145.2916.13012X246.7887210X2012X28211311.2952X314836X34X31211X421261120X41X1721X25121Xn1512 Xn5追赶法:用来求对角方程组；将系数矩阵A转化为A=L*U,L为普通下n-1对角矩阵，U为单位上n-1对角矩阵，然后通过用来求对角方程组；将系数矩阵程组l*y=b,u*x=y,来求解x.、主要仪器设备c++C-freeCodeBlocks、结构程序设计Gauss消元法：#include<iostream>#include<cmath>usingnamespacestd;intmain(){intn,i,j,k;doublea[100][100],b[100],o;cout<<"输入未知数个数："<<endl;cin>>n;cout<<"输入数列:"<<endl;for(i=1;i<=n;i++)for(j=1;j<=n+1;j++)cin>>a皿];for(i=1;i<=n;i++)for(j=i+1;j<=n;j++)if(fabs(a[j][i])>1e-7){o=a[i][i]/a皿i];for(k=i;k<=n+1;k++)a[j][k]=a[j][k]*o-a[i][k];}for(i=n;i>0;i--){b[i]=a皿n+1]/a[i][i]；for(j=i-1;j>0;j--)a皿n+1]=a[j][n+1]-b[i]*a皿i]；}cout<<"解得:"<<endl;for(i=1;i<=n;i++)cout<<b[i]<<endl;lf\n",i,d[i])；}、结果讨论和分析消元法： 列主元素消元法:解法： 追赶法:

讨论从消元过程可以看出，对于n阶线性方程组，只要各步主元素不为零，经过n-1步消元，就可以得到一个等价的系数三角形阵的方程组，然后再利用回代过程可求得原方程组的解.由于列主元素法相似且优于完全主元素法 所以省北消元过程相当于分解A为单位下三角阵L与上三角阵U的乘积，解方程组Ly=b回代过程就是解方程组Ux=y。其中白单位下三角阵、U为上三角阵.在A的LU分解中,L取下三角阵,U取单位上三角阵,这样求解方程组Ax=d的方?法。另外是追赶法和其他方法求同一方程结果不一样，我多次修改源程序，也不知道原因。再就是追赶法有很大的局改良。图：实验地点综合楼五层506室指导教师王峥太原理工大学学生实验报告学院名称计算机科学与技术专业班级计科1002学号14学生姓名xxx实验日期成绩课程名称计算数值方法实验题目实验三线性方程组的迭代解法、课题名称线性方程组的迭代解法使用雅可比迭代法或高斯-赛德尔迭代法对下取」方程组进行求解。10xiX22x3 7.2Xi10X22x38.3x1x25x34.2二、目的和意义学习使用雅可比迭代法或高斯-赛德尔迭代法三、计算公式雅克比迭代法：设线性方程组Ax=b的系数矩阵A可逆且主对角元素aii,a22,…,ann均不为零，令D=diag(aii,a22,…，ann)并将A分解成A=(A-D)+D从而线性方程组可写成Dx=(D-A)x+b则有迭代公式x(k+1)=Bx(k)+fi其中Bi=I-D-1A,fi=D1b。四、主要仪器设备Vc++C-freeCodeBlocks五、结构程序设计雅克比迭代法:#include<>#include<>main(){inti;doublex1[20],x2[20],x3[20];doublex10,x20,x30;printf("请输入x1,x2,x3的初值：\n");scanf("%lf%lf%lf",&x10,&x20,&x30);printf("nx1[n] x2[n] x3[n]\n");for(i=0;i<18;i++){x1[0]=x10;x2[0]=x20;x3[0]=x30;x1[i+1]=*x2[i]+*x3[i]+;x2[i+1]=*x1[i]+*x3[i]+;x3[i+1]=*x1[i]+*x2[i]+;printf("%5d%51f%51f%5lf\n",i,x1[i],x2[i],x3[i]);}}六、实验结果与分析:雅克比迭代法:

鼾通Axl・,乙-四初值？d6£1n xl[n1x2Cn]x3tn]n 叼《目h目眄④0.0000000.0B00001 0,7200080.B30000@-84000021.0700001.1500003 1.0B70901.1571001.24820U4 1,0853501,1853401.2828365 1.0950981.1?50991.294138£ 1.0983381.1?8337t-2980397 1.099442L.l?94421-2992358 1.099B111-2997779 1.0999361.1999361.29992410 1.0999791.199979i-2999^11 1.0999991-1999931-29999112 1.0999981-29999713 1.0999??L-1.29299214 1.100039j,[Gm1.3Q000015 1.1008301.30000016 1.1010。11.30000017 1<1000981.300000分析讨论：其实，这两个迭代法是之前迭代法的升级，多了几个迭代式子而已，而且两者相差不大比较简单，所以选择了雅克比迭代法进行求解， 但是没有与另一种方法高斯赛德尔迭代法进行实质性的比较。流程图：实验地点综合楼五层506室 指导教师王峥太原理工大学学生实验报告学院名称计算机科学与技术专业班级计科1002学号14学生姓名xxx实验日期成绩课程名称计算数值方法实验题目实验四矩阵特征值与特征向量问题、课题名称使用幕法求A模为最大的特征值及其相应的特征向量。TOC\o"1-5"\h\z2 1 0A12 10 1 2二、目的和意义(1)了解矩阵特征值与特征向量问题解法，掌握幕法。(2)加深对矩阵特征值与特征向量问题求解方法的认识，掌握算法。三、计算公式幕法：由已知的非零向量x0和矩阵A的乘幕构造向量序列｛xn｝以计算矩阵A的按模最大特征值及其特征向量的方法，称为幕法。迭代公式：ykAxk1mkmax(yk),k1,2,...Xkyk

mXk结果可取1 mk1 yk或1 xk四、主要仪器设备Vc++C-freeCodeBlocks五、结构程序设计源代码：#include<>#include<>#defineN3#defineeps1e-6#defineKM30floatMaxValue(floatx[],intn){floatMax=x[0];inti;for(i=1;i<n;i++)if(fabs(x[i])>fabs(Max))Max=x[i];returnMax;}voidPowerMethod(float*A){floatU[N],V[N],r1,r2,temp;inti,j,k=0;while(k<KM){k++;for(i=0;i<N;i++){temp=0;for(j=0;j<N;j++)temp+=*(A+i*N+j)*U[j];V[i]=temp;}for(i=0;i<N;i++)U[i]=V[i]/MaxValue(V,N);if(k==1)r1=MaxValue(V,N);elser2=MaxValue(V,N);if(fabs(r2-r1)<eps)break;r1=r2;}printf("r=%f\n",r2);for(i=0;i<N;i++)printf("y[%d]=%f\n",i+1,U[i]);}voidmain(){floatA[N][N尸{{2,-1,0},{-12-1},{0,-1,2}};floatU[N];//A 的值U[0]=1;U[1]=1;U[2]=1;//x0的值PowerMethod(A[0]);}六、结果讨论和分析分析讨论由于该程序将A矩阵和x0的值编写在程序中，所以要想修改成其他矩阵时比较麻烦，所以也有一定的局限性。幕法是一种求任意矩阵A的按模最大特征值及其对应特征向量的迭代算法。该方法的最大优点是计算简单，容易在计算机上实现，对稀疏矩阵较为适合，但有时收敛速度很慢。流程图：实验地点综合楼五层506实验地点综合楼五层506室 指导教师王峥太原理工大学学生实验报告学院名称计算机科学与技术专业班级计科1002学号14学生姓名XXX实验日期成绩课程名称计算数值方法实验题目实验五代数插值一、课题名称(使用拉格朗日插值法或牛顿插值法求解：已知f(x)在6个点的函数值如下表所示，运用插值方法，求f的近似值。Xf(x)、目的和意义学习使用拉格朗日插值法或牛顿插值法求解三、计算公式设函数在区间［a,b］上n+1互异节点Xo,X1,…,xn上的函数值分别为yo,yi,…,yn,求n次插值多项式Pn(x),满足条件Pn(Xj)=yj,j=0,1,…,n令Ln(x)=yolo(x)+yili(x)+…+ynln(x)=/y/i(x)其中lo(x),ll(x),…，ln(x)为以X0,X1，…,Xn为节点的n次插值基函

…，n数，则Ln(x)是一次数不超过n…，nLn(Xj)=yj,L=0,1,再由插值多项式的唯一性，得Pn(X)三Ln(X)四、主要仪器设备Vc++C-freeCodeBlocks五、结构程序设计#include<>#include<>main(){intI;charL;doubleM[100][100];doublex[100],y[100];doubleX=1,xx=0,w=1,N=0,P,R=1;intn=5;//cout<<”请输入所求均差阶数：";//求所有阶差〃cin>>n;//for(inti=0;i<=n;i++)//{/*cout<<"请输入x"<<i<<"的值："<<endl;cin>>x[i];cout<<"请输入y"<<i<<"的值："<<endl;cin>>y[i];M[i][0]=x[i];M[i][1]=y[i];*/// 用二维保存所有数据M[0][0]=;M[0][1]=;M[1][0]=;M[1][1]=;M[2][0]=;M[2][1]=;M[3][0]=;M[3][1]=;M[4][0]=;M[4][1]=;M[5][0]=;M[5][1]=;for(intj=2;j<=n+1;j++)(for(inti=1;i<=n;i++)(M[i][j]=(M[i][j-1]-M[i-1][j-1])/(M[i][0]-M[i-j+1][0]);}}for(inti=1;i<=n;i++)(cout<<"其"<<i<<”阶均差为:"<<M[i][i+1]<<endl;}cout<<"请输入x的值：x=";cin>>xx;for(inti=0;i<n;i++)(X*=xx-M[i][0];N+=M[i+1][i+2]*X;P=M[0][1]+N;}cout<<"其函数值：y="<<P<<endl;//})分析讨论拉格朗日插值的优点是插值多项式特别容易建立，缺点是增加节点是原有多项式不能利用，必须重新建立，即所有基函数都要重新计算，这就造成计算量的浪费。所以该程序选择了牛顿法，由于在输入数据的时候比较麻烦 所以将数据存在了程序里面，只能在程序里修改数据。 另外由于精度问题计算结果有一定的误差。流程图：实验地点综合楼五层506室指导教师王峥太原理工大学学生实验报告学院名称计算机科学与技术学生姓名xxx课程名称计算数值方法专业班级学院名称计算机科学与技术学生姓名xxx课程名称计算数值方法专业班级计科1002学号14实验日期成绩实验题目实验六最小二乘法拟合多项式、课题名称给定数据点（Xi,yi）,用最小二乘法拟合数据的多项式，并求平方误差。Xi0V、1、目的和意义.熟练运用已学计算方法求解方程组.加深对计算方法技巧，选择正确的计算方法来求解各种方程组.培养使用电子计算机进行科学计算和解决问题的能力三、计算公式建立正规方程组：■(二Xij+k)ak=EXijyi,j=0,1，…,n平方误差：=.(.akXi-yi)四、主要仪器设备Vc++C-freeCodeBlocks五、结构程序设计源代码：#include<>#include<>#defineN15doublepower(double&a,intn){doubleb=1;for(inti=0;i<n;i++)b*=a;returnb;}voidGauss();doubleX[N],Y[N],sumX[N],sumY[N],a[N][N],b[N],l[N][N],x[N];

voidmain()(ofstreamoutdata;ifstreamindata;doubles;inti,j,k,n,index;//cout<<"请输入已知点的个数n=";〃cin>>n;n=7;cout<<endl;输入给定数据循环输入任何数据//cout<<"请输入X和Y:"<<endl;//输入给定数据循环输入任何数据X[0]=;Y[0]二;X[1]=;Y[1]=;X[2]=;Y[2]二;X[3]=;Y[3]二;X[4]=;Y[4]二;X[5]=;Y[5]二;X[6]=;Y[6]=;〃 绑定数据〃可以解绑由下forfor(i=0;i<n;i++)(//cout<<"X["<<i<<"]=";〃cin>>X[i];sumX[1]+=X[i];//cout<<"Y["<<i<<"]=";〃cin>>Y[i];sumY[1]+=Y[i];//cout<<endl;)cout<<"sumX[1]="<<sumX[1]<<"\t"<<"sumY[1]="<<sumY[1]<<endl;cout<<"请输入拟合次数index=";cin>>index;cout<<endl;i=n;sumX[0]=i;for(i=2;i<=2*index;i++){sumX[i]=0;for(j=0;j<n;j++)sumX[i]+=power(X[j],i);cout<<"sumX["<<i<<"]="<<sumX[i]<<endl;)for(i=2;i<=index+1;i++){sumY[i]=0;for(j=0;j<n;j++)sumY[i]+=power(X[j],i-1)*Y[j];cout<<"sumY["<<i<<"]="<<sumY[i]<<endl;)for(i=1;i<=index+1;i++){ // 建立正规方程组for(j=1;j<=index+1;j++)a[i][j]=sumX[i+j-2];b[i]=sumY[i];)k=1; // 用高斯消元法解方程组do{for(j=k+1;j<=index+1;j++)l[j][k]=a[j][k]/a[k][k];for(i=k+1;i<=index+1;i++){for(j=k+1;j<=index+1;j++)a[i][j]=a[i][j]-l[i][k]*a[k][j];b[i]=b[i]-l