《MATLAB7 基礎教程》課件第4章

4.1矩陣與線性代數(shù)

4.2多項式與插值

4.3快速傅里葉變換

4.4函數(shù)的函數(shù)

4.5求解微分方程

4.6稀疏矩陣

4.1.1矩陣分析

MATLAB提供的矩陣分析函數(shù)如表4-1所示。4.1矩陣與線性代數(shù)表4-1矩陣分析函數(shù)

1．向量和矩陣的范數(shù)

對于一個n維向量，可以用一個數(shù)?||x||?來度量該n維向量的大小。如果?||x||?滿足以下三個性質(zhì)，則稱?||x||?為向量x的范數(shù)。

●非負性：對于一切x都有

||x||?>?0，且?||x||?=?0的充要條件是x=0(4-1)

●正齊性：對任何實數(shù)α和向量x，有

||αx||?=?|α|?||x||(4-2)

●三角不等式：對任何向量x和y，有

||x+y||≤||x||?+?||y||?4-3)顯然，滿足以上三個條件的函數(shù)有很多種，如果定義：

(4-4)

(4-5)

(4-6)

由于式(4-4)、(4-5)和(4-6)都滿足范數(shù)的三個性質(zhì)，因此它們都是n維向量x的范數(shù)。將、和分別稱為向量x的1-范數(shù)、2-范數(shù)和∞-范數(shù)，并分別用記號?||x||1、||x||2和||x||∞來表示。另外，記號?||x||?泛指向量的范數(shù)。對于一個矩陣A，相應的矩陣范數(shù)?||A||?應當對一切n階矩陣A和一切n維向量x滿足

|(4-7)

滿足式(4-7)的矩陣范數(shù)稱為與向量范數(shù)相容的范數(shù)。如果以?||A||1、||A||2和?||A||∞分別表示與向量的1-范數(shù)、2-范數(shù)和∞-范數(shù)相容的矩陣范數(shù)，那么對于n階矩陣A=(aij)，有(4-8)(4-9)(4-10)

MATLAB中，常用norm函數(shù)計算向量或矩陣的1-范數(shù)、2-范數(shù)和?∞-范數(shù)，其調(diào)用格式為norm(x,p)或norm(A,p)，作用是計算向量x或矩陣A的p-范數(shù)。其中，p為1、2或inf。命令norm(x)或norm(A)相當于命令norm(x,2)或norm(A,2)，而normest(A)實現(xiàn)矩陣A的

2-范數(shù)的快速估算。

【例】求向量的1-范數(shù)、2-范數(shù)和∞-范數(shù)。

在命令窗輸入：

>>x=[123];

>>[norm(x,1),norm(x,2),norm(x,inf)]運行結果：

ans=

6.00003.74173.0000

【例】求矩陣的1-范數(shù)、2-范數(shù)和∞-范數(shù)。

在命令窗輸入：

>>A=[123

456

789];

>>[norm(A,1),norm(A,2),norm(A,inf)]

運行結果：

ans=

18.000016.848124.0000

2．矩陣的秩

對于一個矩陣A，如果A=0，則A的秩為零；如果A≠0，則稱A中非零子式的最高階數(shù)為A的秩。MATLAB中常用命令rank(A)來以默認允許誤差計算矩陣A的秩，而用命令rank(A,tol)來以給定容許誤差tol計算矩陣的秩。

【例】求矩陣的秩。

在命令窗輸入：

>>A=[123

456

789];

>>rank(A)運行結果：

ans=

2

3．矩陣的行列式

當矩陣A為方陣時，可以利用命令det(A)來計算A的行列式。MATLAB中行列式的計算是通過高斯消去法先實現(xiàn)矩陣的LU分解，然后再計算下三角矩陣L和上三角矩陣U的行列式之積(對角線元素之積)，即A的行列式。

【例】求矩陣的行列式。

在命令窗輸入：

>>A=[132

654

789];

>>det(A)

運行結果：

ans=

-39

4．矩陣的跡

矩陣的跡定義為主對角線的元素之和。無論矩陣是否為方陣，MATLAB中都可以用命令trace(A)來計算矩陣A的跡。

【例】求矩陣的跡。

在命令窗輸入：

>>A=[132

654];

>>trace(A)

運行結果：

ans=

6

5．化零矩陣

矩陣A的化零矩陣Z滿足A*Z的元素近似為零，并且Z‘*Z=I。MATLAB中可以用命令Z=null(A)求矩陣A的化零矩陣，用命令Z=null(A,’r‘)求矩陣A的有理形式的化零矩陣，有理形式的化零矩陣Z’*Z≠I。

【例】求矩陣的化零矩陣。

在命令窗輸入：

>>A=[123

321

123];

>>Z=null(A);

>>A*Z運行結果：

ans=

1.0e-015*

0

-0.6661

0

【例】求矩陣的有理形式的化零矩陣。

在命令窗輸入：

>>A=[123

321

123];

>>z=null(A,‘r’);

>>A*z

運行結果：

ans=

0

0

0

6．正交化

矩陣A的標準正交基B的列向量構成的線性空間與矩陣A的列向量構成的線性空間相同，并且標準正交基B滿足B‘*B=eye(rank(A))。MATLAB提供了命令orth(A)來計算矩陣A的標準正交基。

【例】求矩陣的標準正交基。

在命令窗輸入：

>>A=[123

321

123];

>>B=orth(A)運行結果：

B=

-0.60080.3730

-0.5274-0.8496

-0.60080.3730

繼續(xù)在命令窗輸入：

>>B‘*B

運行結果：

ans=

1.00000.0000

0.00001.0000

7．矩陣的約化行階梯形

MATLAB提供了通過使用不完全主元高斯-約當消去法計算矩陣A的約化行階梯形R的命令R=rref(A)。其默認允許誤差為max(size(A))*eps*norm(A,inf)。

【例】求矩陣的約化行階梯形。

在命令窗輸入：

>>A=[123

321

123];

>>R=rref(A)運行結果：

