计算数值试验报告太原理工大学

1、v1.0可编辑可修改本科实验报告课程名称：计算数值方法实验地点：综合楼五层506室专业班级：计科 1002学号：14学生姓名：XXX指导教师：王峥2012 年 6 月 20太原理工大学学生实验报告学院名称计算机科学与技术专业班级计科1002学号14学生姓名xxx实验日期成绩课程名称计算数值方法实验题目实验一方程求根、课题名称方程求根：熟悉使用、迭代法、牛顿法、割线法等方法对给定的方程进行根的求解。选择上述方法中的两种方法求方程：二分法f(x)=x 3+4x2-10=0在1,2内的一个实根，且要求满足精度|x *-x n|< X 10-5迭代法：用迭代公式 x=f(x)进行迭代计算，直到满

2、足|x*-Xn|v X 10-5为止。二分法：设 f (x)在a,b上连续，且 f (al) *f (xl) <0,记(a2, b2) = (x1, bl) 带入计算式进行计算 直到|x *-Xn|V X 10-5为止。、目的和意义(1) 了解非线性方程求根的常见方法，如二分法、迭代法、牛顿法、割线法。(2) 加深对方程求根方法的认识，掌握算法。会进行误差分析，并能对不同方法进行 比较。三、计算公式(1)迭代法1 ).首先对给定的计算公式进行变形使其能够迭代或者找出相应迭代速度较快的 式子。2 ).带入求好的式子到循环中去比如：Xk i(xj(k 0,1, 2，l )(2)二分法：f(x

3、)在区间(x，y)上连续1) .先找到a、b属于区间(x，y)，使f(a)，f(b)异号，说明在区间(a,b)内一 定有零点，然后求 f(a+b)/2,2) .如果f(a+b)/2=0，该点就是零点，如果f(a+b)/2<0,则在区间(a+b)/2 ，b)内有零点，反之在(a，(a+b)/2 )内有零点带入1 )中继续。四、主要仪器设备Vc+ C-free CodeBlocks五、结构程序设计迭代法：#in clude<>#in clude<>mai n()int i;double xn 15,y,x1,x2,m ;printf("请输入 x1,x2 的

