別及jordan數(shù)值陷阱)
簡介這份資源面向學習MATLAB數(shù)值計算與線性代數(shù)的學生和工程人員聚焦矩陣處理中的特征值、特征向量與對角化等核心問題適合作為課堂實驗或自學練習的配套材料。壓縮包內(nèi)共1個doc文檔約98KB內(nèi)容以實驗講義形式組織涵蓋特征多項式求解、eig函數(shù)求特征值與特征向量、實對稱矩陣正交對角化、balance函數(shù)提高計算精度以及jordan函數(shù)求若當標準形等知識點并配有可直接運行的MATLAB代碼示例與運算結果對照。文檔還整理了trace、sum、prod、cdf2rdf等常用矩陣函數(shù)的用法說明便于讀者在練習中隨時查閱。目前已有727人學習下載可作為矩陣分析入門與實驗報告的參考幫助讀者理解特征值分解的數(shù)學原理并掌握相應MATLAB命令的調(diào)用方式。1. 矩陣處理到底在算什么從“找方向”說起很多人第一次聽到“特征值與特征向量”腦子里冒出的是一堆符號和公式完全不知道它在生活里對應什么。其實可以這樣理解你把一個矩陣想象成一臺“變形機器”任何向量丟進去都會被拉伸、旋轉、翻轉。但總有那么幾個特殊方向丟進去之后只被拉長或壓短方向紋絲不動——這幾個方向就是特征向量拉伸的倍數(shù)就是特征值。Matlab 矩陣處理里最核心的一塊就是幫你把這幾條“不動方向”和“拉伸倍數(shù)”找出來。搞懂這件事的價值在于矩陣對角化、若當標準形、實對稱矩陣的正交分解全都建立在這個概念之上。你不需要手算行列式Matlab 的eig、jordan、balance這些函數(shù)已經(jīng)把底層運算封裝好了。但前提是你得知道每個函數(shù)在做什么、結果怎么讀、什么時候會翻車。這篇就把這些概念用大白話拆開讓你看完能分清eig和eigs的區(qū)別知道jordan為什么有時候“不靈”也明白balance到底在平衡什么。2. 特征值與特征向量Matlab 里到底怎么算、結果怎么讀2.1 從特征多項式到poly和roots在 Matlab 里求特征值最“笨”但也最直觀的路子是先算特征多項式。特征多項式的定義是 det(λI ? A)其中 I 是單位陣λ 是符號變量。Matlab 提供了poly(A)直接返回多項式系數(shù)的向量形式。比如一個 2×2 矩陣poly會給你三個系數(shù)分別對應 λ2、λ1 和常數(shù)項。A [1 0; 2 3]; p poly(A); % 返回特征多項式系數(shù)向量 f poly2str(p, x); % 轉成可讀字符串形式poly返回的是數(shù)值向量poly2str只是把它排版成人類能讀的表達式。如果你想知道特征值具體是多少還得對多項式求根用roots函數(shù)。這條路適合教學演示因為你能親眼看到特征多項式長什么樣。但實際工程里沒人這么干——高階矩陣的多項式求根數(shù)值穩(wěn)定性很差Matlab 內(nèi)部也不會用這條路去算eig。另一種更“數(shù)學”的寫法是用符號變量直接算行列式syms x; E eye(2); f det(x*E - A); % 符號形式的特征多項式這樣得到的是符號表達式適合推導和驗證但不適合大規(guī)模數(shù)值計算。符號運算慢而且當矩陣階數(shù)上去之后表達式會膨脹到?jīng)]法看。2.2eig的三種返回形式與eigs的區(qū)別真正干活用的是eig。它有三種常見調(diào)用方式返回的東西不一樣用途也不一樣。第一種d eig(A)只返回一個列向量里面裝著所有特征值。你只關心特征值大小、不關心特征向量時用這個。第二種[V, D] eig(A)返回兩個矩陣。D 是對角矩陣主對角線上的元素就是特征值V 的每一列是對應的特征向量。滿足的關系是 A·V V·D等價于 D V?1AV。這就是矩陣對角化的核心等式。第三種[V, D] eigs(A)注意這里函數(shù)名多了個 s。eigs不是返回全部特征值而是按大小次序返回一部分通常用于大型稀疏矩陣。它底層用的是迭代法只算你指定的前 k 個。如果你對一個 10000×10000 的矩陣用eig它會試圖算全部特征值內(nèi)存和時間都可能爆掉用eigs指定算前 6 個就輕松很多。A [0 1 0 0; 1 0 0 0; 0 0 0 1; 0 0 1 0]; d eig(A); % 全部特征值列向量 [V, D] eig(A); % 特征向量矩陣 V 和特征值對角陣 D運行后你會發(fā)現(xiàn)特征值是 1、-1、1、-1有重根。D 的主對角線上按順序排列V 的列向量對應每個特征值。這里有個細節(jié)當特征值有重根時eig返回的特征向量不一定是正交的甚至可能線性相關。如果你需要正交基得用實對稱矩陣那條路后面會講。2.3 結果怎么驗證V^-1*A*V和符號化輸出算出[V, D]之后怎么確認沒算錯最直接的辦法是驗算 V?1AV 是否等于 D。在 Matlab 里可以寫成V_inv inv(V); check V_inv * A * V; % 應該接近對角陣 D如果 A 可對角化check的非對角元素應該接近 0。但數(shù)值計算有浮點誤差你可能會看到 10?1? 量級的小數(shù)那是正常的。有時候你希望看到精確的分數(shù)或根號形式而不是小數(shù)。這時可以把 V 轉成符號類型V_sym sym(V); V_sym^-1 * A * V_sym這樣輸出就是分數(shù)和根號方便對照手算結果。但注意符號運算只適合小矩陣驗證大矩陣會非常慢。提示eig返回的特征向量默認是歸一化的但歸一化方式不唯一。不同版本的 Matlab 可能在符號上差一個負號這不影響正確性。3. 對角化與若當標準形能對角化最好不能也別慌3.1 可對角化的條件與實對稱矩陣的特殊待遇一個 n 階方陣能不能對角化取決于它有沒有 n 個線性無關的特征向量。如果特征值都是單根那一定可對角化如果有重根就要看重根對應的特征向量個數(shù)夠不夠。不夠的話就只能退而求其次求若當標準形。實對稱矩陣是個特例它一定可以對角化而且特征向量可以選成兩兩正交的。Matlab 的eig對實對稱矩陣會自動返回正交矩陣 P滿足 P?1AP D而且 P?1 就等于 P?。這意味著你不需要求逆轉置一下就行數(shù)值上更穩(wěn)定。A [2 2 -2; 2 5 -4; -2 -4 5]; % 實對稱矩陣 [P, D] eig(A); % P 是正交陣D 是對角陣這個性質在二次型化簡里特別有用。二次型矩陣是實對稱的對角化之后就得到了標準形P 就是變換矩陣。你不需要手動配方eig一步到位。3.2jordan函數(shù)與若當塊的上三角定義當矩陣不能對角化時jordan(A)返回它的若當標準形 J。若當標準形長什么樣對角線上是特征值對角線上方可能有 1形成一個個“塊”。Matlab 里若當塊是按上三角形定義的也就是說 1 出現(xiàn)在對角線上方而不是下方。A [-1 -2 6; -1 0 3; -1 -1 4]; J jordan(A);運行后 J 是一個上三角矩陣對角線上是特征值某些位置有 1。這些 1 的位置和個數(shù)反映了矩陣的“缺陷程度”——也就是特征向量不夠用的程度。但jordan有個血淚經(jīng)驗它對數(shù)值誤差極其敏感。稍微擾動一下矩陣元素若當塊結構就可能完全改變。所以實際工程里除非你確定矩陣是精確的符號矩陣或者小整數(shù)矩陣否則不要輕易依賴jordan的數(shù)值結果。它更適合教學和理論驗證不適合處理實驗測得的浮點數(shù)據(jù)。3.3 復對角矩陣轉實對角cdf2rdf在做什么實矩陣的特征值可能是復數(shù)而且成共軛對出現(xiàn)。eig返回的 D 可能是復對角陣對角線上有 abi 和 a?bi。如果你希望保持實數(shù)運算可以用cdf2rdf把復對角陣轉成實塊對角陣。A [0 2 1; -2 0 3; -1 -3 0]; [v, d] eig(A); % d 是復對角陣 [V, D] cdf2rdf(v, d); % D 是實塊對角陣轉換之后D 的對角線上不再是復數(shù)而是 2×2 的實數(shù)塊。每個 2×2 塊對應一對共軛特征值。這樣做的好處是后續(xù)計算全在實數(shù)域進行避免復數(shù)運算帶來的額外開銷和精度問題。注意cdf2rdf只對實矩陣有意義。如果原矩陣本身就是復矩陣這個轉換不適用。4. 常見誤區(qū)與易混概念別把這幾組函數(shù)搞混了4.1eig和eigs一個算全部一個算部分最容易混的就是這兩個。名字只差一個字母行為差很遠。eig算全部特征值適合中小型稠密矩陣eigs算部分特征值適合大型稀疏矩陣。如果你對一個 5000×5000 的稀疏矩陣用eigMatlab 會試圖算 5000 個特征值內(nèi)存直接吃滿。用eigs(A, 6)只算前 6 個幾秒鐘就出結果。區(qū)分方法很簡單看矩陣大小和稀疏性。小矩陣、需要全部特征值用eig大矩陣、只關心最大或最小的幾個用eigs。4.2poly和polyfit名字像用途完全不同poly是求特征多項式系數(shù)的輸入是方陣輸出是系數(shù)向量。polyfit是多項式擬合輸入是數(shù)據(jù)點輸出是擬合系數(shù)。兩者都和“多項式”有關但一個服務于矩陣特征值一個服務于數(shù)據(jù)擬合。初學者看到poly容易聯(lián)想到擬合這是典型的望文生義。區(qū)分方法看輸入。輸入是方陣用poly輸入是 x 和 y 數(shù)據(jù)向量用polyfit。4.3jordan的數(shù)值陷阱理論成立不等于數(shù)值可靠理論上任何方陣都有若當標準形但數(shù)值上jordan可能給你一個完全錯誤的結果。原因在于若當塊結構對擾動極其敏感。一個特征值接近重根的矩陣在浮點誤差下可能被判定為可對角化也可能被判定為有若當塊結果不穩(wěn)定。區(qū)分方法如果矩陣元素是精確整數(shù)或符號jordan可信如果矩陣來自實驗測量或浮點計算不要依賴jordan的結構只看特征值就夠了。4.4balance不是“平衡方程”而是縮放矩陣balance(A)返回一個相似變換矩陣 T 和平衡矩陣 B使得 B T?1AT。它的作用是縮放矩陣的行和列讓各行各列的范數(shù)盡量接近從而提高特征值計算的精度。它不是用來解平衡方程的也不是用來做物理平衡的。區(qū)分方法看到balance出現(xiàn)在特征值計算上下文里就知道它在做矩陣縮放目的是數(shù)值穩(wěn)定。4.5trace和sum跡是特征值之和但sum不一定是trace(A)返回矩陣的跡等于主對角線元素之和也等于所有特征值之和。sum(A, dim)是元素求和dim1 按列求和dim2 按行求和。sum(sum(A,1),2)才是所有元素之和。初學者容易把trace和sum混為一談尤其是在驗證特征值之和等于跡的時候。區(qū)分方法trace只針對方陣結果是標量sum對任意矩陣都行結果可能是向量。5. 什么時候用、什么時候別用一個具體例子和邊界判斷假設你手頭有一個 3×3 的實對稱矩陣來自某個物理系統(tǒng)的剛度矩陣你想知道它的特征值和特征向量用來判斷系統(tǒng)的振動模態(tài)。這時候eig是最合適的選擇因為實對稱矩陣保證正交對角化數(shù)值穩(wěn)定結果可靠。A [2 2 -2; 2 5 -4; -2 -4 5]; [P, D] eig(A); % P 的列是振動模態(tài)D 的對角線是頻率平方但如果你面對的是一個 20000×20000 的稀疏矩陣來自有限元分析你只關心前 10 個最小特征值。這時候絕對不能用eig必須用eigs(A, 10, sm)其中sm表示 smallest magnitude。否則內(nèi)存和時間都扛不住。再比如你有一個矩陣特征值有重根你想知道它能不能對角化。這時候不要直接上jordan因為數(shù)值結果可能騙你。更穩(wěn)妥的做法是檢查特征向量的秩用[V, D] eig(A)然后看rank(V)是否等于矩陣階數(shù)。如果等于說明有足夠多的線性無關特征向量可對角化如果小于說明不能對角化若當塊存在。還有一個場景你算完特征值發(fā)現(xiàn)結果里有很小的虛部比如 1.0000 0.0000i。這通常是浮點誤差不是真的復數(shù)特征值??梢杂胷eal函數(shù)取實部或者用cdf2rdf轉成實塊對角陣。提示判斷特征值是否為實數(shù)不要看虛部是否嚴格為零而是看虛部是否小于某個容差比如abs(imag(d)) 1e-10。從那以后我每次用eig之前都會先問自己三個問題矩陣是稠密還是稀疏需要全部特征值還是前幾個矩陣有沒有特殊結構對稱、三角、稀疏這三個問題決定了用eig還是eigs用不用balance以及結果怎么驗證。希望幫到你。本文還有配套的精品資源點擊獲取