線性代數應用實例,取自線性代數機算與應用指導（MATLAB）版 2010.12,例 1 平板穩態溫度的計算,為了計算平板形導熱體的溫度分布，將平板劃分為許 多方格，每一個節點上的穩態溫度將等于其周圍四個 節點溫度的平均值。由此可得出階數與節點數相同的 線性方程組，方程的解將取決于平板的邊界條件。 這個方法可以用來計算飛行器的蒙皮溫度等。,平板溫度計算的模型,整理為,A=4,-1,-1,0; -1,4,0,-1; -1,0,4,-1; 0,-1,-1,4; b=30; 50; 60; 80; U=rref(A,b),MATLAB 程序(ma1),運行結果為： U = 1.0000 0 0 0 21.2500 0 1.0000 0 0 26.2500 0 0 1.0000 0 28.7500 0 0 0 1.0000 33.7500,向高階系統擴展,則要解 25 階的線性方程組。 運行書上的程序得溫度分布 如下,將平板分割得愈細，求出的解就愈精確。如果把上 述區域分成 25 個點如右,MATLAB 程序ma2,例 2 交通流的建模,對于一個有雙向車流的十 字路口，根據流出流入車 數相等的規則，可以列出 下列方程組:,節點A：x1360x2260 節點B：x2220x3292 節點C：x3320x4357 節點D：x4260x1251 相應的矩陣方程為：,A=1,-1,0,0;0,1,-1,0;0,0,1,-1;-1,0,0,1; b=-100;72;37;-9; U=rref(A,b),MATLAB 程序(ma3),運行結果為： U = 1 0 0 -1 9 0 1 0 -1 109 0 0 1 -1 37 0 0 0 0 0,由于 U 的最后一行為全零，也就是說，四個方程中實 際上只有三個獨立。x4 可以任設，因為如果有一些車沿 此路口環行，對方程無影響，故方程組的解可如上表示.,把上述模型擴展到多個十字路，乃至整個城市，就構成高階的線性代數方程組。例如下面的 6 節點交通流圖，它就要由 6 個方程和 7 個變量來描述。用行最簡型方法可以知道，它的解將包括兩個自由變量。其物理意義類推。,向高階系統擴展,左圖描述了四個城市之間的航空 航線圖，其中1、2、3、4 表示四 個城市；帶箭頭線段表示兩個城 市之間的航線。設行號表示起點 城市，列號為到達城市，則 定義鄰接矩陣 A 為：,例 3 飛機航線問題,轉機航線的數學模型,不難證明：矩陣 A2=A*A 表示一個人連續坐兩次航班可以到達的城市，矩陣 A3=A*A*A 表示連續坐三次航班可以到達的城市：,其中，第 i 行描述從城市 i 出發，可以到達各個城市的 情況，若能到達第 j 個城市，記 A(i,j)=1，否則 A(i,j)=0， 規定 A(i,i)=0 (其中 i=1,2,3,4)。如第 2 行表示：從城市 2 出發可以到達城市 3 和城市 4 而不能到達城市 1 和 2。,多次轉機到達的城市,分析矩陣 A3 的第二行，可以得出： 某人從城市 2 出發，連續坐三次航班 可以到達城市 2、3 和城市4，不能到 達城市 1，而到達城市 3 和城市 4 的 方法各有兩種。,不難看出，轉機兩次以下的航線的航路矩陣為 At2= A+ A2 + A3 程序為（ma4） A=0,1,1,1; 0,0,1,1; 0,0,0,0; 1,1,0,0; At2=A+A2+A3,例 4 行列式的幾何應用,二階行列式的幾何意義是兩個二維向量構成的平行四邊形的面積，三階行列式的幾何意義是三個 3 維向量構成的平行六面體的體積。如下圖所示，用 MATLAB 軟件來實現面積和體積的運算。,實例 （I）已知三角形ABC三個頂點的坐標分別為：(1,2),(3,3),(4,1)，計算該三角形的面積； （II）已知凸九邊形九個頂點的坐標分別為：(0,8.5),(3,7),(6,0),(3,-4),(1,-5),(-5,-3), (-7,0),(-5,6),(-3,8),計算該九邊形的面積。 （III）在平面坐標系中畫出以上三角形和九邊形。,平行四邊形面積計算,解：(I)如圖所示，三角形 ABC 的面積就等于向量AB和向量AC所構成平行四邊形面積的一半。其中：,由向量 和 所構成的平行四 邊形的面積為行列式 的絕對值。 計算的MATLAB語句為： S=abs(a1*b2-a2*b1) 實例給出的是三角形三個頂點坐標a1,b1， a2,b2， a3,b3，求該三角形面積，則有： MATLAB寫成S=abs(det(a2-a1,b2-b1; a3-a1,b3-b1),多邊形可以劃分為多個三角形來計算。 先對三角形面積計算構成一個函數程序； 這個子程序名為：cal_area3(A,B,C) A,B,C為三個頂點的二維坐標向量 凸多邊形面積只需多次調用這個函數程序； 例如五邊形ABCDE，可由 S5= cal_area3(A,B,C)+ cal_area3(A,C,D)+ cal_area3(A,D,E) 求得。（MATLAB程序ma4） 也可由多邊形面積子程序cal_arean(A)計算。,擴展至多邊形面積計算,解：(II)如圖所示，凸九邊形面積是由9-2=7個三角形面積組成。,在MATLAB命令窗口運行程序ma5.m，即可以算出三角形和九邊形面積，同時可以得到圖形：,MATLAB 程序,function s=cal_area3(a,b,c) % a,b,c 應為同形的 2 維行向量或列向量， % 格式檢驗語句略去 ab=b-a; % 計算向量AB ac=c-a; % 計算向量AC if size(ab)=1,2 % 判讀向量AB是否為行向量 A=ab;ac; % 構造矩陣A else A=ab,ac; end s=abs(det(A)/2; % 根據公式計算三角形面積,例 5 藥方配置問題,（1）某醫院要購買這 7 種特效藥，但藥廠的第 3號和第 6 號特效藥已經賣完，請問能否用其它特效藥配制出這兩種脫銷的藥品。 分析：即 3, 6 向量與其他向量是否線性相關 （2）現在該醫院想用這 7 種中草藥配制三種新的特效藥，下表為新藥所需的成分質量 (單位: 克) 。請問如何配制。 分析：這是新藥向量與原來藥向量是否線性相關的問 題。,問題及分析思路,新藥的成分要求,u1=10;12;5;7;0;25;9;6;8; u2=2;0;3;9;1;5;4;5;2; u3=14;12;11;25;2;35;17;16;12; u4=12;25;0;5;25;5;25;10;0; u5=20;35;5;15;5;35;2;10;0; u6=38;60;14;47;33;55;39;35;6; u7=100;55;0;35;6;50;25;10;20; U1=u1,u2,u3,u4,u5,u6,u7 V1,r=rref(U1),問題 (1) 的 MATLAB 程序（ma6),運行結果,V1 = 1 0 1 0 0 0 0 0 1 2 0 0 3 0 0 0 0 1 0 1 0 0 0 0 0 1 1 0 0 0 0 0 0 0 1 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 r = 1 2 4 5 7,可見這七種特效藥是“相關的”, 3、6 兩種藥可用其它 5種藥線性配制出來, 但第1 、2 、 4、5 、7 種藥“無關”。,因此，8，9 兩種藥可以配出，第 10 種藥則不能配出。,V2 = 1 0 1 0 0 0 0 1 3 0 0 1 2 0 0 3 0 3 4 0 0 0 0 1 0 1 0 2 2 0 0 0 0 0 1 1 0 0 0 0 0 0 0 0 0 0 1 0 1 0 0 0 0 0 0 0 0 0 0 1 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 s = 1 2 4 5 7 10,為求第二個問題, 把 3 種新藥與 7 種原藥組成矩陣 U2， 求 rref，得：,s1=40 62 14 44 53 50 71 41 14 s2=162 141 27 102 60 155 118 68 52 s3=88 67 8 51 7 80 38 21 30 U2=U1,s1,s2,s3 V2 r=rref(U2),問題 (2) 的 MATLAB 程序,假設一個城市的總人口數是固定不變，但人口的分布 情況變化如下：每年都有 5 的市區居民搬到郊區； 而有 15 的郊區居民搬到市區。若開始有 700000 人 口居住在市區，300000 人口居住在郊區。請分析： （1）10 年后市區和郊區的人口各是多少？ （2）30 年后、50 年后市區和郊區的人口各是多少？ （3）分析（2）中數據相似的原因。,例 6 人口遷徙問題,解 這個問題可以用矩陣乘法來描述。令人口變量,其中 xn 為市區人口所占比例，yn 為郊區人口所占比例。 在 n+1年的人口分布狀態為：,用矩陣乘法可寫成：,A=0.95,0.15;0.05,0.85; X0=700000;300000; X10=A10*X0,開始市區和郊區的人口數為,可以得到 n 年后市區和郊區的人口分布：,因此 (1) 10 年后的人口可用程序計算如下：,運行結果為：,故市區和郊區人口數約為：744630和255370。,無限增加時間 n，市區和郊區人口之比將趨向常數 0.75/0.25。為了弄清為什么它趨向于一個穩態值，需要求 Ak, 為此可先將 A 對角化, 再求其 k 次冪。,對角矩陣的冪次可以化為元素的冪次,余下很容易計算。,% 分析 n 年后城市人口分布(ma7) clear A=0.95,0.15; 0.05,0.85; X0=700000; 300000; P,lambda=eig(A); syms n % 定義符號變量 n Xn=P*lamda.n*inv(P)*X0,MATLAB 程序,顯然, 隨 n 增大 (4/5)n 趨近于零, 而 Xn 趨于,運行結果為：,例 7 多項式插值與擬合,求: (1) 過這五個點作一個四次多項式函數,(2) 請根據這五個點，擬合一個二次多項式函數,下表給出了平面坐標系中五個點的坐標。,并求 x=5 時的函數值 p4(5)。用 MATLAB 繪制多項式函 數 p4(x) 的曲線、已知點及插值點 (5, p4(5)。,并用 MATLAB 繪制 p2(x) 的曲線及已知的五個點。,其中矩陣：,解：(1) 根據已知條件，把五個點的坐標值分別代入 四次多項式函數，可以得到如下線性方程組：,系數矩陣 A 的行列式為范德蒙 (Vandermonde) 行列式， 且五個坐標點的橫坐標各不相同，則該行列式不等于零， 所以方程組有唯一解。,MATLAB 程序:(ma8) x=0;1;2;3;4; % 輸入已知點坐標 y=-27;0;21;0;-75; A=x.0, x.1, x.2, x.3, x.4; % 構造 vandermonde 矩陣 a=Ay; % 得到適定方程組的唯一解 a,運行程序，得到 a(1)=-27, a(2)=12, a(3)=26, a(4)=-12, a(5)= 1.,把五個點的坐標值分別代入二次多項式函數，可以得 到如下線性方程組：,其中，,(2) 多項式擬合要解一個超定方程,該方程組有三個未知數，但有五個方程，進一步分析 可以得到該方程組無解，即不存在一個二次多項式曲 線剛好能過已知的五個點。MATLAB 軟件提供了一個 利用最小二乘法解決超定方程組解的方法。求系數的 公式也是 a = Ay，以找到一條二次曲線來近似地描述 已知 5 個點的變化情況。,對比插值和擬合的曲線如下圖,用平面坐標系中的一個閉合圖形來描述剛體，用一個矩 陣 X 來表示它。X 的一列表示剛體一個頂點的坐標。為 了使圖形閉合，X 的最后一列和第一列相同；為了實現 剛體的平移運算，給矩陣 X 添加元素值都為 1 的一行， 使 X 為 3n 矩陣。,若有矩陣：,則可以證明，矩陣 Y1 是剛體 X 沿 x 軸正方向平移 c1， 沿 y 軸正方向平移 c2 后的結果；矩陣 Y2 是剛體 X 以坐 標原點為中心逆時針轉動 t 弧度的結果。,例 8 剛體的平面運動,實 例,用下列數據表示字母 A:,對 A 進行以下平面運動, 并繪制移動前后的圖形。,(1) 向上移動 15, 向左移動 30; (2) 逆時針轉動 /3; (3) 先逆時針轉動3 /4, 然后向上平移 30, 向右平移 20。,解 構造剛體矩陣 X，平移矩陣及轉動矩陣。,MATLAB 程序（ma9), X=0,4,6,10,8,5,3.5,6.1,6.5,3.2,2,0; 0,14,14,0,0,11,6,6,4.5,4.5,0,0; ones(1,12); % 構造剛體矩陣 X M1=1,0,-30; 0,1,15; 0,0,1; % 構造平移矩陣 M1 Y1=M1*X; % 計算平移結果, fill(Y1(1, :), Y1(2, :), red); % 繪制平移后剛體 Y1, plot(X(1, :), X(2, :); % 繪制原來剛體 X hold on axis equal, R1=cos(pi/3), -sin(pi/3), 0; sin(pi/3), cos(pi/3), 0; 0,0,1; % 構造轉動矩陣 R1 Y2=R1*X; % 計算旋轉結果 fill(Y2(1, :), Y2(2, :), blue); % 繪制轉動后剛體 Y2, fill(Y3(1, :), Y3(2, :), black); % 繪制轉動及平移后剛體 grid on hold off, M2=1,0,20; 0,1,30; 0,0,1; % 構造轉動矩陣 M2 R2=cos(3*pi/4), -sin(3*pi/4), 0; % 構