R=

10-1

012

000

8．兩個子空間的夾角

MATLAB提供了計算矩陣A和B的列子空間夾角的命令subspace(A,B)。如果子空間夾角很小，則表明兩個子空間幾乎是線性相關的。

【例】求矩陣的子空間夾角。在命令窗輸入：

>>A=[123

321

123];

>>B=[125

322

124];

>>subspace(A,B)

運行結果：

ans=

1.57084.1.2求解線性方程組

通常遇到的線性方程組都具有如Ax=b的形式。本小節(jié)將對恰定系統(tǒng)、超定系統(tǒng)、欠定系統(tǒng)的情況分別進行討論。

線性代數(shù)中，求解形如Ax=b的線性方程組需要求出所有可能存在的解，即通解。求線性方程組通解分為三個步驟：

(1)利用命令null(A)求相應齊次系統(tǒng)Ax=0的通解，得到方程組Ax=0的基礎解向量，基礎解向量的任意線性組合均為方程組Ax=b的解。

(2)尋找非齊次系統(tǒng)Ax=b的特解。

(3)將步驟(1)與步驟(2)的結果相加，即為非齊次系統(tǒng)Ax=b的通解。

下面將分別討論恰定系統(tǒng)、超定系統(tǒng)、欠定系統(tǒng)特解的求法。

1．恰定系統(tǒng)

恰定系統(tǒng)Ax=b的矩陣A為方陣，但是A可能是奇異矩陣或非奇異矩陣。當A為非奇異矩陣時，可直接利用命令x=A\b或x=inv(A)*b來求系統(tǒng)的特解。

【例】求非奇異矩陣恰定系統(tǒng)的特解。在命令窗輸入：

>>A=[124

321

123];

>>b=[1;2;3];

>>x=A\b

運行結果：

x=

-2.5000

5.7500

-2.0000當A為奇異矩陣時，Ax=b的特解不存在或者不唯一。這種情況下用求非奇異矩陣恰定系統(tǒng)特解的命令不再適用。對于這種系統(tǒng)，可以先通過命令rref([Ab])將增廣矩陣[Ab]約化為行階梯形式并查看結果。如果A的行秩小于增廣矩陣的行秩，則方程組無解；反之，可以利用命令x=pinv(A)*b來求解。

【例】求奇異矩陣恰定系統(tǒng)的特解。在命令窗輸入：

>>A=[137

-144

11018];

>>b=[5;2;12];

>>rref([Ab]),x=pinv(A)*b運行結果：

ans=

1.00000

2.28572.0000

0

1.00001.57141.0000

0

0

0

0

x=

0.3850

-0.1103

0.7066這種情況下，A的行秩等于增廣矩陣的行秩，求得的解即精確解。如果A不變，b=[5;2;11]，在命令窗輸入：

>>rref([Ab])

運行結果：

ans=

1.00000?

2.28570

0?1.0000?1.57140

00

0

1.0000

顯然，A的行秩小于增廣矩陣的行秩，第三行表示方程“0=1”，因而無解。

2．超定系統(tǒng)

超定系統(tǒng)的未知量個數(shù)少于方程組的個數(shù)。超定系統(tǒng)常用于實驗數(shù)據(jù)的曲線擬合問題。

【例】

采集一個頻率為1、直流分量為0的正弦信號，橫軸(時間軸)與縱軸上的點分別為

那么基函數(shù)應該選取sin(t)和cos(t)，問題等價于求解滿足a?sin(t)+b?cos(t)=y的a和b。在M文件中輸入：圖4-1原始數(shù)據(jù)與擬合曲線

3．欠定系統(tǒng)

欠定系統(tǒng)的未知量個數(shù)多于方程組的個數(shù)。如果A的行秩不小于其增廣矩陣的行秩，其解存在無窮多個?？梢杂妹顇=A\b求其特解。

【例】求欠定系統(tǒng)的特解。

在命令窗輸入：

>>A=[6873

6874];

>>b=[2;2];

>>x=R\b運行結果：

x=

0

0.2500

0

04.1.3逆矩陣與偽逆矩陣

對于一個非奇異方陣A，其行列式不為零，因而存在逆矩陣。使用MATLAB命令inv(A)可以求出非奇異方陣A的逆矩陣。如果A的行列式等于零，或者A不是方陣，那么A的逆矩陣將不存在。這時，只有通過計算A的偽逆矩陣來滿足某些特殊的需求，通過命令B=pinv(A)可以求得A的偽逆矩陣B，并且A與B滿足ABA=A、BAB=B。4.1.4矩陣的分解

MATLAB中的線性方程組求解都是基于以下三種矩陣分解，對應函數(shù)分別為chol、lu以及qr：

●?Cholesky分解：分解為對稱正定上下三角矩陣。

●?LU分解：分解為下三角矩陣和上三角矩陣。

●?QR分解：分解為正交矩陣和非奇異上三角矩陣。

1．Cholesky分解

實際應用中，有一類常見且極重要的矩陣，即對稱正定矩陣。如果對于任何n維非零向量x，對稱矩陣A都滿足xTAx?>?0，那么A是對稱正定矩陣?？梢宰C明，對稱正定矩陣的各階順序主子行列式都大于零，并且存在單位上三角矩陣G，使得A=GTG。這種分解稱為Cholesky分解。

MATLAB提供的chol函數(shù)可以直接實現(xiàn)對稱正定矩陣的Cholesky分解，調(diào)用格式為G=chol(A)。

【例】Cholesky分解。

在命令窗輸入：

>>A=pascal(3),G=chol(A)

運行結果：

A=

111

123

136

G=

111

012

001

2．LU分解

MATLAB中的許多矩陣運算都以LU分解為基礎，如方陣的求逆操作、求行列式等。LU分解將方陣A分解為一個下三角矩陣L和一個上三角矩陣U的乘積，常用的調(diào)用格式有兩種：

