西南交通大学数值分析上机实习

PAGEPAGE10序言数值分析是计算数学的一个主要部分，它不仅要研究各种数学问题的数值解法，同时也要分析所用的数值解法在理论上的合理性，如解法所产生的误差能否满足精度要求：解法是否稳定、是否收敛及熟练的速度等。由于求解的迭代过程很复杂，计算量很大，所以需要借助其他辅助手段来求解。本文中采用了软件C++来编写求解程序，C++语言简洁紧凑，使用方便灵活，其程序书写形式自由，压缩了一切不必要的编程工作并且控制语句非常丰富（如for循环，while循环，break语句和if语句），对编写程序提供了很大的便利。在本文中使用C++编写了牛顿法、牛顿-Steffensen法方程求解的程序和雅格比法、高斯－赛德尔迭代法求解方程组的程序及Runge-Kutta4阶算法，并通过实例求解验证了其可行性，比较了求解同一种问题时不同方法之间的优缺性，其中包含解的精确度和解的收敛速度两个重要指标。选用该软件编写主要原因是大学时期学习了C++程序设计，因此对于改软件相对而言要熟悉些。目录第一章牛顿法和牛顿-Steffensen法迭代求解的比较 31.1计算结果 31.2结果分析 4第二章Jacobi和Causs-Seidel迭代法比较 42.1计算结果 42.2结果分析 4第三章观察Runge-Kutta4阶算法稳定区间的作用 53.1计算结果 63.2结果分析 6总结 6附录 7第一章牛顿法和牛顿-Steffensen法迭代求解的比较分别用牛顿法，及基于牛顿算法下的Steffensen加速法求ln(x+sinx)=0的根。初值x0分别取0.1,1,1.5,2,4进行计算。求sinx=0的根。初值x0分别取1，1.4，1.6,1.8，3进行计算。分析其中遇到的现象与问题。1.1计算结果(1)方程可变形为x+注：x0为初始值，x为收敛解，N为迭代次数牛顿法x00.111.524x0.5109730.5109730.5109730.5109730.510973N445640牛顿—Steffensen加速法x0.5109730.510973无收敛解0.5109730.510973N3348(2)牛顿法x011.41.61.83x03.1415931.41596.283193.14159N57843牛顿—Steffensen加速法x无收敛解-3.1415925.13276.28319无收敛解N4631.2结果分析从结果对比我们可发现牛顿—Steffensen加速法比牛顿法要收敛的快，牛顿法对于初值的选取特别重要，比如第（1）问中的初值为4的情况，迭代次数算了40次，远大于其余初值的情况；在第（2）问中的初值为1.6的情况，收敛解得31.4159，分析其原因应该是，x0=1.6,；在牛顿—Steffensen加速法第(1)问中x=1.5、第二问x=0和3的情况，无收敛解。第二章Jacobi和Causs-Seidel迭代法比较用雅格比法与高斯－赛德尔迭代法解下列方程组Ax=b，研究其收敛性，上机验证理论分析是否正确，比较它们的收敛速度，观察右端项对迭代收敛有无影响。(1)A行分别为A1=[6,2,-1],A2=[1,4,-2],A3=[-3,1,4]；b1=[-3,2,4]T,b2=[100,-200,345]T,(2)A行分别为A1=[1,0,8,0.8],A2=[0.8,1,0.8],A3=[0.8,0.8,1]；b1=[3,2,1]T,b2=[5,0,-10]T,(3)A行分别为A1=[1,3],A2=[-7,1]；b=[4,6]T,2.1计算结果（1）x0=[0,0,0]T为初始值,N为迭代次数，0.00001为误差精度,X为收敛解Jacobi迭代Causs-Seidel迭代Ax=b1X[-0.727271,0.808083,0.252524][-0.727275,0.808082,0.252523]N2012Ax=b2X[36.3636,-2.0707,114.04][36.3636,-2.07071,114.04]N2617（2）Jacobi迭代Causs-Seidel迭代Ax=b1X无收敛解[5.76922,0.76922,-4.23075]N37Ax=b2X无收敛解[36.6923,7.69229,-42.3077]N44（3）Jacobi迭代Causs-Seidel迭代Ax=bX无收敛解无收敛解N2.2结果分析（1）第一小题的雅可比迭代矩阵为0-0.33330.1667-0.250000.50000.7500-0.25000其特征值分别为-0.5427，0.2713+0.3708i，0.2713-0.3708i，经计算谱半径为小于1,故方程组雅可比迭代收敛。而高斯－赛德尔迭代矩阵为0-0.33330.166700.08330.45830-0.27080.0104其特征值为0，0.0469+0.3504i，0.0469-0.3504i经计算谱半径为0.3535小于1,故原方程组高斯－赛德尔迭代矩阵收敛。第二小题的雅可比迭代矩阵为0-0.8000-0.8000-0.80000-0.8000-0.8000-0.80000征值分别为-1.6，0.8，0.8，经计算谱半径为1.6大于1,故方程组雅可比迭代发散。而高斯－赛德尔迭代矩阵为0-0.8000-0.800000.6400-0.160000.12800.7680其特征值为0，0.7040+0.1280i，0.704280i经计算谱半径为0.7155小于1,故原方程组高斯－赛德尔迭代矩阵收敛。第三小题的雅可比迭代矩阵为0-370其特征值分别为4.5826i，－4.5826i，谱半径为4.5826大于1,故方程组雅可比迭代发散。而高斯－赛德尔迭代矩阵为0-30-21其特征值为0，－21，谱半径为21大于1,故原方程组高斯－赛德尔迭代矩阵发散。（2）从2.1中的结果列表中可以看到，Causs-Seidel迭代法比Jacobi迭代法收敛速度要快，并且方程组右端项对迭代收敛是无影响的。第三章观察Runge-Kutta4阶算法稳定区间的作用用Runge-Kutta4阶算法对初值问题y/=-20*y，y(0)=1按不同步长求解，用于观察稳定区间的作用，推荐两种步长h=0.1,0.2。注：此方程的精确解为：3.1计算结果=1.0h=0.1h=0.21.6935131251.693331253.2结果分析,h=0.2时，h=-4,而Runge-Kutta4阶算法的绝对稳定区间是[-2.78,0],故h=0.2时计算不稳定；而h=0.1时，h=-2，在绝对稳定区间内，计算稳定，结果可靠。总结通过这次上机练习，让我对数值分析所介绍的迭代求解方法及其理论有了更深层次的理解，了解了各种方法之间的优缺点，并且认识到了自己在以前的学习中所存在的问题，及时的修补了自己知识上的漏洞。同时也提高了我在编写程序上的熟练程度，所以，我认为这次上机实习是非常有收获的，给予我学习数值分析的帮助也是非常大的。附录1.1第(1)问牛顿法#include<iostream>#include<cmath>usingnamespacestd;intmain(){ doublex,x0; inti=0; cout<<"请输入初值:"; cin>>x0; do{x=x0; x0=x0-(x0+sin(x0)-1)/(1+cos(x0)); i++;}while(fabs(x-x0)>=0.000001);cout<<"x="<<x0<<endl;cout<<"N="<<i<<endl;return0;}1.2第(1)问牛顿—Steffensen加速法#include<iostream>#include<cmath>usingnamespacestd;intmain(){ doublex,x0,x1,x2; inti=0; cout<<"请输入初值:"; cin>>x0; do{ x=x0;x1=x0-(sin(x0)+x0-1)/(cos(x0)+1); x2=x1-(sin(x1)+x1-1)/(cos(x1)+1); x0=x0-(x1-x0)*(x1-x0)/(x2-2*x1+x0); i++;}while(fabs(x-x0)>=0.000001);cout<<"x="<<x0<<endl;cout<<"N="<<i<<endl;return0;}1.3第(2)问牛顿法#include<iostream>#include<cmath>usingnamespacestd;intmain(){ doublex,x0; inti=0; cout<<"请输入初值:"; cin>>x0; do{x=x0; x0=x0-sin(x0)/cos(x0); i++;}while(fabs(x-x0)>=0.000001);cout<<"x="<<x0<<endl;cout<<"N="<<i<<endl;return0;}1.4第(2)问牛顿—Steffensen加速法#include<iostream>#include<cmath>usingnamespacestd;intmain(){ doublex,x0,x1,x2; inti=0; cout<<"请输入初值:"; cin>>x0; do{ x=x0;x1=x0-sin(x0)/cos(x0); x2=x1-sin(x1)/cos(x1); x0=x0-(x1-x0)*(x1-x0)/(x2-2*x1+x0); i++;}while(fabs(x-x0)>=0.000001);cout<<"x="<<x0<<endl;cout<<"N="<<i<<endl;return0;}2.1雅格比法#include<iostream>#include<cmath>usingnamespacestd;intmain(){ doublea[10][10],b[10],sum=0.0,x[10],y[10],c,s=0.0; inti,k=0,n,j; cout<<"请输入维数:"; cin>>n; cout<<"请输入数组A:"; for(i=0;i<n;i++) for(j=0;j<n;j++) cin>>a[i][j];cout<<"请输入数组B:"; for(i=0;i<n;i++) cin>>b[i];cout<<"请输入初始X:"; for(i=0;i<n;i++) cin>>x[i]; do{for(i=0;i<n;i++) {for(j=0;j<i;j++) sum=sum+a[i][j]*x[j];for(j=i+1;j<n;j++) sum=sum+a[i][j]*x[j]; y[i]=(b[i]-sum)/a[i][i]; sum=0; } c=0; for(i=0;i<n;i++) if(c<fabs(x[i]-y[i])) c=fabs(x[i]-y[i]); if(c>0.00001) c=0; elses=1;for(i=0;i<n;i++) x[i]=y[i];k++;}while(s!=1);for(i=0;i<n;i++) cout<<"解为："<<x[i]<<endl;cout<<"N="<<k<<endl;return0;}2.2高斯－赛德尔迭代法#include<iostream>#include<cmath>usingnamespacestd;intmain(){ doublea[10][10],b[10],sum=0.0,x[10],y[10],c,s=0.0,z[10]; inti,k=0,n,j; cout<<"请输入维数:"; cin>>n; cout<<"请输入数组A:"; for(i=0;i<n;i++) for(j=0;j<n;j++) cin>>a[i][j];cout<<"请输入数组B:"; for(i=0;i<n;i++) cin>>b[i];cout<<"请输入初始X:"; for(i=0;i<n;i++) cin>>x[i]; do{for(i=1;i<n;i++) z[i]=x[i]; for(i=0;i<n;i++) {for(j=0;j<i;j++) sum=sum+a[i][j]*x[j];for(j=i+1;j<n;j++) sum=sum+a[i][j]*x[j]; y[i]=(b[i]-sum)/a[i][i]; sum=0; x[i]=y[i]; } c=0; for(i=0;i<n;i++) if(c<fabs(z[i]-y