1.5 MATLAB的數值運算

MATLAB具有強大的數值運算能力，它不僅能對矩陣和向量進行相應的運算，而且也可處理多項式的解、數據分析、函數的極值、線性方程組的解、函數的微積分和函數繪圖等問題。

1.5.1 矩陣運算

MATLAB的基本數據單元是不需要指定維數的復數矩陣，它提供了各種矩陣的運算與操作。它既可以對矩陣進行整體的處理，也可以對矩陣的某個或某些元素進行單獨的處理，所以在MATLAB環境下矩陣的操作同數的操作一樣簡單。

1.矩陣的實現

在MATLAB語言中不必描述矩陣的維數和類型，它們是由輸入的格式和內容來確定的。例如，

當時，把A當做一個2 ×2維的矩陣；

當A=[12]時，把A當做一個2維行向量；

當時，把A當做一個2維列向量；

當A=5時，把A當做一個標量；

當A=1+2j時，把A當做一個復數。

（1）矩陣的賦值

在MATLAB中，矩陣可以用以下幾種方式進行賦值。

●直接列出元素的形式；

●通過命令和函數產生；

●建立在文件中；

●從外部的數據文件中裝入。

①簡單矩陣的輸入

對于比較小的簡單矩陣可以使用直接排列的形式輸入，把矩陣的元素直接排列到方括號中，每行內的元素間用空格或逗號分開，行與行的內容用分號隔開。例如，矩陣

在MATLAB下的輸入方式為

      ＞＞A=[1，2，3；4，5，6；7，8，9]

或

      ＞＞A=[123；456；789]

都將得到相同的結果：

      A=
        1  2  3
        4  5  6
        7  8  9

對于比較大的矩陣，可以用回車鍵代替分號，對同一行的內容也可利用續行符號“…”，把一行的內容分兩行來輸入。例如，前面的矩陣還可以等價地由下面兩種方式來輸入：

      ＞＞A=[1  2  3；4  5  6
                  7      8  9]

或

      ＞＞A=[1  2  3；4  5…
                   6；7  8  9]

輸入后A矩陣將一直保存在工作空間中，除非被替代或清除，在MATLAB的命令窗口中可隨時輸入矩陣名查看其內容。

②利用語句或函數產生矩陣

在MATLAB中，向量可利用下面的語句來產生：

s1：s2：s3

其中，s1為起始值，s3為終止值，s2為步矩。

使用這樣的命令就可以產生一個由s1開始，以步距s2自增，并終止于s3的行向量。例如

      ＞＞y=[0：pi/4：pi；0：10/4：10]

結果顯示：

      y=
          0  0.7854  1.5708  2.3562  3.1416
          0  2.5000  5.0000  7.5000  10.0000

如果s2省略，則認為自增步距為1。例如，利用以下命令可產生一個行向量：

      ＞＞x=1：5

結果顯示：

      x=
          1  2  3  4  5

利用linspace（）函數也可產生一個行向量，該函數的調用格式為

x=linspace（n，m，z）

其中，x為產生的等間隔行向量；參數n和m分別為行向量中的起始和終止元素的值；z為需要產生行向量的元素數。例如，對于以下命令：

      ＞＞x=linspace（0，2*pi，5）

結果顯示：

      x=
        0  1.5708  3.1416  4.7124  6.2832

利用size（）函數可測取一個矩陣的維數，該函數的調用格式為

[n，m] =size（A）

其中，A為要測試的矩陣名；而返回的兩個參數n和m分別為A矩陣的行數和列數。

當要測試的變量是一個向量時，仍可由size（）函數來得出其大小；更簡潔地，用戶可以使用length（）函數來求出，該函數的調用格式為

n=length（x）

其中，x為要測試的向量名；而返回的n值為向量x的元素個數。

如果對一個矩陣A用length（A）函數測試，則返回該矩陣行列的最大值，即該函數等效于max（size（A））。

（2）矩陣的元素

MATLAB的矩陣元素可用任何表達式來描述，它既可以是實數，也可以是復數。例如

      ＞＞B=[ -1/3  1.3； sqrt（3）（1+2+3）*j]

結果顯示：

      B=
            -0.3333  1.3000
            1.7321   0+6.0000i

MATLAB允許把矩陣作為元素來建立新的矩陣。例如

      ＞＞A=[1  2  3；4  5  6；7  8  9]；C=[A；[10  11  12]]

結果顯示：

      C=
            1    2    3
            4    5    6
            7    8    9
           10    11   12

MATLAB還允許對一個矩陣的單個元素進行賦值和操作。例如，如果想將A矩陣的第2行第3列的元素賦值為100，則可通過下面的語句來完成：

      ＞＞A=[1  2  3；4  5  6；7  8  9]；A（2，3）=100

這時將只改變此元素的值，而不影響其他元素的值。

如果給出的行數或列數大于原來矩陣的范圍，則MATLAB將自動擴展原來的矩陣，并將擴展后未賦值的矩陣元素置為0。例如，想把以上矩陣A的第4行第5列元素的值定義為8，可以通過下面語句來完成：

      ＞＞A（4，5）=8

結果顯示：

      A=
        1    2     3     0               0
        4    5    100    0               0
        7    8     9     0               0
        0    0     0     0               8

利用上面的語句除了對單個矩陣元素進行定義之外，MATLAB還允許對子矩陣進行定義和處理。例如：

      A（1：3，1：2：5）                %表示取A矩陣的第1行到第3行內，且位于第1，3，5
                                         列上的所有元素構成的子矩陣；

      A（2：3，：）                     %表示取A矩陣的第2行和第3行的所有元素構成的子矩陣；
      A（：，j）                        %表示取A矩陣第j列的全部元素構成的子矩陣；
      B（：，[3，5，10]）=A（：，1：3） %表示將A矩陣的前3列，賦值給B矩陣的第3，5和第10列；
      A（：，n：-1：1）                 %表示由A矩陣中取n至1反增長的列元素組成一個新的
                                         矩陣。

特別地，當A（：）在賦值語句的右邊時，表示將A的所有元素按列在一個長的列向量中展成串。例如：

      ＞＞A=[12；34]，B=A（：）

