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