MATLAB方程式求解

最新推荐文章于 2024-07-29 00:30:00 发布

Augusdi

最新推荐文章于 2024-07-29 00:30:00 发布

阅读量2.2w

点赞数

文章标签： matlab input matrix function c each

本文链接：https://blog.csdn.net/Augusdi/article/details/4063353

版权

本文详细介绍了如何使用MATLAB解决线性联立方程组，涵盖了矩阵运算、线性方程组的解法，包括左除法、反矩阵法，并通过多个实际例子，如桁架负重分析、电电路分析和热传导问题，来展示MATLAB在工程问题中的应用。此外，还探讨了稀疏矩阵和相关指令在数值分析中的重要性。

摘要由CSDN通过智能技术生成

第十章聯立方程式解

解聯立方程式為線性代數之第一步，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) ：簡化方程式型式

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

其使用範圍可以參考線性代數的內容。

所謂線性方程式係指方程式中之變數僅屬一階者，其幕次不能大於一，且任何項中不得有兩變數相乘或相除之情形。在矩陣表示法上，通常採用[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，其結果也正確。

10.2 範例一：桁架負重分析

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

解：

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


x方向:　　T₁cos(a₁)-T₂cos(b₂)=0
y方向:　　T₁sin(a₁)+T₂sin(b₂)=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公斤。

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

這是分別在各分路上之電流值，其中有負號者表示與原先設定之方向相反。

10.2 範例三：穀倉通風系統分析

範例三：

最低0.47元/天解锁文章