結果顯示：

      B=
        1
        3
        2
        4

（3）特殊矩陣的實現

在MATLAB中特殊矩陣可以利用以下函數來建立。

①單位矩陣函數eye（）

      基本格式：A=eye（n）      %產生一個n階的單位矩陣A
              A=eye（size（B））%產生與B矩陣同階的單位矩陣A
              A=eye（n，m）     %產生一個主對角線的元素為1，其余元素全為0的n ×m矩陣

②零矩陣函數zeros（）

      基本格式：A=zeros（n，m）  %產生一個n ×m零矩陣A
              A=zeros（n）       %產生一個n ×n零矩陣A
              A=zeros（size（B））%產生一個與B矩陣同階的零矩陣A

③1矩陣函數ones（）

      基本格式：A=ones（n，m）
              A=ones（n）
              A=ones（size（B））

④隨機元素矩陣函數rand（）

隨機元素矩陣的各個元素是隨機產生的。如果矩陣的隨機元素滿足[0，1]區間上的均勻分布，則可以由MATLAB函數rand（）來生成。該函數的調用格式為

              A=rand（n，m）
              A=rand（n）
              A=rand（size（B））

⑤對角矩陣函數diag（）

如果用MATLAB提供的方法建立一個向量V=[a1，a2，…，an]，則可利用函數diag（V）來建立一個對角矩陣。例如：

      ＞＞V=[1，2，3，4]；A=diag（V）

如果矩陣A為一個方陣，則調用V=diag（A）將提取出A矩陣的對角元素來構成向量V，而不管矩陣的非對角元素是何值。

⑥伴隨矩陣函數compan（）

假設有一個多項式

p（s）=sn+a1 sn-1 +a2 sn-2 +…+an-1 s+an

則可寫出一個伴隨矩陣

生成伴隨矩陣的函數的調用格式為

A=compan（p）

其中，p=[1，a1，a2，…，an]為一個多項式向量。

例如有一個向量p=[1 2 3 4 5]，則可通過下面的命令構成一個伴隨矩陣。

      ＞＞p=[1  2  3  4  5]；A=compan（p）

⑦上三角矩陣函數triu（）和下三角矩陣函數tril（）

調用格式為

A=triu（B）和A=tril（B）

其中，B為一個矩陣。例如：

      ＞＞B=[1  2  3；4  5  6；7  8  9]；A=tril（B）

2.矩陣的基本運算

矩陣運算是MATLAB的基礎，MATLAB的矩陣運算功能十分強大，并且運算的形式和一般的數學表示十分相似。

（1）矩陣的轉置

矩陣轉置的運算符為“′”。例如：

      ＞＞A=[1  2  3；4  5  6]；B=A′

如果A為復數矩陣，則A′為它的復數共軛轉置，非共軛轉置使用A.′，或者用conj（A）實現。

（2）矩陣的加和減

矩陣的加減法的運算符為“+”和“-”。只有同階矩陣方可進行加減運算，標量可以和矩陣進行加減運算，但應對矩陣的每個元素施加運算。例如：

      ＞＞A=[1  2  3；4  5  6；7  8  9]；B=A+1

（3）矩陣的乘法

矩陣乘法的運算符為“*”。當兩個矩陣中前一矩陣的列數和后一矩陣的行數相同時，可以進行乘法運算，這與數學上的形式是一致的。例如：

      ＞＞A=[1  2  3；4  5  6；7  8  9]；B=[1  1；1  1；1  1]；C=A*B

在MATLAB中還可將矩陣和標量進行相乘，其結果為標量與矩陣中的每個元素分別相乘。

（4）矩陣的除法

矩陣的除法有兩種運算符“\”和“/”，分別表示左除和右除。一般地講，x=A\B是A*x=B的解，x=B/A是x*A=B的解。通常A\B≠B/A，而A\B=inv（A）*B，B/A=B*inv（A）。

（5）矩陣的乘方

矩陣乘方的運算符為“^”。一個方陣的乘方運算可以用A^P來表示。P為正整數，則A的P次冪即為A矩陣自乘P次。如果P為負整數，則可以將A自乘P次，然后對結果進行求逆運算，就可得出該乘方結果。如果P是一個分數，如P=m\n，其中n和m均為整數，則首先應該將A矩陣自乘n次，然后對結果再開m次方。例如：

      ＞＞A=[1  2  3；4  5  6；7  8  9]；B=A^2，C=A^0.1

（6）矩陣的翻轉

MATLAB還提供了一些矩陣翻轉處理的特殊命令，對n ×m維矩陣A，如

            B=fliplr（A）  %命令將矩陣A進行左右翻轉再賦給矩陣B，即bij =ai，m+1-j
            C=flipud（A）  %命令將矩陣A進行上下翻轉再賦給矩陣C，即cij =an+1-i，j
            D=rot90（A）   %命令將矩陣A進行旋轉90度后賦給矩陣D，即dij =aj，m+1-i

例如：

      ＞＞A=[1  2  3；4  5  6；7  8  9]；B=fliplr（A），C=flipud（A）

（7）矩陣的超越函數

MATLAB中exp（），sqrt（），sin（），cos（）等基本函數命令可以直接使用在矩陣上，這種運算只定義在矩陣的單個元素上，即分別對矩陣的每個元素進行運算。超越數學函數，可以在函數后加上m而成為矩陣的超越函數，例如expm（A），sqrtm（A），logm（A）分別為矩陣指數、矩陣開方和矩陣對數。矩陣的超越函數要求運算的矩陣必須為方陣。例如

      ＞＞A=[1  2  3；4  5  6；7  8  9]；B=expm（A），C=sqrtm（A）

3.矩陣的特殊運算

矩陣的特殊運算包括以下內容。

（1）矩陣行列式

矩陣A={aij}的行列式定義為

|A |=det（A）∑（-1）ka1k1a2k2…ankn

式中，k1，k2，…，kn是將序列1，2，3，…，n的元素交換k次所得出的一個序列，每個這樣的序列稱為一個置換，而∑表示對k1，k2，…，kn取遍1，2，3，…，n所有的排列的求和。

MATLAB求矩陣行列式函數的調用格式為

det（A）

