全拆解)
簡介東北大學理學院數(shù)學系數(shù)值分析課程實驗報告以解線性方程組的迭代法與迭代格式比較為主線系統(tǒng)介紹了Jacobi迭代法、Gauss-Seidel迭代法與SOR松弛迭代法。報告嚴格按學校規(guī)范撰寫包含實驗目的、原理、內(nèi)容、數(shù)據(jù)記錄和處理、結(jié)果分析等完整環(huán)節(jié)方便高校學生對照學習或完善課程實驗報告。壓縮包僅含1個doc文檔大小186KB結(jié)構(gòu)緊湊。文中給出了三種算法的MATLAB實現(xiàn)從矩陣分解、迭代公式到誤差閾值控制均有對應代碼并記錄了不同初始向量與松弛因子ω下的迭代次數(shù)和收斂誤差例如ω從1.10到1.90時迭代次數(shù)由11次增至11828次能直觀展示松弛參數(shù)對收斂速度的影響。已有2614人學習下載可作為理解迭代法原理、查看算法實現(xiàn)細節(jié)及撰寫實驗報告的重要參考。1. 一份讓你少熬夜的數(shù)值分析實驗報告三種迭代法從原理到 MATLAB 代碼全拆解如果你正在數(shù)值分析課上被“解線性方程組的迭代法”折騰大概率經(jīng)歷過這樣的場景講義上的公式一看就會自己寫 MATLAB 一跑就死循環(huán)或者迭代次數(shù)忽大忽小根本不知道問題出在哪。這份實驗報告正好是解藥——它把 Jacobi 迭代法、Gauss-Seidel 迭代法和 SOR 松弛迭代法從原理、代碼到實驗數(shù)據(jù)完整整理了一遍連不同松弛因子下迭代次數(shù)從 11 次到 11828 次的對比都給你算好了。適合三類人正在做類似課程實驗、需要復現(xiàn)數(shù)據(jù)寫報告的同學準備面試或考研復試、想快速理解迭代法實現(xiàn)細節(jié)的人以及在實際工程里需要手寫迭代求解器、想找個靠譜參考的開發(fā)者。這份資源不只是一堆代碼更是一份能直接跑通、能照著改的完整實驗樣本。2. 實驗設(shè)計解析線性方程組迭代法的收斂條件與矩陣分裂2.1 迭代法為什么可以解線性方程組直接法用高斯消元或 LU 分解一次性求出精確解迭代法則把一個Ax b的問題改寫成固定點形式x Bx f然后從一個初始向量出發(fā)反復計算x_{k1} Bx_k f直到前后兩步的差足夠小。這里最關(guān)鍵的是收斂條件迭代矩陣B的譜半徑必須小于 1也就是所有特征值的模都小于 1。如果譜半徑大于等于 1迭代就會發(fā)散或者原地打轉(zhuǎn)這也是很多同學實驗報告里出現(xiàn)“死循環(huán)”的根本原因。這份實驗報告選用的測試矩陣是三階的A [4, -1, 1; 4, -8, 1; -2, 1, 5]; b [7; -21; 15];它不是嚴格對角占優(yōu)矩陣第二行對角線是 -8但該行其他元素的絕對值之和是 5滿足對角占優(yōu)第一行對角線 4 比 2 大第三行 5 比 3 大實際上整體是嚴格對角占優(yōu)的所以 Jacobi 和 Gauss-Seidel 都收斂。實驗報告特意選了這樣一個行列式不為零、迭代性質(zhì)又不太極端的例子非常合適初學者觀察不同方法的收斂速度差異。2.2 矩陣分解D、L、U 在 MATLAB 里怎么取迭代法把矩陣A拆成三部分對角部分D、嚴格下三角部分L和嚴格上三角部分U。注意這里L和U的方向取決于你自己的定義報告里的定義是D diag(diag(A)); % 對角矩陣 L -tril(A, -1); % 下三角矩陣取負 U -triu(A, 1); % 上三角矩陣取負為什么取負因為后續(xù)迭代公式里需要的是A D - L - U。tril(A, -1)返回 A 中主對角線以下的部分但不含主對角線triu(A, 1)返回主對角線以上的部分。如果你直接用tril(A, -1)得到的是一個“帶正號的下三角”那么A D L U就會寫錯正負號。這個細節(jié)是很多實現(xiàn)翻車的起點后面避坑章節(jié)還會再提到。以這個測試矩陣為例D是對角線[4, -8, 5]L是[0,0,0; -4,0,0; 2,-1,0]的相反數(shù)U是[0,1,-1; 0,0,-1; 0,0,0]的相反數(shù)。寫出迭代公式后就明白為什么 MATLAB 代碼里到處是inv(D)、inv(D-L)。2.3 三種迭代格式的矩陣表達與收斂性Jacobi 迭代公式的核心是“用舊值計算新值”X_1 inv(D) * (L U) * X_0 inv(D) * b;把A D - L - U代入Ax b得到Dx (LU)x b所以迭代矩陣是inv(D)*(LU)。每次迭代都基于上一輪的整個向量X_0所以各個分量可以并行計算。Gauss-Seidel 迭代的關(guān)鍵是“用最新值計算下一個分量”X_1 inv(D - L) * U * X_0 inv(D - L) * b;這里把下三角部分并入左側(cè)求逆等價于在計算第 i 個分量時已經(jīng)使用本輪更新過的前 i-1 個分量。Gauss-Seidel 通常比 Jacobi 收斂快實驗數(shù)據(jù)也印證了初始向量為零時Jacobi 需要 14 次Gauss-Seidel 只要 9 次。SOR 是 Gauss-Seidel 的加速版公式里多了一個松弛因子ω代碼是k 1.75; % 松弛因子 omega X_1 inv(D - k*L) * ((1-k)*D k*U) * X_0 k * inv(D - k*L) * b;當ω1時 SOR 退化成 Gauss-Seidelω1叫超松弛通常用于加速收斂ω1叫低松弛有時用于改善穩(wěn)定性。理論要求0 ω 2超出這個區(qū)間大概率發(fā)散。2.4 SOR 松弛因子的“玄學”取值不同迭代次數(shù)從 11 次變成 11828 次實驗報告里最有意思的部分是 SOR 松弛因子掃描。初始向量為零向量時ω 值迭代次數(shù)1.10111.20141.30161.40191.50211.60321.70551.801531.9011828這個數(shù)據(jù)我第一次看到時也嚇了一跳。1.9 和 2.0 只差 0.1迭代次數(shù)卻從一百多次暴漲到上萬次。原因在于 SOR 的最優(yōu)松弛因子和系數(shù)矩陣的特征值分布強相關(guān)當ω接近上界 2 時迭代矩陣的譜半徑會無限接近 1收斂速度變得極慢。這不是代碼 bug而是數(shù)值方法本身的特性。這也是為什么實驗報告的數(shù)據(jù)比單純跑一次代碼有價值——它幫你提前劃出了“危險區(qū)間”。3. 核心代碼拆解一個迭代循環(huán)從初始化到誤差收斂的完整鏈路3.1 原始代碼逐行解讀矩陣、初值、迭代步報告的 Jacobi 部分雖然不長但每一行都值得細看A [4,-1,1; 4,-8,1; -2,1,5]; b [7;-21;15]; D diag(diag(A)); % 對角矩陣 L -tril(A,-1); % 下三角矩陣取負 U -triu(A,1); % 上三角矩陣取負 X_0 [0;0;0]; % 初始向量 X_1 inv(D) * (L U) * X_0 inv(D) * b; % 初次迭代 dis abs(X_1 - X_0); % 絕對誤差向量 p 1; for i 1:3 while(dis(i) 10^(-5)) X_1 inv(D) * (L U) * X_0 inv(D) * b; dis abs(X_1 - X_0); X_0 X_1; p p 1; end end這里p是迭代次數(shù)計數(shù)器。初始迭代執(zhí)行了一次所以p從 1 開始。for i 1:3加while的組合很有意思它并不是一次性判斷所有分量是否滿足誤差而是分別檢查每個分量。我們稍后會說這種寫法的隱患但先理解執(zhí)行流程當?shù)?1 個分量誤差大于閾值時不斷迭代直到它小于閾值然后跳出 while進入 i2再讓第 2 個分量滿足最后 i3。因為每輪迭代都會更新整個向量所以這個流程實際上會讓所有分量最終都滿足但可能比統(tǒng)一的范數(shù)判斷多算一些迭代。3.2 誤差判斷的寫法按分量判斷和按范數(shù)判斷的差別實驗報告里的誤差判斷是while(dis(i) 10^(-5))只觀察某個分量的絕對誤差。常見做法是計算整個向量的無窮范數(shù)或 2 范數(shù)while max(abs(X_1 - X_0)) 1e-5 X_1 inv(D) * (L U) * X_0 inv(D) * b; dis abs(X_1 - X_0); X_0 X_1; p p 1; end兩種寫法在大多數(shù)情況下結(jié)果一致但按分量判斷存在一個潛在問題如果某個分量在迭代中波動比如先收斂又變壞外層 for 循環(huán)已經(jīng)越過了它就不會再檢查可能最終這個分量的誤差其實略超閾值。對這個三階測試矩陣而言不會出問題但換成更復雜的矩陣時就可能“表面收斂、實際沒收斂”。所以我通常建議改用max(abs(...))或norm(...)作為退出條件。3.3 實驗數(shù)據(jù)三種方法在不同初始向量下的迭代次數(shù)這份實驗報告記錄了三組初始向量的結(jié)果這部分數(shù)據(jù)可以直接用到你自己的報告里初始向量[0;0;0]Jacobi 14 次Gauss-Seidel 9 次SORω1.10 時 11 次ω1.50 時 21 次ω1.75 時 80 次初始向量[1;1;1]Jacobi 13 次Gauss-Seidel 9 次SORω1.25 時 14 次ω1.50 時 20 次ω1.75 時 79 次初始向量[-1;-1;-1]Jacobi 14 次Gauss-Seidel 9 次SORω1.25 時 15 次ω1.50 時 21 次ω1.75 時 87 次可以看出Gauss-Seidel 對初始向量不敏感始終 9 次收斂Jacobi 基本穩(wěn)定在 13~14 次SOR 則對 ω 和初值都很敏感。這個數(shù)據(jù)告訴我們一個實際經(jīng)驗如果你的問題是線性的Gauss-Seidel 是性價比很高的選擇SOR 的加速效果需要謹慎調(diào)參調(diào)不好反而更慢。3.4 一個隱藏問題inv 求逆的效率與數(shù)值穩(wěn)定性報告代碼里反復出現(xiàn)inv(D)、inv(D-L)、inv(D-k*L)。在小規(guī)模三階矩陣下毫無問題但如果你把這個代碼改造成大規(guī)模稀疏矩陣求解器inv會先把矩陣變成稠密矩陣內(nèi)存和計算量瞬間爆炸。更穩(wěn)的做法是改用左除運算符X_1 (D - L) \ (U * X_0) (D - L) \ b;\對稀疏矩陣會走 LU 分解或迭代求解數(shù)值穩(wěn)定性也更好。這是從“能跑”到“能用在工程里”的關(guān)鍵一步二次開發(fā)章節(jié)會繼續(xù)展開。4. 運行與復現(xiàn)在 MATLAB 里把三種解法跑通并畫出收斂曲線4.1 環(huán)境準備MATLAB 或 Octave腳本文件怎么建最簡單的方式是打開 MATLAB新建一個腳本命名為iterative_demo.m把三種方法的代碼分三段貼進去然后點運行。如果沒有 MATLAB 授權(quán)用 GNU Octave 也能跑通兩者語法在這份報告里幾乎完全兼容。需要注意的一點是文件名不要用中文也不要和內(nèi)置函數(shù)重名比如別叫sum.m或plot.m。復制代碼時最容易踩的坑是MATLAB 的矩陣換行用分號但如果你從網(wǎng)頁復制可能把全角分號或空格粘進去。建議先檢查矩陣行的分隔符是不是半角的;或換行。4.2 跑通 Jacobi 與 Gauss-Seidel注意看迭代次數(shù) p 的累加邏輯運行 Jacobi 部分后命令窗口不會直接顯示迭代次數(shù)需要自己加一行輸出disp([Jacobi 迭代次數(shù): , num2str(p)]);我在復現(xiàn)時習慣把p的初值改成 0然后讓循環(huán)里的p p 1統(tǒng)一計數(shù)。因為原報告在 while 之前已經(jīng)做了一次迭代p初始為 1這樣最終得到的數(shù)字才和實驗數(shù)據(jù)對得上。Gauss-Seidel 部分的思路完全相同只是迭代矩陣變成了inv(D-L)*U。跑完你應該看到p14和p9。4.3 批量測試不同松弛因子用 for 循環(huán)復現(xiàn) SOR 實驗數(shù)據(jù)實驗報告里的 ω 掃描是手寫一個個改的但復現(xiàn)時更推薦用循環(huán)。下面這段代碼可以直接復現(xiàn)表里的數(shù)據(jù)A [4,-1,1; 4,-8,1; -2,1,5]; b [7;-21;15]; D diag(diag(A)); L -tril(A,-1); U -triu(A,1); omegas [1.10,1.20,1.30,1.40,1.50,1.60,1.70,1.80,1.90]; for k omegas X_0 [0;0;0]; X_1 inv(D - k*L) * ((1-k)*D k*U) * X_0 k * inv(D - k*L) * b; dis abs(X_1 - X_0); p 1; while(dis(1) 1e-5 || dis(2) 1e-5 || dis(3) 1e-5) X_1 inv(D - k*L) * ((1-k)*D k*U) * X_0 k * inv(D - k*L) * b; dis abs(X_1 - X_0); X_0 X_1; p p 1; end fprintf(omega%.2f, 迭代次數(shù)%d\n, k, p); end注意這里我把 while 條件改成“任一分量不滿足就繼續(xù)”避免原始代碼里for while結(jié)構(gòu)可能帶來的漏檢問題。輸出結(jié)果會和實驗報告基本一致ω1.10 時 11 次ω1.90 時 11828 次。跑最后一條時要有耐心可能等幾秒。4.4 修改初始向量和誤差閾值觀察收斂速度的變化把初始向量改成[1;1;1]只需要改X_0那一行。改成[-1;-1;-1]同理。誤差閾值10^(-5)也可以改成10^(-8)你會看到迭代次數(shù)普遍增加。如果改得太小比如10^(-15)三階矩陣用 double 浮點數(shù)就可能無法達到循環(huán)會進入無限震蕩。這個現(xiàn)象可以當作實驗報告里的“延伸思考”來寫。如果想要更直觀可以加一個誤差記錄數(shù)組把每次迭代的max(abs(X_1-X_0))存下來最后用semilogy畫收斂曲線。這個技巧放到第六章。5. 避坑指南迭代法實驗里最容易翻車的 5 個細節(jié)5.1 現(xiàn)象while 循環(huán)直接卡死跑幾分鐘都不停原因迭代矩陣譜半徑大于等于 1或者誤差閾值設(shè)成了不可達的數(shù)。最常見的是代碼里L和U的正負號寫反導致迭代矩陣根本不是預期那個。解決先打印D、L、U看符號再計算迭代矩陣的特征值比如max(abs(eig(inv(D)*(LU))))如果譜半徑大于等于 1就不要再跑了。這是我做迭代法實驗最深的血淚教訓——不查譜半徑就盲目跑循環(huán)等于讓電腦陪你死等。5.2 現(xiàn)象Gauss-Seidel 結(jié)果和 Jacobi 幾乎一樣原因矩陣分解時沒有取負導致D-L和Dtril(A,-1)混用。很多同學抄代碼時把L -tril(A,-1)寫成L tril(A,-1)那么D-L實際上變成了D 下三角根本不符合推導。解決統(tǒng)一用A D - L - U的約定打印A - (D - L - U)看是否為零矩陣。這個檢查成本極低能避免一整晚的困惑。5.3 現(xiàn)象命令窗口報錯“Matrix must be square”或“奇異矩陣”原因inv(D-k*L)在某個松弛因子下矩陣不可逆。理論上只要 A 非奇異D-kL通常也可逆但如果你把 ω 設(shè)成負數(shù)或大于 2這個矩陣可能變得病態(tài)。解決在循環(huán)里用 try-catch 捕獲奇異矩陣錯誤或者改用左除\左除對近奇異矩陣會給出警告但至少能返回一個結(jié)果。另外檢查一下是不是把矩陣的一行數(shù)據(jù)輸錯了比如[4,-1,1;4,-8,1;-2,1,5]中間的分號是否被替換成了逗號。5.4 現(xiàn)象松弛因子 1.9 時迭代次數(shù)跳到 11828 次這是 bug 嗎不是 bug是 SOR 方法的固有現(xiàn)象。當 ω 接近上界 2 時迭代矩陣譜半徑趨近 1收斂速度趨近于零。遇到這種情況不要慌張把它寫進實驗分析里反而是加分項。解決如果你想看到完整的收斂過程可以把誤差閾值放大到10^(-3)迭代次數(shù)會顯著減少也能觀察到誤差緩慢下降的曲線。如果只是驗證代碼正確性建議只測到 1.8。5.5 現(xiàn)象用 Octave 運行結(jié)果和 MATLAB 不一樣原因Octave 的浮點數(shù)運算和 MATLAB 基本一致但某些內(nèi)置函數(shù)如inv、eig底層實現(xiàn)有細微差異在接近收斂邊界時可能表現(xiàn)出不同迭代次數(shù)。另外循環(huán)里的變量作用域和顯示精度也可能不同。解決關(guān)鍵結(jié)果不要直接用比較而是用abs(p1 - p2) 1容忍。如果迭代次數(shù)差很多檢查是不是用了不同的數(shù)據(jù)文件或沒清空工作區(qū)。6. 二次開發(fā)把實驗代碼改造成通用的迭代法求解器6.1 封裝成函數(shù)支持 method 參數(shù)和松弛因子原報告的代碼是平鋪的每換一個方法就要復制改一遍。更好的做法是封裝成一個函數(shù)function [x, iter, err_history] iterative_solve(A, b, method, omega, tol, x0) % method: jacobi, gs, sor % omega: 僅 SOR 使用 % tol: 迭代終止閾值 % x0: 初始向量默認零向量 D diag(diag(A)); L -tril(A, -1); U -triu(A, 1); n size(A, 1); x x0; err_history []; iter 0; while true switch method case jacobi x_new D \ ((L U) * x) D \ b; case gs x_new (D - L) \ (U * x) (D - L) \ b; case sor x_new (D - omega*L) \ (((1-omega)*D omega*U) * x) ... omega * ((D - omega*L) \ b); otherwise error(未知方法); end err max(abs(x_new - x)); err_history [err_history; err]; x x_new; iter iter 1; if err tol || iter 10000 break; end end end這個函數(shù)把三種方法統(tǒng)一到了一套switch結(jié)構(gòu)里最大的好處是測試不同方法時不需要反復改代碼。我把原來的 while 條件改成了while true靠迭代上限和誤差閾值雙重退出避免死循環(huán)。6.2 增加誤差記錄與收斂曲線在函數(shù)里我把每次迭代的誤差存進err_history這樣調(diào)用后可以直接畫收斂曲線[x, iter, his] iterative_solve(A, b, sor, 1.1, 1e-5, [0;0;0]); semilogy(1:length(his), his, -o); xlabel(迭代次數(shù)); ylabel(誤差對數(shù)坐標);用semilogy畫誤差下降曲線能直觀看到 SOR 在不同 ω 下收斂速度的差異。你會發(fā)現(xiàn) ω1.9 的曲線下降得非常平緩這就是 11828 次的來源。6.3 驗證方法用嚴格對角占優(yōu)矩陣和譜半徑判斷改造完之后必須驗證函數(shù)沒有跑偏。我一般會做兩步第一步構(gòu)造一個嚴格對角占優(yōu)矩陣比如三對角矩陣用直接法A\b的結(jié)果和迭代解對比第二步計算迭代矩陣的譜半徑確認小于 1。另外再跑一遍原始實驗數(shù)據(jù)確保 Jacobi 返回 14 次、Gauss-Seidel 返回 9 次。驗證通過后再去調(diào) ω就心中有數(shù)了。從那以后我每次做迭代法實驗都會先寫一個iterative_solve這樣的通用函數(shù)再批量跑參數(shù)而不是復制粘貼三次代碼。這個習慣幫我少走了很多彎路希望你也能直接用上。希望幫到你。本文還有配套的精品資源點擊獲取 解線性方程組的三種迭代法原理與MATLAB實現(xiàn)全拆解
簡介東北大學理學院數(shù)學系數(shù)值分析課程實驗報告以解線性方程組的迭代法與迭代格式比較為主線系統(tǒng)介紹了Jacobi迭代法、Gauss-Seidel迭代法與SOR松弛迭代法。報告嚴格按學校規(guī)范撰寫包含實驗目的、原理、內(nèi)容、數(shù)據(jù)記錄和處理、結(jié)果分析等完整環(huán)節(jié)方便高校學生對照學習或完善課程實驗報告。壓縮包僅含1個doc文檔大小186KB結(jié)構(gòu)緊湊。文中給出了三種算法的MATLAB實現(xiàn)從矩陣分解、迭代公式到誤差閾值控制均有對應代碼并記錄了不同初始向量與松弛因子ω下的迭代次數(shù)和收斂誤差例如ω從1.10到1.90時迭代次數(shù)由11次增至11828次能直觀展示松弛參數(shù)對收斂速度的影響。已有2614人學習下載可作為理解迭代法原理、查看算法實現(xiàn)細節(jié)及撰寫實驗報告的重要參考。1. 一份讓你少熬夜的數(shù)值分析實驗報告三種迭代法從原理到 MATLAB 代碼全拆解如果你正在數(shù)值分析課上被“解線性方程組的迭代法”折騰大概率經(jīng)歷過這樣的場景講義上的公式一看就會自己寫 MATLAB 一跑就死循環(huán)或者迭代次數(shù)忽大忽小根本不知道問題出在哪。這份實驗報告正好是解藥——它把 Jacobi 迭代法、Gauss-Seidel 迭代法和 SOR 松弛迭代法從原理、代碼到實驗數(shù)據(jù)完整整理了一遍連不同松弛因子下迭代次數(shù)從 11 次到 11828 次的對比都給你算好了。適合三類人正在做類似課程實驗、需要復現(xiàn)數(shù)據(jù)寫報告的同學準備面試或考研復試、想快速理解迭代法實現(xiàn)細節(jié)的人以及在實際工程里需要手寫迭代求解器、想找個靠譜參考的開發(fā)者。這份資源不只是一堆代碼更是一份能直接跑通、能照著改的完整實驗樣本。2. 實驗設(shè)計解析線性方程組迭代法的收斂條件與矩陣分裂2.1 迭代法為什么可以解線性方程組直接法用高斯消元或 LU 分解一次性求出精確解迭代法則把一個Ax b的問題改寫成固定點形式x Bx f然后從一個初始向量出發(fā)反復計算x_{k1} Bx_k f直到前后兩步的差足夠小。這里最關(guān)鍵的是收斂條件迭代矩陣B的譜半徑必須小于 1也就是所有特征值的模都小于 1。如果譜半徑大于等于 1迭代就會發(fā)散或者原地打轉(zhuǎn)這也是很多同學實驗報告里出現(xiàn)“死循環(huán)”的根本原因。這份實驗報告選用的測試矩陣是三階的A [4, -1, 1; 4, -8, 1; -2, 1, 5]; b [7; -21; 15];它不是嚴格對角占優(yōu)矩陣第二行對角線是 -8但該行其他元素的絕對值之和是 5滿足對角占優(yōu)第一行對角線 4 比 2 大第三行 5 比 3 大實際上整體是嚴格對角占優(yōu)的所以 Jacobi 和 Gauss-Seidel 都收斂。實驗報告特意選了這樣一個行列式不為零、迭代性質(zhì)又不太極端的例子非常合適初學者觀察不同方法的收斂速度差異。2.2 矩陣分解D、L、U 在 MATLAB 里怎么取迭代法把矩陣A拆成三部分對角部分D、嚴格下三角部分L和嚴格上三角部分U。注意這里L和U的方向取決于你自己的定義報告里的定義是D diag(diag(A)); % 對角矩陣 L -tril(A, -1); % 下三角矩陣取負 U -triu(A, 1); % 上三角矩陣取負為什么取負因為后續(xù)迭代公式里需要的是A D - L - U。tril(A, -1)返回 A 中主對角線以下的部分但不含主對角線triu(A, 1)返回主對角線以上的部分。如果你直接用tril(A, -1)得到的是一個“帶正號的下三角”那么A D L U就會寫錯正負號。這個細節(jié)是很多實現(xiàn)翻車的起點后面避坑章節(jié)還會再提到。以這個測試矩陣為例D是對角線[4, -8, 5]L是[0,0,0; -4,0,0; 2,-1,0]的相反數(shù)U是[0,1,-1; 0,0,-1; 0,0,0]的相反數(shù)。寫出迭代公式后就明白為什么 MATLAB 代碼里到處是inv(D)、inv(D-L)。2.3 三種迭代格式的矩陣表達與收斂性Jacobi 迭代公式的核心是“用舊值計算新值”X_1 inv(D) * (L U) * X_0 inv(D) * b;把A D - L - U代入Ax b得到Dx (LU)x b所以迭代矩陣是inv(D)*(LU)。每次迭代都基于上一輪的整個向量X_0所以各個分量可以并行計算。Gauss-Seidel 迭代的關(guān)鍵是“用最新值計算下一個分量”X_1 inv(D - L) * U * X_0 inv(D - L) * b;這里把下三角部分并入左側(cè)求逆等價于在計算第 i 個分量時已經(jīng)使用本輪更新過的前 i-1 個分量。Gauss-Seidel 通常比 Jacobi 收斂快實驗數(shù)據(jù)也印證了初始向量為零時Jacobi 需要 14 次Gauss-Seidel 只要 9 次。SOR 是 Gauss-Seidel 的加速版公式里多了一個松弛因子ω代碼是k 1.75; % 松弛因子 omega X_1 inv(D - k*L) * ((1-k)*D k*U) * X_0 k * inv(D - k*L) * b;當ω1時 SOR 退化成 Gauss-Seidelω1叫超松弛通常用于加速收斂ω1叫低松弛有時用于改善穩(wěn)定性。理論要求0 ω 2超出這個區(qū)間大概率發(fā)散。2.4 SOR 松弛因子的“玄學”取值不同迭代次數(shù)從 11 次變成 11828 次實驗報告里最有意思的部分是 SOR 松弛因子掃描。初始向量為零向量時ω 值迭代次數(shù)1.10111.20141.30161.40191.50211.60321.70551.801531.9011828這個數(shù)據(jù)我第一次看到時也嚇了一跳。1.9 和 2.0 只差 0.1迭代次數(shù)卻從一百多次暴漲到上萬次。原因在于 SOR 的最優(yōu)松弛因子和系數(shù)矩陣的特征值分布強相關(guān)當ω接近上界 2 時迭代矩陣的譜半徑會無限接近 1收斂速度變得極慢。這不是代碼 bug而是數(shù)值方法本身的特性。這也是為什么實驗報告的數(shù)據(jù)比單純跑一次代碼有價值——它幫你提前劃出了“危險區(qū)間”。3. 核心代碼拆解一個迭代循環(huán)從初始化到誤差收斂的完整鏈路3.1 原始代碼逐行解讀矩陣、初值、迭代步報告的 Jacobi 部分雖然不長但每一行都值得細看A [4,-1,1; 4,-8,1; -2,1,5]; b [7;-21;15]; D diag(diag(A)); % 對角矩陣 L -tril(A,-1); % 下三角矩陣取負 U -triu(A,1); % 上三角矩陣取負 X_0 [0;0;0]; % 初始向量 X_1 inv(D) * (L U) * X_0 inv(D) * b; % 初次迭代 dis abs(X_1 - X_0); % 絕對誤差向量 p 1; for i 1:3 while(dis(i) 10^(-5)) X_1 inv(D) * (L U) * X_0 inv(D) * b; dis abs(X_1 - X_0); X_0 X_1; p p 1; end end這里p是迭代次數(shù)計數(shù)器。初始迭代執(zhí)行了一次所以p從 1 開始。for i 1:3加while的組合很有意思它并不是一次性判斷所有分量是否滿足誤差而是分別檢查每個分量。我們稍后會說這種寫法的隱患但先理解執(zhí)行流程當?shù)?1 個分量誤差大于閾值時不斷迭代直到它小于閾值然后跳出 while進入 i2再讓第 2 個分量滿足最后 i3。因為每輪迭代都會更新整個向量所以這個流程實際上會讓所有分量最終都滿足但可能比統(tǒng)一的范數(shù)判斷多算一些迭代。3.2 誤差判斷的寫法按分量判斷和按范數(shù)判斷的差別實驗報告里的誤差判斷是while(dis(i) 10^(-5))只觀察某個分量的絕對誤差。常見做法是計算整個向量的無窮范數(shù)或 2 范數(shù)while max(abs(X_1 - X_0)) 1e-5 X_1 inv(D) * (L U) * X_0 inv(D) * b; dis abs(X_1 - X_0); X_0 X_1; p p 1; end兩種寫法在大多數(shù)情況下結(jié)果一致但按分量判斷存在一個潛在問題如果某個分量在迭代中波動比如先收斂又變壞外層 for 循環(huán)已經(jīng)越過了它就不會再檢查可能最終這個分量的誤差其實略超閾值。對這個三階測試矩陣而言不會出問題但換成更復雜的矩陣時就可能“表面收斂、實際沒收斂”。所以我通常建議改用max(abs(...))或norm(...)作為退出條件。3.3 實驗數(shù)據(jù)三種方法在不同初始向量下的迭代次數(shù)這份實驗報告記錄了三組初始向量的結(jié)果這部分數(shù)據(jù)可以直接用到你自己的報告里初始向量[0;0;0]Jacobi 14 次Gauss-Seidel 9 次SORω1.10 時 11 次ω1.50 時 21 次ω1.75 時 80 次初始向量[1;1;1]Jacobi 13 次Gauss-Seidel 9 次SORω1.25 時 14 次ω1.50 時 20 次ω1.75 時 79 次初始向量[-1;-1;-1]Jacobi 14 次Gauss-Seidel 9 次SORω1.25 時 15 次ω1.50 時 21 次ω1.75 時 87 次可以看出Gauss-Seidel 對初始向量不敏感始終 9 次收斂Jacobi 基本穩(wěn)定在 13~14 次SOR 則對 ω 和初值都很敏感。這個數(shù)據(jù)告訴我們一個實際經(jīng)驗如果你的問題是線性的Gauss-Seidel 是性價比很高的選擇SOR 的加速效果需要謹慎調(diào)參調(diào)不好反而更慢。3.4 一個隱藏問題inv 求逆的效率與數(shù)值穩(wěn)定性報告代碼里反復出現(xiàn)inv(D)、inv(D-L)、inv(D-k*L)。在小規(guī)模三階矩陣下毫無問題但如果你把這個代碼改造成大規(guī)模稀疏矩陣求解器inv會先把矩陣變成稠密矩陣內(nèi)存和計算量瞬間爆炸。更穩(wěn)的做法是改用左除運算符X_1 (D - L) \ (U * X_0) (D - L) \ b;\對稀疏矩陣會走 LU 分解或迭代求解數(shù)值穩(wěn)定性也更好。這是從“能跑”到“能用在工程里”的關(guān)鍵一步二次開發(fā)章節(jié)會繼續(xù)展開。4. 運行與復現(xiàn)在 MATLAB 里把三種解法跑通并畫出收斂曲線4.1 環(huán)境準備MATLAB 或 Octave腳本文件怎么建最簡單的方式是打開 MATLAB新建一個腳本命名為iterative_demo.m把三種方法的代碼分三段貼進去然后點運行。如果沒有 MATLAB 授權(quán)用 GNU Octave 也能跑通兩者語法在這份報告里幾乎完全兼容。需要注意的一點是文件名不要用中文也不要和內(nèi)置函數(shù)重名比如別叫sum.m或plot.m。復制代碼時最容易踩的坑是MATLAB 的矩陣換行用分號但如果你從網(wǎng)頁復制可能把全角分號或空格粘進去。建議先檢查矩陣行的分隔符是不是半角的;或換行。4.2 跑通 Jacobi 與 Gauss-Seidel注意看迭代次數(shù) p 的累加邏輯運行 Jacobi 部分后命令窗口不會直接顯示迭代次數(shù)需要自己加一行輸出disp([Jacobi 迭代次數(shù): , num2str(p)]);我在復現(xiàn)時習慣把p的初值改成 0然后讓循環(huán)里的p p 1統(tǒng)一計數(shù)。因為原報告在 while 之前已經(jīng)做了一次迭代p初始為 1這樣最終得到的數(shù)字才和實驗數(shù)據(jù)對得上。Gauss-Seidel 部分的思路完全相同只是迭代矩陣變成了inv(D-L)*U。跑完你應該看到p14和p9。4.3 批量測試不同松弛因子用 for 循環(huán)復現(xiàn) SOR 實驗數(shù)據(jù)實驗報告里的 ω 掃描是手寫一個個改的但復現(xiàn)時更推薦用循環(huán)。下面這段代碼可以直接復現(xiàn)表里的數(shù)據(jù)A [4,-1,1; 4,-8,1; -2,1,5]; b [7;-21;15]; D diag(diag(A)); L -tril(A,-1); U -triu(A,1); omegas [1.10,1.20,1.30,1.40,1.50,1.60,1.70,1.80,1.90]; for k omegas X_0 [0;0;0]; X_1 inv(D - k*L) * ((1-k)*D k*U) * X_0 k * inv(D - k*L) * b; dis abs(X_1 - X_0); p 1; while(dis(1) 1e-5 || dis(2) 1e-5 || dis(3) 1e-5) X_1 inv(D - k*L) * ((1-k)*D k*U) * X_0 k * inv(D - k*L) * b; dis abs(X_1 - X_0); X_0 X_1; p p 1; end fprintf(omega%.2f, 迭代次數(shù)%d\n, k, p); end注意這里我把 while 條件改成“任一分量不滿足就繼續(xù)”避免原始代碼里for while結(jié)構(gòu)可能帶來的漏檢問題。輸出結(jié)果會和實驗報告基本一致ω1.10 時 11 次ω1.90 時 11828 次。跑最后一條時要有耐心可能等幾秒。4.4 修改初始向量和誤差閾值觀察收斂速度的變化把初始向量改成[1;1;1]只需要改X_0那一行。改成[-1;-1;-1]同理。誤差閾值10^(-5)也可以改成10^(-8)你會看到迭代次數(shù)普遍增加。如果改得太小比如10^(-15)三階矩陣用 double 浮點數(shù)就可能無法達到循環(huán)會進入無限震蕩。這個現(xiàn)象可以當作實驗報告里的“延伸思考”來寫。如果想要更直觀可以加一個誤差記錄數(shù)組把每次迭代的max(abs(X_1-X_0))存下來最后用semilogy畫收斂曲線。這個技巧放到第六章。5. 避坑指南迭代法實驗里最容易翻車的 5 個細節(jié)5.1 現(xiàn)象while 循環(huán)直接卡死跑幾分鐘都不停原因迭代矩陣譜半徑大于等于 1或者誤差閾值設(shè)成了不可達的數(shù)。最常見的是代碼里L和U的正負號寫反導致迭代矩陣根本不是預期那個。解決先打印D、L、U看符號再計算迭代矩陣的特征值比如max(abs(eig(inv(D)*(LU))))如果譜半徑大于等于 1就不要再跑了。這是我做迭代法實驗最深的血淚教訓——不查譜半徑就盲目跑循環(huán)等于讓電腦陪你死等。5.2 現(xiàn)象Gauss-Seidel 結(jié)果和 Jacobi 幾乎一樣原因矩陣分解時沒有取負導致D-L和Dtril(A,-1)混用。很多同學抄代碼時把L -tril(A,-1)寫成L tril(A,-1)那么D-L實際上變成了D 下三角根本不符合推導。解決統(tǒng)一用A D - L - U的約定打印A - (D - L - U)看是否為零矩陣。這個檢查成本極低能避免一整晚的困惑。5.3 現(xiàn)象命令窗口報錯“Matrix must be square”或“奇異矩陣”原因inv(D-k*L)在某個松弛因子下矩陣不可逆。理論上只要 A 非奇異D-kL通常也可逆但如果你把 ω 設(shè)成負數(shù)或大于 2這個矩陣可能變得病態(tài)。解決在循環(huán)里用 try-catch 捕獲奇異矩陣錯誤或者改用左除\左除對近奇異矩陣會給出警告但至少能返回一個結(jié)果。另外檢查一下是不是把矩陣的一行數(shù)據(jù)輸錯了比如[4,-1,1;4,-8,1;-2,1,5]中間的分號是否被替換成了逗號。5.4 現(xiàn)象松弛因子 1.9 時迭代次數(shù)跳到 11828 次這是 bug 嗎不是 bug是 SOR 方法的固有現(xiàn)象。當 ω 接近上界 2 時迭代矩陣譜半徑趨近 1收斂速度趨近于零。遇到這種情況不要慌張把它寫進實驗分析里反而是加分項。解決如果你想看到完整的收斂過程可以把誤差閾值放大到10^(-3)迭代次數(shù)會顯著減少也能觀察到誤差緩慢下降的曲線。如果只是驗證代碼正確性建議只測到 1.8。5.5 現(xiàn)象用 Octave 運行結(jié)果和 MATLAB 不一樣原因Octave 的浮點數(shù)運算和 MATLAB 基本一致但某些內(nèi)置函數(shù)如inv、eig底層實現(xiàn)有細微差異在接近收斂邊界時可能表現(xiàn)出不同迭代次數(shù)。另外循環(huán)里的變量作用域和顯示精度也可能不同。解決關(guān)鍵結(jié)果不要直接用比較而是用abs(p1 - p2) 1容忍。如果迭代次數(shù)差很多檢查是不是用了不同的數(shù)據(jù)文件或沒清空工作區(qū)。6. 二次開發(fā)把實驗代碼改造成通用的迭代法求解器6.1 封裝成函數(shù)支持 method 參數(shù)和松弛因子原報告的代碼是平鋪的每換一個方法就要復制改一遍。更好的做法是封裝成一個函數(shù)function [x, iter, err_history] iterative_solve(A, b, method, omega, tol, x0) % method: jacobi, gs, sor % omega: 僅 SOR 使用 % tol: 迭代終止閾值 % x0: 初始向量默認零向量 D diag(diag(A)); L -tril(A, -1); U -triu(A, 1); n size(A, 1); x x0; err_history []; iter 0; while true switch method case jacobi x_new D \ ((L U) * x) D \ b; case gs x_new (D - L) \ (U * x) (D - L) \ b; case sor x_new (D - omega*L) \ (((1-omega)*D omega*U) * x) ... omega * ((D - omega*L) \ b); otherwise error(未知方法); end err max(abs(x_new - x)); err_history [err_history; err]; x x_new; iter iter 1; if err tol || iter 10000 break; end end end這個函數(shù)把三種方法統(tǒng)一到了一套switch結(jié)構(gòu)里最大的好處是測試不同方法時不需要反復改代碼。我把原來的 while 條件改成了while true靠迭代上限和誤差閾值雙重退出避免死循環(huán)。6.2 增加誤差記錄與收斂曲線在函數(shù)里我把每次迭代的誤差存進err_history這樣調(diào)用后可以直接畫收斂曲線[x, iter, his] iterative_solve(A, b, sor, 1.1, 1e-5, [0;0;0]); semilogy(1:length(his), his, -o); xlabel(迭代次數(shù)); ylabel(誤差對數(shù)坐標);用semilogy畫誤差下降曲線能直觀看到 SOR 在不同 ω 下收斂速度的差異。你會發(fā)現(xiàn) ω1.9 的曲線下降得非常平緩這就是 11828 次的來源。6.3 驗證方法用嚴格對角占優(yōu)矩陣和譜半徑判斷改造完之后必須驗證函數(shù)沒有跑偏。我一般會做兩步第一步構(gòu)造一個嚴格對角占優(yōu)矩陣比如三對角矩陣用直接法A\b的結(jié)果和迭代解對比第二步計算迭代矩陣的譜半徑確認小于 1。另外再跑一遍原始實驗數(shù)據(jù)確保 Jacobi 返回 14 次、Gauss-Seidel 返回 9 次。驗證通過后再去調(diào) ω就心中有數(shù)了。從那以后我每次做迭代法實驗都會先寫一個iterative_solve這樣的通用函數(shù)再批量跑參數(shù)而不是復制粘貼三次代碼。這個習慣幫我少走了很多彎路希望你也能直接用上。希望幫到你。本文還有配套的精品資源點擊獲取