迭代法解线性方程组作业
沈欢00986096
北京大学工学院,北京100871
2011年10月12日
摘要
由所给矩阵生成系数矩阵A和右端项b,分析系数矩阵A,并用Jacobi迭代法、GS迭代法、SOR(逐步松弛迭代法)解方程组Ax=b
1生成系数矩阵A、右端项b,并分析矩阵A
由文件”gr900900c rg.mm”得到了以.mm格式描述的系数矩阵A。
A矩阵是900∗900的大型稀
疏对称矩阵。
于是,在matlaB中,使用”A=zeros(900,900)”语句生成900∗900的零矩阵。
再
按照.mm文件中的描述,分别对第i行、第j列的元素赋对应的值,就生成了系数矩阵A,并
将A存为.mat文件以便之后应用。
由于右端项是全为1的列向量,所以由语句”b=ones(900,1)”生成。
得到了矩阵A后,求其行列式,使用函数”det(A)”,求得结果为”Inf”,证明行列式太大,matlaB无法显示。
由此证明,矩阵A可逆,线性方程组
Ax=b
有唯一解。
接着,判断A矩阵是否是对称矩阵(其实,这步是没有必要的,因为A矩阵本身是对称矩阵,是.mm格式中的矩阵按对称阵生成的)。
如果A是对称矩阵,那么
A−A T=0。
于是,令B=A−A T,并对B求∞范数。
结果显示: B ∞=0,所以,B是零矩阵,也就是:A是对称矩阵。
然后,求A的三个条件数:
Cond(A)= A ∗ A−1
所求结果是,对应于1范数的条件数为:377.2334;对应于2范数的条件数为:194.5739;对应
于3范数的条件数为:377.2334;
1
从以上结果我们看出,A是可逆矩阵,但是A的条件数很大,所以,Ax=b有唯一解并且矩阵A相对不稳定。
所以,我们可以用迭代方法来求解该线性方程组,但是由于A的条件数太大迭代次数一般而言会比较多。
2Jacobi迭代法
Jacobi迭代方法的程序流程图如图所示:
图1:Jacobi迭代方法程序流程图
在上述流程中,取x0=[1,1,...,1]T将精度设为accuracy=10−3,需要误差满足:
error= x k+1−x k
x k+1
<accuracy
时迭代结束。
在上述Jacobi迭代中,只有求取x k+1一步需要用到Jacobi迭代的知识,Jacobi迭代的方法简单。
可以并行计算,其分量形式为:
x k+1 1=(b1−
j=1
a1j x k j)/a11;
x k+1 2=(b2−
j=2
a2j x k j)/a22;
2
.........
x k+1 n =(b n−
j=n
a nj x k j)/a nn.
于是,按上述分量形式可以直接编制程序求取x k+1.
最终,程序中设置的迭代次数标志flag显示,经过36次迭代,计算出了满足精度要求的解x,并附于本文最后。
3GS迭代法
GS迭代法思路基本和Jacobi迭代法相同。
只是在求x k+1
i 时,同时使用x k+1
1
、x k+1
2
...x k+1
i−1
和x k
i+1
、x k
i+1
...x k n的
信息。
所以,GS迭代法的程序流程图与Jacobi迭代法的相同,如下:
图2:GS迭代方法程序流程图
对于GS迭代求解x k+1的方法与Jacobi迭代有所不同,其分量形式为:
x k+1 i =(b i−
i−1
j=1
a ij x k+1
j
−
n
j=i+1
a ij x k j)/a ii
i=1,2......n
3
根据分量形式可以编制求x k+1的程序。
其他步骤和设置,均与Jacobi迭代法相同。
最终,程序中设置的迭代次数标志flag显示,经过39次迭代,GS迭代法计算出了满足和Jacobi迭代相同精度要求的解x。
4SOR迭代
SOR迭代是GS迭代的扩展,其x k+1
i 是GS迭代方法得到的x k+1
i
和x k
i
的加权平均。
即:
x k+1 i =w∗((b i−
i−1
j=1
a ij x k+1
j
−
n
j=i+1
a ij x k j)/a ii)+(1−w)∗x k i
其中,w是松弛因子,当0<w<2时,SOR方法收敛。
所以,选择合适的w就是一个关键的问题。
SOR
图3:SOR迭代方法程序流程图
在所有其余步骤和设置与GS方法相同的情况下,本文笔者随意设置了的w值为0.88,最终需要的迭代次数为40,得到与GS方法具有相同精度的解。
笔者又尝试将w设置为更小的0.60,需要迭代34次.笔者又将w设置为较大的1.50,需要迭代30次。
以上数据说明,SOR方法对w的敏感性很高,选择较为适当的w对计算才有帮助,否则反而可能回增加迭代次数。
4
结论和说明
本文用三种迭代解法求取线性方程组Ax =b 的解。
解向量如下(由于三种解精度相同,所以就只附上用Jacobi 解法的解):0100200300400500600700800900
02
4
6
8
10
12
图4:解的图示:图中点的意义为(k,x(k)=y)
5。