計算矩陣的行列式有多種算法，在MATLAB中采用的方法為LU分解法。

（2）矩陣求逆

對于一個已知的n ×n維非奇異方陣A來說，如果有一個同樣大小的C矩陣滿足

AC =CA=I

式中，I為單位陣，則稱C矩陣為A矩陣的逆矩陣，并記為C=A-1。

矩陣求逆的算法是多種多樣的，比較常用的有行（列）主元素Gauss消去法、LU分解法和基于奇異值分解的方法等，MATLAB提供了一個求取逆矩陣的函數inv（），其調用格式為

inv（A）

如果A矩陣為奇異的或接近奇異，則利用此函數有可能產生錯誤的結果。

（3）矩陣的跡

假設一個方陣為A={aij}，（i，j=1，2，…，n），則矩陣A的跡定義為

亦即矩陣的跡為該矩陣對角線上各個元素之和。由代數理論可知矩陣的跡和該矩陣的特征值之和是相同的。在MATLAB中提供了求取矩陣跡的函數trace（），其調用格式為

trace（A）

（4）矩陣的秩

對于n ×m維的矩陣A，若矩陣所有的列向量中共有rc 個線性無關，則稱矩陣的列秩為rc；如果rc =m，則稱A為列滿秩矩陣。相應地，若矩陣A的行向量中有rr個是線性無關的，則稱矩陣A的行秩為rr；如果rr=n，則稱A為行滿秩矩陣。可以證明，矩陣的行秩和列秩是相等的，記

rank{A} =rc =rr

這時矩陣的秩為rank{A}。矩陣的秩也表示該矩陣中行列式不等于0的子式的最大階次。所謂子式，即為從原矩陣中任取k行及k列所構成的子矩陣。

矩陣求秩的算法也是多種多樣的，有的算法是穩定的，而有的算法可能因矩陣的條件數變化而不是很穩定。MATLAB中采用的算法是基于矩陣的奇異值分解的算法，首先對矩陣做奇異值分解，得出矩陣A的n個奇異值σi（i=1，2，…，n），在這n個奇異值中找出大于給定誤差限ε的個數r，這時r就可以認為是矩陣A的秩。

MATLAB提供了一個內部函數rank（）來用數值方法求取一個已知矩陣的秩，其調用格式為

k=rank（A，tol）

其中，A為要求秩的矩陣；k為所求矩陣A的秩；而tol=ε為判0用的誤差限，一般可以取默認值eps。這樣調用格式就可以簡化為

k=rank（A）

這里的判0用誤差限就取做機器中的默認值eps。如果eps取的不合適，則求出的數值秩可能和原矩陣的秩不同。所以在使用數值秩時應當引起注意。

（5）矩陣的三角分解

矩陣的三角分解又稱為LU分解，它的目的是將一個矩陣A分解成一個下三角矩陣L和一個上三角矩陣U的乘積，亦即可以寫成A=LU，其中L和U矩陣分別可以寫成

這樣產生的矩陣與原來的矩陣A的關系可以寫成

a11 =u11，a12 =u12，…，a1n =u1n

a21 =l21 u11，a22 =l21 u12 +u22，…，a2n =l21 u1n +u2n

?

其中

由上式可以立即得出求lij和uij的遞推計算公式

該公式的遞推初值為u1i =a1i，i=1，2，…，n。

注意，在上述的算法中并未對主元素進行任何選取，因此該算法并不一定數值穩定。

在MATLAB下也給出了矩陣的LU分解函數lu（），該函數的調用格式為

[L，U] =lu（A）

其中，L和U分別為變換后的下三角矩陣和上三角矩陣。在MATLAB的lu（）函數中考慮了主元素選取的問題，所以該函數一般會給出可靠的結果。由該函數得出的下三角矩陣L并不一定是一個真正的下三角矩陣。因為選取它可能進行了一些元素行的交換，這樣主對角線的元素可能不是1。例如：

      ＞＞A=[1  2  3；4  5  6；7  8  0]；[L，U] =lu（A）

LU分解使用的算法是高斯變量消去法，這種分解是矩陣求逆、矩陣求行列式和線性方程求解的基礎，也是方陣“/”或“\”兩種矩陣除法的基礎。

（6）矩陣的奇異值分解

矩陣的奇異值也可以看成是矩陣的一種測度，對任意的n ×m矩陣A來說，總有ATA≥0，AAT≥0，且有

rank{ATA} =rank{AAT} =rank（A）

進一步可以證明，ATA與AAT有相同的非零特征值λi，且相同的非零特征值總是為正數。在數學上把這些非零的特征值的平方根稱做矩陣A的奇異值，記為

矩陣的奇異值大小通常決定矩陣的性態。如果矩陣的奇異值變化特別大，則矩陣中某個參數有一個微小的變化將嚴重影響到原矩陣的參數，如其特征值的大小，這樣的矩陣又稱為病態矩陣，有時也稱為奇異矩陣。

假設矩陣A為n ×m矩陣，且rank{A} =r，則矩陣A可以分解為

其中，L和M為正交矩陣；Δ=diag{σ1，…，σr}為對角矩陣，且其對角元素均不為0。

MATLAB提供了直接求取矩陣奇異值分解的函數，其調用格式為

[U，S，V] =svd（A）

其中，A為原始矩陣；S為對角矩陣，其對角元素就是A的奇異值；而U和V均為正交矩陣，并滿足

A=USVT

矩陣最大奇異值σmax和最小奇異值σmin的比值又稱為該矩陣的條件數，記為cond{A}，即cond{A} =σmax/σmin。矩陣的條件數越大，則對參數變化越敏感。在MATLAB下也提供了求取矩陣條件數的函數cond（），其調用格式為

cond（A）

（7）矩陣的范數

矩陣的范數也是對矩陣的一種測度。在介紹矩陣的范數之前，首先要介紹向量范數的基本概念。

如果對線性空間中的一個向量x存在一個函數ρ（x）滿足下面三個條件：

① ρ（x）≤0且ρ（x）=0的充要條件是x=0；

② ρ（ax）=|a|ρ（x），a為任意標量；

③對向量x和y有ρ（x+y）≤ρ（x）+ρ（y）。