●?[L,U]=lu(A)：得到一個上三角矩陣U和一個準下三角矩陣L，使得A=LU。準下三角矩陣L實際上是下三角矩陣的置換形式。

●?[L,U,P]=lu(A)：得到一個單位下三角矩陣L、上三角矩陣U和一個由0、1組成的置換矩陣P，使得PA=LU。

【例】LU分解。

在命令窗輸入：

>>A=[102030

204580

3080171];

>>[L1U1]=lu(A),[L2U2P2]=lu(A)運行結果：

L1=

0.3333?0.8000?1.0000

0.6667?1.0000

0

1.0000

0

0

U1=

30.000080.0000

171.0000

0

-8.3333

-34.0000

0

0

0.2000L2=

1.00000

0

0.6667?1.00000

0.3333?0.80001.0000

U2=

30.000080.0000?171.0000

0

-8.3333

-34.0000

0

0

0.2000

P2=

001

010

100

LU分解可用于求矩陣A的逆和行列式：

det(A)=det(L)*det(U)

inv(A)=inv(U)*inv(L)

3．QR分解

矩陣的正交分解是將矩陣分解為一個正交矩陣Q和一個上三角矩陣R的乘積。此種分解適用于任意矩陣，是非常重要的分解形式。其調(diào)用格式如下：

●?[Q,R]=qr(A)：生成一個與A同階的上三角矩陣R和一個正交矩陣Q，使得A=QR。

●?[Q,R,P]=qr(A)：得到一個置換矩陣P、上三角矩陣R和正交矩陣Q，使得AP=QR，選擇置換矩陣P使得abs(diag(R))按降序排列?！?[Q,R]=qr(A,0)：生成一種“經(jīng)濟”的分解。如果矩陣A是m×n階，并且m>n，則只計算Q的前n列。

●?[Q,R,P]=qr(A,0)：生成一種“經(jīng)濟”的分解。其中P是一個置換向量，使得Q*R=A(:,P)，選擇置換列向量使得abs(diag(R))按降序排列。

【例】QR分解。在命令窗輸入：

>>A=?[111

222

333];

>>[Q,R]=qr(A)

運行結果：

Q=

-0.26730.95620.1195

-0.5345?-0.0439-0.8440

-0.8018?-0.28950.5228R=

?-3.7417-3.7417?-3.7417

?0

0.0000

0.0000

?0

0

-0.00004.1.5矩陣的非線性運算

MATLAB支持矩陣的幾種非線性運算，分別為求矩陣的正整數(shù)次冪、負冪、分數(shù)次冪、矩陣元素的冪以及指數(shù)矩陣。

1．矩陣的正整數(shù)冪

利用命令A^p可以求矩陣A的p次冪，其中A為方陣，p為正整數(shù)。

【例】矩陣的正整數(shù)次冪。

在命令窗輸入：

>>A=[111;123;136];

>>X=A^2

運行結果：

X=

3610

61425

102546

2．矩陣負冪、分數(shù)次冪

如果方陣A為非奇異矩陣，則A的逆存在，可以用A^(-p)求A的逆矩陣A-1的p次冪。當p為分數(shù)時，則可以利用A^p求A的分數(shù)次冪。

【例】矩陣的負冪和分數(shù)次冪。

繼續(xù)在命令窗輸入：

>>Y=A^(-2)運行結果：

Y=

19.0000?-26.000010.0000

-26.000038.0000?-15.0000

10.0000?-15.00006.0000

再輸入：

>>Z=A^(1/2)運行結果：

Z=

0.87750.43870.1937

0.43871.00990.8874

0.19370.88742.2749

3．矩陣元素的冪

利用?.^?運算符可以求由矩陣每個元素的冪組成的矩陣。

【例】矩陣元素的冪。

繼續(xù)在命令窗輸入：

>>A.^2

運行結果：

ans=

111

149

1936

4．指數(shù)矩陣

對于一個線性系統(tǒng)，其狀態(tài)方程通?？梢詫懽鳎?/p>

其中，x=x(t)為關于t的函數(shù)向量，A為與t無關的矩陣，則系統(tǒng)的解可以表示為

【例】非時變線性系統(tǒng)的解。

在命令窗輸入：

>>A=?[0-5-1

51-9

-48-6];

>>x0=[1;0;1];

>>X=[];

fort=0:.01:1

X=[Xexpm(t*A)*x0];

end

>>plot3(X(1,:),X(2,:),X(3,:),‘-.’)

運行結果如圖4-2所示。圖4-2非時變線性系統(tǒng)的解曲線4.1.6特征值與特征向量

對于一個方陣A，它的特征值l和特征向量x滿足Ax=lx。如果將所有特征值l組成對角矩陣，將所有特征向量x組成矩陣V，則有AV=V。根據(jù)V是否奇異，可將A化為對角陣或約當陣。

1．化為對角陣

當V為非奇異矩陣時，A的每個特征值都對應一個特征向量，可以寫成A=VV-1的形式。MATLAB中，通常使用命令[V,D]=eig(A)求出矩陣A的所有特征值組成的對角矩陣D和A的所有特征向量組成的矩陣V。

【例】化對角陣。

在命令窗輸入：

>>A=[61116

-9-10-13

5712];

>>[V,D]=eig(A)運行結果：可以驗證，V*D*inv(V)的值與A的值相等。

2．化為約當陣

當V為奇異矩陣時，A的獨立特征向量數(shù)小于A的階數(shù)，因而不可以寫成A=VV-1的形式對其進行對角化，但是可以將A化為與對角型相似的約當型。MATLAB中，通常使用命令[X,J]=jordan(A)求出矩陣A的約當矩陣J和廣義特征向量組X，使得A=XJX-1。

【例】化對角陣。在命令窗輸入：

>>A=[61219

-9-20-33

4915];

>>[X,J]=jordan(A)運行結果：

X=

-1.75001.50002.7500

3.0000-3.0000?-3.0000

-1.25001.50001.2500

J=

-100

011

001

可以驗證，X*J*inv(X)的值與A的值相等。