4、值：n");sca nf("%lf%lf", &x1, &x2);prin tf（"请输入精度要求：n"）;sca nf("%lf",&m);printf(" n xnn”);i=0;doxn0=(x1+x2)/2 ;xn i+1=sqrt(10/(4+x ni);370.3 10 1559.1431X159.17123X1145.2916.13012X246.78012X2811.29521X3241X3131211X42x174215x,87210x24836x31261120x4X25

5、21 Xn 112 Xn追赶法:用来求对角方程组；将系数矩阵A转化为A=L*U, L为普通下n-1对角矩阵，U为单位上n-1对角矩阵，然后通过程组l*y=b,u*x=y,来求解x.、主要仪器设备C+ C-free CodeBlocks、结构程序设计Gauss消肖元法：#i nclude <iostream>#in clude <cmath>using n amespace std;int mai n()int n,i,j,k;double a100100,b100,o;cout<<"输入未知数个数:"<<endl;cin

6、87;n;cout<<"输入数列:"<<endl;for (i=1;i<=n ;i+)for (j=1;j<=n +1;j+)cin >>aij;for (i=1;i<=n ;i+)for (j=i+1;j<=n ;j+)if (fabs(aji)>1e-7)o=aii/aji;for (k=i;k<=n+1;k+)ajk=ajk*o-aik;for (i=n; i>0;i-)bi=ai n+1/aii;for (j=i-1;j>0;j-)aj n+1=aj n+1-bi*aji;cout&

7、lt;<"解得:"<<endl;for (i=1;i<=n ;i+) cout<<bi<<e ndl;lfn",i,di);、结果讨论和分析消元法：列主元素消元法:解法:追赶法:讨论从消元过程可以看出，对于n阶线性方程组，只要各步主元素不为零，经过n-1步消元，就可以得到一个等价的系数三角形阵的方程组，然后再利用回代过程可求得原方程组的解由于列主元素法相似且优于完全主元素法所以省略消元过程相当于分解A为单位下三角阵L与上三角阵U的乘积，解方程组Ly=b回代过程就是解方程组Ux=y。其中白 单位下三角阵、U为上三角阵.

8、在A的LU分解中,L取下三角阵,U取单位上三角阵,这样求解方程组Ax=d的方法 法。另外是追赶法和其他方法求同一方程结果不一样，我多次修改源程序，也不知道原因。再就是追赶法有很大的局 改良。图：实验地点综合楼五层506室指导教师王峥太原理工大学学生实验报告学院名称计算机科学与技术专业班级计科1002学号14学生姓名xxx实验日期成绩课程名称计算数值方法实验题目实验三线性方程组的迭代解法、课题名称线性方程组的迭代解法使用雅可比迭代法或高斯-赛德尔迭代法对下列方程组进行求解。10禺x2 2x37.2禺 10x2 2x38.3XiX2 5X34.2-、目的和意义学习使用雅可比迭代法或高斯-赛德尔迭代

9、法三、计算公式雅克比迭代法：设线性方程组Ax=b的系数矩阵A可逆且主对角元素aii,a22,ann均不为零，令D=diag(an,a 22,a nn)并将A分解成A=(A-D)+D从而线性方程组可写成Dx=(D-A) x+b则有迭代公式(k+1)(k)-x =Bx +f 1-1-1其中，Bi=l-D A,f i=D bo四、主要仪器设备Vc+ C-free CodeBlocks五、结构程序设计雅克比迭代法：#in clude<>#in clude<>main ()int i;double x120 ,x220,x320;double x10, x20, x30;prin

10、tf(”请输入 x1,x2,x3 的初值：n");sca nf("%lf%lf%lf", &x10, &x20, &x30);prin tf(" n x1 nx2 nx3 n n");for(i=0;i<18;i+)x10=x10;x20=x20;x30=x30;x1i+1=*x2i+*x3i+;x2i+1=*x1i+*x3i+;x3i+1=*x1i+*x2i+;prin tf("%5d %5lf %5lf %5lfn",i,x1i,x2i,x3i); 六、实验结果与分析：雅克比迭代法：a时訂n

11、cl n Jx2£nx3tnRta.aonaoo队 00000010.7200000.B30000.84000020.971000丄.070眄®01-15000031.0570901.1571001-24820041.QS535&1.1S53401.28282051.095098li950991-294138£1.0982381.19833?1-29803971.0994421.1994421-29922581.0996111.198111.29977291.0999361.1999361-299924101.0999791.1999791-29997S1

12、11-0999931-1999931-299991121.0999981A999981.299997131.0999991.1?1.299999141.1000991 9AAAAA1 3S0000IS1A00&&9*.300000161.100890.30000Bj?j lf10000.300600分析讨论：其实，这两个迭代法是之前迭代法的升级，多了几个迭代式子而已，而且两者相差不大比较简单，所以选择了雅克比迭代法进行求解， 但是没有与另一种方法 高斯赛德尔迭代法进行实质 性的比较。流程图：实验地点综合楼五层506室指导教师王峥太原理工大学学生实验报告学院名称计算机科学与技术专

13、业班级计科1002学号 14学生姓名xxx实验日期成绩课程名称计算数值方法实验题目实验四矩阵特征值与特征向量问题、课题名称使用幕法求A模为最大的特征值及其相应的特征向量。、目的和意义（1）了解矩阵特征值与特征向量问题解法，掌握幕法。（2）加深对矩阵特征值与特征向量问题求解方法的认识，掌握算法。二、计算公式幕法：由已知的非零向量x0和矩阵A的乘幕构造向量序列xn以计算矩阵 A的按模最大特征值及其特征向量的方法，称为幕法。迭代公式:y k Ax k 1m k max( yQ, k1,2，Xkyk mk结果可取1mkiyk或 iXk四、主要仪器设备Vc+ C-free CodeBlocks五、结构程

14、序设计源代码：#in clude<>#in clude<>#define N 3#defi ne eps 1e-6#define KM 30float MaxValue(float x,int n)float Max=x0;int i;for (i=1;i <n ;i+) if(fabs(xi)>fabs(Max)Max=xi; return Ma x;void PowerMethod(float *A)float UN,VN,r1,r2,temp;int i,j,k=0;while(kvKM)k+;for(i=0;i<N;i+)temp=0;for(

15、j=0;j<N;j+)temp+=*(A+i*N+j)*Uj;Vi=temp;for(i=0;i<N;i+)Ui=Vi|/MaxValue(V,N);if(k=1)r 1=MaxValue(V,N);else r2=MaxValue(V,N);if(fabs(r2-r1)<eps)break;r1=r2;prin tf("r=%fn",r2);for(i=0;i<N;i+)pri ntf("y%d=%fn",i+1,Ui);void mai n()float ANN=2,-1,0,-1,2,-1,0,-1,2 ;float UN;

16、 /A的值U0=1;U1=1;U2=1;/x0 的值PowerMethod(A0);六、结果讨论和分析分析讨论由于该程序将A矩阵和x0的值编写在程序中，所以要想修改成其他矩阵时比较麻烦，所以也有一定的局限性。幕法是一种求任意矩阵A的按模最大特征值及其对应特征向量的迭代算法。该方法的最大优点是计算简单，容易在计 算机上实现，对稀疏矩阵较为适合，但有时收敛速度很慢。流程图：实验地点综合楼五层5。6室指导教师王峥太原理工大学学生实验报告学院名称计算机科学与技术专业班级计科1002学号 14学生姓名xxx实验日期成绩课程名称计算数值方法实验题目实验五代数插值一、课题名称(使用拉格朗日插值法或牛顿插值法

17、求解：已知f(x)在6个点的函数值如下表所示，运用插值方法，求f的近似值。xf(x)、目的和意义学习使用拉格朗日插值法或牛顿插值法求解三、计算公式设函数在区间a,b上n+1互异节点Xo,X i,Xn上的函数值分别为yo,y i,y n,求n次插值多项式Pn(x),满足条件Pn(Xj)=yj, j=0,1,n令Ln(x)=y ol o(x)+y il i(x)+ +yn(x)= 刀 yd i(x)其中l o(x),l l(x),l n(x)为以xo,x 1,x n为节点的n次插值基函，n数，则Ln(x)是一次数不超过n的多项式，且满足Ln(Xj)=yj, L=0,1,再由插值多项式的唯一性，得P

18、n(X)三 Ln(X)四、主要仪器设备Vc+ C-free CodeBlocks五、结构程序设计#in clude <>#in clude <>mai n()int I;char L;double M100100;double x100,y100;double X=1,xx=0,w=1,N=0,P,R=1;int n=5;cout<<"请输入所求均差阶数："II求所有阶差cin»n;IIfor(i nt i=0;i<=n ;i+)III*cout<<"请输入 x"<<i<&

19、lt;"的值："<<endl;cin> >xi;cout<<"请输入 y"<<i<<"的值："<<endl;cin> >yi;MiO=xi;Mi1=yi; */用二维保存所有数据M00=;M01=;M10=;M11=;M20=;M21=;M30=;M31=;M40=;M41=;M50=;M51=;for( int j=2;j<=n+1;j+)for(i nt i=1;i<=n ;i+)Mij=(Mij-1-Mi-1j-1)/(Mi0-Mi-

20、j+10);for(i nt i=1;i<=n ;i+)cout<<"其"<<i<<"阶均差为："<<Mii+1<<endl;cout<<"请输入x的值：x="cin»xx;for(i nt i=0;i< n; i+)X*=xx-Mi0;N+=Mi+1i+2*X;P=M01+N;cout<<"其函数值：y="<<p<<endl; / 六、结果讨论和分析耳邛fi均差为谕丼童为;6觸I勢0苴

21、齒;_ 请按任意犍继统H2 甘3 -d116 翻 197333 012381 0084981?0y-0.631717分析讨论拉格朗日插值的优点是插值多项式特别容易建立，缺点是增加节点是原有多 项式不能利用，必须重新建立，即所有基函数都要重新计算，这就造成计算量 的浪费。所以该程序 选择了牛顿法，由于在输入数据的时候比较麻烦所以将数据存在了程序里面，只能在程序里修改数据。另外由于精度问题计算结果有一定的误差。流程图：实验地点综合楼五层506室指导教师王峥太原理工大学学生实验报告学院名称计算机科学与技术专业班级计科1002学号 14学生姓名xxx实验日期成绩课程名称计算数值方法实验题目实验六最小二

22、乘法拟合多项式、课题名称给定数据点（Xi , yi），用最小二乘法拟合数据的多项式，并求平方误差。Xi0yi1、目的和意义1 熟练运用已学计算方法求解方程组2 加深对计算方法技巧，选择正确的计算方法来求解各种方程组3 培养使用电子计算机进行科学计算和解决问题的能力三、计算公式建立正规方程组：i+kjE(E Xi) ak=EXi yi , j=0,1,n平方误差：-k2I =E(E akXi -y i)四、主要仪器设备Vc+ C-free CodeBlocks五、结构程序设计源代码：#in clude<>#in clude<>#define N 15double powe

23、r(double &a,i nt n)double b=1;for(i nt i=0;i <n ;i+)b*=a;return b;void Gauss();double XN,YN,sumXN,sumYN,aNN,bN,lNN,xN;void main()ofstream outdata;ifstream in data;double s;int i,j,k,n,index;cout<<"请输入已知点的个数n="cin»n;n=7;cout<<e ndl;输入给定数据循环输入任何数据cout<<"请输入

24、 X 和 Y:"<<endl;/X0=;Y0=;X1=;Y1=;X2=;Y2=;X3=;Y3=;X4=;Y4=;X5=;Y5=;X6=;Y6=; II绑定数据II可以解绑由下forfor(i=0;i <n ;i+)cout<<"X"<<i<<"="I/cin >>Xi;sumX1+=Xi;cout<<"Y"<<i<<"="cin >>Yi;sumY1+=Yi;cout<<e ndl

25、;cout<<"sumX1="<<sumX1<<"t"<<"sumY1="<<sumY1<<e ndl;cout<<"请输入拟合次数index="cin»in dex;cout<<e ndl;i=n;sumX0=i;for(i=2;i<=2* in dex;i+)sumXi=0;for(j=0;j< n;j+)sumXi+=power(Xj,i);cout<<"sumX&quo

26、t;<<i<<"="<<sumXi<<e ndl;for(i=2;i<=i ndex+1;i+)sumYi=0;for(j=0;j< n;j+)sumYi+=power(Xj,i-1)*Yj;cout<<"sumY"<<i<<"="<<sumYi<<e ndl;for(i=1;i<=i ndex+1;i+)/建立正规方程组for(j=1;j<=i ndex+1;j+)aij=sumXi+j-2;bi=sumYi;k=1;/用高斯消元法解方程组dofor(j=k+1;j<=i ndex+1;j+) ljk=ajk/akk;for(i=k+1;i<=i ndex+1;i+)for(j=k+1;j<=i ndex+1;j+)aij=aij-lik*akj;bi=bi-lik*bk;if(k=in dex+1) break;k+;while(1);xin dex+1=bi ndex+1/ai ndex+1 in dex+1;for(i=i ndex;i>=1;i-)s=0;for(j=i+1;j<=i ndex+1;j+