則稱ρ（x）為x向量范數。

范數的形式是多種多樣的。可以證明，下面給出的一族式子都滿足上述的三個條件：

且

這里用到了向量范數的記號‖x‖p。

對于任意的非零向量x，矩陣A的范數為

和向量的范數一樣，矩陣的范數定義如下：

其中，s{A}為矩陣A的特征值，而smax{ATA}即為矩陣A的最大奇異值的平方。換句話說，‖A‖2 為矩陣A的最大奇異值。

MATLAB提供了求取矩陣范數的函數norm（），它允許求各種意義下的矩陣范數，該函數的調用格式為

N=norm（A，選項）

其中，選項如表1-8所示。

（8）矩陣的特征值與特征向量

對一個矩陣A來說，如果存在一個非零的向量x，且有一個標量λ滿足

Ax =λx

則稱λ為矩陣A的一個特征值，而x為對應于特征值λ的特征向量。嚴格說來，x應該稱為A的右特征向量。如果矩陣A的特征值不包含重復的值，則對應的各個

特征向量為線性獨立的，這樣由各個特征向量可以構成一個非奇異的矩陣，如果用它對原始矩陣做相似變換，則可以得出一個對角矩陣。

矩陣的特征值與特征向量可以由MATLAB提供的函數eig（）容易地求出，該函數的調用格式為

[V，D] =eig（A）

其中，A為要處理的矩陣；D為一個對角矩陣，其對角線上的元素為矩陣A的特征值，而每個特征值對應的V矩陣的列為該特征值的特征向量。該矩陣是一個滿秩矩陣，它滿足AV=VD，且每個特征向量各元素的平方和（即2范數）均為1。如果調用該函數時只給出一個返回變量，則將只返回矩陣A的特征值。即使A為復數矩陣，也同樣可以由eig（）函數得出其特征值與特征向量矩陣。

（9）矩陣的特征多項式、特征方程和特征根

對于給定的n ×n階矩陣A，稱多項式

f（s）=det（sI-A）=a0 sn+a1 sn-1 +…+an-1 s+an

為矩陣A的特征多項式，其中系數a0，a1，…，an稱為矩陣的特征多項式系數。

MATLAB提供了求取矩陣特征多項式系數的函數poly（），其調用格式為

p=poly（A）

其中，A為給定的矩陣；返回值p為一個行向量，其各個分量為矩陣A的降冪排列的特征多項式系數。即p=[a0 a1 … an]。

令特征多項式等于零所構成的方程稱為該矩陣的特征方程，而特征方程的根稱為該矩陣的特征根。MATLAB中根據矩陣特征多項式求特征根的函數為roots（），其調用格式為

V=roots（p）

其中，p為特征多項式的系數向量，而V為特征多項式的解，即原始矩陣的特征根。

【例1-11】 求的特征多項式、特征根及其模值。

解 MATLAB命令如下

      ＞＞A=[1  2  3；4  5  6；7  8  9]；p=poly（A），r=roots（p）′，abs（r）

結果顯示：

      p=
          1.0000  -15.0000     -18.0000   -0.0000
      r =
          16.1168    -1.1168   -0.0000
      ans =
          16.1168    1.1168    0.0000

1.5.2 向量運算

雖然在MATLAB中向量和矩陣在形式上有很多的一致性，但它們實際上遵循著不同的運算規則。MATLAB向量運算符由矩陣運算符前面加一點“.”來表示，如“.*”、“./”和“.^”等。

1.向量的加減

向量的加減運算與矩陣的運算相同，所以“+”和“-”既可被向量接受又可被矩陣接受。

2.向量的乘法

向量乘法的操作符為“.*”。如果x和y兩向量具有相同的維數，則x.*y表示向量x和y單個對應元素之間的相乘。例如：

      ＞＞x=[1  2  3]；y=[4  5  6]；z=x.*y

結果顯示：

      z=
        4  10  18

可見向量的輸入和輸出與矩陣具有相同的格式，但它們的運算規則不同。例如，如果x是一個向量，則求取向量x平方時不能直接寫成x*x，而必須寫成x.*x，否則將給出錯誤信息。

但是對于矩陣可以使用向量運算符號，這時實際上就相當于把矩陣看成了向量來進行運算。例如，對于兩個維數相同矩陣A和B，C=A.*B表示A和B矩陣的對應元素之間直接進行乘法運算，然后將結果賦給C矩陣，把這種運算稱為矩陣的點積運算。兩個矩陣之間的點積是它們對應元素的直接運算，它與矩陣的乘法是不同的。例如：

      ＞＞A=[1  2  3；4  5  6；7  8  9]；B=[2  3  4；5  6  7；8  9  0]；C=A.*B

結果顯示：

      C=
            2   6  12
        20  30     42
        56  72   0

3.向量的除法

向量除法的操作符為“./”或“.\”。它們的運算結果一樣。例如：

      ＞＞x=[123]；y=[456]；z=y./x

結果顯示：

      z=
        4.0000  2.5000  2.0000

對于向量除法運算x.\y和y./x一樣，將得到相同的結果，這與矩陣的左、右除是不一樣的，因向量的運算是它們對應元素間的運算。

對于矩陣也可使用向量的除法操作符，這時就相當于把矩陣看成向量進行運算。

4.向量的乘方

向量乘方的運算符為“.^”。向量的乘方是對應元素的乘方，在這種底與指數均為向量的情況下，要求它們的維數必須相同。例如：

      ＞＞x=[1  2  3]；y=[4  5  6]；z=x.^y

結果顯示：

      z=
        1  32  729

它相當于：z=[1 2 3].^[4 5 6] =[14 25 36]。

若指數為標量時，會得到如下結果。例如：

      ＞＞x=[1  2  3]；z=x.^2

結果顯示：

      z=
        1  4  9

以上運算相當于：z=[1 2 3].^2=[12 22 32]。

若底為標量時，則會得到如下結果。例如：

      ＞＞x=[1  2  3]；y=[4  5  6]；z=2.^[x  y]

結果顯示：

      z=
        2  4  8  16  32  64

以上運算相當于：z=2.^[1 2 3 4 5 6] =[21 22 23 24 25 26]。

