前面我們已經知道對于線性方程組，一般有兩種數值解法：直接法和迭代法。直接法前面已經寫過了，沒看的同學可以移步閱讀：直接法。本次主要講述迭代法及其相應的MATLAB代碼。

考慮線性方程組

當

為低階稠密陣時一般采取直接法進行求解。但是我們在工程技術中經常遇到的是一些大型稀疏矩陣方程組( 的階數很大，但 元素較多，解PDE(Partial Differential Equation)時常遇到)，此時利用迭代法進行求解是合適的。在計算機內存和運算兩方面，迭代法通常都可利用 中有大量 元素的特點。

本文主要介紹迭代法的基本思路和雅可比(Jacobi)迭代法、高斯-塞德爾(Gauss-Seidel)迭代法、(逐次)超松弛(SOR, Successie over Relaxation)迭代法、共軛梯度(CG, Conjugate Gradient)法。

1. 迭代法的基本思路

對于線性方程組

為了使用迭代法對其進行求解，我們首先需要構造迭代序列，我們將上式改寫為

其中

和 都是已知的，從而我們就可以構造出一個迭代序列為

對向量序列

一組初始值 ，如果上述迭代序列收斂(有課本定理知其充要條件為 )，最終向量序列 將逐步逼近原方程組的精確解 。迭代次數越多，計算誤差就越小，因此，任意給定一個精度，都可以通過迭代計算出滿足精度要求的解(理論上，實際情況還需考慮計算機存儲字長)。

這就是我們構造迭代法的基本思路。MATLAB代碼為

function [x,t,it] = Iteration(A,b,I,eps) % 輸入： % A: 系數矩陣 % b: 載荷矩陣 % I: 最大迭代次數 % 輸出： % x: 解矩陣 % t: 時間 % it: 迭代次數 % 迭代初值默認為 0if nargin < 4eps = 1e-6; endtic [n,~] = size(A); B = zeros(size(A)); x = zeros(n,1); f = zeros(n,1);for r = 1:nfor l = 1:nif l == rB(r,r) = 0;elseB(r,l) = -A(r,l)/A(r,r);endendf(r) = b(r)/A(r,r); endx_exact = Ab; it = 1; for k = 1:I-1x = B*x+f;if norm(x-x_exact)>epsit = it+1;end end t = toc; end

示例：

clear A = [8 -3 2; 4 11 -1; 6 3 12]; b = [20; 33; 36]; I = 1000; [x,t,it] = Iteration(A,b,I)

運行結果

x =321t =0.016490086000000it =16

可見只是迭代

次，就滿足默認精度( )。

2. 雅可比(Jacobi)迭代法

其實上面第1節的代碼使用的方法就是最原始的雅可比迭代法，只是在形式上看起來不那么規范。下面我們來看規范的解

的雅可比迭代法的形式：

其中

， ，這里 分別為對角陣、下三角陣、上三角陣。并稱 為解 的雅可比迭代法的迭代矩陣。

其分量計算公式為

但實際編程計算一般采取簡潔且高效的矩陣形式，其MATLAB代碼為

function [x,t,it] = Jacobi(A,b,I,eps) % 雅可比迭代法 % 輸入： % A: 系數矩陣 % b: 載荷矩陣 % I: 最大迭代次數 % 輸出： % x: 解矩陣 % t: 時間 % it: 迭代次數 % 迭代初值默認為 0if nargin < 4eps = 1e-6; endtic [n,~] = size(A); x = zeros(n,1);D = diag(diag(A)); %求 A 的對角矩陣 L = -tril(A,-1); %求 A 的下三角矩陣,不帶對角線 U = -triu(A,1); %求 A 的上三角矩陣 J = D(L+U); f = Db;x_exact = Ab; it = 1; for k = 1:I-1x = J*x+f;if norm(x-x_exact)>epsit = it+1;end endt = toc; end

仍然計算第一節的例子，只需把函數名改一下：

clear A = [8 -3 2; 4 11 -1; 6 3 12]; b = [20; 33; 36]; I = 1000; [x,t,it] = Jacobi(A,b,I)

運行結果也與第1節一致：

x =321t =0.037686548000000it =16

3. 高斯-塞德爾(Gauss-Seidel)迭代法

高斯-塞德爾迭代法與雅可比迭代法很相似，可以看成是其的一種改進。其形式為：

其中

， ，這里 分別為對角陣、下三角陣、上三角陣。并稱 為解 的高斯-塞德爾迭代法的迭代矩陣。

其分量計算公式為

其MATLAB代碼亦采取矩陣形式：

function [x,t,it] = GaussSeidel(A,b,I,eps) % 高斯-賽德爾(Gauss-Seidel)迭代法 % 輸入： % A: 系數矩陣 % b: 載荷矩陣 % I: 最大迭代次數 % 輸出： % x: 解矩陣 % t: 時間 % it: 迭代次數 % 迭代初值默認為 0if nargin < 4eps = 1e-6; endtic [n,~] = size(A); x = zeros(n,1);D = diag(diag(A)); %求 A 的對角矩陣 L = -tril(A,-1); %求 A 的下三角矩陣,不帶對角線 U = -triu(A,1); %求 A 的上三角矩陣 G = (D-L)U; f = (D-L)b;x_exact = Ab; it = 1; for k = 1:I-1x = G*x+f;if norm(x_exact-x)>epsit = it+1;end endt = toc; end