3．Schur分解

MATLAB中還提供了無需特征值分解的Schur分解，即將A分解為正交矩陣U和分塊上三角矩陣S(由1×1、2×2塊組成)，且滿足A=USUT。調(diào)用格式為[U,S]=schur(A)。

【例】Schur分解。

在命令窗輸入：

>>A=[61219

-9-20-33

4?915];

>>[U,S]=schur(A)運行結果：

U=

-0.47410.66480.5774

0.81270.07820.5774

-0.3386-0.74300.5774

S=

-1.000020.7846-44.6948

0

1.0000-0.6096

0?0

1.0000

可以驗證，U*S*U'?的值與A的值相等。4.1.7奇異值分解

奇異值分解在矩陣分析中占有極為重要的位置。對于一個矩陣A，它的一個奇異值和對應的奇異向量u、v滿足Av=u，ATu=v。如果將所有奇異值組成對角矩陣S，將奇異列向量u、v組成正交矩陣U、V，則AV=US，ATU=VS，A=USVT。

MATLAB中提供的奇異值分解函數(shù)svd的調(diào)用格式如下：

●?[U,S,V]=svd(A)：得到一個與A同階的奇異對角矩陣S，兩個酉矩陣U和V。若A為m×n陣，則U為m×m陣，V為n×n陣。S中奇異值非負且按降序排列。

●?[U,S,V]=svd(A,0)：得到一種“經(jīng)濟”的分解，只計算出矩陣U的前n列，矩陣S的大小為n×n。

【例】奇異值分解。

在命令窗輸入：

>>A=[94

68

27];

>>[U,S,V]=svd(A)運行結果：

U=

-0.61050.71740.3355

-0.6646-0.2336?-0.7098

-0.4308-0.65630.6194

S=

14.93590

0

5.1883

0

0

V=

-0.69250.7214

-0.7214-0.6925繼續(xù)在命令窗輸入：

>>[U,S,V]=svd(A,0)

運行結果：

U=

-0.61050.7174

-0.6646-0.2336

-0.4308-0.6563

S=

14.93590

0

5.1883

V=

-0.69250.7214

-0.7214-0.6925

可以驗證，U*S*V'?的值與A的值相等。4.2.1多項式

MATLAB提供了豐富的多項式運算函數(shù)，包括計算多項式的根、矩陣特征多項式的系數(shù)、估計多項式的值、計算多項式的乘除法、多項式的求導、多項式擬合以及多項式之比與部分分式展開式之間的相互轉(zhuǎn)化。表4-2給出了有關多項式運算的函數(shù)。4.2多項式與插值表4-2多項式運算函數(shù)

1．多項式的表示方法

MATLAB中的多項式可表示為按降冪排列的多項式系數(shù)組成的行向量。例如，多項式p(x)=x3-2x+3可表示為如下向量形式：

p=[10-23];

2．多項式的根

利用roots函數(shù)可以計算多項式的根。例如求方程x3-1=0的根。

【例】求多項式的根。

在命令窗輸入：

>>p=[100-1];

>>r=roots(p)

運行結果：

r=

-0.5000+0.8660i

-0.5000-0.8660i

1.0000

由多項式的根可以反推出多項式，仍以方程x3-1=0的根為例。

【例】由多項式的根反推多項式。

繼續(xù)在命令窗輸入：

>>poly(r)

運行結果：

ans=

1.00000.00000.0000-1.0000

3．特征多項式

利用poly函數(shù)也可以計算矩陣特征多項式的系數(shù)，且系數(shù)按降冪排列。

【例】

求矩陣的特征多項式。

在命令窗輸入：

>>A=[13?-1

526

901];

>>poly(A)

運行結果：

ans=

1.0000-4.0000-1.0000-167.0000

4．多項式估值

對于一個多項式，可以利用polyval函數(shù)估計指定自變量時的函數(shù)值。例如估計多項式f(x)=x2-2在x=3處的值。

【例】

多項式的單點估值。

在命令窗輸入：

>>f=[10-2];

>>polyval(f,3)

運行結果：

ans=

7用戶也可以在自變量x取多個值時對多項式進行估值，此時polyval函數(shù)的第二個輸入?yún)⒘靠梢允窍蛄炕蚓仃嚒?/p>

【例】

多項式的矩陣估值。

繼續(xù)在命令窗輸入：

>>A=[12

34];

>>polyval(f,A)運行結果：

ans=

-12

714

但是很多情況下需要估計自變量為矩陣的多項式的值，如f(X)=X2-2I的估值。這時需要用到矩陣多項式估值函數(shù)polyvalm。

【例】

自變量為矩陣的多項式估值。

繼續(xù)在命令窗輸入：

>>X=[12

34];

>>polyvalm(f,X)

運行結果：

ans=

510

1520

5．多項式的乘除法

利用卷積函數(shù)conv和反卷積函數(shù)deconv可以分別實現(xiàn)多項式的乘法和除法。例如對多項式a(x)?=x2+2x+3和多項式b(x)=3x2+2x+1進行乘法和除法運算。

【例】多項式的乘法。

在命令窗輸入：

>>a=[123];

>>b=[321];

>>c=conv(a,b)運行結果：

c=

381483

【例】多項式的除法。

繼續(xù)在命令窗輸入：

>>[q,r]=deconv(c,a)運行結果：

q=

321

r=

00000

其中，q代表商式的系數(shù)，r代表余式的系數(shù)。

6．多項式求導

利用polyder函數(shù)可以對多項式求導，還可以對兩個多項式相乘或相除之后得到的多項式求導。

【例】多項式求導。

在命令窗輸入：

>>p=[237];

>>q=polyder(p)運行結果：

q=

43

當polyder函數(shù)有兩個輸入?yún)⒘?，一個輸出參量時，實現(xiàn)輸入多項式相乘后求導。

【例】多項式相乘后求導。

在命令窗輸入：

>>a=[135];

>>b=[246];

>>c=polyder(a,b)運行結果：

c=