同樣對于矩陣也可以采用運算符“.^”。例如：

      ＞＞A=[1  2  3；4  5  6；7  8  0]；B=A.^A

結果顯示：

      B=
              1          4      27
            256       3125    46656
        823543    16777216       1

即矩陣B中的每個元素都是矩陣A元素的相應乘方。例如，3125=5^5。

可見如果對矩陣使用向量運算符號，實際上就相當于把矩陣看成了向量的運算，否則將作為矩陣運算。

1.5.3 關系和邏輯運算

1.關系運算

MATLAB常用的關系操作符如表1-9所示。

MATLAB的關系操作符可以用來比較兩個大小相同的矩陣，或者比較一個矩陣和一個標量。比較兩個元素大小時，結果是1表明為真，結果是0表明為假。函數find（）在關系運算中很有用，它可以在矩陣中找出一些滿足一定關系的數據元素。例如：

      ＞＞A=1：9；B=A＞4

結果顯示：

      B=
        0  0  0  0  1  1  1  1  1
      ＞＞A=1：9；C=A（A＞4）

或

      ＞＞A=1：9；C=find（A＞4）

結果顯示：

      C=
        5  6  7  8  9

2.邏輯運算

MATLAB的邏輯操作符有 &（與）、|（或）和～（非）。它們通常用于元素或0-1矩陣的邏輯運算。

與運算符和或運算符可比較兩個標量或兩個同階矩陣。對于矩陣，邏輯運算符是作用于矩陣中的元素。邏輯運算結果信息也用“0”和“1”表示，邏輯操作符認定任何非零元素都表示為真。給出1為真，0為假。

非是一個元操作符，當A非零時，～A返回的信息為0；當A為零時，～A返回的信息為1。因而就有：P|（～P）返回值為1，P&（～P）返回值為0。例如：

      ＞＞A=1：9；C=～（A＞4）

結果顯示：

      C=
        1  1  1  1  0  0  0  0  0
      ＞＞A=1：9；C=（A＞4）&（A＜7）

結果顯示：

      C=
        0  0  0  0  1  1  0  0  0

3.關系和邏輯運算函數

除了上面介紹的關系和邏輯運算符外，MATLAB中還提供了一些關系和邏輯運算函數，如表1-10所示。

對于矩陣，any（）和all（）命令按列對其進行處理，并返回帶有處理列所得結果的一個行向量。

1.5.4 多項式運算

多項式運算是數學中最基本的運算之一。在MATLAB中同樣可對多項式進行相應的一系列運算。

1.多項式的表示

在高等數學中，多項式一般可表示成以下形式：

f（x）=a0 xn+a1 xn-1 +…+an-1 x+an

其中，a0，a1，…，an稱為多項式的系數。

所以多項式很容易用其系數組成的行向量p=[a0 a1 … an]來表示，其中行向量是按其系數降冪排列組成的系數向量。在MATLAB中，構造多項式正是采用了把多項式的各項系數依降冪次序排放在行向量的對應元素位置，直接輸入其系數向量的方法來實現的。對于缺項的系數一定要進行補零。例如，對于多項式

f（x）=x4 +5x3 +3x+2

可用以下MATLAB命令來表示：

      ＞＞p=[1  5  0  3  2]

在MATLAB中，利用函數poly2str（）可將多項式的系數向量表示成相應多項式的習慣表示形式，該函數的調用格式為

f=poly2str（p，′s′）

其中，p為多項式的系數向量，s為多項式的變量名，f為相應的多項式。例如：

      ＞＞p=[1  5  0  3  2]；f=poly2str（p，′x′）

結果顯示：

      f=
        x^4+5 x^3+3x+2

2.多項式的四則運算

多項式的四則運算是指多項式的加、減、乘和除運算。其中多項式的加、減運算要求兩個相加、減的多項式的系數向量維數的大小必須相等。

（1）多項式的加減

在MATLAB中，當兩個相加、減的多項式階次不同時，低階多項式的系數向量必須用首零填補，使其與高階多項式的系數向量有相同維數。

【例1-12】 求以下兩個多項式的和。

f1（x）=x4 +5x3 +3x+2，f2（x）=x2 +6x+5

解 MATLAB命令如下：

      ＞＞p1=[1  5  0  3  2]；p2=[0  0  1  6  5]；p=p1+p2

結果顯示：

      p=
        1  5  1  9  7

（2）多項式的乘法

在MATLAB中，多項式的乘法運算，利用函數conv（）來實現。函數conv（）相當于執行兩個數組的卷積，其調用格式為

p=conv（p1，p2）

其中，p1和p2為多項式的系數按降冪排列構成的系數向量；p為多項式p1和p2的乘積多項式，按其系數降冪排列構成的多項式積的系數向量。

【例1-13】 求以下兩個多項式的乘積。

f1（x）=x4 +5x3 +3x+2，f2（x）=x2 +6x+5

解 MATLAB命令如下：

      ＞＞p1=[1  5  0  3  2]；p2=[1  6  5]；p=conv（p1，p2）

結果顯示：

      p=
        1  11  35  28  20  27  10

需要說明的是，當對多個多項式執行乘法運算時，可重復使用conv（）函數。

【例1-14】 求多項式f（x）=（x+1）2（x2 +6x+5）的展開式。

解 MATLAB命令如下：

      ＞＞p=conv（[1  1]，conv（[1  1]，[1  6  5]））

結果顯示：

      p =
          1  8  18  16  5

（3）多項式的除法

在MATLAB中，多項式的除法運算是利用函數deconv（）來實現的，其調用格式為

[p，r] =deconv（p1，p2）

其中，p1，p2為多項式的系數按降冪排列構成的系數向量；p為多項式p1被p2除的商多項式，按其系數降冪排列構成的多項式商的系數向量，而余多項式為r。函數deconv（）相當于兩個數組的解卷運算，使p1=conv（p2，p）+r成立。

3.多項式的值及其導數

如果多項式函數為

f（x）=a0 xn+a1 x（n-1）+…+an-1 x+an

則可以求出該函數的導數函數為

f′（x）=na0 xn-1 +（n-1）a1 xn-2 +…+an-1

