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