分享一下我老師大神的人工智慧教程！零基礎，通俗易懂！http://blog.csdn.net/jiangjunshow

也歡迎大家轉載本篇文章。分享知識，造福人民，實現我們中華民族偉大復興！

               

第十章 聯立方程式解

解聯立方程式為線性代數之第一步，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=A/C ：使用左除法求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.3769
17.9531  -22.8718    8.6384    0.7080

根據定義，原矩陣與其反矩陣相乘應等於單位矩陣，且相乘之順序不拘，亦即：Ａ*B與B*A均應等於Ｉ：

B*A
ans =
 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     1

d=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.1987
V =
 0.4053    0.6812    0.0680
 0.9136   -0.7124   -0.4004
 0.0319   -0.1688    0.9138
D =
 1.2167         0         0
      0    0.0068         0
      0         0    0.1987
A*V

ans =
 0.4931    0.0046    0.0135
 1.1116   -0.0048   -0.0796
 0.0388   -0.0011    0.1816
V*D

ans =
 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 線性聯立方程式矩陣解法

所謂線性方程式係指方程式中之變數僅屬一階者，其幕次不能大於一，且任何項中不得有兩變數相乘或相除之情形。在矩陣表示法上，通常採用[A][x]=[b]之型式，其解為[x]=[A]/[b]或[x]=[b]/[A]。

對於線性聯立方程式目前已經發展出許多種解法，其中包括變數迭代消去法及克雷蒙法(Creamer's 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求解時，只要用矩陣倒除即可，或稱為左除法。例如一組聯立方程式以矩陣表示為[A][X]=[C]時，其未知數項[X]可以利用MATLAB倒除的指令求得，即[X]＝[A]/[C]：

A=[3 4;5 -2];
C=[10;8];
X=A/C 

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=A/C
C=A*x 

x =
   0.0479
   2.1737
   3.6467
C =
  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 範例一：桁架負重分析

範例一：有一重物Ｗ架在天花板之ＡＢ兩點上，其剛性支架ＡＣ與ＢＣ與點Ｃ處以梢聯結並由此吊掛重物。支架ＡＣ與ＢＣ與天花板間之角度分別為a₁與b₂。試求其間受力之關係。


解：

設支架ＡＣ與ＢＣ上之支撐力分別為T1、T2，則依直角座標分析，其作用在點Ｃ之X與Y方向力的關係如下：

x方向:　　T1cos(a1)-T2cos(b2)=0
y方向:　　T1sin(a1)+T2sin(b2)=W

以矩陣表示，其AT=W之型式如下：

　　
因此，由T=A/W之左除法可以得到答案。其相關程式如下：

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=A/W;T1=T(1);T2=T(2);

執行例：

設夾角a1=30度, b2=60度, W=100公斤，求張力T1及T2。
>> [T1,T2]=findT(30,60,100)
T1 =           50
T2 =       86.603

得到AC、BC桿上之張力分別為50公斤與86.603公斤。

 

張貼者： Martin Fon 位於 11/22/2006 11:20:00 下午 0 意見   此文章的連結

標籤： Chap10

 

10.2 範例二：電路分析

範例二：有一電路如下圖，試寫一MATLAB程式，求解各電阻器上之電流。

圖3.電路圖


就克希荷夫(Kirchhoff)定律，可以在電路上任選定三個迴圈建立電阻與電壓降三項聯立方程式，及電流經Ａ、Ｂ兩點分流時建立之二項聯立方程式，其結果如下：




整理此組聯立方程式，其變數Ｉ有五項，其方程式有五組，故應可得解。茲以矩陣[R][I]=[V]表示，利用左除法即為[I]=[R]/[V]。其各變數矩陣分別為：



　　　　　

程式內容：

function Current=findC(v1,v2,r)
% Prog using Kirchhoff's 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=R/V;

執行例：設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水柱高表示其範圍100mmH ₂O；風量以cmm表示。依風阻之定義，其關係式為P=QR，其單位應為mmH ₂O/cmm。實際上之風阻可參考穀物乾燥之相關書籍。其值域約為3,000-5,000mmH ₂O/cmm之間，榖物之風阻則此值之數十倍。故若以１風歐姆（ＷΩ）＝１mmH ₂O/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=R/P;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。熱阻的觀念可以用來類比電路，並進行解析。圖為一道牆由木板、玻璃纖維、水泥及磚牆等材料。這些材質之熱導係數可參閱此網站 http://hyperphysics.phy-astr.gsu.edu/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

熱阻與電阻之類比關係


將各層產生之熱阻，以電路的方式進行串聯，即可形成類似網路，並計算各層之溫度。其各層間之關係式如下：



　　 　　
以矩陣表示，即為[A][T]=[C]。其溫度T可解如後：

    [T]=[A]/[C]

程式說明

程式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, W[m.C]
%   thick:thickness of each layer, mm
%   area:the crosssection area, m^2
% Outputs:
%   temp:[heatflow & temperatures at each layer, W/m^2,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, 2006
if nargin<5, area=1;end
R=thick./k/area/1000;
nn=length(thick);C=zeros(1,nn);
CC=-eye(nn)+[[zeros(nn-1,1) eye(nn-1)];C];
A=[R',CC(:,2:end)]
C=[ti C(2:end-1) -to]';
x=[0 cumsum(thick)];
TT=A/C;q=TT(1);T=[ti TT(2:end)' to];
line(x,T,'marker','s')
xlabel('distance,mm')
ylabel('Temperature,C')

執行例：

設室內溫度為攝氏25度，室外為-10度。牆壁由內而外為木板(50mm)、玻璃纖維(100mm)、水泥面(20mm)、磚牆(200mm)，其對應之熱導係數分別為[0.12 0.04 0.8 0.6]。以此代入heat_wall程式中執行，得其結果：

>> [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.0000

T =
20.0000   19.0217   -5.4348   -9.3478  -10.0000

x =
  0    10   135   195   245

q =
 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.0000

T =

20.0000   13.4307   -2.9927   -5.6204   -6.0584  -10.0000

x =
  0   100   225   285   335   395

q =

 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=A/C;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=A/C;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 =     0
ans =     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 =     3
　R2 =     3

由於兩者之階值相同，表示其解存在且為唯一。下面就一個多數解的例子進一步說明：

　A=[3 2 1; 10 -25 5];C=[5000 2000]';
　R1=rank(A),R2=rank([A C])
　R1 =     2
　R2 =     2

顯然這組方程式應有解，但由於階數2小於變數個數3，故其解為多數解。利用左除法亦可得到其中一解，即令T3等於零時，所得之答案。其過程如下：

　T=A/C
　T =
     1357.9
     463.16
          0

範例六：

有一電路如圖，接於一電源Ｖ上。試求通過各電阻器上之電流。



根據克希荷夫Kirchhoff定律，可以選定三個迴圈建立屬於電壓降之三項聯立方程式，及電流經由Ａ、Ｂ、Ｃ三點產生分流時所能建立之四個關係方程式得之：



整理此組聯立方程式，其電流變數Ｉ有五項，其方程式有五組，故應可得解。茲以矩陣[R][I]=[V]表示，利用左除法即為[I]=[R]/[V]。其各變數矩陣分別為：



　 　

程式內容

function Curr=ele_amp(v,r)
% Prog using Kirchihoff's 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=R/V;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 ==> inv
Matrix 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.9
　T2-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 3]
A =
     3     5     8
     2     4     3
>> b=[100 150]'
b =
   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 x₁+b x₂+c x₃想適配下列資料：

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=A/C仍可以得到相關係數，因此得到