数值分析大作业

数值分析大作业学号：*********专业：机械工程学生姓名：***2014年10月摘要：在自然科学与工程技术中，很多问题的解决常常归结为求解线性方程组Ax=b。随着计算机的发展，利用计算机这个强有力的计算工具去求解线性方程组是一个非常实用的问题。在求解大型线性方程组时，直接法在多次消元，回代的过程中，四则运算的误差累计与传播无法控制，致使计算结果的精度就无法保证，特别是求解大型稀松矩阵时，还要对系数矩阵进行分解。而迭代法相对于直接法而言，具有保持迭代矩阵不变的特点，计算程序一般也比较简单，且对于许多问题收敛速度比较快。比较常用的迭代法有雅克比迭代法、高斯一塞德尔迭代法和逐次超松弛迭代法等，本次研究目的是通过求解一个线性方程组来比较它们的迭代效果，验证一些已有的结论。1．数学原理1.1雅可比迭代法将线性方程组的系数矩分解为A=D+L+U，其中D是由A的主对角元素构成的对角矩阵，L是由A的严格下三角部分构成的严格下三角矩阵，U是由A的严格上三角部分构成的严格上三角矩阵，即若系数矩阵A的对角元素，则矩阵D非奇异，取M=D，N=-（L+U)，则，因而，构造的迭代法为：1.2高斯-赛得尔迭代法将线性方程组的系数矩分解为A=D+L+U。若系数矩阵A的对角元素不等于0，则矩阵D非奇异，取M=L+D，N=-U，则因而，构造的迭代法为：1.3逐次超松弛迭代法线性方程组的系数矩分解为A=D+L+U。取则因而，构造的迭代法为：2.程序设计2.1雅可比迭代法求解雅可比迭代法MATLAB程序如下：%majacobi.mfunctionx=majacobi(A,b,x0,ep,N)%用途：用Jacobi迭代法解线性方程组Ax=b%格式;x=majacobi(A,b,x0,ep,N)A为系数矩阵，b为右端向量，%x0为初始向量（默认零向量），ep为精度（默认1e-6),%N为最大迭代次数（默认500次），x返回近似解向量n=length(b);ifnargin<5,N=500;endifnargin<4,ep=1e-6;endifnargin<3,x0=zeros(n,1);endx=zeros(n,1);k=0;whilek<Nfori=1:nx(i)=(b(i)-A(i,[1:i-1,i+1:n])*x0([1:i-1,i+1:n]))/A(i,i);endifnorm(x-x0,inf)<ep,break;endx0=x;k=k+1;endifk==N,Warning('100');enddisp(['k=',num2str(k)])在MATLAB命令窗口执行创建的m文件>>editmajacobi>>A=[10.9,1.2,2.1,0.9;1.2,11.2,1.5,2.5;2.1,1.5,9.8,1.3;0.9,2.5,1.3,12.3]A=10.90001.20002.10000.90001.200011.20001.50002.50002.10001.50009.80001.30000.90002.50001.300012.3000>>b=[-7.0,5.3,10.3,24.6]';>>x=majacobi(A,b)k=17x=-0.99860.00711.00321.96562.2高斯-赛得尔迭代法求解高斯-赛得尔迭代法MATLAB程序如下：%maseidel.mfunctionx=maseidel(A,b,x0,ep,N)%用途：用Gauss-Seidel迭代法解线性方程组Ax=b%格式：x=maseidel（A,b,x0,ep,N)A为系数矩阵，b为右端向量，%x0为初始向量（默认零向量），ep为精度（默认1e-6),%N为最大迭代次数（默认500次），x返回近似解向量n=length(b);ifnargin<5,N=500;endifnargin<4,ep=1e-6;endifnargin<3,x0=zeros(n,1);endx=zeros(n,1);k=0;whilek<Nfori=1:nifi==1x(1)=(b(1)-A(1,2:n)*x0(2:n))/A(1,1);elseifi==nx(n)=(b(n)-A(n,1:n-1)*x(1:n-1))/A(n,n);elsex(i)=(b(i)-A(i,1:i-1)*x(1:i-1)-A(i,i+1:n)*x0(i+1:n))/A(i,i);endendendifnorm(x-x0,inf)<ep,break;endx0=x;k=k+1;endifk==N,Warning('100');enddisp(['k=',num2str(k)])在MATLAB命令窗口执行创建的m文件：>>editmaseidel>>A=[10.9,1.2,2.1,0.9;1.2,11.2,1.5,2.5;2.1,1.5,9.8,1.3;0.9,2.5,1.3,12.3]A=10.90001.20002.10000.90001.200011.20001.50002.50002.10001.50009.80001.30000.90002.50001.300012.3000>>b=[-7.0,5.3,10.3,24.6]';>>x=maseidel(A,b)k=7x=-0.99860.00711.00321.96562.3逐次超松弛迭代法求解逐次超松弛迭代法MATLAB程序如下：%masor.mfunctionx=masor(A,b,omega,x0,ep,N)%用途：用SOR迭代法解线性方程组Ax=b%格式：x=maseidel（A,b,x0,ep,N)A为系数矩阵，b为右端向量，%omega为松弛因子（默认1.5),x0为初始向量（默认零向量），ep为精度（默认1e-6),%N为最大迭代次数（默认500次），x返回近似解向量n=length(b);ifnargin<6,N=500;endifnargin<5,ep=1e-6;endifnargin<4,x0=zeros(n,1);endifnargin<3,omega=1.5;endx=zeros(n,1);k=0;whilek<Nfori=1:nifi==1x1(1)=(b(1)-A(1,2:n)*x0(2:n))/A(1,1);elseifi==nx1(n)=(b(n)-A(n,1:n-1)*x(1:n-1))/A(n,n);elsex1(i)=(b(i)-A(i,1:i-1)*x(1:i-1)-A(i,i+1:n)*x0(i+1:n))/A(i,i);endendx(i)=(1-omega)*x0(i)+omega*x1(i);endifnorm(x0-x,inf)<ep,break;endk=k+1;x0=x;endifk==N,Warning('100');enddisp(['k=',num2str(k)])取不同的松弛因子进行试验，在MATLAB命令窗口执行创建的m文件：>>editmasor>>A=[10.9,1.2,2.1,0.9;1.2,11.2,1.5,2.5;2.1,1.5,9.8,1.3;0.9,2.5,1.3,12.3]A=10.90001.20002.10000.90001.200011.20001.50002.50002.10001.50009.80001.30000.90002.50001.300012.3000>>b=[-7.0,5.3,10.3,24.6]';>>x=masor(A,b,0.4)k=33x=-0.99860.00711.00321.9656取其它的松弛因子得：取松弛因子为0.8>>x=masor(A,b,0.8)得到k=12；x=-0.9986，0.0071，1.0032，1.9656。取松弛因子为1.2>>x=masor(A,b,1.2)得到k=11；x=-0.9986，0.0071，1.0032，1.9656。取松弛因子为1.6>>x=masor(A,b,1.6)得到k=33；x=-0.9986，0.0071，1.0032，1.9656。三．结果分析和讨论由程序解得的结果可知：雅可比迭代法在第17次迭代达到收敛，高斯-塞德尔迭代法在第7次迭代达到收敛。从此例可以看出，高斯-塞德尔迭代法比雅可比迭代法收敛速度快(达到同样的精度所需迭代次数少)。由雅可比迭代公式可知，在迭代的每一步计算过程中使用的是旧的分量来计算，没有用到迭代出的新的分量。从直观上看,最新计算出的分量可能比旧的分量接近真实值，要更好一些。而高斯-赛德尔迭代法就在计算过程中加入了新值的因素，在每次迭代时将前面所得到的新值迅速地被利用到计算下一个值上，从而可以减少计算的次数，节约时间。由本次计算和迭代公式可以得出，高斯-赛德尔迭代法在收敛速度上要优于雅可比迭代法。而迭代序列的收敛取决于迭代矩阵的谱半径，而与初值向量的选取无关。因为在求解同一个方程组时，雅可比迭代法和高斯-赛德尔迭代法的谱半径不一定相同，而且并无包含关系，因此，有时雅可比迭代法收敛，而高斯-赛德尔不收敛。当然，也有雅可比迭代法不收敛，而高斯-赛德尔迭代法收敛的情形。从超松弛迭代法的结果来看，松弛因子分别取0.4,0.8,1.2,1.6时，达到最终精度