在MATLAB中提供了多項式求值函數polyval（）和多項式求導函數polyder（），它們的調用格式分別為

f0=polyval（p，x0）和dp=polyder（p）

其中，p為由多項式系數按降冪排列構成的系數向量；x0為求值點的x值。該函數將返回f（x）多項式在x=x0的值f0。函數polyder（p）返回多項式導數的系數向量，亦即向量

dp=[na0（n-1）a1 … an-1]

同樣，MATLAB也提供了多項式矩陣的求值函數polyvalm（），其調用格式為

fA=polyvalm（p，A）

其中，p為矩陣多項式函數按降冪排列構成的向量，即

p=[a0 a1 … an]

而A為一個給定矩陣，返回值fA為下面的矩陣多項式的值：

f（A）=a0An+a1An-1 +…+an-1A+anI

4.多項式的求解

MATLAB中多項式的求解運算同樣可利用前面介紹的roots（）函數來實現，其調用格式為

r=roots（p）

其中，p為多項式的系數向量；r為多項式的解。

【例1-15】 求方程f（x）=x2 +5x+6=0的解。

解 MATLAB命令如下：

      ＞＞p=[1  5  6]；x=roots（p）

結果顯示：

      x=
        -3.0000

        -2.0000

5.多項式的擬合

在MATLAB中，多項式的擬合可用polyfit（）函數來實現，其調用格式為

[p，s] =polyfit（x，y，n）

其中，x，y為利用最小二乘法進行擬合的數據；n為要擬合的多項式的階次；p為要擬合的多項式的系數向量；s為使用該函數獲得的錯誤預估計值。一般來說，多項式擬合的階數越高，擬合的精度就越高。

【例1-16】 用6階多項式對[0，2*pi]區間上的sint函數進行擬合。

解 MATLAB命令如下：

      ＞＞x=linspace（0，2*pi，80）；y=sin（x）；p=polyfit（x，y，6），y1=polyval（p，x）；
      ＞＞plot（x，y，′ro′，x，y1）

擬合結果如圖1-10所示。

6.多項式的插值

在MATLAB中，多項式的插值方法包括一維線性插值、二維線性插值、三維線性插值和三次樣條插值等。

（1）一維線性插值

對于一維線性數據，可以通過插值或查表來求得離散點之間的數據值。在MATLAB中，一維線性插值可用函數interp1（）來實現，其調用格式為

yi=interp1（x，y，xi，′method′）

其中，返回值yi為在插值向量xi處的函數值向量，它是根據向量x與y插值而來的。如果y是一個矩陣，那么對y的每一列進行插值，返回的矩陣yi的大小為length（xi）×length（y，2）。如果xi中有元素不在x的范圍內，則與之相對應的yi返回NaN。如果x省略，表示x=1：n，此處n為向量y的長度或為矩陣y的行數，即size（y，1）。參數′method′表示插值方法，它可采用以下方法：最近插值（′nearest′）、線性插值（′linear′）、三次插值（′cubic′）和三次樣條插值（′spline′）。′method′的默認值為線性插值。

需要注意，上述函數interp1（）的所有插值調用格式，都要求向量x為單調。當x為單調且等間距時，可以使用快速插值法，此時可將′method′參數的值設置為′nearest′，′linear′或′cubic′。例如：

      ＞＞x=linspace（0，2*pi，80）；y=sin（x）；
      ＞＞x1=[0.5 1.4 2.6 4.2]，y1=interp1（x，y，x1，′cubic′）

結果顯示：

      x1=
          0.5000    1.4000    2.6000    4.2000
      y1=
          0.4794    0.9855    0.5155    -0.8716

（2）二維線性插值和三維線性插值

與一維線性插值一樣，二維線性插值也可以通過插值或查表來求得離散點之間的數據值。在MATLAB中，二維線性插值可用函數interp2（）來實現，其調用格式為

zi=interp2（x，y，z，xi，yi，′method′）

其中，返回值zi為在插值向量xi，yi處的函數值向量，它是根據向量x，y與z插值而來的。x，y和z也可以是矩陣。如果xi，yi中有元素不在x，y的范圍內，則與之相對應的zi返回NaN。如果x，y省略，表示x=1：m，y=1：n，此處n與m為矩陣z的行數和列數，即[n，m] =size（z）。參數′method′的定義同上。

需要注意，上述的所有插值調用格式，都要求向量x與y為單調且具有相同的格式。當x與y為單調且等間距時，同樣可以使用快速插值法。例如：

      ＞＞z=[3  5  7  9  11  10  4  15；1  3  2  6  11  5  7  13]；
      ＞＞z1=interp2（z，[1，3]，[1，2]）

結果顯示：

      z1=
          3  2

與一維線性插值和二維線性插值一樣，三維線性插值也可以通過插值或查表來求得離散點之間的數據值。在MATLAB中，三維線性插值可用函數interp3（）來實現，其使用方法與in-terp1（）函數和interp2（）函數基本類似，這里就不詳細介紹了。

（3）三次樣條插值

使用高階多項式的插值常常會產生病變的結果，而三次樣條插值就能消除這種病變。在三次樣條插值中，要尋找三次多項式，以逼近每對數據點之間的曲線。在MATLAB中，三次樣條插值利用spline（）函數來實現，其調用格式為

yi=spline（x，y，xi）

其中，返回值yi為在插值向量xi處的函數值向量，它是根據向量x與y插值而來的。此函數的作用等同于interp1（x，y，xi，′spline′）。

1.5.5 數據分析

MATLAB強大的數組運算功能，決定了它很容易對一大批數據進行一般的數據分析，如求數組的極值、平均值、中值、和、積、標準差、方差、協方差和排序等。

1.隨機數

在MATLAB中提供了兩個用來產生隨機數組的函數rand（）和randn（）。其調用格式為

A=rand（n，m）和A=randn（n，m）

其中，n和m分別為將要產生隨機數組的行和列；A為n ×m維的隨機數組。函數rand（n，m）用來產生一個n ×m維在[0，1]區間上均勻分布的隨機數組；函數randn（n，m）用來產生一個n ×m維的均值為0、標準差為1的正態分布的隨機數組。

