簡(jiǎn)介這份資源面向?qū)W習(xí)MATLAB數(shù)值計(jì)算與線性代數(shù)的學(xué)生和工程人員聚焦矩陣處理中的特征值、特征向量與對(duì)角化等核心問(wèn)題適合作為課堂實(shí)驗(yàn)或自學(xué)練習(xí)的配套材料。壓縮包內(nèi)共1個(gè)doc文檔約98KB內(nèi)容以實(shí)驗(yàn)講義形式組織涵蓋特征多項(xiàng)式求解、eig函數(shù)求特征值與特征向量、實(shí)對(duì)稱矩陣的正交對(duì)角化、balance函數(shù)提高計(jì)算精度以及jordan函數(shù)求若當(dāng)標(biāo)準(zhǔn)形等知識(shí)點(diǎn)并配有可直接運(yùn)行的MATLAB代碼示例與運(yùn)算結(jié)果。文檔還整理了trace、sum、prod、cdf2rdf等常用矩陣函數(shù)的用法說(shuō)明便于讀者對(duì)照實(shí)驗(yàn)步驟復(fù)現(xiàn)結(jié)果、理解矩陣分析的基本流程。目前已有727人學(xué)習(xí)下載適合需要系統(tǒng)梳理MATLAB矩陣運(yùn)算命令、完成相關(guān)實(shí)驗(yàn)報(bào)告或查漏補(bǔ)缺的讀者參考使用。1. 從一次翻車的對(duì)角化說(shuō)起這套 MATLAB 矩陣處理實(shí)驗(yàn)到底值不值得跑帶過(guò)幾屆做數(shù)值計(jì)算方向的學(xué)生之后我發(fā)現(xiàn)一個(gè)規(guī)律凡是自稱學(xué)過(guò)線性代數(shù)的人第一次用 MATLAB 處理特征值問(wèn)題十有八九會(huì)在對(duì)角化驗(yàn)證那一步翻車。不是概念不懂而是從紙面推導(dǎo)到工程實(shí)現(xiàn)之間隔著一層?xùn)|西——數(shù)值精度、符號(hào)運(yùn)算、若當(dāng)塊的排列方向每一個(gè)都能讓結(jié)果對(duì)不上。這份實(shí)驗(yàn)六的資料恰好就是沖著這層隔閡來(lái)的。它圍繞特征值、特征向量、若當(dāng)標(biāo)準(zhǔn)形這幾個(gè)矩陣分析里的硬骨頭把 MATLAB 里poly、eig、jordan、balance、cdf2rdf這些函數(shù)串成了一條完整的操作鏈。適合誰(shuí)正在做線性代數(shù)實(shí)驗(yàn)、需要把矩陣對(duì)角化落到代碼上的學(xué)生以及工作中偶爾要處理模態(tài)分析、系統(tǒng)穩(wěn)定性判據(jù)、二次型標(biāo)準(zhǔn)化的工程師。它不教你線性代數(shù)它教你怎么讓 MATLAB 把線性代數(shù)算對(duì)。2. 特征多項(xiàng)式與特征值poly、eig 和符號(hào)求解三條路怎么選2.1 特征多項(xiàng)式poly 的向量表示和 poly2str 的可讀化MATLAB 里求特征多項(xiàng)式最直接的是poly(A)它返回的是多項(xiàng)式系數(shù)的行向量按降冪排列。這個(gè)向量本身不好讀得配合poly2str轉(zhuǎn)成人類能看的字符串形式。% 定義方陣 A A [1 0; 2 3]; % poly 返回特征多項(xiàng)式的系數(shù)向量按降冪排列 p poly(A) % 把系數(shù)向量轉(zhuǎn)成可讀的多項(xiàng)式字符串 f poly2str(p, x)運(yùn)行后p是[1 -4 3]對(duì)應(yīng)x^2 - 4x 3。這里有個(gè)容易忽略的點(diǎn)poly(A)返回的系數(shù)向量長(zhǎng)度是 n1首項(xiàng)系數(shù)恒為 1因?yàn)樘卣鞫囗?xiàng)式是首一的。如果你拿這個(gè)向量去做roots得到的就是特征值。參數(shù)上沒(méi)什么可調(diào)的poly只接受方陣非方陣會(huì)直接報(bào)錯(cuò)。另一種路子是從定義出發(fā)用符號(hào)變量硬算行列式A [1 0; 2 3]; E eye(2); % 2 階單位陣 syms x % 聲明符號(hào)變量 f det(x*E - A) % 直接按定義算特征多項(xiàng)式這條路的好處是結(jié)果以符號(hào)形式呈現(xiàn)(x-1)*(x-3)一眼就能看出根。代價(jià)是符號(hào)運(yùn)算慢矩陣一大就卡。我一般只在矩陣階數(shù)小于 5 且需要精確表達(dá)式時(shí)才用符號(hào)法其余情況一律走poly加roots。2.2 eig 的四種調(diào)用形式與特征向量的歸一化問(wèn)題eig是這套流程的核心。它至少有四種調(diào)用形式每種返回的東西不一樣混用是新手最常見(jiàn)的翻車點(diǎn)。A [0 1 0 0; 1 0 0 0; 0 0 0 1; 0 0 1 0]; % 形式一只返回特征值列向量 d eig(A) % 形式二返回特征向量矩陣 V 和特征值對(duì)角陣 D滿足 A*V V*D [V, D] eig(A) % 形式三eigs 按大小次序返回適合大型稀疏矩陣 [V, D] eigs(A) % 形式四廣義特征值問(wèn)題 A*x lambda*B*x % [V, D] eig(A, B)形式二里有個(gè)關(guān)鍵細(xì)節(jié)D的主對(duì)角線元素是特征值但順序不保證按大小排。上面這個(gè)例子里D的對(duì)角線是-1, -1, 1, 1特征值 -1 是二重的對(duì)應(yīng)V的第 1、2 列。這里V的列向量是歸一化的模長(zhǎng)為 1但符號(hào)可能和你手算的不一樣——這是數(shù)值算法的正常行為不是錯(cuò)誤。驗(yàn)證對(duì)角化是否成立標(biāo)準(zhǔn)操作是% 用符號(hào)形式輸出 V避免浮點(diǎn)誤差干擾判斷 V_sym sym(V) % 驗(yàn)證 V^-1 * A * V 是否等于 D V_sym^-1 * A * V_sym結(jié)果應(yīng)該是一個(gè)對(duì)角陣對(duì)角線上就是特征值。如果這里算出來(lái)非對(duì)角元有微小數(shù)值別慌那是浮點(diǎn)誤差用sym轉(zhuǎn)成符號(hào)再算就干凈了。2.3 提高精度balance 什么時(shí)候該用balance這個(gè)函數(shù)很多人沒(méi)用過(guò)但它在特征值計(jì)算精度上確實(shí)能救場(chǎng)。原理是對(duì)矩陣做相似變換把行范數(shù)和列范數(shù)拉到接近的量級(jí)減少數(shù)值計(jì)算中的舍入誤差。% 返回相似變換矩陣 T 和平衡矩陣 B滿足 B T^-1 * A * T [T, B] balance(A); % 只要平衡矩陣 B balance(A);什么時(shí)候該用當(dāng)你發(fā)現(xiàn)eig算出來(lái)的特征值對(duì)擾動(dòng)特別敏感或者矩陣元素量級(jí)差異很大比如有的元素是 1e-8有的是 1e8先balance一下再求特征值精度會(huì)明顯改善。注意balance不改變特征值只改變矩陣的數(shù)值條件。我一般在對(duì)精度有要求且矩陣條件數(shù)偏大時(shí)才會(huì)加這一步日常小矩陣直接eig就夠了。3. 對(duì)角化實(shí)戰(zhàn)實(shí)對(duì)稱矩陣、復(fù)對(duì)角轉(zhuǎn)實(shí)對(duì)角與若當(dāng)標(biāo)準(zhǔn)形3.1 實(shí)對(duì)稱矩陣的對(duì)角化eig 直接給出正交陣實(shí)對(duì)稱矩陣有個(gè)漂亮性質(zhì)一定可以對(duì)角化而且特征向量可以選成正交的。MATLAB 的eig對(duì)實(shí)對(duì)稱矩陣會(huì)自動(dòng)返回正交矩陣P滿足P * A * P D。A [2 2 -2; 2 5 -4; -2 -4 5]; % 實(shí)對(duì)稱矩陣 [P, D] eig(A) % P 是正交陣D 是對(duì)角陣 % 驗(yàn)證正交性P * P 應(yīng)該是單位陣 P * P % 驗(yàn)證對(duì)角化P * A * P 應(yīng)該等于 D P * A * P運(yùn)行后D的對(duì)角線是1, 1, 10P的列向量?jī)蓛烧磺夷ｉL(zhǎng)為 1。這里P * P會(huì)得到一個(gè)接近單位陣的結(jié)果非對(duì)角元是 1e-16 量級(jí)的浮點(diǎn)噪聲。如果你在做二次型標(biāo)準(zhǔn)化P就是變換矩陣D就是標(biāo)準(zhǔn)形系數(shù)。注意實(shí)對(duì)稱矩陣不會(huì)出現(xiàn)復(fù)特征值所以不需要cdf2rdf處理。3.2 復(fù)對(duì)角矩陣轉(zhuǎn)實(shí)對(duì)角cdf2rdf 的 2x2 實(shí)數(shù)塊非對(duì)稱實(shí)矩陣的特征值可能是復(fù)數(shù)而且成共軛對(duì)出現(xiàn)。eig返回的D是對(duì)角陣但對(duì)角線上是復(fù)數(shù)。如果你需要實(shí)數(shù)形式的塊對(duì)角陣就得用cdf2rdf。A [0 2 1; -2 0 3; -1 -3 0]; [v, d] eig(A) % v 和 d 都是復(fù)數(shù)形式 [V, D] cdf2rdf(v, d) % 轉(zhuǎn)成實(shí)數(shù)塊對(duì)角形式轉(zhuǎn)換后D不再是對(duì)角陣而是在對(duì)角線上用 2x2 實(shí)數(shù)塊代替共軛復(fù)數(shù)對(duì)。上面這個(gè)例子里D會(huì)變成一個(gè) 3x3 的塊對(duì)角陣其中有一個(gè) 2x2 塊對(duì)應(yīng)特征值0 ± 3.7417i。這個(gè)函數(shù)在控制系統(tǒng)里做模態(tài)分析時(shí)特別有用因?yàn)閷?shí)矩陣形式更方便后續(xù)計(jì)算。參數(shù)上cdf2rdf要求輸入的v和d必須是eig輸出的復(fù)數(shù)形式不能先手動(dòng)取實(shí)部。3.3 若當(dāng)標(biāo)準(zhǔn)形jordan 函數(shù)與上三角塊的排列規(guī)則若當(dāng)標(biāo)準(zhǔn)形是矩陣相似理論的終極形式MATLAB 里用jordan求。A [-1 -2 6; -1 0 3; -1 -1 4]; J jordan(A) % 返回 A 的若當(dāng)標(biāo)準(zhǔn)形運(yùn)行結(jié)果是J 1 1 0 0 1 0 0 0 1這里有個(gè)必須記住的坑MATLAB 的若當(dāng)塊是按上三角形定義的也就是說(shuō)非對(duì)角元 1 出現(xiàn)在對(duì)角元的上方。這和某些教材里下三角的定義相反。如果你拿教材上的結(jié)果直接對(duì)比發(fā)現(xiàn) 1 的位置不對(duì)先確認(rèn)定義方向。另外jordan對(duì)數(shù)值矩陣的魯棒性一般矩陣接近虧損時(shí)結(jié)果可能不穩(wěn)定符號(hào)矩陣會(huì)好很多。3.4 跡與行列式的驗(yàn)證trace、sum、prod 的配合特征值有兩個(gè)基本恒等式所有特征值之和等于矩陣的跡所有特征值之積等于矩陣的行列式。用 MATLAB 驗(yàn)證這兩個(gè)性質(zhì)能快速判斷特征值算得對(duì)不對(duì)。A [1 2 2; 2 1 2; 2 2 1]; a trace(A) % 矩陣的跡 d eig(A) % 特征值 b sum(d, 1) % 特征值求和dim1 按列求和 e det(A) % 矩陣行列式 f prod(d, 1) % 特征值求積dim1 按列求積sum和prod的dim參數(shù)容易搞混dim1按列操作dim2按行操作。對(duì)于列向量dsum(d, 1)和sum(d, 2)結(jié)果一樣但習(xí)慣上求全部元素和用sum(sum(A, 1), 2)更穩(wěn)妥。驗(yàn)證通過(guò)的標(biāo)準(zhǔn)是a和b相等、e和f相等浮點(diǎn)誤差在 1e-10 量級(jí)以內(nèi)都算正常。4. 避坑與排查特征值計(jì)算里那些讓人懷疑人生的報(bào)錯(cuò)4.1 現(xiàn)象eig 返回的特征向量驗(yàn)證不通過(guò)V^-1AV 不是對(duì)角陣原因通常有兩個(gè)。一是矩陣本身不可對(duì)角化虧損矩陣特征向量不夠 n 個(gè)線性無(wú)關(guān)的這時(shí)候V是奇異的V^-1根本不存在或者數(shù)值上極不穩(wěn)定。二是浮點(diǎn)誤差累積尤其是矩陣條件數(shù)很大時(shí)V^-1 * A * V的非對(duì)角元可能跑到 1e-8 量級(jí)看起來(lái)不像對(duì)角陣。解決辦法先用rank(V)檢查特征向量矩陣是否滿秩。如果不滿秩說(shuō)明矩陣虧損得走若當(dāng)標(biāo)準(zhǔn)形路線。如果是浮點(diǎn)誤差把V轉(zhuǎn)成符號(hào)形式sym(V)再算V^-1 * A * V結(jié)果就干凈了。我一般會(huì)在驗(yàn)證前加一句format short控制顯示精度避免被一堆小數(shù)位干擾判斷。4.2 現(xiàn)象jordan 函數(shù)報(bào)錯(cuò)或返回結(jié)果和教材對(duì)不上jordan對(duì)數(shù)值矩陣的容差很敏感。如果矩陣元素是浮點(diǎn)數(shù)且接近虧損MATLAB 可能判定為非虧損返回一個(gè)對(duì)角陣而不是帶若當(dāng)塊的矩陣。另外前面提過(guò)若當(dāng)塊的方向是上三角教材可能是下三角。解決辦法把矩陣轉(zhuǎn)成符號(hào)形式再調(diào)jordanjordan(sym(A))會(huì)走精確運(yùn)算結(jié)果穩(wěn)定得多。方向問(wèn)題只能靠對(duì)照定義解決記住 MATLAB 是上三角塊。如果矩陣規(guī)模大且必須用數(shù)值方法考慮先用balance改善條件數(shù)再算。4.3 現(xiàn)象cdf2rdf 報(bào)錯(cuò)說(shuō)輸入維度不匹配cdf2rdf要求輸入的v和d必須來(lái)自eig的復(fù)數(shù)輸出且d必須是對(duì)角陣形式。如果你手動(dòng)構(gòu)造了一個(gè)實(shí)數(shù)矩陣傳進(jìn)去或者d不是對(duì)角陣就會(huì)報(bào)錯(cuò)。解決辦法確保調(diào)用鏈?zhǔn)荹v, d] eig(A)之后直接[V, D] cdf2rdf(v, d)中間不要對(duì)v或d做任何修改。如果A是實(shí)對(duì)稱矩陣eig返回的本來(lái)就是實(shí)數(shù)不需要cdf2rdf強(qiáng)行調(diào)用反而會(huì)出錯(cuò)。4.4 現(xiàn)象符號(hào)運(yùn)算 solve 求特征值結(jié)果里有 root 表達(dá)式或復(fù)數(shù)用solve(f)求符號(hào)多項(xiàng)式的根時(shí)如果多項(xiàng)式次數(shù)大于 4MATLAB 可能返回root形式的符號(hào)解或者帶復(fù)數(shù)。這不是錯(cuò)誤是代數(shù)基本定理的必然結(jié)果——五次以上沒(méi)有通用根式解。解決辦法對(duì)高次多項(xiàng)式直接用roots求數(shù)值解或者用vpa把符號(hào)解轉(zhuǎn)成高精度數(shù)值。如果只關(guān)心實(shí)特征值用real取實(shí)部后再判斷虛部是否可忽略。我一般對(duì)超過(guò) 4 階的矩陣直接走eig不碰符號(hào)求解省得跟root表達(dá)式較勁。4.5 現(xiàn)象balance 之后特征值變了balance做的是相似變換理論上不改變特征值。如果你發(fā)現(xiàn)變了大概率是變換后的矩陣條件數(shù)改善數(shù)值計(jì)算更準(zhǔn)了所以結(jié)果和沒(méi)平衡之前有微小差異。這不是 bug是精度提升的表現(xiàn)。解決辦法把balance前后的特征值都算出來(lái)對(duì)比差異應(yīng)該在合理范圍內(nèi)。如果差異巨大檢查T矩陣是否可逆balance返回的T在極端情況下可能接近奇異。我一般只在確認(rèn)矩陣條件數(shù)很差時(shí)才用balance日常計(jì)算不加這一步。5. 進(jìn)階技巧用符號(hào)工具箱做精確對(duì)角化與二次型標(biāo)準(zhǔn)化前面幾章走的是數(shù)值路線快但有浮點(diǎn)誤差。如果你需要精確結(jié)果——比如做二次型標(biāo)準(zhǔn)化、推導(dǎo)變換矩陣的解析形式——符號(hào)工具箱是更好的選擇。下面這套流程是我在需要干凈結(jié)果時(shí)的固定操作。先看一個(gè)完整的符號(hào)對(duì)角化例子A [1 2 2; 2 1 2; 2 2 1]; E eye(3); syms x f det(x*E - A) % 特征多項(xiàng)式符號(hào)形式 eigenvalues solve(f) % 精確特征值f會(huì)返回x^3 - 3*x^2 - 9*x - 5solve給出5, -1, -1。接下來(lái)對(duì)每個(gè)特征值求特征向量% 特征值 5 對(duì)應(yīng)的特征向量 y 5; B y*E - A; b1 sym(null(B)) % 基礎(chǔ)解系符號(hào)形式 % 特征值 -1 對(duì)應(yīng)的特征向量 y -1; B y*E - A; b2 sym(null(B)) % 返回兩個(gè)基向量null(B)返回的是齊次方程組B*Z 0的基礎(chǔ)解系。對(duì)于二重特征值 -1b2是一個(gè) 3x2 的矩陣兩列分別是兩個(gè)線性無(wú)關(guān)的特征向量。把b1和b2拼起來(lái)得到變換矩陣TT [b1, b2] % 所有特征向量拼成的矩陣 D T^-1 * A * T % 驗(yàn)證對(duì)角化結(jié)果應(yīng)該是精確對(duì)角陣D會(huì)返回一個(gè)精確的對(duì)角陣對(duì)角線上是5, -1, -1非對(duì)角元嚴(yán)格為零。這就是符號(hào)運(yùn)算的價(jià)值——沒(méi)有浮點(diǎn)噪聲結(jié)果可以直接寫(xiě)進(jìn)推導(dǎo)過(guò)程。再看二次型標(biāo)準(zhǔn)化的場(chǎng)景。給定二次型矩陣A要把它化成標(biāo)準(zhǔn)形步驟是步驟操作對(duì)應(yīng)代碼1確認(rèn)矩陣對(duì)稱isequal(A, A)2求特征值和正交特征向量[P, D] eig(A)3構(gòu)造正交變換x P * y4代入二次型得到標(biāo)準(zhǔn)形y * D * y5驗(yàn)證變換矩陣正交性P * P接近單位陣這套流程的關(guān)鍵在第 2 步實(shí)對(duì)稱矩陣的eig返回的P一定是正交陣所以P^-1就是P不用真的求逆。如果你拿到的P不正交檢查矩陣是否真的對(duì)稱——數(shù)值上A - A的范數(shù)應(yīng)該在 1e-12 以下才算對(duì)稱。最后一個(gè)技巧是關(guān)于精度控制的。當(dāng)你需要高精度數(shù)值結(jié)果時(shí)vpa配合eig的符號(hào)輸出能派上用場(chǎng)A [1 2 2; 2 1 2; 2 2 1]; [V, D] eig(sym(A)); % 符號(hào)形式求特征值 V_num double(V); % 轉(zhuǎn)回雙精度 D_num double(D); % 需要更高精度時(shí)用 vpa V_high vpa(V, 32); % 32 位有效數(shù)字vpa的第二個(gè)參數(shù)控制有效數(shù)字位數(shù)默認(rèn)是 32。位數(shù)越高計(jì)算越慢一般 32 位足夠工程使用。我通常在驗(yàn)證理論推導(dǎo)時(shí)才用vpa日常計(jì)算double就夠了。從那以后我每次做對(duì)角化驗(yàn)證都強(qiáng)制走一遍符號(hào)轉(zhuǎn)數(shù)值再對(duì)比的流程——先用sym算精確結(jié)果再用double轉(zhuǎn)回來(lái)對(duì)比數(shù)值解差異超過(guò) 1e-10 就回頭查矩陣條件數(shù)。這個(gè)習(xí)慣幫我逮住過(guò)好幾次隱藏的虧損矩陣問(wèn)題。希望幫到你。本文還有配套的精品資源點(diǎn)擊獲取