8305638

當polyder函數(shù)有兩個輸入?yún)⒘浚瑑蓚€輸出參量時，實現(xiàn)輸入多項式相除后求導。第一個輸入?yún)⒘繛楸怀龜?shù)，第二個輸入?yún)⒘繛槌龜?shù)；第一個輸出參量為分子，第二個輸出參量為分母。

【例】多項式相除后求導。

繼續(xù)在命令窗輸入：

>>[q,d]=polyder(a,b)

運行結果：

q=

-2-8-2

d=

416404836

7．多項式擬合

利用polyfit函數(shù)可以尋找最小二乘法擬合的多項式系數(shù)。調(diào)用格式為p=polyfit(x,y,n)。其中，x、y為需要擬合的原始數(shù)據(jù)向量，n為多項式的次數(shù)。

【例】多項式擬合。在命令窗輸入：

>>x=[12345];

>>y=[5.544123.2291467];

>>p=polyfit(x,y,4);

>>x2=1:.1:5;

>>y2=polyval(p,x2);

>>plot(x,y,‘*’,x2,y2)

運行結果如圖4-3所示。圖4-3原始數(shù)據(jù)與擬合曲線

8．多項式的部分分式展開

residue函數(shù)用于實現(xiàn)兩個多項式之比與部分分式展開式之間的相互轉(zhuǎn)化。對于多項式b(s)和a(s)，有命令[r,p,k]=residue(b,a)將多項式之比轉(zhuǎn)化為部分分式展開式；命令[b,a]=residue(r,p,k)將部分分式展開式轉(zhuǎn)化為多項式之比。其中，r為由rn組成的向量，p為由pn組成的向量，k為ks，b為b(s)，a為a(s)。下面以傳遞函數(shù)為例進行轉(zhuǎn)換。

【例】多項式之比轉(zhuǎn)化為部分分式展開式。

在命令窗輸入：

>>b=[2-1];

>>a=[12-3];

>>[r,p,k]=residue(b,a)

運行結果：

r=

1.7500

0.2500

p=

-3.0000

1.0000

k=

[?]

【例】部分分式展開式轉(zhuǎn)化為多項式之比。

繼續(xù)在命令窗輸入：

>>[b1,a1]=residue(r,p,k)

運行結果：

b1=

2-1

a1=

1.00002.0000-3.00004.2.2插值

插值是利用已知數(shù)據(jù)的規(guī)律，計算未測的中間點值或預測數(shù)據(jù)所在范圍以外的值的過程。MATLAB提供了豐富的插值方法，使用戶可以在光滑性、執(zhí)行速度與占用存儲空間之間進行選擇。

MATLAB的polyfun路徑下的插值函數(shù)如表4-3所示。表4-3polyfun路徑下的插值函數(shù)

1．一維插值

MATLAB中的一維插值包括多項式插值和基于FFT的插值。函數(shù)interp1常用于實現(xiàn)一維多項式插值，這在數(shù)據(jù)分析和曲線擬合中非常有用。該函數(shù)由已知數(shù)據(jù)(xi,yi)構造一個多項式P(x)，使得在xi點處的函數(shù)值P(xi)等于yi，其余x點處的值由P(x)決定。

函數(shù)interp1的調(diào)用格式為y=interp1(xi,yi,x,method)。其中，xi、yi為已知數(shù)據(jù)點組成的向量；x為插值點組成的向量；y為插值點向量處的函數(shù)值向量；method為可選字符串，用于指定插值方法：

●?method='nearest'：最鄰近插值。該方法將插值點x處的函數(shù)值y置為與之最接近的已知點xi處的值，即yi。●?method=‘linear’：線性插值。該方法用直線連接相鄰的已知數(shù)據(jù)點(xi,yi)，插值點x處的函數(shù)值y由該直線確定。

●?method=‘spline’：三次樣條插值。用三次函數(shù)連接相鄰的已知數(shù)據(jù)點(xi,yi)，插值函數(shù)在已知點處連續(xù)，并且它的一階導數(shù)和二階導數(shù)也連續(xù)。

●?method='pchip'OR'cubic'：三次插值。這兩種方法等價，用于實現(xiàn)分段光滑Hermite插值，并保留數(shù)據(jù)的單調(diào)性和形狀。選擇使用插值方法時，還應該考慮到各種方法的光滑性、執(zhí)行速度與占用存儲空間：

●最鄰近插值：速度最快，但是光滑性最差。

●線性插值：比最鄰近插值占用的內(nèi)存更多，執(zhí)行速度稍慢，插值結果是連續(xù)的，但是端點處斜率不連續(xù)。

●三次樣條插值：最長的相對運行時間，比三次插值占用的內(nèi)存要少。產(chǎn)生最光滑的曲線。但是如果輸入數(shù)據(jù)不均勻或者間距太小，則會得到不希望得到的結果?！袢尾逯担罕茸钹徑逯捣ɑ蚓€性插值法需要占用更多的內(nèi)存和執(zhí)行時間，但是得到的函數(shù)及其一階導數(shù)都是連續(xù)的。

如果插值點x超出了已知點xi的取值范圍，則默認情況下返回函數(shù)值為NaN。下面的例子將對以上幾種插值方法進行比較：

【例】多項式插值。創(chuàng)建一個M文件，輸入如下代碼：

%OriginalData

x=1:10;

y=[?0.0701-0.06890.0505-0.03260.0196...

-0.01130.0062-0.00330.0018-0.0009];

%Interpolation

xi=1:0.1:10;yi1=interp1(x,y,xi,'nearest');

yi2=interp1(x,y,xi,'linear');

yi3=interp1(x,y,xi,'spline');

yi4=interp1(x,y,xi,'pchip');

figure;

subplot(221);

plot(x,y,'*',xi,yi1);

title('最鄰近插值');subplot(222);

plot(x,y,'*',xi,yi2);

title('線性插值');

subplot(223);

plot(x,y,'*',xi,yi3);

title('三次樣條插值');