2.最大值和最小值

如果給定一組數據{xi}，i=1，2，…，n，則可利用MATLAB將這些數據用一個向量表示出來，即

x=[x1，x2，…，xn]

利用MATLAB中的函數max（）和min（）便可求出這組數據的最大值和最小值，調用格式為

[xmax，i] =max（x）和 [xmin，i] =min（x）

其中，返回的xmax及xmin分別為向量x的最大值和最小值，而i為最大值或最小值所在的位置。當然這兩個函數均可以只返回一個參數，而不返回i。如果給出的x不是向量而是矩陣，則采用min（）和max（）函數得出的結果將不是數值，而是一個向量。它的含義是得出由每一列構成的向量的最大值或最小值所構成的行向量。當x為復數時，通過計算max（abs（x））返回結果。

對于函數max（）和min（）也可以采用以下格式：

z=max（x，y）和z=min（x，y）

其中，x與y為向量或矩陣，它們的大小一樣；返回結果z是一個包含最大或最小元素，且與它們大小一樣的向量或矩陣。

3.平均值

利用MATLAB中的mean（）函數可求出向量或矩陣的平均值，調用格式為

y=mean（x）

其中，x為向量或矩陣；y為向量或矩陣x的平均值。若x為向量，則返回向量元素的平均值。若x為矩陣，則返回一行向量，包括矩陣每列元素的平均值。

4.中值

利用MATLAB中的median（）函數可求出向量或矩陣的中值，調用格式為

y=median（x）

其中，x為向量或矩陣；y為向量或矩陣x的中值。若x為向量，則返回向量元素的中值；若x為矩陣，則返回一行向量，包括矩陣每列元素的中值。

5.求和

利用MATLAB中的sum（）函數可求出向量或矩陣中元素的和，調用格式為

y=sum（x）

其中，x為向量或矩陣；y為向量或矩陣x中元素的和。若x為向量，則返回向量元素的總和；若x為矩陣，則返回一行向量，包括矩陣每列元素的和。

6.求積

利用MATLAB中的prod（）函數可求出向量或矩陣中元素的積，調用格式為

y=prod（x）

其中，x為向量或矩陣；y為向量或矩陣x中元素的積。若x為向量，則返回向量元素的積；若x為矩陣，則返回一行向量，包括矩陣每列元素的積。

7.標準差

利用MATLAB中的std（）函數可求出向量或矩陣中元素的標準差，調用格式為

y=std（x）

其中，x為向量或矩陣；y為向量或矩陣x中元素的標準差。若x為向量，則返回向量元素的標準差；若x為矩陣，則返回一行向量，包括矩陣每列元素的標準差。

8.方差

利用MATLAB中的var（）函數可求出向量或矩陣中元素的方差，調用格式為

y=var（x）

其中，x為向量或矩陣；y為向量或矩陣x中元素的方差。若x為向量，則返回向量元素的方差；若x為矩陣，則返回一行向量，包括矩陣每列元素的方差。

9.協方差

利用MATLAB中的cov（）函數可求出向量或矩陣中元素的協方差，調用格式為

y=cov（x）

其中，x為矩陣；y為矩陣x中各列間的協方差陣。函數cov（x，y）相當于

cov（x（：），y（：））

10.按實部或幅值對特征值進行排序

利用MATLAB中的esort（）和dsort（）函數，可對特征值按實部或幅值進行排序，函數的調用格式為

[s，ndx] =esort（p）和 [s，ndx] =dsort（p）

其中，esort（p）針對連續系統，根據實部按遞減順序對向量p中的復特征值進行排序；ndx為索引矢量，對于連續特征值，先列出不穩定特征值。dsort（p）針對離散系統，根據幅值按遞減順序對向量p中的復特征值進行排序；ndx為索引矢量。對于離散特征值，先列出不穩定特征值。

另外，sort（）函數用于對元素按升序進行排序。

1.5.6 函數極值

在MATLAB中，提供了兩個基于單純形算法求解多元函數極小值的函數fmins（）（僅適用于MATLAB 6.5及以前的版本）和fminsearch（），以及一個基于擬牛頓法求解多元函數極小值的函數fminunc（），其調用格式分別為

                          x=fmins（fun，x0，options）
                          x=fminsearch（fun，x0，options）
                          x=fminunc（fun，x0，options）

其中，fun表示函數名；x0表示函數的初值，其大小往往能決定最后解的精度和收斂速度；x為所求極小值點；options表示選項，它是由一些控制變量構成的向量，比如它的第一個分量不為0表示在求解時顯示整個動態過程（其默認值為0），第二個分量表示求解的精度（默認值為1e-4）。可以指定這些參數來控制求解的條件。在調用以上函數時，首先應該寫一個描述f（·）的函數，其格式為

y=fun（x）

【例1-17】 求函數f（x）=3x-x3的最小值。

解 首先根據方程編寫一個函數ex1 -17.m。

            %ex1  17.m
            function y=ex1  17（x）
            y=3*x-x^3；

然后利用下面的命令調用fminsearch（）函數求方程的解。

      ＞＞x=fminsearch（′ex1  17′，0），ymin=ex1  17（x）

運行結果：

      x=
        -1.0000
      ymin=
              -2

結果表明，函數在x=-1時有最小值ymin=-2。

以上函數的極值問題，也可直接利用以下命令得到與上面相同的結果。

＞＞f=′3*x-x^3′，x=fminsearch（f，0），ymin=3*x-x^3

因為函數f（x）的極小值問題等價于函數-f（x）的極大值問題，所以函數fminsearch（）也可用來求解函數 f（x）的極大值，這時只需在所求函數前加一負號即可。例如，求函數f（x）=3x-x3 的最大值，可利用以下命令：

      ＞＞f=′-（3*x-x^3）′，x=fminsearch（f，0），ymax=3*x-x^3

運行結果：

      x =
            1.0000
      ymax=
            2

結果表明，函數在x=1時有最大值ymax=2。

1.5.7 代數方程求解

利用MATLAB中求函數f（.）零點的函數fzero（）和fsolve（），可以方便地求得非線性方程f（·）=0的解，它們的調用格式分別為