同樣計算前面的例子，其結果為：

x =321t =0.054796674000000it =8

從結果可以看出，對于同樣的精度要求，高斯-塞德爾迭代法只需要

次迭代(雅可比迭代法需要 次)，收斂明顯要更快。

3. (逐次)超松弛(SOR)迭代法

簡單來說，逐次超松弛迭代法就是在前面高斯-塞德爾迭代法的基礎上加了松弛因子

，所以解 的SOR方法為

，初始向量，

其中

， 。

其分量計算公式為

顯然，當

時，SOR方法即為高斯-塞德爾迭代法。當 時，稱為超松弛法；當 時，稱為低松弛法。

其MATLAB代碼為

function [x,t,it,w] = SOR(A,b,I,eps,w) % 逐次超松馳迭代法(successive over relaxation method)迭代法 % 輸入： % A: 系數矩陣 % b: 載荷矩陣 % I: 最大迭代次數 % w: 松弛因子(w=1 時即為 Gauss-Seidel 迭代法) % 輸出： % x: 解矩陣 % t: 時間 % it: 迭代次數 % w: 松弛因子 % 迭代初值默認為 0tic [n,~] = size(A); x = zeros(n,1);D = diag(diag(A)); %求 A 的對角矩陣 L = -tril(A,-1); %求 A 的下三角矩陣,不帶對角線 U = -triu(A,1); %求 A 的上三角矩陣 w_opt = 2/(1+sqrt(1-(vrho(D(L+U)))^2)); % 最佳松弛因子 if nargin < 4eps = 1e-6;w = w_opt; end if nargin < 5w = w_opt; end Lw = (D-w*L)((1-w)*D+w*U); f = w*((D-w*L)b);x_exact = Ab; it = 1; for k = 1:I-1x = Lw*x+f;if norm(x-x_exact)>epsit = it+1;end endt = toc; end

還是計算前面那個例子，選擇精度為

，松弛因子 ：clear A = [8 -3 2; 4 11 -1; 6 3 12]; b = [20; 33; 36]; I = 1000; [x,t,it,w] = SOR(A,b,I,1e-6,1.2)

運行結果為

x =3.0000000000000002.0000000000000001.000000000000000t =0.003814113000000it =13w =1.200000000000000

如果采取默認松弛因子(自動選擇最佳松弛因子)：

clear A = [8 -3 2; 4 11 -1; 6 3 12]; b = [20; 33; 36]; I = 1000; [x,t,it,w] = SOR(A,b,I)

其結果為

x =3.0000000000000002.0000000000000001.000000000000000t =0.009697537000000it =8w =1.034531942537068

此時迭代次數減少為

次，松弛因子也選擇為接近 的數，與高斯-塞德爾迭代法效果基本一致。

4. 共軛梯度(CG)法

也稱共軛斜量法，它是一種變分方法，對應于求一個二次函數的極值。CG方法是一種求解大型稀疏對稱正定方程組十分有效的方法。

其算法描述為：

（1）任取

，計算 ，取 。

（2）對

，計算

（3）若

，或 ，計算停止，則 。由于 正定，故當 時， ，而 ，也即 。

寫成MATLAB代碼為

function [x,t,it] = CG(A,b) % 共軛梯度法 CG(conjugate gradient) % 針對大型稀疏對稱正定矩陣方程組 % 輸入： % A: 系數矩陣 % b: 載荷矩陣 % 輸出： % x: 解矩陣 % t: 時間 % it: 迭代次數tic [n,~] = size(A); x = zeros(n,1); r = b - A*x; p = r; it = 0;while (sum(r.*r) ~= 0) && (sum(p.*(A*p)) ~= 0)it = it+1;rr = sum(r.*r);alpha = rr/sum(p.*(A*p));x = x + alpha*p;r = r - alpha*A*p;beta = sum(r.*r)/rr;p = r + beta*p; end t = toc; end

用此CG方法求解一個

階Hilbert矩陣：clear n = 6; A = zeros(n,n); % Hilbert 矩陣 for i = 1:nfor j = 1:nA(i,j) = 1/(i+j-1);end end x_star = ones(n,1); b = A*x_star; [x,t,it] = CG(A,b)

運行結果為：

x =0.9999999999998971.0000000000040020.9999999999680481.0000000000921610.9999999998907641.000000000045352t =0.006673698000000it =350

而如果采用Gauss-Seidel迭代法，精度要求為

，其結果為x =0.9999999999998711.0000000000043290.9999999999676931.0000000000897200.9999999998961901.000000000042388t =2.914480589000000it =5438741

對比結果可以發現，達到類似精度，G-S迭代法需要迭代

次迭代，是CG方法的近 倍，足以見得在處理大型稀疏正定矩陣時CG方法的優勢之明顯。

總結

以上是生活随笔為你收集整理的matlab 矩阵jocobi迭代_第6章 解线性方程组的迭代法(基于MATLAB)的全部內容，希望文章能夠幫你解決所遇到的問題。

如果覺得生活随笔網站內容還不錯，歡迎將生活随笔推薦給好友。