簡介這是東北大學(xué)數(shù)值分析課程的一份完整實驗報告面向數(shù)學(xué)系及軟件等專業(yè)學(xué)生可用于學(xué)習(xí)線性方程組迭代法、撰寫實驗報告及參考MATLAB代碼。資源包內(nèi)僅含1個doc文檔大小約186KB內(nèi)容涵蓋Jacobi迭代法、Gauss-Seidel迭代法與SOR松弛迭代法的原理、迭代公式、數(shù)據(jù)記錄、結(jié)果分析及可直接運行的MATLAB實現(xiàn)。報告以典型三階方程組為例在相同精度下對比三種方法的迭代次數(shù)與誤差并針對SOR法給出ω從1.10到1.90的完整實驗數(shù)據(jù)直觀反映松弛因子對收斂速度的影響同時記錄了不同初始向量下的迭代表現(xiàn)可幫助理解初值對迭代法收斂過程的作用。內(nèi)容遵循東北大學(xué)理學(xué)院實驗報告規(guī)范包含實驗?zāi)康摹⒃?、步驟、數(shù)據(jù)處理、結(jié)果分析等要素并預(yù)留指導(dǎo)教師簽字欄可直接作為課程實驗報告模板參考。目前已有2614人學(xué)習(xí)下載適合正在準(zhǔn)備數(shù)值分析實驗、復(fù)習(xí)迭代算法或需要MATLAB實現(xiàn)參考的同學(xué)。1. 一個ω從1.80調(diào)到1.90迭代次數(shù)為什么從153暴漲到11828SOR迭代法里松弛因子ω從1.80調(diào)到1.90表面看只是把“步子邁大”了一點點可同一份代碼、同一個10^(-5)的誤差閾值下迭代次數(shù)從153次暴漲到11828次。反過來ω1.10時這個數(shù)字只有11次。也就是說一個看起來無足輕重的參數(shù)能把同一個線性方程組的求解時間從“瞬間完成”拖到“肉眼可見卡住”。這正是數(shù)值分析里迭代法最容易被當(dāng)成玄學(xué)的地方Jacobi、Gauss-Seidel、SOR三種迭代法拆開看不過是一行矩陣遞推合起來卻牽涉矩陣分裂、譜半徑、收斂性三大底層概念。本文要拆的這份某高校數(shù)值分析課程的實驗報告恰好用一套完整的MATLAB代碼和三組實驗數(shù)據(jù)把這三件事串了起來值得逐行讀懂。2. 一個矩陣分裂三個遞推式Jacobi、Gauss-Seidel、SOR的數(shù)學(xué)同構(gòu)2.1 把A拆成D、L、U負(fù)號約定不是玄學(xué)實驗報告的代碼第一屏就是在拆矩陣。A是那個3階方陣b是右端項然后一行一行取出D、L、UA [4,-1,1; 4,-8,1; -2,1,5]; b [7; -21; 15]; D diag(diag(A)); % 對角矩陣 L -tril(A,-1); % 嚴(yán)格下三角取負(fù)號 U -triu(A,1); % 嚴(yán)格上三角取負(fù)號這段代碼有三個細(xì)節(jié)值得停下來看。第一tril(A,-1)和triu(A,1)都跳過了對角線取出的是“嚴(yán)格”下三角和嚴(yán)格上三角第二取出來之后立刻乘了負(fù)號第三D、L、U 合起來滿足 A D - L - U而不是直覺上的 A D L U。為什么要用減號因為后面三種迭代法都可以寫成“矩陣分裂”的形式把 A 拆成 M - N其中 M 要求可逆并且容易求逆。Jacobi 取 M D、N -(LU)Gauss-Seidel 取 M D - L、N USOR 取 M (1/ω)(D - ωL)、N (1/ω)((1-ω)D ωU)。一旦按這個約定拆三種迭代公式的右端項都干凈得像教材原樣迭代矩陣的譜半徑也直接落在我們關(guān)心的位置。如果反過來用 A D L U 的拆法代碼里會到處冒出負(fù)號稍不留神就把迭代矩陣寫錯。常見的實現(xiàn)思路就是像實驗報告這樣tril/triu取完之后統(tǒng)一乘 -1讓公式和教材保持完全對齊。2.2 Jacobi三個分量同一輪齊步走Jacobi 的推導(dǎo)很直接由 A D - L - U 得 DX b - (LU)X于是 X D^{-1}(LU)X D^{-1}b。寫成迭代式就是X_{k1} D^{-1}(LU)X_k D^{-1}b實驗報告里的代碼完全對應(yīng)這個式子D diag(diag(A)); % 對角矩陣 L -tril(A,-1); % 嚴(yán)格下三角 U -triu(A,1); % 嚴(yán)格上三角 X_0 [0;0;0]; % 初始向量 X_1 inv(D)*(LU)*X_0 inv(D)*b; % 第一步迭代 dis abs(X_1 - X_0); % 相鄰兩步解的絕對差 p 1; % 已迭代1次 for i 1:3 while dis(i) 1e-5 X_1 inv(D)*(LU)*X_0 inv(D)*b; dis abs(X_1 - X_0); X_0 X_1; p p 1; end end注意D^{-1}(LU)和D^{-1}b其實跟迭代步無關(guān)理想實現(xiàn)應(yīng)該在循環(huán)外只算一次但對3階矩陣這種開銷無所謂。真正重要的行為是每一步迭代X 的三個新分量都完全由 X_0 的舊值計算出來誰也不占誰的便宜。這就是“同時更新”天然可并行。壞處也在這里第一個分量 x1 的新值已經(jīng)算出來了第二個方程里卻還拿著舊值在算等于每輪都浪費了一次“信息復(fù)用”的機會。2.3 Gauss-Seidel新值剛算出來就被下一個方程拿去用Gauss-Seidel 的分量形式是算第 i 個新分量時小于 i 的角標(biāo)全部用本輪新值大于 i 的角標(biāo)才用舊值。寫成方程組再整理就是X_{k1} (D-L)^{-1}U X_k (D-L)^{-1}b實驗報告里 GS 的循環(huán)體跟 Jacobi 幾乎長得一模一樣只是把inv(D)換成了inv(D-L)。這里有個容易誤讀的點矩陣形式的代碼依然是“X_0 全部舊值算出 X_1”看起來不像就地更新。但實際上(D-L)^{-1}右乘一個向量等價于一次自下而上的三角回代新算出的分量在同一個表達式的后續(xù)方程里立刻生效。換句話說矩陣形式把分量式里“邊算邊用”的過程壓縮進了一次回代里。為什么 GS 通常更快因為每一輪迭代里已經(jīng)更新的信息被后面的方程及時消化掉了不像 Jacobi 要等到下一輪才用上。代價是這種更新天然串行x2 必須等 x1x3 必須等 x2沒法直接并行。對3階小矩陣這無所謂但對大規(guī)模并行求解這個差異會直接影響算法選型。2.4 SOR給Gauss-Seidel裝上一個可調(diào)“油門”SOR 的矩陣形式是把 GS 的更新量再乘一個 ω然后加回舊值X_{k1} (D-ωL)^{-1}((1-ω)D ωU)X_k ω(D-ωL)^{-1}b實驗代碼里直接用 k 表示 ωk 1.75; % 松弛因子 ω D diag(diag(A)); % 對角矩陣 L -tril(A,-1); % 嚴(yán)格下三角 U -triu(A,1); % 嚴(yán)格上三角 X_0 [0;0;0]; % 初始向量 X_1 inv(D-k*L)*((1-k)*Dk*U)*X_0 k*inv(D-k*L)*b;當(dāng) ω1 時(1-ω)D ωU退化成 U公式就完整退回 Gauss-Seidelω1 叫欠松弛ω1 叫超松弛。所謂“松弛”字面意思是把這一步的更新量朝目標(biāo)方向再推一把或收一把。推過頭了不會立刻發(fā)散但會在目標(biāo)附近來回震蕩白白多跑很多輪——實驗里 ω1.90 那 11828 次就是這么來的。三種方法寫全了就是一張對照表方法迭代矩陣每步主要成本JacobiD^{-1}(LU)一次矩陣乘向量Gauss-Seidel(D-L)^{-1}U一次下三角回代SOR(D-ωL)^{-1}((1-ω)DωU)一次下三角回代三種迭代矩陣同構(gòu)于“A M - N 分裂”這一個框架區(qū)別只在 M 拿走多少、N 剩下多少。這也是為什么同一份實驗報告能一口氣對比三種方法它們不是三個獨立算法而是一個分裂框架下的三個參數(shù)化實例。3. 收斂速度的底層機制譜半徑、漸近收斂因子與那11828次3.1 收斂的充要條件與“嚴(yán)格對角占優(yōu)”這個安全墊把三種迭代都寫成 X_{k1} B X_k f 之后誤差服從 e_{k1} B e_k于是 e_k B^k e_0。要讓誤差對任意初值都趨于0核心條件只有一個迭代矩陣 B 的譜半徑 ρ(B) 1。譜半徑是 B 所有特征值模長的最大值它直接決定誤差每步被壓縮多少。譜半徑不好算所以實驗課都會給一個更省事的充分條件嚴(yán)格對角占優(yōu)。也就是每一行對角元絕對值的和必須大于該行所有非對角元絕對值之和。對報告里這個矩陣逐行檢查|4| |-1| |1| 2|-8| |4| |1| 5|5| |-2| |1| 3。三行全滿足所以 Jacobi 和 Gauss-Seidel 的收斂性在理論上就有保證。這也解釋了實驗結(jié)論里“三種迭代法均收斂”不是運氣而是矩陣本身選得好一個嚴(yán)格對角占優(yōu)矩陣天然就是演示迭代法收斂性最安全的教學(xué)道具。SOR 的收斂性比前兩者更挑矩陣。嚴(yán)格對角占優(yōu)并不自動保證任意 ω∈(0,2) 的 SOR 都收斂本實驗里 ω 從 1.10 一路測到 1.90 確實都收斂但不能把這個結(jié)論推廣到所有矩陣。能做的最穩(wěn)妥的驗證就是像實驗報告這樣逐點實測。順便驗證一下解的精確值。報告表格里寫的“精確解243”其實是 [2;4;3] 的連寫。帶回原方程4×2 - 1×4 1×3 74×2 - 8×4 1×3 -21-2×2 1×4 5×3 15。右端項完全吻合說明這個例子選的是有唯一解且解很規(guī)則的3階系統(tǒng)便于學(xué)生在迭代次數(shù)上做對比。3.2 為什么Gauss-Seidel比Jacobi快信息復(fù)用同一初值 [0;0;0]、同一閾值 1e-5 下Jacobi 用了 14 次Gauss-Seidel 用了 9 次。差的這 5 次就是“信息復(fù)用”的收益。用分量視角看Jacobi 每一輪的三次更新都讀舊值GS 每算出一個新 x1立刻被第二個方程拿去用新 x1、x2 又立刻被第三個方程拿去用。相當(dāng)于每一輪迭代里GS 比 Jacobi 多消化了半輪的新信息所以同樣走一步GS 的“有效信息量”更高。需要注意這并不絕對教材上有 Jacobi 收斂而 Gauss-Seidel 不收斂的反例但在嚴(yán)格對角占優(yōu)的矩陣上GS 的優(yōu)勢是典型表現(xiàn)。3.3 ω的譜半徑軌跡誤差按ρ^k衰減尾部必然爆炸實驗報告里 [0;0;0] 初值下有一組完整的 ω 掃描數(shù)據(jù)ω迭代次數(shù)1.10111.20141.30161.40191.50211.60321.70551.801531.9011828這組數(shù)據(jù)有兩個特征值得琢磨前半個區(qū)間幾乎是平緩的線性增長到 1.80 之后開始失控1.90 直接給出一根刺穿圖表的天線。原因是 SOR 迭代矩陣的譜半徑隨 ω 變化不是線性的。已知 Jacobi 迭代矩陣的譜半徑可以算一下對這個 Aρ(B_J) ≈ 0.335。代入最優(yōu)松弛因子公式 ω_opt 2 / (1 √(1 - ρ2))得到 ω_opt ≈ 1.03。實驗數(shù)據(jù)里 1.10 最快11次和理論最優(yōu)位置完全對得上一旦越過最優(yōu)點譜半徑開始回升ω→2 時譜半徑趨近于1。誤差衰減是 e_k ≈ ρ^k e_0想達到 1e-5需要的步數(shù)粗略估計是 k ≈ ln(ε) / ln(ρ)。ρ0.335 時這個值只有約 10 步出頭對應(yīng) Jacobi 的 14 次非常合理而當(dāng) ρ 逼近 0.99 甚至 0.999 時k 會從幾百跳到上萬正是 153 次到 11828 次之間發(fā)生的躍變。這也是為什么迭代法實驗里“跑不完”往往不是死循環(huán)而是譜半徑已經(jīng)貼到了1。4. 避坑改這份MATLAB代碼之前先看這四個細(xì)節(jié)4.1 停機判斷用的是相鄰迭代差不是殘差現(xiàn)象代碼里dis abs(X_1 - X_0)閾值 1e-5看起來“誤差”達標(biāo)了但如果你另外算一下殘差可能還沒到 1e-5。原因dis衡量的是兩次迭代之間解的變化量不是方程本身的殘差 ‖b - AX_k‖。收斂時相鄰迭代差確實趨近于0但反過來不一定成立。矩陣條件數(shù)大時可能出現(xiàn)“解已經(jīng)幾乎不動但殘差還明顯偏大”的情況這時候按相鄰差停機就會提前收工。解決做實驗時至少加一行驗證norm(b - A*X_1, inf)讓殘差也低于閾值再停止循環(huán)。更嚴(yán)謹(jǐn)?shù)膶懛ㄊ前淹C條件直接改成殘差范數(shù)而不是相鄰解差。4.2 for套while的怪異結(jié)構(gòu)等價于while all不是提前停機現(xiàn)象初次讀這段代碼很容易以為某個分量先達標(biāo)循環(huán)就停了甚至擔(dān)心它會漏算其他分量。原因仔細(xì)看for i 1:3里套的是三個串行的while dis(i) 1e-5。第一個 while 退出后第二個、第三個 while 還會繼續(xù)執(zhí)行而每一個 while 內(nèi)部一旦進入都對全部三個分量做一輪完整迭代。所以整體行為等價于“只要還有分量沒達標(biāo)就繼續(xù)迭代”也就是while any(dis 1e-5)的逐段展開。解決想要一眼可讀的版本直接寫while max(dis) 1e-5效果完全一樣。實驗報告里這個寫法在3階問題上沒問題但它的可讀性確實容易被誤解。4.3 小矩陣?yán)锔矣胕nv規(guī)模一上來就翻車現(xiàn)象3階矩陣用inv(D)、inv(D-L)看起來很直觀算得也快把同樣的寫法搬到一個幾百階的稀疏矩陣上速度瞬間崩掉還可能算出明顯不對的解。原因inv是顯式求逆。稀疏矩陣的逆通常不再稀疏顯式求逆會破壞稀疏性存儲和乘法成本都失控此外求逆本身的數(shù)值誤差會被放大直接乘到右端項上。解決循環(huán)外先把迭代矩陣算好一次循環(huán)體內(nèi)只做矩陣乘向量需要解(D-L)\v時用反斜杠運算符(D-L)\v做回代而不是inv(D-L)*v。這在小矩陣上結(jié)果一致在大矩陣上是能不能跑完的區(qū)別。4.4 ω1.90已經(jīng)踩在懸崖邊上再往2走就是死循環(huán)現(xiàn)象1.80 是 153 次1.90 是 11828 次。如果試 ω1.99這個實驗基本跑不完真取 ω2代碼會永遠(yuǎn)不退出。原因SOR 的迭代矩陣在 ω→2 時譜半徑趨近于1誤差每步幾乎不衰減1e-5 可能要百萬步才能摸到ω2 時譜半徑達到1誤差不再被壓縮循環(huán)條件永遠(yuǎn)滿足。解決選 ω 之前先估算最優(yōu)點比如用 Jacobi 譜半徑代入公式 ω_opt 2/(1√(1-ρ2))把實驗范圍縮在最優(yōu)值附近?？拷?的區(qū)域不是“更快的捷徑”而是“收斂性正在消失的懸崖”。5. 數(shù)據(jù)背后的設(shè)計三組初值和一串ω到底在驗證什么5.1 三組初值的迭代次數(shù)幾乎一樣線性收斂的初值不敏感性實驗報告特意測了三組初始向量[0;0;0]、[1;1;1]、[-1;-1;-1]但把三組結(jié)果放在一起看會發(fā)現(xiàn)迭代次數(shù)幾乎沒怎么動初始向量JacobiGauss-SeidelSOR(ω1.25)SOR(ω1.50)SOR(ω1.75)[1;1;1]139142079[-1;-1;-1]149152187初值翻來覆去GS 始終是 9 次Jacobi 在 13 和 14 之間浮動SOR 的差別也在個位數(shù)。這不是巧合而是線性迭代的固有性質(zhì)誤差 e_k B^k e_0初值只影響前面的常數(shù)系數(shù)真正決定每步壓縮多少的是譜半徑 ρ。初值差幾步本質(zhì)是前幾輪各自的路徑差異迭代次數(shù)的大頭完全由 ρ 決定。實驗故意設(shè)計三組初值恰恰是想讓觀察者自己看出“初值對迭代次數(shù)影響很小”從而把注意力從“調(diào)初值”移到“調(diào)迭代格式和松弛因子”上來。理解了這一點以后遇到“換個初值迭代次數(shù)暴漲”的情況第一反應(yīng)就不該是調(diào)初值而是檢查迭代矩陣的譜半徑是不是已經(jīng)逼近1。5.2 ω序列從1.10排到1.90用一次實驗“畫”出譜半徑曲線1.10 到 1.90 每隔 0.1 取一個點九個點連起來就是 SOR 譜半徑右半支的采樣。前段 11、14、16、19、21 幾乎是線性爬坡說明譜半徑還在溫和回升到 1.70 的 55、1.80 的 153斜率突然變陡1.90 的 11828 則直接掉進“譜半徑無限逼近1”的區(qū)域。這種排點方式本身就是實驗設(shè)計的核心先放一個接近理論最優(yōu)的 1.10再逐步向危險區(qū) 2 靠攏用迭代次數(shù)當(dāng)“探針”去反推譜半徑的形狀。所以數(shù)據(jù)表不是記錄一堆無意義的數(shù)字而是在做一次繞開特征值計算的數(shù)值實驗。5.3 這份文檔資料的真正用法可復(fù)跑、可改造的迭代法基準(zhǔn)代碼這份材料既是規(guī)范的實驗報告模板也是一份能直接跑的 MATLAB 迭代法基準(zhǔn)代碼。拿到手之后最值得做的不是照著抄一遍而是改三個地方觀察行為變化把閾值從 1e-5 改成 1e-8三種方法的迭代次數(shù)會同步上升但 J、GS、SOR 之間的相對差距基本不變可以驗證“閾值只改變截止位置不改變收斂速度的排序”。把矩陣 A 的某個對角元改小破壞嚴(yán)格對角占優(yōu)條件比如把 A(1,1) 從 4 改成 1立刻能看到 Jacobi 和 GS 不再保證收斂實驗會翻車或死循環(huán)。固定 ω1.90 完整跑一次用 11828 次循環(huán)直觀感受“譜半徑接近1”意味著什么。這三組改造跑完收斂性理論就不再是書上的定理而是自己手里的實驗現(xiàn)象。6. 設(shè)計取舍迭代法不是萬能藥譜半徑才是總開關(guān)6.1 直接法與迭代法的分界線回到最根本的問題這個3階方程組用A\b一行就能解出來為什么還要研究迭代法答案是規(guī)模。直接法的復(fù)雜度是 O(n^3)幾百階就明顯吃力迭代法每步只是矩陣乘向量或三角回代配合稀疏存儲單步成本是 O(nnz)對大規(guī)模稀疏系統(tǒng)優(yōu)勢巨大。但迭代法的總成本不是固定值它取決于譜半徑——碰到 ρ0.999 這種病態(tài)情況一萬步疊加起來反而比直接法還貴。所以工程上的取舍很清楚稠密小矩陣直接法省心稀疏大矩陣迭代法才有意義而迭代法能不能用先看譜半徑不看迭代格式的“高級感”。6.2 我做迭代實驗的固定動作從那以后我拿到一個線性方程組第一件事永遠(yuǎn)是先算迭代矩陣的譜半徑再決定要不要用迭代法然后才輪得到選 GS 還是 SOR、ω 取多少。跑完循環(huán)也一定會加一行殘差驗證而不是看到相鄰誤差達標(biāo)就收工。有一次我圖省事沒檢查譜半徑就把 ω 設(shè)成 1.9 掛機跑實驗回來發(fā)現(xiàn)循環(huán)還在轉(zhuǎn)白白等了半小時——這種黑匣子式的試錯本質(zhì)上是把譜半徑這個總開關(guān)誤當(dāng)成“調(diào)參游戲”。希望幫到你。本文還有配套的精品資源點擊獲取