subplot(224);

plot(x,y,'*',xi,yi4);

title('三次插值');

程序運行結果如圖4-4所示。圖4-4各種插值方法的對比通過函數(shù)interpft可以實現(xiàn)基于FFT的一維插值。該方法首先計算周期函數(shù)值向量的傅里葉變換，再計算更多點數(shù)的傅里葉反變換，調(diào)用格式為y=interpft(x,n)。其中，x為等間隔周期函數(shù)采樣值，n為需要返回的等間隔點數(shù)。

【例】多項式插值。在命令窗輸入：

>>y=sin(0:9);

>>N=length(y);

>>L=5;

>>M=N*L;

>>x=0:L:L*N-1;

>>xi=0:M-1;

>>yi=interpft(y,M);

>>plot(x,y,'o',xi,yi);

>>legend('原始數(shù)據(jù)','插值函數(shù)')

運行結果如圖4-5所示。圖4-5原始數(shù)據(jù)與插值函數(shù)

2．二維插值

二維插值在圖像處理和數(shù)據(jù)可視化中非常重要，可以利用函數(shù)interp2來實現(xiàn)。常用調(diào)用格式為ZI=interp2(X,Y,Z,XI,YI,method)，其中，X、Y為已知數(shù)據(jù)點，Z為已知點處的函數(shù)值，XI、YI為插值點，ZI為求得的插值點的函數(shù)值。method為可選字符串，用于指定插值方法：

●?method='nearest'：最鄰近插值。該方法通過已知數(shù)據(jù)生成分段光滑的常數(shù)曲面，插值點處的函數(shù)值即最鄰近點處的值?！駇ethod=‘linear’：雙線性插值。該方法由已知點生成一個雙線性曲面，插值點處的函數(shù)值由最近的4個已知點確定。

●?method=‘cubic’：雙三次插值。該方法由已知點生成一個雙三次曲面，插值點處的函數(shù)值由最近的16個已知點確定。這種方法產(chǎn)生的曲面最光滑，因而在圖像處理中有重要的應用。在要求插值函數(shù)及其導數(shù)連續(xù)時，這種方法非常必要。

下面的例子將對以上幾種二維插值方法進行比較。

【例】二維插值。

創(chuàng)建一個M文件，輸入如下代碼：

z=peaks(x,y);

[xi,yi]=meshgrid(-3:0.25:3);

zi1=interp2(x,y,z,xi,yi,‘nearest’);

zi2=interp2(x,y,z,xi,yi,‘bilinear’);

zi3=interp2(x,y,z,xi,yi,‘bicubic’);

figure;

subplot(221);

surf(x,y,z);

title('原始數(shù)據(jù)');subplot(222);

surf(xi,yi,zi1);

title('最鄰近插值');

subplot(223);

surf(xi,yi,zi2);

title('雙線性插值');

subplot(224);

surf(xi,yi,zi3);

title('雙三次插值');

figure;

subplot(221);

contour(x,y,z);

title(‘原始數(shù)據(jù)’);

subplot(222);

contour(xi,yi,zi1);

title(‘最鄰近插值’);

subplot(223);

contour(xi,yi,zi2);

title(‘雙線性插值’);

subplot(224);

contour(xi,yi,zi3);

title(‘雙三次插值’);

運行這段程序，得到插值曲面圖和等高線圖如圖4-6和圖4-7所示。圖4-6插值曲面圖圖4-7等高線圖4.3.1快速傅里葉變換的概念

序列x(n)的離散傅里葉變換(DFT)及離散傅里葉反變換(IDFT)定義如下式所示：4.3快速傅里葉變換通常FFT算法采取的是按時間抽取(DIT)或按頻率抽取(DIF)基2算法，即把一個N(N=2m)點序列DFT分裂為2個N/2點序列的DFT，再把這2個N/2的DFT分裂為4個N/4點序列的DFT，如此不斷分裂(共分裂成lb?N級)，直至分為N/2個2點序列的DFT。不難看出，N點序列的DFT計算量為復乘N2次，復加N(N-1)次；而基2FFT的計算量為復乘N?lb?N次，復加(N?-1)lb?N次，大大節(jié)省了運算量。4.3.2快速傅里葉變換的應用

有關FFT的函數(shù)如表4-4所示。表4-4FFT的相關函數(shù)通常用fft(x)或fft(x,n)命令來求序列的DFT。若序列x的長度為N，fft(x)則實現(xiàn)計算x的N點DFT的功能；而fft(x,n)則實現(xiàn)計算x的n點DFT，n>N時對序列補零，n<N時對序列進行截斷。建議使用2m點序列的DFT，這樣能夠?qū)崿F(xiàn)FFT快速算法。

利用fft函數(shù)分析序列的頻率分量非常有效。下面給出一個計算周期正弦信號諧波的實例，設N點序列為

【例】周期信號的頻譜分析。

創(chuàng)建一個M文件，輸入如下代碼：

%SineWave

N=50;

n=0:N-1;

w1=2*pi/N;

w2=3*2*pi/N;

x=2*sin(w1*n)+0.5*sin(w2*n);

plot(n,x);

xlabel(‘n’);

ylabel('x');

%FrequencySpectrum

H=abs(fft(x,N))/N*2;

figure;

stem(0:N/2-1,H(1:N/2),‘markersize’,1);

xlabel(‘k’);

ylabel(‘Magnitude’);

運行這段程序，得到序列x的波形圖和幅度頻率圖分別如圖4-8和圖4-9所示。

可以看出，直流分量(k=0)以及其他頻率分量的幅值為0，基波(k=1)幅值為2，3次諧波(k=3)幅值為0.5，符合構成序列x的各頻率分量幅值。圖4-8x的波形圖圖4-9x的幅度頻率圖函數(shù)的函數(shù)是指一個函數(shù)以其他函數(shù)作為輸入?yún)⒘?。例如函?shù)的畫圖函數(shù)fplot，其調(diào)用格式為fplot(@fun,limits)。其中，@fun是需要繪制的函數(shù)句柄，limits為指定的繪制范圍。表4-5給出了本節(jié)所要討論的函數(shù)。4.4函?數(shù)?的?函?數(shù)表4-5函?數(shù)?的?函?數(shù)4.4.1函數(shù)的表示方法

