线性方程组6种数值解法的对比研究
线性方程组数值解法实验研究
一、实验目的
熟悉求解线性方程组的有关理论和方法；会编写Gauss消去法、LU分解法、Jacobi迭代法、Gauss-Seidel迭代法、超松弛(SOR)迭代法及共轭梯度法的程序；通过实际计算，进一步了解各种方法的优缺点，选择合适的数值方法。
二、实验内容
- 编程实现Gauss消去法、LU分解法，并用这两种方法求解线性方程组：
⎧⎩⎨⎪⎪x1+2x2−2x3=7x1+x2+x3=22x1+2x2+x3=5(1)
\left\{ \begin{array}{l} {x_1} + 2{x_2} - 2{x_3} = 7\\ {x_1} + {x_2} + {x_3} = 2\\ 2{x_1} + 2{x_2} + {x_3} = 5 \end{array} \right. \tag 1
- 编写Jacobi迭代法、高斯-赛德尔迭代法程序，并用其求解以下84阶线性方程组，并指出各迭代法是否收敛？实验结果说明什么问题？
⎛⎝⎜⎜⎜⎜⎜⎜⎜⎜⎜⎜⎜⎜6816816⋱1⋱8⋱6816816⎞⎠⎟⎟⎟⎟⎟⎟⎟⎟⎟⎟⎟⎟⎛⎝⎜⎜⎜⎜⎜⎜⎜⎜⎜⎜⎜⎜x1x2x3⋮x82x83x84⎞⎠⎟⎟⎟⎟⎟⎟⎟⎟⎟⎟⎟⎟=⎛⎝⎜⎜⎜⎜⎜⎜⎜⎜⎜⎜⎜⎜71515⋮151514⎞⎠⎟⎟⎟⎟⎟⎟⎟⎟⎟⎟⎟⎟(2)
\left( {\begin{array}{*{20}{c}} 6&1\\ 8&6&1\\ {}&8&6&1\\ {}& \ddots & \ddots & \ddots \\ {}&8&6&1\\ {}&8&6&1\\ {}&8&6 \end{array}} \right)\left( {\begin{array}{*{20}{c}} {{x_1}}\\ {{x_2}}\\ {{x_3}}\\\vdots \\ {{x_{82}}}\\ {{x_{83}}}\\ {{x_{84}}} \end{array}} \right) = \left( {\begin{array}{*{20}{c}} 7\\ {15}\\ {15}\\\vdots \\ {15}\\ {15}\\ {14} \end{array}} \right)\tag 2
- 编写超松弛(SOR)迭代法程序，并用其求解以下4阶方程组，并分析超松弛因子 对SOR收敛性的影响。实验结果说明了什么问题？
⎛⎝⎜⎜⎜13.422−13.42200012.252−12.25200012.377−12.37700011.797⎞⎠⎟⎟⎟⎛⎝⎜⎜⎜x1x2x3x4⎞⎠⎟⎟⎟=⎛⎝⎜⎜⎜750.530010230⎞⎠⎟⎟⎟(3)
\left( {\begin{array}{*{20}{c}} {13.422}&0&0&0\\ { - 13.422}&0&0\\ 0&0\\ 0&0\end{array}} \right)\left( {\begin{array}{*{20}{c}} {{x_1}}\\ {{x_2}}\\ {{x_3}}\\ {{x_4}} \end{array}} \right) = \left( {\begin{array}{*{20}{c}} {750.5}\\ {300}\\ {102}\\ {30} \end{array}} \right)\tag 3
- 编写共轭梯度法程序 ，并用其求解问题(1)(2)(3)中的方程组，并求出共轭梯度法在这三种问题下的收敛速度，列出详细的分析过程。
三、实验结果分析
1. Gauss消去法与LU分解法求解线性方程组
根据线性方程组公式(1)可得，系数矩阵 AA 和常数向量 bb
A = \left( {\begin{array}{*{20}{c}} 1&2\\ 1&1&1\\ 2&2&1 \end{array}} \right),b = \left( {\begin{array}{*{20}{c}} 7\\ 2\\ 5 \end{array}} \right)
假设
x = {\left[ {\begin{array}{*{20}{c}} {{x_1}}}} \end{array}} \right]^{\rm{T}}}
由此可得
Ax = b
其中，根据矩阵 AA 的初等行变换，可以得到矩阵 AA 的秩为3，因此，该线性方程组具有唯一确定解。Gauss消去法和LU分解法都可以得到解。其结果如表1所示。
表1 实验1的结果
| 方法 | xx | 计算误差(>2.22E−16>2.22E-16) |
|---|---|---|
| Gauss消去法 | [12−1]T{\left[ {\begin{array}{*{20}{c}}1&2}1}\end{array}} \right]^{\rm{T}}} | 0 |
| LU分解法 | [12−1]T{\left[ {\begin{array}{*{20}{c}}1&2}1}\end{array}} \right]^{\rm{T}}} | 0 |
2. Jacobi迭代法与Gauss-Seidel迭代法求解84阶线性方程组
根据公式(2)，建立线性方程组标准方程
Ax = b
式中
A = \left( {\begin{array}{*{20}{c}} 6&1\\ 8&6&1\\ {}&8&6&1\\ {}& \ddots & \ddots & \ddots \\ {}&8&6&1\\ {}&8&6&1\\ {}&8&6 \end{array}} \right),b = \left( {\begin{array}{*{20}{c}} 7\\ {15}\\ {15}\\\vdots \\ {15}\\ {15}\\ {14} \end{array}} \right)
x = {\left( {\begin{array}{*{20}{c}} {{x_1}}}}& \cdots }}}}}} \end{array}} \right)^{\rm{T}}}
基于迭代公式
{x^{\left( {k + 1} \right)}} = G{x^{\left( k \right)}} + f\tag 4
将系数矩阵 AA 分解成对角矩阵和上(下)三角矩阵的组合，即
A = D - L - U
其中
D = {\rm{diag}}\left( {\begin{array}{*{20}{c}} 6&6&6& \cdots &6&6&6 \end{array}} \right)
L = - \left( {\begin{array}{*{20}{c}} 0\\ 8&0\\ {}&8&0\\ {}& \ddots & \ddots \\ {}&8&0\\ {}&8&0\\ {}&8&0 \end{array}} \right)
U = - \left( {\begin{array}{*{20}{c}} 0&1\\ {}&0&1\\ {}&0&1\\ {}& \ddots & \ddots \\ {}&0&1\\ {}&0&1\\ {}&0 \end{array}} \right)
根据问题的特殊性，很容易发现该问题的解为
{x^*} = \left( {\begin{array}{*{20}{c}} 1&1&1& \cdots &1&1&1 \end{array}} \right)
由于常系数矩阵 AA 不是满秩矩阵，因此该问题存在多个解。
2.1 Jacobi迭代法
Jacobi迭代法的迭代公式为
D{x^{\left( {k + 1} \right)}} = \left( {L + U} \right){x^{\left( k \right)}} + b\tag 5
将公式(5)转化为公式(4)的形式，则可得
G = {D^{ - 1}}\left( {L + U} \right)
f = {D^{ - 1}}b
接下来分析该迭代法时候收敛，计算矩阵 GG 的谱半径 ρ(G)\rho \left( G \right)
\rho \left( G \right) = \max \left\{ {\left| {{\lambda _i}} \right|{\rm{|}}{\lambda _i} \in \lambda ,A\varepsilon = \lambda \varepsilon } \right\}{\rm{ = }}1.1195 > 1
因此，Jacobi迭代法对该问题不收敛，无法求解。
2.2 Gauss-Seidel迭代法
Gauss-Seidel迭代法的迭代公式为
\left( {D - L} \right){x^{\left( {k + 1} \right)}} = U{x^{\left( k \right)}} + b\tag 6
将公式(6)转化为公式(4)的形式，则可得
G = {\left( {D - L} \right)^{ - 1}}U
f = {\left( {D - L} \right)^{ - 1}}b
接下来分析该迭代法时候收敛，计算矩阵 GG 的谱半径 ρ(G)\rho \left( G \right)
\rho \left( G \right) = \max \left\{ {\left| {{\lambda _i}} \right|{\rm{|}}{\lambda _i} \in \lambda ,A\varepsilon = \lambda \varepsilon } \right\}{\rm{ = }}0.88747
因此，Gauss-Seidel迭代法可以求解该问题，计算结果如表2所示。
表2 问题2的Gauss-Seidel迭代法计算结果
| 迭代次数 | 允许误差 | 实际误差 |
|---|---|---|
| 598 | 1E-4 | 9.4250E-5 |
以上2种方法求解该问题时，发现Jacobi迭代法不收敛，Gauss-Seidel迭代法收敛，在允许精度为1E-4的情况下，迭代次数为598次，可见Gauss-Seidel迭代法在求解该线性方程组时，收敛速度不佳。
3. 超松弛(SOR)迭代法
在Gauss-Seidel迭代法的基础上，为获得更快的收敛效果，在修正量前乘以一个修正系数 ω\omega，即得到逐次超松弛迭代法，简称SOR迭代法，其中 ω\omega 为松弛因子。SOR迭代法的迭代格式为
{x^{\left( {k + 1} \right)}} = {x^{\left( k \right)}} + \omega {D^{ - 1}}\left( {L{x^{\left( {k + 1} \right)}} + U{x^{\left( k \right)}} - D{x^{\left( k \right)}} + b} \right)\tag 7
将公式(7)转化为公式(4)的形式，即
G = {\left( {D - \omega L} \right)^{ - 1}}\left[ {\left( {1 - \omega } \right)D + \omega U} \right]{x^{\left( k \right)}}
f = \omega {\left( {D - \omega L} \right)^{ - 1}}
接下来分析该迭代法时候收敛，计算矩阵 GG 的谱半径 ρ(G)\rho \left( G \right)
\rho \left( G \right) = \max \left\{ {\left| {{\lambda _i}} \right|{\rm{|}}{\lambda _i} \in \lambda ,A\varepsilon = \lambda \varepsilon } \right\}
则SOR迭代法的谱半径 ρ(G)\rho \left( G \right)和松弛因子 ω\omega 相关，因此，我们讨论 ω\omega 和谱半径及迭代次数的关系。采用问题(3)的线性方程组，分析松弛因子 ω\omega 对收敛速度即谱半径 ρ(G)\rho \left( G \right)的影响，其中 ω∈(0,2)\omega \in \left( {0,2} \right)，结果如图1所示。
系数矩阵A和常数向量b分别为
A = \left( {\begin{array}{*{20}{c}} {13.422}&0&0&0\\ { - 13.422}&0&0\\ 0&0\\ 0&0\end{array}} \right),b = \left( {\begin{array}{*{20}{c}} {750.5}\\ {300}\\ {102}\\ {30} \end{array}} \right)
由计算结果可得，当 ω=1\omega = 1 时，迭代次数最小，其迭代次数为1。图1描述了随松弛因子变化时，曲率半径和迭代次数的变化。由此，我们发现，在这个算例下，迭代次数和曲率半径呈正相关，曲率半径最小的时候，迭代次数最小。
图1 问题3的SOR迭代法计算结果(计算精度为1E-4)
4. 共轭梯度法
共轭梯度法与之前的迭代法不同，它属于不定常迭代法，用于求解对称正对线性方程组，其本质上是一种变分方法。
共轭梯度法的迭代方程为
{x^{\left( {k + 1} \right)}} = {x^{\left( k \right)}} + {\alpha _k}{p^{\left( k \right)}}
其中
αk=∥∥r(k)∥∥22/(Ap(k),p(k)){\alpha _k}=\left\| r^\left( k \right)\right\|_2^2 / \left( Ap^{(k)},p^{(k)}\right), p(k){p^{\left( k \right)}}为第 kk次迭代过程中的搜索方向，它与梯度方向 r(k){r^{\left( k \right)}} 相关，即
{p^{\left( {k + 1} \right)}} = {r^{\left( {k + 1} \right)}} + {\beta _k}{p^{\left( k \right)}}
其中
{r^{\left( {k + 1} \right)}} = {r^{\left( k \right)}} - \alpha A{p^{\left( k \right)}}
{\beta _k} = \left\| r^\left( k+1 \right)\right\|_2^2/\left\| r^\left( k \right)\right\|_2^2
当 ∥βk∥<TOL\left\| {{\beta _k}} \right\| 时，则迭代结束，其中 TOLTOL为允许误差。
对于初始解，共轭迭代法要求其为最速下降方向，这可以保证所有的迭代过程中的 x(k){x^{\left( k \right)}}都是共轭的，且具有较快地收敛速度。
针对问题(1)(2)(3)中的算例，共轭迭代法首先要求对系数矩阵作对角化处理，即
b = {A^{\rm{T}}}b,A = {A^{\rm{T}}}A
其计算结果如表3所示。
表3 问题1,2,3的共轭迭代法计算结果
| 问题序号 | 迭代次数 | 计算误差 | 允许误差 |
|---|---|---|---|
| 1 | 2 | 7.9422E-5 | |
| 2 | 11 | 8.8953E-5 | 1E-4 |
| 3 | 4 | 1.9327E-5 |
由表3所得，共轭迭代法在问题(1)(2)(3)中的线性方程组求解上具有较好的收敛效果，尤其对问题(2)中的大型稀疏矩阵上，其收敛效果明显高于Gauss-Seidel迭代法。
四、实验结论
线性方程组的直接求法，Gauss消去法和LU分解法在维度较小的非病态问题上，具有良好的效果。这两种算法理论上的复杂度为 O(n3)O\left( {{n^3}} \right)，使其能有效地求出低维非病态线性方程组的精确解，但在髙维度或病态的线性方程组问题上出现无法在有限时间内收敛或无法求解的情况。
线性方程组的迭代求法，Jacobi迭代法、Gauss-Seidel迭代法、超松弛(SOR)迭代法和共轭梯度法针对大型稀疏矩阵（如问题(2)的84阶线性方程组）表现出不同的效果。Jacobi迭代法的不收敛、Gauss-Seidel迭代法的极大的迭代次数以及共轭梯度法较优的计算结果，表明了这些迭代法不同的计算特性。SOR迭代法的松弛因子 ω\omega 对算法的收敛性影响较大，同时，在问题(3)的分析与求解过程中，我们发现谱半径与SOR迭代法呈正相关特性。
基于以上实验结果，我们得到以下结论：
1. 直接求法在一般问题，维数较小下(系数矩阵的秩小于100)，非常实用。
2. 迭代法的收敛性在不同问题反映出的现象差异很大。
3. 采用Jacobi迭代法、Gauss-Seidel迭代法求解大型稀疏矩阵，收敛性难以保证，而共轭迭代法收敛极好。
4. 超松弛迭代法的收敛性受松弛因子影响，并可能存在最优的松弛因子，使其收敛速度最快。
五、实验程序
所有程序都在MATLAB(R2012b)平台上实现。
1. Gauss 消去法
| 源代码 |
|---|
function [rankA, rankB, N, X] = CNumbericGauss(A, b)
% 调用格式: [rankA, rankB, N, X] = CNumbericGauss(A, b)
%
% 作者: 王瑞
% 时间: 2015.10.27 17:12 - 20.12
% 版本: Version 1.0
%
% 任务: 高斯消去法求解线性方程组的解 Ax = b, 换主元的高斯消去法, 取最邻近非 0 首项
%
% 输入: A = 系数矩阵, 方阵
% b = 常系数向量, 行向量
%
% 输出: rankA = 系数矩阵 A 的秩
% rankB = 增广矩阵 B 的秩, 其中 B = [A|b]
% N = 方程组未知量个数
% X = 方程组的解向量
% B = [A b];
rankA = rank(A);
rankB = rank(B);
N = length(b);
X = zeros(size(b));
if rankA ~= rankBdisp('None: 矩阵 A 的秩不等于 [A b] 的秩， 无解！');return ;
end
if rankA < Ndisp(['Waring: 矩阵的秩小于' num2str(N) '，存在无穷多解。']);return;
end
if N == 1disp('别闹，用手算的。');return;
end
%% Gauss 消去法
if rankA == Ndisp(['Good: 矩阵的秩为' num2str(N) '，存在唯一解。']);for i = 1 : Nif abs(B(i, i)) <= epsB(i, i) = 0;[flag, j] = max(abs(B(i:N, i)));j = j + i -1;temp = B(i, :);B(i, :) = B(j, :);B(j, :) = temp;if abs(flag) <= eps % 无非 0 首项判定disp('这里出现了一个 Error，首项极小！');return;endendB(i, :) = B(i, :)/B(i, i);for j = i+1 : N % Gauss 消去B(j, :) = B(j, :) - B(j, i)*B(i, :);endendfor i = N : -1 : 1 % Gauss 回代求解if i == NX(i) = B(i, N+1)/B(i, N);elseX(i) = (B(i, N+1) - B(i, i+1:N)*X(i+1:N)) / B(i, i);endend
end
end2. LU分解法
| 源代码 |
|---|
function [rankA, rankB, N, X] = CNumbericLU(A, b)
% 调用格式: [rankA, rankB, N, X] = CNumbericLU(A, b)
%
% 作者: 王瑞
% 时间: 2015.10.27 20:14 - 21:37
% 版本: Version 1.0
%
% 任务: 选主元的三角分解法求解线性方程组的解 Ax = b, LU 法
%
% 输入: A = 系数矩阵, 方阵
% b = 常系数向量, 行向量
%
% 输出: rankA = 系数矩阵 A 的秩
% rankB = 增广矩阵 B 的秩, 其中 B = [A|b]
% N = 方程组未知量个数
% X = 方程组的解向量
%
B = [A b];
rankA = rank(A);
rankB = rank(B);
N = length(b);
X = zeros(size(b));
if rankA ~= rankBdisp('None: 矩阵 A 的秩不等于 [A b] 的秩， 无解！');return ;
elseif rankA ~= Ndisp(['Waring: 矩阵的秩小于' num2str(N) '，存在无穷多解。']);return;
end
if rankA == 1disp('别闹，用手算的。');return;
end
disp(['Good: 矩阵的秩为' num2str(N) '，存在唯一解。']);
%% 求 L U P [核心]
[L, U, P] = lu(A);
b = P*b;
%% 求 y
y = zeros(size(b));
for i = 1 : Nif i == 1y(1) = b(1);elsey(i) = b(i) - L(i, 1:i-1)*y(1:i-1);end
end
%% 求 X
for i = N : -1 : 1if i == NX(i) = y(i) / U(i, i);elseX(i) = (y(i) - U(i, i+1:N)*X(i+1:N)) / U(i, i);end
end
end
3. Jacobi迭代法
| 源代码 |
|---|
function [rankA, rankB, N, X, ite, tol] = CNumbericJacobiIteration(A, b, TOL, ITE, initX)
% 调用格式: [rankA, rankB, N, X, ite, tol] = CNumbericJacobiIteration(A, b, initX, TOL, ITE, initX)
% [rankA, rankB, N, X, ite, tol] = CNumbericJacobiIteration(A, b, initX, TOL, ITE)
%
% 作者: 王瑞
% 时间: 2015.10.27 21:40; 2015.10.28 13:37 - 15:00
% 版本: Version 1.0
%
% 任务: Jacobi迭代法求解线性方程组的解 Ax = b
% 构建 x(k+1) = Bx(k) + f
% B = E - D^-1*A = D^-1*(L+U), f = D^-1*b;
%
% 输入: A = 系数矩阵, 方阵
% b = 常系数向量, 行向量
% initX = 初始解
% ITE = 迭代次数上限
% TOL = 解的精度(范数)
%
% 输出: rankA = 系数矩阵 A 的秩
% rankB = 增广矩阵 B 的秩, 其中 B = [A|b]
% N = 方程组未知量个数
% X = 方程组的解向量
% ite = 求解的迭代次数
% tol = 实际误差
%
if nargin == 5X = initX;
elseif nargin == 4X = zeros(size(b));
endB = [A, b];
rankA = rank(A);
rankB = rank(B);
N = length(b);
ite = 0;
tol = inf;
if rankA ~= rankBdisp('None: 矩阵 A 的秩不等于 [A b] 的秩， 无解！');return ;
elseif rankA ~= Ndisp(['Waring: 矩阵的秩小于' num2str(N) '[' num2str(rankA) ']''，存在无穷多解。']);
end
if rankA == 1disp('别闹，用手算的。');return;
end
if rankA == Ndisp(['Good: 矩阵的秩为' num2str(N) '，存在唯一解。']);
end
%% B0, f
D = diag(A);
L = (-1) * tril(A, -1);
U = (-1) * triu(A, 1);
invD = diag(1./D);
B0 = invD*(L+U);
f = invD*b;
R = max(abs(eig(B0))); % 谱半径，收敛判定
if R >= 1disp(['Error: 谱半径 R = ' num2str(R) '，大于等于 1，算法不收敛。']);return;
end
%% Jacobi 迭代
while ite < ITEite = ite + 1;X = B0*X + f;tol = norm(b - A*X) / norm(b);if tol < TOLdisp('Excatly: 求得解。');break;end
end
if ite > ITEdisp(['Message: 在 ' num2str(ITE) ' 次迭代过程中，无法求得解，'...'建议增大总的迭代次数 或者 分析算法的收敛性。']);
end
end
4. Gauss-Seidel迭代法
| 源代码 |
|---|
function [rankA, rankB, N, X, ite, tol] = CNumbericGaussSeidelIteration(A, b, TOL, ITE, initX)
% 调用格式: [rankA, rankB, N, X, ite, tol] = CNumbericGaussSeidelIteration(A, b, TOL, ITE, initX)
% [rankA, rankB, N, X, ite, tol] = CNumbericGaussSeidelIteration(A, b, TOL, ITE)
%
% 作者: 王瑞
% 时间: 2015.10.28 14:15 - 15:00
% 版本: Version 1.0
%
% 任务: Gauss-Seidel迭代法求解线性方程组的解 Ax = b
% 构建 x(k+1) = Gx(k) + f
% G = (D - L)^-1*U, f = (D - L)^-1*b
%
% 输入: A = 系数矩阵, 方阵
% b = 常系数向量, 行向量
% initX = 初始解
% ITE = 迭代次数上限
% TOL = 解的精度(范数)
%
% 输出: rankA = 系数矩阵 A 的秩
% rankB = 增广矩阵 B 的秩, 其中 B = [A|b]
% N = 方程组未知量个数
% X = 方程组的解向量
% ite = 求解的迭代次数
% tol = 实际误差
%
if nargin == 5X = initX;
elseif nargin == 4X = zeros(size(b));
endB = [A, b];
rankA = rank(A);
rankB = rank(B);
N = length(b);
ite = 0;
tol = inf;
if rankA ~= rankBdisp('None: 矩阵 A 的秩不等于 [A b] 的秩， 无解！');return ;
elseif rankA ~= Ndisp(['Waring: 矩阵的秩小于' num2str(N) '[' num2str(rankA) ']''，存在无穷多解。']);
end
if rankA == 1disp('别闹，用手算的。');return;
end
if rankA == Ndisp(['Good: 矩阵的秩为' num2str(N) '，存在唯一解。']);
end
%% G, f
D = diag(diag(A));
L = (-1) * tril(A, -1);
U = (-1) * triu(A, 1);
invDL = eye(length(b))/(D - L);
G = invDL*U;
f = invDL*b;
R = max(abs(eig(G))); % 谱半径，收敛判定
if R >= 1disp(['Error: 谱半径 R = ' num2str(R) '，大于等于 1，算法不收敛。']);return;
end
%% Gauss-Seidel 迭代
while ite < ITEite = ite + 1;Xk = X;X = G*Xk + f;tol = norm(b - A*X) / norm(b);if tol < TOLdisp('Excatly: 求得解。');break;end
end
if ite > ITEdisp(['Message: 在 ' num2str(ITE) ' 次迭代过程中，无法求得解，'...'建议增大总的迭代次数 或者 分析算法的收敛性。']);
end
end
5. 超松弛(SOR)迭代法
| 源代码 |
|---|
function [rankA, rankB, N, X, ite, tol] = CNumbericSORIteration(A, b, Omega, TOL, ITE, initX)
% 调用格式: [rankA, rankB, N, X, ite, tol] = CNumbericSORIteration(A, b, Omega, TOL, ITE, initX)
% [rankA, rankB, N, X, ite, tol] = CNumbericSORIteration(A, b, Omega, TOL, ITE)
%
% 作者: 王瑞
% 时间: 2015.10.28 16:00 - 16:38
% 版本: Version 1.0
%
% 任务: 超松弛(SOR)迭代法求解线性方程组的解 Ax = b
% 构建 x(k+1) = Gx(k) + f
% G = (D-Omege*L)^-1*[(1-Omega)*D+Omega*U],
% f = Omega*(D-Omega*L)^-1*b
%
% 输入: A = 系数矩阵, 方阵
% b = 常系数向量, 行向量
% initX = 初始解
% Omege = 松弛因子 (0, 2)
% ITE = 迭代次数上限
% TOL = 解的精度(范数)
%
% 输出: rankA = 系数矩阵 A 的秩
% rankB = 增广矩阵 B 的秩, 其中 B = [A|b]
% N = 方程组未知量个数
% X = 方程组的解向量
% ite = 求解的迭代次数
% tol = 实际误差
%
if nargin == 5X = zeros(size(b));
elseif nargin == 6X = initX;
endB = [A, b];
rankA = rank(A);
rankB = rank(B);
N = length(b);
ite = 0;
tol = inf;
if rankA ~= rankBdisp('None: 矩阵 A 的秩不等于 [A b] 的秩， 无解！');return ;
elseif rankA ~= Ndisp(['Waring: 矩阵的秩小于' num2str(N) '[' num2str(rankA) ']''，存在无穷多解。']);
end
if rankA == 1disp('别闹，用手算的。');return;
end
if rankA == Ndisp(['Good: 矩阵的秩为' num2str(N) '，存在唯一解。']);
end
%% D, L, U
D = diag(diag(A));
L = (-1) * tril(A, -1);
U = (-1) * triu(A, 1);
%% G, f
G = (D - Omega*L)\((1-Omega)*D + Omega*U);
f = (D - Omega*L)\(Omega*b);
R = max(abs(eig(G))); % 谱半径，收敛判定
if R >= 1disp(['Error: 谱半径 R = ' num2str(R) '，大于等于 1，算法不收敛。']);return;
end
disp(['R = ' num2str(R)]);
%% SOR 迭代
while ite < ITEite = ite + 1;X = G*X + f;tol = norm(b - A*X) / norm(b);if tol < TOLdisp('Excatly: 求得解。');break;end
end
if ite > ITEdisp(['Message: 在 ' num2str(ITE) ' 次迭代过程中，无法求得解，'...'建议增大总的迭代次数 或者 分析算法的收敛性。']);
end
end
6. 共轭梯度法
| 源代码 |
|---|
function [rankA, rankB, N, X, ite, tol] = CNumbericCGMIteration(A, b, TOL, ITE, initX)
% 调用格式: [rankA, rankB, N, X, ite, tol] = CNumbericCGMIteration(A, b, TOL, ITE, initX)
% [rankA, rankB, N, X, ite, tol] = CNumbericCGMIteration(A, b, TOL, ITE)
%
% 作者: 王瑞
% 时间: 2015.10.28 16:54 - 19:21
% 版本: Version 1.0
%
% 任务: 共轭梯度法(Conjugate Gradient Method, CGM)迭代法求解线性方程组的解 Ax = b
% 适用于系数矩阵为对称阵的线性方程组(函数内包含矩阵对称化 b = A'*b, A = A'*A)
% 构建 x(k+1) = x(k) + alpha(k)*p(k)
% 等同 MATLAB 内置函数 cgs
%
% 输入: A = 系数矩阵, 方阵
% b = 常系数向量, 行向量
% ITE = 迭代次数上限
% TOL = 解的精度(范数)
% initX = 初始解
%
% 输出: rankA = 系数矩阵 A 的秩
% rankB = 增广矩阵 B 的秩, 其中 B = [A|b]
% N = 方程组未知量个数
% X = 方程组的解向量
% ite = 求解的迭代次数
% tol = 实际误差
%
if nargin == 5X = initX;
elseif nargin == 4X = zeros(size(b));
end
%% 解的判定
B = [A, b];
rankA = rank(A);
rankB = rank(B);
N = length(b);
ite = 0;
tol = inf;
if rankA ~= rankBdisp('None: 矩阵 A 的秩不等于 [A b] 的秩， 无解！');return ;
elseif rankA ~= Ndisp(['Waring: 矩阵的秩小于' num2str(N) '[' num2str(rankA) ']''，存在无穷多解。']);
end
if rankA == 1disp('别闹，用手算的。');return;
end
if rankA == Ndisp(['Good: 矩阵的秩为' num2str(N) '，存在唯一解。']);
end
% 系数矩阵对称，放大 A' 倍
if rank(A-A') ~= 0b = A'*b;A = A'*A;
end
% GCM 迭代
r = b - A*X;
while ite < ITEerr = r'*r;ite = ite + 1;if ite == 1p = r;elsebeta = err / errold;p = r + beta*p;endAp = A*p;alpha = err / ((Ap)'*p);X = X + alpha*p;r = r - alpha*Ap;errold = err;tol = norm(r) / norm(b);if tol < TOLdisp('Excatly: 求得解。');break;end
end
if ite > ITEdisp(['Message: 在 ' num2str(ITE) ' 次迭代过程中，无法求得解，'...'建议增大总的迭代次数 或者 分析算法的收敛性。']);
end
end线性方程组6种数值解法的对比研究相关推荐
- matlab 非线性方程数值解法,非线性方程组的几种数值解法+matlab源代码
摘要很多领域都有涉及到非线性方程组,例如天气预报,石油地质勘探,电力系统计算等,甚至商业领域也有非线性优化问题,这些问题要从本质上解决就是求出非线性方程组的解.但是目前已知的数值解法并不完善,选择不同 ...
- 常微分方程数值解的c语言程序,常微分方程的数值解法 一阶常微分方程数值解的C语言编程实现.doc...
常微分方程的数值解法 一阶常微分方程数值解的C语言编程实现 导读:就爱阅读网友为您分享以下"一阶常微分方程数值解的C语言编程实现"资讯,希望对您有所帮助,感谢您对92的支持! 一阶 ...
- 数值计算方法第三章—线性方程组的数值解法知识点总结
线性方程组的数值解法 本文参考书为马东升著<数值计算方法> 高斯消去法 顺序高斯消去法 通过初等变换消去方程组系数矩阵主对角线以下的元素,而使方程组化为等价的上三角形方程组 列主元高斯消去 ...
- matlab二阶迎风差分格式,热传导方程几种差分格式的MATLAB数值解法比较
第2 5卷 第 2期 沈 阳 化 工 大 学 学 报 V0 . 5 N . 1 o2 2 J n. 01l u 2 2 011 0 . 6 J OURNAL HENYANG OF S UNI VERS ...
- 二维非稳态导热微分方程_一种二维非稳态导热问题的数值解法.pdf
一种二维非稳态导热问题的数值解法 维普资讯 第 9卷 第 2f/I 石 油 化 工 高 等 学 校 学 报 Vo1 ． 9No．2 1996~ 6月 JOURNAL OF PETROCHEM ICAL ...
- 通用求根算法zeroin_Modern Robotics运动学数值解法及SVD算法（C matlab）
前言 原著之前CSDN已经注销,新CSDN Galaxy_Robot的博客_CSDN博客-机器人,C语言,我是谁?领域博主blog.csdn.net 这半个月的业余时间研究了机器人逆运动学的解析解法 ...
- [常微分方程的数值解法系列五] 龙格-库塔(RK4)法
龙格-库塔法 简介 基本思想 具体方法 一阶 二阶 求解参数 特殊二阶 三阶 高阶 步长选择 例子 在惯性导航以及VIO等实际问题中利用IMU求解位姿需要对IMU测量值进行积分得到需要的位置和姿态,其 ...
- python计算机器人运动学分析_V-rep学习笔记：机器人逆运动学数值解法（The Jacobian Transpose Method）...
机器人运动学逆解的问题经常出现在动画仿真和工业机器人的轨迹规划中:We want to know how the upper joints of the hierarchy would rotate ...
- 二阶边值问题的数值解matlab,《二阶常微分方程边值问题的数值解法》-毕业论文.doc...
w 摘 要 本文主要研究二阶常微分方程边值问题的数值解法.对线性边值问题,我们总结了两类常用的数值方法,即打靶法和有限差分方法,对每种方法都列出了详细的计算步骤和Matlab程序代码,通过具体的算例对 ...
- 一阶常微分方程的数值解法（二阶显式、隐式 Adams 公式及 Milne 方法）
一阶常微分方程的数值解法 这里我们介绍二阶显式.隐式 Adams 公式及 Milne 方法求解方程. 题目: 对初值问题 u ′ = u − t 2 , 0 ≤ t ≤ 1 , u ( 0 ) = 0 ...
最新文章
- python键_在Python中创建键命令
- Lucene系列-facet--转
- python(matplotlib5)——Contours 等高线图
- 用于稠密检索的无监督领域适应方法—Generative Pseudo Labeling (GPL)
- LeetCode 460. LFU缓存（哈希双链表）
- C语言课后习题(60)
- Python线程类首先是一个类
- codeforces 1221 A B C D
- [转载] Python中Numpy基础
- error: 'vector' does not name a type
- js constructor 和 instanceof
- 白话CSS3的新特性
- 【新手教程】从零搭建php动态网站
- python控制电机转动_Micropython TurnipBit 旋转按钮控制直流电机转速（儿时记忆中的吊扇）...
- java十进制二进制之间的互相转换
- Java中String使用及分析（UTF-8简单编码/解码器实现）
- 现代战争——僵尸网络的历史 上篇
- 测试人收入情况大曝光，你的收入在什么水平
- Eclipse官网地址
- 计算机台式机硬盘,台式电脑硬盘和笔记本硬盘有什么区别【详解】