MATLAB方程式求解.doc

matlab方程式求解 第十章 聯立方程式解 解聯立方程式為線性代數之第一步，matlab 在這方面具有強大的解題功能。利用矩陣除法可以在解題技術上發揮相當大的功效。在工程上，一般線性聯立方程式常用作結構分析或機構設計，計算桁架之應力等。其他如化學程序中之質能平衡、電子學中之電路分析等等均可利用線性方程式得解。在前述之二維繪圖指令中，部份可配合統計上的需要加以應用，對於某些資料之說明甚有幫助。有關線性聯立方程式解題之指令可以參考：inv(a) ：矩陣之反矩陣，又稱為 ，需為方矩陣 det(a) ：矩陣之行列式值，需為方矩陣 =pinv(a) ：虛擬反矩陣，為與a同大小之矩陣，且a*x*a=a及x*a*x=x rank(a) ：矩陣之階數 x=inv(a)*c ：使用反矩陣求ax=c之解 x=ac ：使用左除法求ax=c之解 x=d/c ：使用右除法求x=之解 b=rref(a) ：簡化方程式型式 張貼者： martin fon 位於 11/22/2006 11:56:00 下午 0 意見 此文章的連結 標籤： chap10 10.1 線性矩陣應用指令 線性代數以矩陣表示，會有諸多特殊名稱，在應用上相當普遍，有必要進行瞭解。有些部份在前面之說明中業已經提到，但在這裡則特別另加說明。行列式與反矩陣在工程數學或線性代數中，行列式與反矩陣是常用的矩陣特性。尤其在討論聯立方程式之解的過程中更常用到。反矩陣有如一個矩陣的倒數，一個方矩陣若為，則其反矩陣可以a-1表示，而其與原矩陣之關係為： aa-1=i a-1a=i其中，稱為單位矩陣(identity matrix)，其大小為方矩陣，而對角線元素值均為1。在matlab中有一個指令稱為eye，可以用以建立這種單位矩陣。例如：eye(4)ans = 1 0 0 0 0 1 0 0 0 0 1 0 0 0 0 1矩陣之反矩陣以inv(a)表示，有如一個數之倒數一樣。例如，設為一魔術方陣magic(4)：a=rand(4)a = 0.9355 0.0579 0.1389 0.2722 0.9169 0.3529 0.2028 0.1988 0.4103 0.8132 0.1987 0.0153 0.8936 0.0099 0.6038 0.7468其反矩陣即為inv(a)表示，其結果設為，則：b=inv(a)b =-1.5054 3.6825 -1.4860 -0.4013 5.3255 -7.0376 3.9063 -0.1473-20.0637 22.9531 -8.5486 1.376917.9531 -22.8718 8.6384 0.7080根據定義，原矩陣與其反矩陣相乘應等於單位矩陣，且相乘之順序不拘，亦即：*b與b*a均應等於：b*aans = 1.0000 0.0000 0.0000 -0.0000 0.0000 1.0000 -0.0000 0.0000-0.0000 0.0000 1.0000 -0.0000-0.0000 0.0000 -0.0000 1.0000上述結果即為單位矩陣，可以試試以a*b，其結果應相同。此處、及均為方矩陣，且矩陣之行列式值或det(a)需不得為零，否則其反矩陣不能存在。此種矩陣不存在的情形，稱為奇異矩陣(singular)。一個矩陣若具奇異特性，則其行列值為零，例如：d=magic(4)d = 16 2 3 13 5 11 10 8 9 7 6 12 4 14 15 1d=det(d)d = 0行值式之值為零，表示其反矩陣不存在，故即使用inv(d)指令也會產生一些數字，但其結果並不可靠，而且會有一些警告訊息出現，例如：dd=inv(d)warning: matrix is close to singular or badly scaled. results may be inaccurate. rcond = 1.306145e-017.dd =1.0e+014 * 0.9382 2.8147 -2.8147 -0.9382 2.8147 8.4442 -8.4442 -2.8147-2.8147 -8.4442 8.4442 2.8147-0.9382 -2.8147 2.8147 0.9382上述結果之值也變成很大，顯然不是正確的值。不信的話可以用aa-1=i， a-1a=i 印證一下（在此至少證實一下垃圾進拉圾出(gabage in, gabage out這句名言）。所以，某矩陣是否為奇異矩陣，可以使用det進行檢驗。其值若為零則屬奇異矩陣。特徵值及特徵向量(eigen value & eigen vector)線性代數中，最重要的工具是利用特徵值與特徵向量（eigenvalues & eigenvectors）來說明方矩陣之多項式特性，此兩者均為方型矩陣。利用矩陣之操作，可以作出相等效果之矩陣等式，但型式更為簡化，或容易得解。設特徵矩陣值與特徵向量其分別為d與v，則其與原矩陣之關係如下： a * v= v * d假設有一個方矩陣係由亂數函數產生，其與可以利用eig()指令如下求得：a=rand(3);v, d= eig(a)a = 0.4103 0.3529 0.1389 0.8936 0.8132 0.2028 0.0579 0.0099 0.1987v = 0.4053 0.6812 0.0680 0.9136 -0.7124 -0.4004 0.0319 -0.1688 0.9138d = 1.2167 0 0 0 0.0068 0 0 0 0.1987a*vans = 0.4931 0.0046 0.0135 1.1116 -0.0048 -0.0796 0.0388 -0.0011 0.1816v*dans = 0.4931 0.0046 0.0135 1.1116 -0.0048 -0.0796 0.0388 -0.0011 0.1816由上述的演算，得知a*及*是一樣的。eig()之指令尚有另外一種型式，即：v,d = eig(a,b)其中，分別為兩矩陣，為對角矩陣特徵值，而為對應之特徵向量。其關係如下： a * v = b * v * d其使用範圍可以參考線性代數的內容。 張貼者： martin fon 位於 11/22/2006 11:40:00 下午 0 意見 此文章的連結 標籤： chap10 10.2 線性聯立方程式矩陣解法 所謂線性方程式係指方程式中之變數僅屬一階者，其幕次不能大於一，且任何項中不得有兩變數相乘或相除之情形。在矩陣表示法上，通常採用ax=b之型式，其解為x=ab或x=b/a。對於線性聯立方程式目前已經發展出許多種解法，其中包括變數迭代消去法及克雷蒙法(creamers method)。目前matlab 所使用之解法則以變數迭代消去法為主，或稱為高斯(gauss limination)消去法。這是利用兩線性方程式分別乘以某特定常數，使其與另一方程式之同一變數係數相同，因而兩式相減得以消除該變數。如此展轉消除，最後可以得聯立方程之解。茲舉例說明如下： 3x +4y =10 5x -2y =8上述聯立方程式中，先將第一式乘兩邊5、第二次兩邊乘3，兩式相減，即可得： 5(4y) - 3(-2y) = 5* 10 - 3*8由上式整理之後，可得果為y = 1；代入原式任一式可得x = 2，終得解。如果利用第九章之繪圖指令可以繪出此兩條曲線，由其交點即可得到相同的解。 ezplot(3*x+4*y-10,-10,10);hold on ezplot(5*x-2*y-8,-10,10) 圖10.1利用這樣的解法並不是解線性方程式之重點，因為線性代數方程式可用矩陣的型式表示，甚至以矩陣之乘除法可以得解。以上面之聯立線性方程式為例，可以化成之型式，設x=x1，y=x2：左除法求解利用matlab求解時，只要用矩陣倒除即可，或稱為左除法。例如一組聯立方程式以矩陣表示為ax=c時，其未知數項x可以利用matlab倒除的指令求得，即xac：a=3 4;5 -2;c=10;8;x=ac x = 2 1 結果立即得到x1=2、x2=1。反矩陣法求解解上式矩陣有時可用反矩陣法。反矩陣有倒數的意義，但在矩陣中必須為方矩陣，且須為非特異矩陣(non-singular)。通常用(-1)次方表示之，如a-1 。其基本特質為與原矩陣相乘後，將得到單位矩陣： a-1a=aa-1=i以此乘於之兩側，可得： a-1ax=a-1c 或a-1ax=ix=x= a-1c在matlab反矩陣a-1 之求法為inv(a)；或稱為單位矩陣，表示為eye(a)。就上面之聯立方程式之解，可以用matlab演算如下：a=3 4;5 -2;c=10;8;x=inv(a)*c x = 2.0000 1.0000 所得的結果與前述之分析相同。理論上雖然如此，在演算時求反矩陣較費時間，使用左除法是為較佳的操作法。範例：有一組聯立方程式如下： 6x - 3y +5z = 12 10x + 4y -8z =-20 -6x + 2y +3z = 15解：先安排及矩陣，再利用左除法求得x。a=6 -3 5;10 4 -8;-6 2 3;c=12 -20 15;x=acc=a*x x = 0.0479 2.1737 3.6467c = 12.0000 -20.0000 15.0000 由上述的計算，可以迅速得到結果，即x= 0.0479, y= 2.1737, z= 3.6467。通常利用matlab解聯立方程式，亦會碰到無解的情形，此時或稱為特異狀況（或 |a| 為零的情形）。前述之左除結果，可用以反求值，x，其結果也正確。 張貼者： martin fon 位於 11/22/2006 11:30:00 下午 0 意見 此文章的連結 標籤： chap10 10.2 範例一：桁架負重分析 範例一：有一重物架在天花板之兩點上，其剛性支架與與點處以梢聯結並由此吊掛重物。支架與與天花板間之角度分別為a1與b2。試求其間受力之關係。解：設支架與上之支撐力分別為t1、t2，則依直角座標分析，其作用在點之x與y方向力的關係如下：x方向:t1cos(a1)-t2cos(b2)=0y方向:t1sin(a1)+t2sin(b2)=w以矩陣表示，其at=w之型式如下：因此，由t=aw之左除法可以得到答案。其相關程式如下：function t1,t2=findt(a1,b2,w)% program to find tensions in supports% inputs:% a1,b2:inclined angles, in degrees% w:weight, in kg% ouputs:t:tensions in supports, kg% example: t=findt(30,60,100)a=cosd(a1) -cosd(b2); sind(a1) sind(b2);w=0 w;t=aw;t1=t(1);t2=t(2);執行例：設夾角a1=30度, b2=60度, w=100公斤，求張力t1及t2。 t1,t2=findt(30,60,100)t1 = 50t2 = 86.603得到ac、bc桿上之張力分別為50公斤與86.603公斤。 張貼者： martin fon 位於 11/22/2006 11:20:00 下午 0 意見 此文章的連結 標籤： chap10 10.2 範例二：電路分析 範例二：有一電路如下圖，試寫一matlab程式，求解各電阻器上之電流。圖3.電路圖就克希荷夫(kirchhoff)定律，可以在電路上任選定三個迴圈建立電阻與電壓降三項聯立方程式，及電流經、兩點分流時建立之二項聯立方程式，其結果如下：整理此組聯立方程式，其變數有五項，其方程式有五組，故應可得解。茲以矩陣ri=v表示，利用左除法即為i=rv。其各變數矩陣分別為：程式內容：function current=findc(v1,v2,r)% prog using kirchhoffs law to find the currents in a circuit.% inputs:% v1, v2: voltage of generators. volts% r: the elements of resistances in kohms. r=r1 r2 r3 r4 r5% outputs:% current: currents through resistors c1 c2 c3 c4 c5, in ma.% example: current=findc(100,50,5 100 200 150 250)v=v1 0 v2 0 0;r=r(1) 0 0 r(4) 0;0 r(2) 0 -r(4) r(5);0 0 -r(3) 0 r(5);.1 -1 0 -1 0;0 1 -1 0 -1;current=rv;執行例：設v1=100v，v2=75v，r=80 120 250 150 200 k，其電流之安培數(ma)分別如下： current=findc(100,75,80 120 250 150 200)current = 0.5082 0.1126 -0.1166 0.3956 0.2292這是分別在各分路上之電流值，其中有負號者表示與原先設定之方向相反。 張貼者： martin fon 位於 11/22/2006 11:00:00 下午 0 意見 此文章的連結 標籤： chap10 10.2 範例三：穀倉通風系統分析 範例三：圖.穀倉通風系統某一組穀倉乾燥系統，其風量係由一台送風機以管路分送二倉，並經過穀倉筒內之穀物以進行乾燥。為簡化整個系統之分析，可以將其類比為一般之串並聯電路。整個系統因此可視為一個電阻分路，各風管對所送之風量均有一定的阻力，而風機所產生之壓力則可類比為電壓。在實際的情況下，其間之關係並不盡然為線性，不過在低風量下，將其關係視為線性應不致產生太大的誤差，其結果應可接受。分析：此種類比係簡化風壓、風量與風阻間之關係為線性，故可以利用電路進行類比模擬。在系統中，送風機係將大氣壓的空氣壓縮，提高其風壓。故若以大氣壓比擬為電路之接地時，高壓空氣經過風管、穀倉中之穀層，其後將再排放至大氣中，其最後風壓又恢復為大氣壓。因此，以電路模擬時，可改用下圖所示。其中r1為主管之風阻，r2與r3則為分岐管之風阻，然後通過穀層，其風阻較大，分別為r4與r5。圖5.通風系統之電路模擬在單位方面，風壓常以mm水柱高表示其範圍100mmh2o；風量以cmm表示。依風阻之定義，其關係式為p=qr，其單位應為mmh2o/cmm。實際上之風阻可參考穀物乾燥之相關書籍。其值域約為3,000-5,000mmh2o/cmm之間，榖物之風阻則此值之數十倍。故若以風歐姆（）mmh2o/cmm，則可使用為單位計算。利用兩迴圈及點之風流量建立聯立方程式之關係式，其相關程式演算如下： 程式內容：function qflow=duct_flow(p,r)% prog using circuit to find airflows in a system.% inputs:% p: prssure of fan. mmh2o% r: wind resistances r=r1 r2 r3 r4 r5, in kwohms.% outputs:% qf: airflow, in cmm q1 q2 q3% example: qf=duct_flow(100,1 2 2 10 10)p=p 0 0;r=r(1) 0 r(2)+r(4);0 r(3)+r(5) -r(2)-r(4);1 -1 -1;q=rp;qflow=q;執行例： qf=duct_flow(100,1 2 2 10 10)qf = 14.286 7.1429 7.1429此處風量之單位為cmm。 張貼者： martin fon 位於 11/22/2006 10:50:00 下午 0 意見 此文章的連結 標籤： chap10 10.2 範例四：穩態熱傳導的問題 穩態熱傳導的問題熱傳導是溫差產生的熱能流動，其流動速率與材料之熱阻係數有關。其公式如下： q = q/t = k a (th-tc)/d = t/(d/ka) = /r其中d為材料之厚度，為截面積，k為熱導係數，其單位為w/mk。r則為熱阻，等於d/ka。熱阻的觀念可以用來模擬電路，並進行解析。圖為一道牆由木板、玻璃纖維、水泥及磚牆等材料。這些材質之熱導係數可參閱此網站/hbase/thermo/heatcond.html#c1。其他相關值如下：材質名稱 熱導係數w/mk 磚牆 0.6 木材 0.12-0.04 石膏 0.08 玻璃纖維 0.04 玻璃 0.8 水泥 0.8 鐵板 50.2 水 0.6 熱阻與電阻之類比關係將各層產生之熱阻，以電路的方式進行串聯，即可形成類似網路，並計算各層之溫度。其各層間之關係式如下： 以矩陣表示，即為at=c。其溫度t可解如後： t=ac程式說明程式heat_wall.m可以執行上述之操作。輸入項包括內外溫度(ti,to)、各層熱導係數(k)、厚度(thick)及截面積(area)等。執行後即可得到溫度分佈及熱流(q)等。程式內容function t,x,q=heat_wall(ti,to,k,thick,area)% prog calculating heat transfer through a wall.% inputs:% ti,to: inside & outside temperature, c% k: thermoconductivities of each layer, wm.c% thick:thickness of each layer, mm% area:the crosssection area, m2% outputs:% temp:heatflow & temperatures at each layer, w/m2,c% example:% t,x,q=heat_wall(20,-10,0.08 0.04 0.12 0.6,.% 10 125 60 50)% designed by d.s. fon. date: nov. 26, 2006if nargin t,x,q=heat_wall(20,-10,0.08 0.04 0.12 0.6, 10 125 60 50)a = 0.1250 1.0000 0 0 3.1250 -1.0000 1.0000 0 0.5000 0 -1.0000 1.0000 0.0833 0 0 -1.0000t =20.0000 19.0217 -5.4348 -9.3478 -10.0000x = 0 10 135 195 245q = 7.8261其溫度分佈則如下圖：沿軸溫度分佈圖若內外側均採用木板，但內側厚度為100mm，外側為60mm，則其溫度分佈及熱流量如下： t,x,q=heat_wall(20,-10,0.08 0.04 0.12 0.6 0.08, 100 125 60 50 60)a = 1.2500 1.0000 0 0 0 3.1250 -1.0000 1.0000 0 0 0.5000 0 -1.0000 1.0000 0 0.0833 0 0 -1.0000 1.0000 0.7500 0 0 0 -1.0000t =20.0000 13.4307 -2.9927 -5.6204 -6.0584 -10.0000x = 0 100 225 285 335 395q = 5.2555其溫度變化雖不大，但熱流量則有顯著的減少。 張貼者： martin fon 位於 11/22/2006 10:40:00 下午 0 意見 此文章的連結 標籤： chap10 10-2 範例五：簡易之熱流分析 範例五：熱流問題熱流的問題有時是一個平面分佈的情形。以下圖為例，在一個侷限的空間中，設將其分為四等分，每等分之接觸面積均相同。進口處之溫度為i，出口溫度為to，其餘外圍部份均為絕熱狀態，熱流停止，故接觸之該面溫度不起變化。同樣，若將熱流比做電流，則可以形成類比電路如圖(b)，設其進出介面時之熱阻以表示。由於材質假設相同，其截面積也相同，故每一介面之熱阻也應相同。由電路之網絡觀之，熱流由外界流入t1節點，再分兩路流入t2與t3兩節點，最後滙流至t4，由此流至外界to。注意t1不能直接流至t4，因為兩者之接觸面積為零。因此由電路網絡很容易解出此並聯電路之解。由於熱阻均相同，在公式中可以不考慮。其熱流之平衡式可表列如下：根據熱流之關係，可以改變為：由上式得到的結論是任何一點之溫度應為其相連週圍之溫度點之平均值。其聯立方程式最後為下列型式：程式內容function t=heat_plate(ti,to)% prog calculating heat transfer through a plate.% inputs:% ti,to: input & output temperature, c% outputs:% temp:temperatures at each layer,c% example:% t=heat_plate(20,-10)c=ti 0 0 to;a=3 -1 -1 0;-1 2 0 -1;-1 0 2 -1;-1 -1 0 3;tt=ac;t=ti tt to;執行範例 t=heat_plate(100,-10)t = 100.0000 68.5714 52.8571 52.8571 37.1429 -10.0000深入分析如果要更確的數據，則可以增加方格，但其運算方程式之數目亦將增加。例如增加為3x3方格，則總格數為九格，必須要有九條方程式才能處理。根據前面所述之原理，某格之溫度等於其鄰近方格溫度之平均值。因此t5之溫度應為t2、t4、t6、t8等四溫度之平均值；而t3與t7則僅為鄰近兩項值之平均，其餘均為三個鄰近值之平均。若將其改為矩陣at=c之型式，則各矩陣之內容變為：程式內容 function t=heat_plate3(ti,to)% prog calculating heat transfer through a plate by 3x3.% inputs:% ti,to: input & output temperature, c% outputs:% temp:temperatures at each layer,c% example:% t=heat_plate3(20,-10)c=ti 0 0 0 0 0 0 0 to;a=3 -1 0 -1 0 0 0 0 0;-1 3 0 -1 -1 0 0 0 0; 0 -1 2 0 0 -1 0 0 0; -1 0 0 3 -1 0 -1 0 0; 0 -1 0 -1 4 -1 0 -1 0; 0 0 -1 0 -1 3 0 0 -1; 0 0 0 -1 0 0 2 -1 0; 0 0 0 0 -1 0 -1 3 -1; 0 0 0 0 0 -1 0 -1 3;tt=ac;t=reshape(tt,3,3);執行例： t=heat_plate3(20,-10)t =12.4272 8.2767 6.0922 9.0049 6.3107 3.9078 6.5291 4.0534 -0.6796張貼者： martin fon 位於 11/22/2006 10:30:00 下午 0 意見 此文章的連結 標籤： chap10 10.3 範例六、det & rank 指令之應用 det & rank 指令利用左除法的求解過程甚為簡單，但有時並不一定有解，有時也可能有多組解。聯立方程式中，型式上其變數與方程式之數目相同，應可以得解。但有些方程式是相依的，或者其中方程可能由其他兩式相互加減而得的結果，此時即無法無法獲得唯一解，因此過程上仍需先行檢驗。行列式值與聯立方程式求解有密切的關係。其值若為零，所得之解可能為非唯一組解；若不為零，則可能得到一組解。在matlab中，可以使用位階指令rank檢驗。舉例：一魔術函數造成之矩陣，其行列值及位階分別由det & rank 指令進行檢驗：a=magic(4);det(a)rank(a)ans = 0ans = 3矩陣雖為4x4，利用行列值檢視，其值為det(a)=0，表示無法獲得單一解。同樣利用rank位階指令檢查，其結果為3階，兩者均顯示與矩陣相關之聯立方程式可能為多組解。要檢測一組聯立方程式是否有解，使用上述rank的指令功能可以研判。例如檢測ax=c這一組聯聯立方程式，先求rank(a)與 rand(a c)之階值，後者a c稱為增廣矩陣(augrmented matrix)。若所得兩個階值相同，且其階值等於變數個數時，應為有解且其解屬唯一。若階值小於變數個數時，其解應為多數解，其變數可用兩者之差數之線性組合表示。範例：a=6 -3 5;10 4 -8;-6 2 3;c=12 -20 15;r1=rank(a)r2=rank(a c)r1 = 3r2 = 3由於兩者之階值相同，表示其解存在且為唯一。下面就一個多數解的例子進一步說明：a=3 2 1; 10 -25 5;c=5000 2000;r1=rank(a),r2=rank(a c)r1 = 2r2 = 2顯然這組方程式應有解，但由於階數2小於變數個數3，故其解為多數解。利用左除法亦可得到其中一解，即令t3等於零時，所得之答案。其過程如下：t=act = 1357.9 463.16 0範例六：有一電路如圖，接於一電源上。試求通過各電阻器上之電流。根據克希荷夫kirchhoff定律，可以選定三個迴圈建立屬於電壓降之三項聯立方程式，及電流經由、三點產生分流時所能建立之四個關係方程式得之：整理此組聯立方程式，其電流變數有五項，其方程式有五組，故應可得解。茲以矩陣ri=v表示，利用左除法即為i=rv。其各變數矩陣分別為： 程式內容function curr=ele_amp(v,r)% prog using kirchihoffs law to find the currents in a circuit.% inputs:% v: voltage of generators. volts% r: the resistances r=r1 r2 r3 r4 r5 r6, in kohms.% outputs:% curr: currents through resistors c1 c2 c3 c4 c5 c6, in ma% example: current=ele_amp(100,1 1 2 5 5 10)v=v 0 0 0 0 0;r=r(1) 0 0 0 r(5) r(6); 0 r(2) r(3) 0 -r(5) 0; 0 0 -r(3) r(4) 0 -r(6); 1 0 -1 0 -1 0; 0 1 -1 -1 0 0; -1 0 0 1 0 -1;current=rv;curr=current;執行結果： current=ele_amp(100,1 1 2 5 5 10)current = 8.0338 17.97 3.1712 14.799 4.8626 6.7653 current=ele_amp(100,0 1 2 5 5 10)current = 8.7356 19.54 3.4483 16.092 5.2874 7.3563張貼者： martin fon 位於 11/22/2006 10:20:00 下午 0 意見 此文章的連結 標籤： chap10 10.3 虛反矩陣指令pinv之應用 pinv指令在多數解的例子中，有時並不是僅要將其中一變數設定為零之解。為使整個系統得到最佳化，亦可利用pinv指令求得最小模組之合理解。pinv(a)又稱為虛反矩陣(pseudoinverse)，其功能與反矩陣之計算相同，但它會基於svd(a)函數(或稱奇異值分解函數)之計算方式，求得一個不是屬於全階之矩陣a之反矩陣。這是長方形矩陣求解時，在多重解中求其反矩陣之折衷方式。故若矩陣a為方矩陣或非零矩陣，則其結果應與inv(a)相同。只是在這樣的狀況，寧可使用inv(a)較為省事。處理這些長方矩陣或特異矩陣時，使用pinv(a)會有意想不到的效果。其解法是根據反矩陣法：a=3 2 1; 10 -25 5;c=5000 2000; t=inv(a)*c? error using = invmatrix must be square.t=pinv(a)*c t = 1203.9 485.16 418 上面之例因為不是方形矩陣，故求其反矩陣時會有錯誤的信息，但若用虛反矩陣指令pinv，反而相安無事，這是將t1、t2以其餘一變數t3表示之情況下，求得其最小平方之組合。其結果是否合用則端視問題之限制與應用而定。 pinv(a,tol) 之指令後面另有參數tol，可以輸入容許值。其預設值為max(size(a) * norm(a) * eps(class(a)，讀者可參考手冊之說明，以瞭解其使用方法。rref指令在處理聯立方程式時，若直接有解，應可利用左除法得到其答案。若為多數解，則需要將其重組並簡化至梯形的表列方式，將某些變數設定為自變數，並將其餘因變數之係數均轉換為1。此時可以利用rref(a c)這個指令達到目的。這個指令執行結果，會產生一個增廣矩陣，其內容為a c。以前述之t組變數為例：m=rref(a c) m = 1 0 0.36842 1357.9 0 1 -0.052632 463.16 最後簡化之梯形聯立方程式為：t1+0.3684t3=1357.9t2-0.0526325t3=463.16此時之聯立方程式經過rref這個指令會將結果簡化。例題：下面為一農場作業之資料，比較重要的作業可以分為整地、噴藥及收割等三項。下表為甲乙兩工人對個別作業每公頃所需之作業時數。設若此兩工人今年總工作時數分別為100及150小時，試求其今年各人總作業面積為多少。小時公頃 整地作業 噴藥作業 收割作業= = = =甲工人 3 5 8乙工人 2 4 3設x, y ,z分別為整地、噴藥及收割作業之作業面積數，則甲工人一年之工作情形可以下列方程式表示：甲工人：3x + 5y + 8z=100乙工人：2x + 4y + 3z =150顯然這是一個多數解的系統。設 a=3 5 8;2 4 3a = 3 5 8 2 4 3 b=100 150b = 100 150利用增廣矩陣可以縮減成簡易之x y z系統： rref(a b)ans = 1.0000 0 8.5000 -175.0000 0 1.0000 -3.5000 125.0000解之，x + 8.5z=-175 y-3.5z =125或x = 8.5z-175, y=125 +3.5z在這種情況下，由於變數在均為正值之限制之下，故應可以得到一個範圍解。 張貼者： martin fon 位於 11/22/2006 10:10:00 下午 0 意見 此文章的連結 標籤： chap10 10.4 最小平方解法 當一個系統之聯立方程式個數大於未知數之數目時，會形成過度確定性系統(overdetermined system)。此時，反矩陣法已無法應用，只能用左除法求解。但這種解如同多數解的過程一樣，必須尋最小平方的技巧,使解能夠以最接近所有點為準，成為一般所謂之回歸公式。茲以下面例子做介紹：設有一方程式的關係式為y= a x1+b x2+c x3想適配下列資料：y x1 x2 x3 3 2 3 1 6 3 -2 4-5 -5 3 6 4 2 7 3 2 4 8 -7依問題之需要，可利用rank指令執行如下：a=2 3 1; 3 -2 4; -5 3 6; 2 7 3; 4 8 -7;c=3;6;-5;4;2;r1=rank(a), r2=rank(a c)由階數r1及r2觀之，兩者不相同，故無法以正常方式得解。然而利用左除法b=ac仍可以得到相關係數，因此得到此函數為：b=acr1 = 3r2 = 4b = 1.442 -0.067601 0.4352y= 1.442 x1 -0.067601x2+0.4352 x3然而，此結果並非其真實解，而是利用最小平方法求得。若以此解作驗證，即y=a*b，其結果為：a*bans = 3.1163 6.2018 -4.8014 3.7163 2.1806顯然與原值 3 6 -5 4 2仍有誤差。只是因採用最小平方的關係，可以求得最接近之數。最小平方的取法是將各數項之平方值之和以fmins求得最小值。即：j=i=15ax1i+bx2i+cx3i-y2最小平方之應用將在第十一章統計與迴歸部份進一步說明。 張貼者： martin fon 位於 11/22/2006 10:00:00 下午 0 意見 此文章的連結 標籤： chap10 10.5 網絡分析之應用 梯形數值分析法(trapezoidal numerical integration)數值分析法常用以解積分的問題，積分之解法依問題的性質有不同的方式。比較簡單的方式是將其分成細段，然後細段之微小面積，加總後成為最後的結果。此種積分稱為梯形法。利用梯形法可以解決複雜函數之積分問題。在matlab之中，也有相同的指令，trapz，其語法如下： z = trapz(y) z = trapz(x,y) z = trapz(.,dim)若參數中僅有項，則其分段值設為1；若區段不為，則可就其結果乘以需要之區段值。值可為複數，可以針對實數及虛數部份分開積分。若有對應之x值，則x與y應為成對的向量，此時區間可不必均等值。若為長方形矩陣，則可利用第三個參數dim決定是進行行向(=1)或列向(=2)積分。下面的程式則是就f(x)=1/x函數，在區間a=1與b=2之間作積分，設此區間分成n=10段，進行比較直接運算與使用trapz指令所得的結果：% trapezoid.m% find the area under a curve.%f=(x) 1./x;a=1;b=2;n=10;x=linspace(a,b,n);y=f(x);ts=0;for i=1:n-1, dx=x(i+1)-x(i); ts=ts+dx*(y(i+1)+y(i)/2;endarea=trapz(x,f(x);disp(using trapezodal law:,num2str(ts);disp(using matlab trapez command:,num2str(area);執行trapezoid.m程式後，結果如下： trapezoidusing trapezodal law:0.69392using matlab trapez command:0.69392兩者之結果完全吻合。讀者可使用不同之函數及區間，作不同的演算，並相印證。三角形面積前面已經討論過如何應用程式指令計算由三個節點構成的三角形面積。但在有限元素法的作業上，通常均採用直接與座標點具有關連的結果，即：這是一個行列式值，也可以看成矩陣。其對應的三角形面積det(a)/2。此時(x1, y1)