MATLAB的數(shù)學函數(shù)可以有兩種表示方法，一種方法是將表達式保存在M文件內(nèi)作為MATLAB函數(shù)，另一種方法是直接作為匿名函數(shù)而表示。例如對于函數(shù)創(chuàng)建一個名為ff的M文件，輸入如下代碼并保存為ff.m：

functiony=ff(x)

y=(2*x-3)./(x.^2-5*x-6);這樣，可以創(chuàng)建一個ff函數(shù)的句柄fh，可以按照直接調(diào)用函數(shù)ff的方法調(diào)用句柄fh。

【例】調(diào)用M文件函數(shù)的句柄。

在命令窗輸入：

>>fh=@ff;

>>y1=ff(3),y2=fh(3)

運行結果：

y1=

-0.2500

y2=

-0.2500

表達數(shù)學函數(shù)的另一種方法是創(chuàng)建匿名函數(shù)，仍以前面的函數(shù)為例進行說明。

【例】創(chuàng)建匿名函數(shù)。

在命令窗輸入：

>>fh=@(x)(2*x-3)./(x.^2-5*x-6);

>>y=fh(3)

運行結果：

y=

-0.2500

MATLAB同樣支持多元函數(shù)的匿名函數(shù)。

【例】創(chuàng)建多元匿名函數(shù)。

在命令窗輸入：

>>fh=@(x,y)x.^2+2*y.^2;

>>y=fh(3,1)

運行結果：

y=

11

利用fplot函數(shù)可以繪制已知函數(shù)的圖像。

【例】繪制數(shù)學函數(shù)圖像。

在命令窗輸入：

>>fplot(@(x)2*sin(3*x),[-11]);

運行結果如圖4-10所示。圖4-10利用fplot函數(shù)繪圖4.4.2函數(shù)的最小值與零點

函數(shù)的最小值與零點問題對于最優(yōu)化問題的求解非常重要。

1．一元函數(shù)的最小值

利用fminbnd函數(shù)可以求出一元函數(shù)在指定區(qū)間內(nèi)的最小值。fminbnd函數(shù)的常用調(diào)用格式為x=fminbnd(fun,x1,x2)，其中，fun為一元函數(shù)的句柄，x1、x2分別為指定區(qū)間的下限和上限。

【例】

計算一元函數(shù)的最小值。

在命令窗輸入：

>>x=fminbnd(@(x)-sin(x),0,2*pi)

運行結果：

x=

1.5708

2．多元函數(shù)的最小值

利用fminsearch函數(shù)可以求出多元函數(shù)在指定向量附近的最小值。fminsearch函數(shù)的常用調(diào)用格式為x=fminsearch(fun,x0)，其中，fun為多元函數(shù)的句柄，x0為指定初始向量。

【例】計算多元函數(shù)的最小值。

創(chuàng)建一個M文件函數(shù)，輸入如下代碼并保存為val3.m：

functionf=val3(v)

x=v(1);

y=v(2);

z=v(3);

f=abs(x)+y.^2+z.^4;在命令窗輸入：

>>x=fminsearch(@val3,[1,1,1])

運行結果：

x=

1.0e-004*

-0.00000.0004-0.3072

3．一元函數(shù)的零點

利用fzero函數(shù)可以求出一元函數(shù)在指定值附近或者指定區(qū)間內(nèi)的最小值。fzero函數(shù)的常用調(diào)用格式為x=fzero(fun,x0)，其中，fun為一元函數(shù)句柄，當x0為標量時，求得與x0最近的零點；當x0為二值向量，且fun(x0(1))與fun(x0(2))異號時(否則會產(chǎn)生出錯信息)，求得該區(qū)間內(nèi)的零點。

【例】計算一元函數(shù)的零點。在命令窗輸入：

>>y1=fzero(@(x)cos(x),0)

運行結果：

y1=

-1.5708

或輸入：

>>y2=fzero(@(x)cos(x),[0pi])

運行結果：

y2=

1.57084.4.3數(shù)值積分

在大多數(shù)情況下很難求得一個函數(shù)的原函數(shù)，或者就沒有有限形式的原函數(shù)解析表達式，因而計算函數(shù)的定積分往往采用數(shù)值積分的方法。例如通過quad(@humps,0,1)命令用Simpson法計算函數(shù)humps在0～1區(qū)間上的積分。實際中往往遇到參數(shù)方程的曲線積分，MATLAB中通常運用quad函數(shù)或quadl函數(shù)計算曲線的長度。例如計算曲線x(t)=sin(t)，y(t)=cos(t)，z(t)=t在t[0,]上的長度，其表達式為

【例】

計算曲線長度。在命令窗輸入：

>>leng=quad(@(t)sqrt(cos(t).^2+sin(t).^2+1),0,3*pi)運行結果：

leng=

13.3286利用dblquad函數(shù)可以計算形如的二重積分，dblquad函數(shù)的調(diào)用格式為result=dblquad(fun2,xmin,xmax,ymin,ymax)，其中，fun2代表二元函數(shù)的句柄，xmax、xmin、ymax、ymin分別代表自變量x、y的上下限。

【例】計算f(x,y)?=?ysin(x)?+?xcos(y)的二重積分。

在命令窗輸入：

>>result=dblquad(@(x,y)(y*sin(x)+x*cos(y)),pi,2*pi,0,pi)

運行結果：

result=

-9.86964.4.4嵌套函數(shù)與匿名函數(shù)

