1、实验报告四线性方程组的求解迭代浙江大学城市学院实验报告课程名称 科学计算 实验项目名称 线性方程组的求解迭代法 实验成绩 指导老师(签名 ) 日期 2012-4-6 一. 实验目的和要求1 掌握Jacobi迭代方法.Gauss-Seidel迭代方法.SOR迭代方法的编程思想.能够分别用分量形式和矩阵形式编写相关程序。2 观察SOR迭代法中松弛因子变化情况对收敛的影响。 3 了解Hilbert矩阵的病态性和作为线性方程组系数矩阵的收敛性。二. 实验内容和原理编程题2-1要求写出Matlab源程序(m文件).并有适当的注释语句;分析应用题2-2.2-3.2-4要求将问题的分析过程、Matlab源程
2、序和运行结果和结果的解释、算法的分析写在实验报告上。2-1 编程注释设 对下述求解线性方程组的Matlab程序添上注释语句.其中和分别为线性方程组的系数矩阵和右端向量;为迭代初始向量;为容许迭代最大次数.为迭代终止条件的精度(容许误差).终止条件为前后两次迭代解的差的向量2-范数。1) Jacobi迭代: 2) Gauss-Seidel迭代:GaussSeidelmethod(A,b,x0,Nmax,eps)3) 松弛迭代:SORmethod(A,b,x0,Nmax,eps,w) 2-2 分析应用题利用2-1中的程序来分析用下列迭代法解线性方程组:的收敛性.并求出使的近似解及相应的迭代次数.其
3、中取迭代初始向量为零向量。1)Jacobi迭代法;2)Gauss-Seidel迭代法;3)松弛迭代法(松弛因子依次取1.334.1.95.0.95)。2-3 分析应用题考虑方程组.其中系数矩阵为Hilbert矩阵. 选择问题的维数分别为2、3、5、10.并通过首先给定解再定出右端的办法确定问题.解的给定可以使用函数定义.并取迭代初始向量为零向量.迭代误差为.编写程序: 其中n为Hilbert矩阵的维数.分别构造求解该问题的Jacobi迭代和Gauss-Seidel迭代.看它们是否收敛。2-4 分析应用题解线性方程组的解向量.取.其中为任一非零的六元向量;编写程序输出结果: 认真观察之.能发现什
4、么有趣的现象?【MATLAB相关函数】 提取(产生)对角阵v=diag(x) 若输入向量x.则输出v是以x为对角元素的对角阵;若输入矩阵x.则输出v是x的对角元素构成的向量v=diag(diag(x) 输入矩阵x.输出v是x的对角元素构成的对角阵.可用于迭代法中从A中提取D。 提取(产生)上(下)三角阵v=triu(x) 输入矩阵x.输出v是x的上三角阵;v=tril(x) 输入矩阵x.输出v是x的下三角阵;v=triu(x,1) 输入矩阵x.输出v是x的上三角阵.但对角元素为0.可用于迭代法中从A中提取U。v=tril(x,-1) 输入矩阵x.输出v是x的下三角阵.但对角元素为0.可用于迭代
5、法中从A中提取L。 矩阵特征值b=eig(A) 输入矩阵A.输出b是A的所有特征值。 范数n=norm(x) 输入x为向量或矩阵.输出为x的2范数;n=norm(x,p) 输入x为向量或矩阵.当p=1, inf时分别输出为x的1.无穷范数; Hilbert矩阵h=hilb(n) 输出h为n阶Hilbert矩阵三. 操作方法与实验步骤(包括实验数据记录和处理)2-1 编程注释设 对下述求解线性方程组的Matlab程序添上注释语句.其中和分别为线性方程组的系数矩阵和右端向量;为迭代初始向量;为容许迭代最大次数.为迭代终止条件的精度(容许误差).终止条件为前后两次迭代解的差的向量2-范数。1)Jac
6、obi迭代: function X = Jacobimethod( A,b,X0,P,Nmax,eps)%A和b分别为线性方程组的系数矩阵和右端向量;x0为迭代初始向量X0;Nmax为容许迭代最大%次数.eps为迭代终止条件的精度(容许误差).终止条件为前后两次迭代解的差的向量2-范数%P=1.2.inf或fro.%输出的量:系数矩阵的的值和有关雅可比迭代收敛性的相关信息及AX=b的精确解jX和近似解X n m=size(A);for j=1:m a(j)=sum(abs(A(:,j)-2*(abs(A(j,j);endfor i=1:n if a(i)=0 disp(请注意:系数矩阵A不是严
7、格对角占优的.此雅可比迭代不一定收敛) return endendif a(i)0 disp(请注意:系数矩阵A是严格对角占优的.此方程组由唯一解.且雅可比迭代收敛)endfor k=1:Nmax k for j=1:m X(j)=(b(j)-A(j,1:j-1,j+1:m)*X0(1:j-1,j+1:m)/A(j,j); end X,djwcX=norm(X-X0,P); xdwcX=djwcX/(norm(X,P)+eps); X0=X; X1=Ab; if (djwcXeps)&(xdwcXeps)&(xdwcXeps) disp(请注意:雅可比迭代次数已经超过最大迭代次数Nmax)en
8、da,X=X;jX=X1,end2)Gauss-Seidel迭代:GaussSeidelmethod(A,b,x0,Nmax,eps)function x = GaussSeidelmethod( A,b,x0,P,Nmax,wucha )% 输入的量:线性方程组AX=b的系数矩阵A和b.初始向量x0.范数的名称P=1.2.inf或fro.% 近似解x的误差(精度)wucha和迭代的最大次数Nmax。 % 输出的量:以系数矩阵的对角元构成对角矩阵D、A的上三角形矩阵U.但对角% 元为0、A的下三角矩阵L.但对角元为0和有关高斯-赛德尔迭代收敛性的相关信息及其AX=b% 的精确解jx和近似解x。
9、D=diag(diag(A);U=-triu(A,1);L=-tril(A,-1);dD=det(D);if dD=0 disp(请注意:因为对角矩阵D奇异.所以此方程组无解.)else disp(请注意:因为对角矩阵D非奇异.所以此方程组有解.) iD=inv(D-L);B2=iD*U;f2=iD*b;jx=Ab; x=x0;n m=size(A); for k=1:Nmax x1=B2*x+f2;djwcx=norm(x1-x,P); xdwcx=djwcx/(norm(x,P)+eps); if (djwcxwucha)|(xdwcxwucha) return else k;x1;k=k
10、+1;x=x1; end end if (djwcxwucha)|(xdwcxwucha) disp(请注意:高斯-赛德尔迭代收敛.此A的分解矩阵D.U,L和方程组的精确解jx和近似解如下:) else disp(请注意:高斯-赛德尔迭代的结果没有达到给定的精度.并且迭代次数已经超过最大迭代次数Nmax.方程组的精确解jx和迭代向量x如下:) x=x;jx=jx endendx=x;D,U,L,jx=jx3)松弛迭代:SORmethod(A,b,x0,Nmax,eps,w)function x = SORmethod( A,b,x,Nmax,wucha,w )% 输入的量:线性方程组AX=b的
11、系数矩阵A和b.初始向量x.范数名称P=1.2.inf.或fro.% 松弛因子w.近似解x的误差(精度)wucha和迭代的最大次数Nmax。% 输出的量:谱半径mH.以系数矩阵A的对角元构成的对角矩阵D、A的上三角矩阵U.但对角% 元为0、A的下三角矩阵L.但对角元为0.迭代次数i.有关超松弛迭代收敛性的相关信息% 及其AX=b的精确解jx和近似解x。D=diag(diag(A);U=-triu(A,1);L=-tril(A,-1);jx=Ab;n m=size(A);iD=inv(D-w*L);B2=iD*(w*U+(1-w)*D);H=eig(B2);mH=norm(H,inf);for
12、k=1:Nmax iD=inv(D-w*L);B2=iD*(w*U+(1-w)*D); f2=w*iD*b;x1=B2*x+f2; x=x1;djwcx=norm(x1-jx,inf);xdwcx=djwcx/(norm(x,inf)+eps); if (djwcxwucha)|(xdwcxNmax disp(迭代次数已经超过最大迭代次数Nmax.谱半径mH.方程组的精确解jx.迭代次数i如下:) mH,D,U,L,jx=jx,i=k-1, end endendif mH=1 disp(请注意:因为谱半径不小于1.所以超松弛迭代序列发散.) disp(谱半径mH.A的分解矩阵D,U,L和方程组
13、的精确解jx.迭代次数i和迭代序列x如下:) i=k-1,mH,D,U,L,jx,else disp(因为谱半径小于1.所以超松弛迭代序列收敛.近似解x如下:)end2-2分析应用题利用2-1中的程序来分析用下列迭代法解线性方程组:的收敛性.并求出使的近似解及相应的迭代次数.其中取迭代初始向量为零向量。1)Jacobi迭代法;2)Gauss-Seidel迭代法;3)松弛迭代法(松弛因子依次取1.334.1.95.0.95)。解:1) A=4 -1 0 -1 0 0;-1 4 -1 0 -1 0;0 -1 4 -1 0 -1;-1 0 -1 4 -1 0;0 -1 0 -1 4 -1;0 0 -
14、1 0 -1 4;b=0;5;-2;5;-2;6;X0=0 0 0 0 0 0;X=Jacobimethod( A,b,X0,inf,0.0001,100)请注意:系数矩阵A是严格对角占优的.此方程组有唯一解.且雅可比迭代收敛k = 1X = 0 1.2500 -0.5000 1.2500 -0.5000 1.5000k = 2X = 0.6250 1.0000 0.5000 1.0000 0.5000 1.2500k = 3X = 0.5000 1.6563 0.3125 1.6563 0.3125 1.7500k = 4X = 0.8281 1.5313 0.7656 1.5313 0.7
15、656 1.6563k = 5X = 0.7656 1.8398 0.6797 1.8398 0.6797 1.8828k = 6X = 0.9199 1.7813 0.8906 1.7813 0.8906 1.8398k = 7X = 0.8906 1.9253 0.8506 1.9253 0.8506 1.9453k = 8X = 0.9626 1.8979 0.9490 1.8979 0.9490 1.9253k = 9X = 0.9490 1.9651 0.9303 1.9651 0.9303 1.9745k = 10X = 0.9826 1.9524 0.9762 1.9524 0.
16、9762 1.9651k = 11X = 0.9762 1.9837 0.9675 1.9837 0.9675 1.9881k = 12X = 0.9919 1.9778 0.9889 1.9778 0.9889 1.9837k = 13X = 0.9889 1.9924 0.9848 1.9924 0.9848 1.9944k = 14X = 0.9962 1.9896 0.9948 1.9896 0.9948 1.9924k = 15X = 0.9948 1.9965 0.9929 1.9965 0.9929 1.9974k = 16X = 0.9982 1.9952 0.9976 1.9
17、952 0.9976 1.9965k = 17X = 0.9976 1.9983 0.9967 1.9983 0.9967 1.9988k = 18X = 0.9992 1.9977 0.9989 1.9977 0.9989 1.9983k = 19X = 0.9989 1.9992 0.9985 1.9992 0.9985 1.9994k = 20X = 0.9996 1.9989 0.9995 1.9989 0.9995 1.9992k = 21X = 0.9995 1.9996 0.9993 1.9996 0.9993 1.9997k = 22X = 0.9998 1.9995 0.99
18、98 1.9995 0.9998 1.9996k = 23X = 0.9998 1.9998 0.9997 1.9998 0.9997 1.9999k = 24X = 0.9999 1.9998 0.9999 1.9998 0.9999 1.9998k = 25X = 0.9999 1.9999 0.9998 1.9999 0.9998 1.9999k = 26X = 1.0000 1.9999 0.9999 1.9999 0.9999 1.9999k = 27X = 0.9999 2.0000 0.9999 2.0000 0.9999 2.0000请注意:雅可比迭代收敛.此方程组的精确解jx
19、和近似解x如下:X = 0.9999 2.0000 0.9999 2.0000 0.9999 2.00002) A=4 -1 0 -1 0 0;-1 4 -1 0 -1 0;0 -1 4 -1 0 -1;-1 0 -1 4 -1 0;0 -1 0 -1 4 -1;0 0 -1 0 -1 4;b=0;5;-2;5;-2;6;x0=0 0 0 0 0 0;x=GaussSeidelmethod( A,b,x0,inf,100,0.0001 )请注意:因为对角矩阵D非奇异.所以此方程组有解.k = 1ans = 0 1.2500 -0.1875 1.2031 0.1133 1.4814k = 2an
20、s = 0.6133 1.3848 0.5173 1.5610 0.6068 1.7810k = 3ans = 0.7364 1.7151 0.7643 1.7769 0.8183 1.8956k = 4ans = 0.8730 1.8639 0.8841 1.8938 0.9133 1.9494k = 5ans = 0.9394 1.9342 0.9444 1.9493 0.9582 1.9756k = 6ans = 0.9709 1.9684 0.9733 1.9756 0.9799 1.9883k = 7ans = 0.9860 1.9848 0.9872 1.9883 0.9903 1
21、.9944k = 8ans = 0.9933 1.9927 0.9938 1.9944 0.9954 1.9973k = 9ans = 0.9968 1.9965 0.9970 1.9973 0.9978 1.9987k = 10ans = 0.9984 1.9983 0.9986 1.9987 0.9989 1.9994k = 11ans = 0.9993 1.9992 0.9993 1.9994 0.9995 1.9997k = 12ans = 0.9996 1.9996 0.9997 1.9997 0.9998 1.9999k = 13ans = 0.9998 1.9998 0.9998
22、 1.9999 0.9999 1.9999x = 0.9998 1.9998 0.9998 1.9999 0.9999 1.99993)当w=1.334:A=4 -1 0 -1 0 0;-1 4 -1 0 -1 0;0 -1 4 -1 0 -1;-1 0 -1 4 -1 0;0 -1 0 -1 4 -1;0 0 -1 0 -1 4;b=0;5;-2;5;-2;6;x0=0 0 0 0 0 0;x=SORmethod( A,b,x,100,0.0001,1.334 )谱半径mH,A的分解矩阵D,U,L和方程组的精确解jx.迭代次数i如下:mH = 0.4305D = 4 0 0 0 0 0 0
23、4 0 0 0 0 0 0 4 0 0 0 0 0 0 4 0 0 0 0 0 0 4 0 0 0 0 0 0 4U = 0 1 0 1 0 0 0 0 1 0 1 0 0 0 0 1 0 1 0 0 0 0 1 0 0 0 0 0 0 1 0 0 0 0 0 0L = 0 0 0 0 0 0 1 0 0 0 0 0 0 1 0 0 0 0 1 0 1 0 0 0 0 1 0 1 0 0 0 0 1 0 1 0jx = 1.0000 2.0000 1.0000 2.0000 1.0000 2.0000i = 0x = 0.9999 2.0000 1.0000 2.0000 1.0000 2.0
24、000当w=1.95: A=4 -1 0 -1 0 0;-1 4 -1 0 -1 0;0 -1 4 -1 0 -1;-1 0 -1 4 -1 0;0 -1 0 -1 4 -1;0 0 -1 0 -1 4;b=0;5;-2;5;-2;6;x0=0 0 0 0 0 0;x=SORmethod( A,b,x,100,0.0001,1.95 )谱半径mH,A的分解矩阵D,U,L和方程组的精确解jx.迭代次数i如下:mH = 0.9589D = 4 0 0 0 0 0 0 4 0 0 0 0 0 0 4 0 0 0 0 0 0 4 0 0 0 0 0 0 4 0 0 0 0 0 0 4U = 0 1 0
25、 1 0 0 0 0 1 0 1 0 0 0 0 1 0 1 0 0 0 0 1 0 0 0 0 0 0 1 0 0 0 0 0 0L = 0 0 0 0 0 0 1 0 0 0 0 0 0 1 0 0 0 0 1 0 1 0 0 0 0 1 0 1 0 0 0 0 1 0 1 0jx = 1.0000 2.0000 1.0000 2.0000 1.0000 2.0000i = 0x = 1.0000 2.0000 1.0000 2.0000 1.0000 2.0000当w=0.95: A=4 -1 0 -1 0 0;-1 4 -1 0 -1 0;0 -1 4 -1 0 -1;-1 0 -1
26、4 -1 0;0 -1 0 -1 4 -1;0 0 -1 0 -1 4;b=0;5;-2;5;-2;6;x0=0 0 0 0 0 0;x=SORmethod( A,b,x,100,0.0001,0.95 )谱半径mH,A的分解矩阵D,U,L和方程组的精确解jx.迭代次数i如下:mH = 0.5260D = 4 0 0 0 0 0 0 4 0 0 0 0 0 0 4 0 0 0 0 0 0 4 0 0 0 0 0 0 4 0 0 0 0 0 0 4U = 0 1 0 1 0 0 0 0 1 0 1 0 0 0 0 1 0 1 0 0 0 0 1 0 0 0 0 0 0 1 0 0 0 0 0 0
27、L = 0 0 0 0 0 0 1 0 0 0 0 0 0 1 0 0 0 0 1 0 1 0 0 0 0 1 0 1 0 0 0 0 1 0 1 0jx = 1.0000 2.0000 1.0000 2.0000 1.0000 2.0000i = 0x = 1.0000 2.0000 1.0000 2.0000 1.0000 2.00002-3分析应用题考虑方程组.其中系数矩阵为Hilbert矩阵. 选择问题的维数分别为2、3、5、10.并通过首先给定解再定出右端的办法确定问题.解的给定可以使用函数定义.并取迭代初始向量为零向量.迭代误差为.编写程序: 其中n为Hilbert矩阵的维数.分别构造求解该问题的Jacobi迭代和Gauss-Seidel迭代.看它们是否收敛。2-4分析应用题解线性方程组的解向量.取.其中为任一非零的六元向量;编写程序输出结果: 认真观察之.能发现什么有趣的现象?四. 实验结果与分析