昆明理工数值分析大作业线性方程组的迭代法

I工科研究生数值分析上机实验实验报告题目：线性方程组的迭代法学院：专业：年级：201级学生姓名：学号：任课教师：1一、实验目的:1、通过上机计算体会迭代法求解线性方程组的特点，并能和消去法比较；2、运用所学的迭代法算法，解决各类线性方程组，编出算法程序；3、体会上机计算时，终止步骤||||)()1(kkxx或k大于所给的迭代次数，对迭代法敛散性的意义；4、体会初始解x0,松弛因子的选取，对计算结果的影响。二、问题简述：分别选用Jacobi迭代法，Gauss-Seidol迭代法和SOR方法计算下列线性方程组的解。方程组一:2方程组二(设对称正定阵系数阵线方程组)：方程组三(三对角形线性方程组)：3三、本题目用到的三种基本解线性方程组的迭代法简述：设有线性方程组：Ax=bx为n维向量A为行列标分别为i,j的元素aij素构成的n阶系数矩阵Jacobi迭代法：将A分解为三部分：A=nnaaa2211-00001,22,111,121nnnnnnaaaaaa-0000,1211,21,112nnnnnnaaaaaa=D-L-U.Jacobi迭代法之迭代公式如下：x(k+1)=D-1(L+U)x(k)+D-1bGauss-Seidol迭代法迭代公式:x(k+1)=(D-L)-1Ux(k)+(D-L)-1bSOR方法迭代公式:()x(k+1)=(D-ωL)-1((1-ω)D+ωU)x(k)+(D-ωL)-1ωb其中ω为SOR方法权因子。根据实验要求设计计算流程图：4Jacobi迭代法：否是否是带入初次迭代向量x0.B=D-1（L+U）将系数矩阵A三角分解为L，U和D，即A=D-L-U设置迭代值x（k-1），x（k）设置初始迭代向量x0量选取精确度ε。定义迭代次数k=1读取b得到下一次迭代的值：x（k）||x（k）-x（k-1）||ε求范数||x（k）-x（k-1）||;记录迭代次数k=k+1;带入上一次次迭代向量x（k-1）.输出最终计算值x*=x（k）;以及终迭代次数k;终止；迭代系数矩阵谱半径小于1x（k）=D-1（L+U）x（k-1）+D-1b5Gauss-Seidol迭代法：否是否是带入初次迭代向量x0.B=(D-L)-1U将系数矩阵A三角分解为L，U和D，即A=D-L-U设置迭代值x（k-1），x（k）设置初始迭代向量x0量选取精确度ε。定义迭代次数k=1读取b得到下一次迭代的值：x（k）||x（k）-x（k-1）||ε求范数||x（k）-x（k-1）||;记录迭代次数k=k+1;带入上一次次迭代向量x（k-1）.输出最终计算值x*=x（k）;以及终迭代次数k;终止；迭代系数矩阵谱半径小于1x（k）=(D-L)-1Ux(k-1)+(D-L)-1b6SOR方法：否是否是带入初次迭代向量x0.B=(D-ωL)-1((1-ω)D+ωU)将系数矩阵A三角分解为L，U和D，即A=D-L-U设置迭代值x（k-1），x（k）设置初始迭代向量x0量选取精确度及权因子ε、ω。定义迭代次数k=1读取b得到下一次迭代的值：x（k）||x（k）-x（k-1）||ε求范数||x（k）-x（k-1）||;记录迭代次数k=k+1;带入上一次次迭代向量x（k-1）.输出最终计算值x*=x（k）;以及终迭代次数k;终止；迭代系数矩阵谱半径小于1x（k）=(D-ωL)-1((1-ω)D+ωU)x(k-1)+(D-ωL)-1ωb7五、实验内容：本次实验共使用三种迭代方法采用不同精确度和相关因子对三种不同类型的线性方程组进行了上机编程计算。除了在不同精度下采用三种迭代方法进行计算不同方法同一方程的迭代次数不同，不同精度同一方程的迭代次数也不同。其中SOR方法因权因子选取不同而迭代次数不同，本实验采用同一精确值分别对三个方程对应选取不同的三组权因子进行了计算。并且采用二分法逼近的方法寻找了相对最佳权因子。六、计算结果：通过编程和迭代得到以下计算结果Jacobi迭代法：方程一：方程一用Jacobi迭代法计算时，迭代系数矩阵谱半径大于1，故迭代法发散。方程二：方程三：分析：Gauss-Seidol迭代法：方程一：方程一用Gauss-Seidol迭代法：算时，迭代系数矩阵谱半径大于1，故迭代法发散。方程二：方程三：分析：8SOR方法：方程一：方程一用SOR迭代法计算时,以0.001，从(0,2)选取权因子ω,迭代系数矩阵谱半径均大于1，选权因子为0，)B(ρ=1，且)B(ρ随着ω增大程增大趋势，故判断迭代法发散。从而方程一用此三种方法均无法计算！！！！方程二：9方程三：选取迭代向量初值：w=0.8w=0.9w=1.0e=e-1e=e-2e=e-3e=e-1e=e-2e=e-3e=e-1e=e-2x12.0082.01842.00182.00152.0122.00122.00091.9995x20.88321.02871.00181.00131.01251.00071.00021.0063x3-3.087-2.9782-2.9996-3-2.9936-3.0008-3.001-3.0052x40.04990.0047-0.0016-0.0016-0.006-0.0021-0.0019-0.0076x51.06340.98020.99710.99750.98320.99760.99820.9909x6-1.9766-2.0371-2.003-2.0024-2.0175-2.0017-2.0011-2.0027x73.00682.96722.99792.99842.99032.99922.99953.0002x8-0.0182-0.0151-0.0011-0.0008-0.0023-0.0003-0.00010.0005x90.95240.99830.99960.99971.00030.999911.0003x10-1.0335-0.9989-1.0001-1.0001-0.9997-1-1-0.9999迭代次数35825726w=1.0w=1.1w=1.2e=e-3e=e-1e=e-2e=e-3e=e-1e=e-2e=e-3x11.99981.99962.00431.99791.99741.99251.9989x20.99940.99940.99411.00061.0011.00651.0002x3-3.0008-3.0007-2.998-3.0004-3.0005-3.0036-3.0004x4-0.0006-0.0003-0.00140.00060.00080.00210.0002x50.999811.00230.99980.99960.99760.9999x6-2-2-2.0006-2-2-1.9989-2x7332.9999333.00013x800000-0.00010x91111111x10-1-1-1-1-1-1-1迭代次数7267279分析：10七、结果分析11八、实验结论12附件：三种迭代法的迭代程序：Jacobi迭代法：A=A;b=ba;x0=x0a;%带已知数据(a:aready)D=diag(diag(A));%求A的对角矩阵L=-tril(A,-1);%求A的下三角阵U=-triu(A,1);%求A的上三角阵B=D\(L+U);f=D\b;x=B*x0+f;n=[];%迭代次数ifvrho(B)1k=[10e-110e-210e-3];foreps=k;i=1;n(i)=1;whilenorm(x-x0)=epsx0=x;x=B*x0+f;n(i)=n+1;endendelsen=0;endGauss-Seidol迭代法：A=A;b=ba;x0=x0a;%带已知数据(a:aready)D=diag(diag(A));%求A的对角矩阵L=-tril(A,-1);%求A的下三角阵U=-triu(A,1);%求A的上三角阵B=(D-L)\U;f=(D-L)\b;x=B*x0+f;n=[];%迭代次数ifvrho(B)1k=[10e-110e-210e-3];foreps=k;i=1;n(i)=1;whilenorm(x-x0)=epsx0=x;x=B*x0+f;n(i)=n+1;endendelsen1=0;end2SOR方法：A=Amn;b=ba;x0=x0a;%带已知数据(a:aleady)D=diag(diag(A));%求A的对角矩阵L=-tril(A,-1);%求A的下三角阵U=-triu(A,1);%求A的上三角阵w=0;wm=(0.8:0.1:1.2)；x=zeros([m15]);i=1;j=1;n=zeros([53]);c=zeros([215]);forw=wmi=i+1;B=(D-w*L)\((1-w)*D+w*U);f=(D-w*L)\(w*b);x(:,j)=B*x0+f;ifvrho(B)1k=[10e-110e-210e-3];foreps=k;whilenorm(x(:,j)-x0)=epsx0=x(:,j);x(:,j)=B*x0+f;n(i,j)=n1+1;endc(2,j)=w;j=j+1;endelseifvrho(B)=1x(:,j)=0;c(1,j)=vrho(B);j=j+3;endendEndclearx=-1*[1,1,1,1,1,1,1,1,1,];c=4*ones(1,10);A=diag(c)+diag(x,1)+diag(x,-1);b=[75-1326-1214-45-5]';x0=[21-301-2301-1]';3D=diag(diag(A));%求A的对角矩阵L=-tril(A,-1);%求A的下三角阵U=-triu(A,1);%求A的上三角阵w=0;wm=(0.8:0.1:1.2);x=zeros([1015]);i=1;j=1;n=zeros([53]);c=zeros([215]);x0=[21-301-2301-1]';forw=wmB=(D-w*L)\((1-w)*D+w*U);f=(D-w*L)\(w*b);x(:,j)=B*x0+f;ifvrho(B)1k=[10e-110e-210e-3];foreps=k;whilenorm(x(:,j)-x0)=epsx0=x(:,j);x(:,j)=B*x0+f;endn(i,j)=n(i,j)+1;c(2,j)=w;j=j+1;endelseifvrho(B)=1x(:,j)=0;c(1,j)=vrho(B);j=j+3;endi=i+1;endend4