x=fzero（fun，x0）和x=fsolve（fun，x0）

其中，fun表示函數名，函數名定義同前；x0表示函數的初值，是為求解過程所假設的起始點，當x0為標量時，該命令將在它兩側尋找一個與之最靠近的解，當x0為二元向量[a，b]時，該命令將在區間[a，b]內尋找一個解；x為返回的解。fzero（）用來對一元方程求解，fsolve（）用來對多元方程求解。當然利用函數fzero（）和fsolve（）也可對線性方程進行求解。

【例1-18】 試求函數f=sinx的零點。

解 MATLAB命令如下：

      ＞＞y=fzero（′sin（x）′，[0，pi]）

結果顯示：

      y=
        0

【例1-19】 試求以下非線性方程組的解。

解 首先根據三元方程編寫一個函數ex1-19.m。

            %ex1  19.m
            function q=ex1  19（p）
            q（1）=sin（p（1））+p（2）^2+log（p（3））-7；
            q（2）=3*p（1）+2*p（2）-p（3）^3+1；
            q（3）=p（1）+p（2）+p（3）-5；

然后利用下面的命令在初值x0 =1，y0 =1，z0 =1下調用fsolve（）函數直接求出方程的解。

      ＞＞x=fsolve（′ex1  19′，[1  1  1]）

結果顯示：

      x=
        0.5991  2.3959  2.0050

【例1-20】 試求以下線性方程組的解。

解 首先根據兩元方程編寫一個函數ex1-20.m。

            %ex1  20.m
            function z=ex1  20（p）
            z（1）=p（1）-p（2）；
            z（2）=p（1）+p（2）-6；

然后利用下面的命令在任意給的初值x0 =0，y0 =0下調用fsolve（）函數直接求出方程的解。

      ＞＞x=fsolve（′ex1  20′，[0 0]）

結果顯示：

      x=
        3.0000  3.0000

特別指出，對線性方程Ax=B，在A的逆存在的條件下，有更簡單的求解方法。在MAT-LAB下對該線性方程，可利用以下命令進行求解。

x=inv（A）*B

例如，對于例1-20也可利用以下命令，得到上述結果。

      ＞＞A=[1  -1；1  1]；B=[0；6]；x=inv（A）*B

1.5.8 微分方程求解

MATLAB中提供了求解常微分方程的函數ode45（），其調用格式為

x=ode45（fun，[t0，tf]，x0，tol）

其中，fun為函數名，其定義同前；[t0，tf]為求解時間區間；x0為微分方程的初值；tol用來指定誤差精度，其默認值為10 -3；x為返回的解。

【例1-21】 求下列微分方程在初始條件x（0）=1，下的解。

解 首先將微分方程寫成一階微分方程組。

令x1 =x，，則可得

然后根據以上微分方程組編寫一個函數ex1-21.m。

            %ex1  21.m
            function dx=ex1  21（t，x）
            dx=[x（2）；（1-x（1）^2）*x（2）-x（1）]；

最后利用以下的MATLAB命令，即可求出微分方程在時間區間[0，30]上的解曲線，如圖1-11所示。

      ＞＞[t，x] =ode45（′ex1  21′，[0，30]，[1；0]）；
      ＞＞plot（t，x（：，1），t，x（：，2））；xlabel（′t′）；ylabel（′x（t）′）

1.5.9 函數積分

對于函數f（x）的定積分

在被積函數f（x）相當復雜時，一般很難采用解析的方法求出定積分的值來，而往往要采用數值方法來求解。求解定積分的數值方法是多種多樣的，如簡單的梯形法、Simpson法和Rom-berg法等，它們的基本思想都是將整個積分空間[a，b]分割成若干個子空間[xi，xi+1]，i=1，2，…，n。其中x1 =a，xn+1 =b。這樣整個積分問題就分解為下面的求和形式：

而在每一個小的子空間上都可以近似地求解出來。例如可以采用下面給出的Simpson方法來求解出[xi，xi+1]上積分近似值：

式中，hi =xi+1 -xi。

MATLAB基于此算法，采用自適應變步長方法用quad（）函數來求取定積分，該函數的調用格式為

y=quad（fun，a，b，tol）

其中，fun為函數名，其定義和其他函數一致；a，b分別為定積分的上下限；tol為變步長用的誤差限，如不給出誤差限，則將自動地假定tol=1 ×10 -3；y為返回的值。

【例1-22】 試求積分的解。

解 這一無窮定積分可以由有限積分來近似，一般情況下，選擇積分的上下限為±15就能保證相當的精度。

首先根據給定的被積函數編寫下面的函數ex1-22.m。

            %ex1  22.m
            function f=ex1  22（x）
            f=1/sqrt（2*pi）*exp（-x.^2/2）；

然后通過下面的MATLAB語句可求出所需函數的定積分。

      ＞＞format long；y=quad（′ex1  22′，-15，15）

結果顯示：

      y=
        1.00000007247356

以上函數的積分問題，也可直接利用以下命令得到與上面相同的結果。

      ＞＞format long；y=quad（′1/sqrt（2*pi）*exp（-x.^2/2）′，-15，15）

或

      ＞＞format long；f=′1/sqrt（2*pi）*exp（-x.^2/2）′；y=quad（f，-15，15）

另外，MATLAB給出一種利用插值運算來更精確更快速地求出所需要的定積分的函數quadl（），該函數的調用格式為

y=quadl（fun，a，b，tol）

其中，tol的默認值為10 -6，其他參數的定義及使用方法和quad（）幾乎一致。該函數可以更準確地求出積分的值，且一般情況下函數調用的步數明顯小于quad（）。

對于上例，使用quadl（）函數可得如下結果：

      ＞＞format long；f=′1/sqrt（2*pi）*exp（-x.^2/2）′；y=quadl（f，-15，15）

結果顯示：

      y=
        1.00000000000338

在二維情況下求積分實質上是求函數曲線與坐標軸之間所夾的封閉圖形的面積，故利用MATLAB中的trapz（）函數，也可求取定積分。例如，對上例有

      ＞＞format long；x=-15：0.01：15；y=ex1  21（x）；area=trapz（x，y）

結果顯示：

      area=
          1.00000000000000