在某些情況下不僅需要調(diào)用函數(shù)，而且還需要與該函數(shù)相關的參數(shù)。例如，如果要使用fzero函數(shù)尋找三次多項式x3+bx+c在不同的系數(shù)b、c情況下的零點，則需要一個函數(shù)能夠在系數(shù)b、c變化時計算該多項式，即能夠在調(diào)用fzero函數(shù)的同時提供多項式函數(shù)中的b、c取值。通常的做法有兩種：使用嵌套函數(shù)或使用匿名函數(shù)。

1．嵌套函數(shù)

以三次多項式x3+bx+c為例，對于不同的b、c取值，可以使用嵌套函數(shù)實現(xiàn)零點求解。

【例】利用嵌套函數(shù)求解零點。

創(chuàng)建一個M文件，輸入如下代碼并保存為findzero.m：

functiony=findzero(b,c,x0)

options=optimset(‘Display’,‘off’);%TurnoffDisplay

y=fzero(@poly,x0,options);

functiony=poly(x)%Computethepolynomial.

y=x^3+b*x+c;

end

end如果要尋找x3+x+1在x=0附近的零點，則在命令窗輸入：

>>findzero(1,1,0)

運行結果：

ans=

-0.6823

2．匿名函數(shù)

仍以三次多項式x3+bx+c為例，對于不同的b、c取值，可以使用匿名函數(shù)實現(xiàn)零點求解。

【例】利用匿名函數(shù)求解零點。

創(chuàng)建一個M文件，輸入如下代碼并保存為poly.m：

functiony=poly(x,b,c)%Computethepolynomial.

y=x^3+b*x+c;

如果要尋找在x=0附近的零點，則在命令窗輸入：

>>b=1;

>>c=1;

>>fzero(@(x)poly(x,b,c),0)

運行結果：

ans=

-0.68234.5.1常微分方程初值問題

本小節(jié)主要介紹常微分方程的初值問題的數(shù)值解法。MATLAB中的ode求解器主要針對以下幾種類型的常微分方程初值問題：

●顯式常微分方程形如y‘?=f(t,y)。

●線性隱式常微分方程形如M(t,y)·y’?=f(t,y)，M(t,y)為矩陣。

●完全隱式常微分方程形如f(t,y,y‘)=0。

常微分方程的初值問題的ode求解器如表4-6所示。4.5求解微分方程表4-6ode求解器

1．求解顯式常微分方程和線性隱式常微分方程

這類問題的求解方法的實質(zhì)是從給定初值開始，計算每一個步長內(nèi)的微分方程的數(shù)值積分。對于非剛性問題，常采用ode45、ode23、ode113求解器求解；而對于剛性問題，常采用ode15s、ode23s、ode23t、ode23tb等求解器求解。這7種ode求解器的調(diào)用格式是相同的，如果需要改變求解的計算方法，只需要改變求解器的函數(shù)名即可。最簡單常用的求解器調(diào)用格式為

[t,y]=solver(odefun,tspan,y0,options)其中，solver為使用的ode求解器函數(shù)；odefun為常微分方程組的函數(shù)句柄，形如dydt=f(t,y)或M(t,y)*dydt=f(t,y)，t為時間標量，y為列向量；tspan為指定求解區(qū)間，tspan=[t0tf]則計算[t0tf]區(qū)間的解，tspan=[t0,t1,t2,…,tf]則計算每個時間點0,t1,t2,…,tf(保證時間點單調(diào))上的解；y0為初值列向量；options為改變默認積分屬性的可選參數(shù)。

對于非剛性問題，常用ode45來求解。下面以一個LC振蕩電路為例，說明求解非剛性問題時ode的用法。

【例】

求解振蕩電路(非剛性問題)。

設L=1，C=1，初始電壓u0=1，初始電流i0=0。電路如圖4-11所示。

則電壓、電流滿足如下微分方程：圖4-11電路連接圖

要求解這個電路方程，創(chuàng)建一個M文件，輸入如下代碼并保存為LC.m：

functiondydt=LC(t,y)

dydt=[y(2)

-y(1)];

在命令窗輸入：

>>[t,y]=ode45(@LC,[02*pi],[0;1]);

>>plot(t,y(:,1),'b-');

>>holdon;

>>plot(t,y(:,2),'r:')

>>holdoff;

>>axis([02*pi-11]);

>>legend('電流','電壓')

運行結果如圖4-12所示。圖4-12電路的電壓和電流波形對于剛性問題，它的解在很小的時間區(qū)間內(nèi)可能有大的跳變。這時不再適合用ode45求解，而通常采用ode15s。下面以vanderPol(范德坡)方程為例，說明求解剛性問題時ode的用法。剛性vanderPol方程的函數(shù)名為vdp1000，方程式為在命令窗輸入：

>>[t,y]=ode15s(@vdp1000,[03000],[2;0]);

>>plot(t,y(:,1),‘-’);

運行結果如圖4-13所示。圖4-13vanderPol方程的波形結果

2．求解完全隱式常微分方程

完全隱式常微分方程通常采用ode15i求解器來求解。該方法采用變階反向微分公式，其基本調(diào)用格式為

[t,y]=ode15i(odefun,tspan,y0,yp0,options)

其中，odefun為完全隱式常微分方程組的函數(shù)句柄，形如f(t,y,y')=0，t為時間標量，y為列向量；tspan為指定求解區(qū)間，tspan=[t0tf]則計算[t0tf]區(qū)間的解，tspan=[t0,t1,t2,…,tf]則計算每個時間點0,t1,t2,…，tf(保證時間點單調(diào))上的解；y0、yp0為代表初值條件y(t0)、y'(t0)的列向量，并且必須滿足f(t0,y0,yp0)=0；options為改變默認積分屬性的可選參數(shù)。下面以Weissinger方程ty2(y')3-y3(y')2+t(t2+1)y'-t2y=0為例，說明ode15i的用法。該方程解析解為在命令窗輸入：

>>t0=1;

>>y0=sqrt(3/2);

>>yp0=0;

>>[t,y]=ode15i(@weissinger,[110],y0,yp0);

>>ytrue=sqrt(t.^2+0.5);

>>plot(t,y,t,ytrue,‘*’);

>>leg