《通信系統(tǒng)中MATLAB基礎與仿真應用》課件-第6章

第6章MATLAB在工程教學中的應用6.1解線性方程組6.2多項式運算6.3曲線擬合6.4插值和樣條6.5數(shù)值積分和微分6.6數(shù)據(jù)分析函數(shù)和傅立葉變換小結(jié)

習題

本章作為本書的基礎知識核心部分，從工程教學的角度，系統(tǒng)地介紹了MATLAB在高等數(shù)學、線性代數(shù)和數(shù)據(jù)處理等方面的應用，同時給出了大量的實例，并結(jié)合圖形圖像處理和數(shù)據(jù)可視化等MATLAB功能進行講述。通過本章的學習，讀者可熟練掌握MATLAB關(guān)于數(shù)學計算方面的內(nèi)容。需要說明的一點是，本章是本著解決問題的目的進行編寫的，因此不再區(qū)分此解決方法是屬于數(shù)值計算還是屬于符號計算，無論采取何種計算方法，都是解決該類問題的一種選擇。在確實需要區(qū)分時，會予以說明。6.1解線性方程組在工程教學和實踐中，對線性方程組的求解是一個重要的內(nèi)容。在實際工程應用、實驗數(shù)據(jù)分析以及進行理論研究等情況下，很多問題都可歸結(jié)為線性方程組的求解。因而，線性方程組求解的應用非常廣泛。求解線性方程組不可避免地要用到矩陣分解的概念，矩陣分解函數(shù)在第3章已有簡單介紹，本節(jié)首先介紹用于解線性方程組的幾個矩陣分解函數(shù)，然后具體講解各類線性方程組的解法。6.1.1矩陣的分解

矩陣分解是一個非常重要的概念，如在求方程組的解、矩陣的特征值、特征向量和矩陣的秩等重要參數(shù)時，都要用到矩陣的分解。在實際工程中的某些場合也要對矩陣進行特定形式的分解?？傊?，無論在理論證明或是科學計算上，對矩陣進行分解操作都十分重要。MATLAB提供了許多矩陣分解函數(shù)，利用這些函數(shù)可以比較容易地進行矩陣分解。下面著重介紹其中幾種矩陣分解的方法，其相關(guān)函數(shù)如表6-1所示。表6-1矩陣分解函數(shù)在MATLAB中，線性方程組的求解主要用到三種基本的矩陣分解，即對稱正定矩陣的cholesky分解、一般方程的gaussian消去法和矩陣的正交分解，這三種分解由函數(shù)chol、lu和qr完成。分解時都使用三角矩陣，其中所有位于對角線以上或以下的元素都為0。矩陣經(jīng)分解成為三角矩陣后，線性方程組不論是用左除還是右除，都可以簡單、快速地求解。

1.對稱正定矩陣的cholesky分解并不是所有的對稱矩陣都可以進行cholesky分解，能進行此種分解的矩陣必須是正定的，即矩陣X的所有對角線元素都是正的，而非對角線元素不會太大。此種分解可用于復矩陣，此時復矩陣必須滿足X′=X是hermitian正定的。設X是一個n×n的對稱正定矩陣，則存在對角線為正的上三角矩陣R，使

X＝R′R。在MATLAB中，cholesky分解的函數(shù)調(diào)用格式有兩種:

·R＝chol(X):X是對稱正定矩陣，R是上三角矩陣，使R′R＝X。如果矩陣X是非正定的，則會給出錯誤信息。

·［R，p］=chol(X):返回兩個參數(shù)，不會給出出錯信息。如果X是正定的，則P等于0，R同上；如果X是非正定的，則P等于正整數(shù)，R是一個階數(shù)為q=p-1的上三角陣，使得R＇*R=X(1:q，1:q)。下面通過例6-1來說明其用法。

【例6-1】cholesky分解。

X＝［22-2；25-4；-2-45］%輸入對稱正定矩陣X

R＝chol(X)%進行cholesky分解

X＝

2

2

-2

2

5

-4

-2

-4

5

R＝

1.4142

1.4142

-1.4142

0

1.7321

-1.1547

001.2910

R＇*R

ans＝

2.00002.0000-2.0000

2.00005.0000-4.0000

-2.0000-4.00005.0000

2.lu分解矩陣的lu分解在工程教學和實踐中非常重要，許多運算都以lu分解為基礎，如方陣的求逆操作inv()、行列式求值det()以及解線性方程組等。lu分解是除法運算的基礎。gaussian消去法或lu分解是將任何方陣表示為一個下三角矩陣L和一個上三角矩陣U的乘積。其調(diào)用格式有兩種:

·［L，U］=lu(X):此調(diào)用格式將得到一個上三角矩陣且存儲在U中和一個準下三角矩陣且存儲在L中，使得X＝LU。準下三角矩陣L實際上是下三角矩陣的轉(zhuǎn)置矩陣。

·［L，U，P］=lu(X):此調(diào)用格式將得到一個主對角元素為1的下三角矩陣L、上三角矩陣U和一個由0和1組成的置換矩陣P，使得PX=LU。注意:進行l(wèi)u分解時，矩陣X必須是方陣。在MATLAB中，lu分解允許線性方程組Ax=b進行如下快速運算:

X=U＼(L＼b)矩陣行列式的值和矩陣求逆也可以利用lu分解進行計算:

det(X)=det(L)*det(U)

inv(X)=inv(U)*inv(L)下面通過例6-2說明如何進行l(wèi)u分解。

【例6-2】lu分解。[HJ2.3mm]

X＝［3-12；12-1；-214］%輸入矩陣X

X＝

3

-1

2

1

2

1

-2

1

4［L，U，P］＝lu(X)%進行l(wèi)u分解

L＝

1.0000

0

0

0.3333

1.0000

0

-0.6667

0.1429

1.0000

U＝

3.0000

-1.0000

2.0000

0

2.3333

-1.6667

0

0

5.5714

P＝

1

0

0

0

1

0

0

0

1

P*X

ans＝

3

-1

2

1

2

-1

-2

1

4

L*U

ans＝

3.0000

-1.0000

2.0000

1.0000

2.0000

-1.0000

-2.0000

1.0000

4.0000

3.qr分解

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

·［Q，R］=qr(X):此調(diào)用格式生成一個與X同階的上三角矩陣R和一個正交矩陣Q，使得X＝QR。

·［Q，R，E］=qr(X):此調(diào)用格式將得到一個置換矩陣E、上三角矩陣R和正交矩陣Q，使得XE=QR，選擇置換矩陣E使得abs(diag(R))遞減。

·［Q，R］=qr(X，0):此調(diào)用格式將生成一種“經(jīng)濟”的分解，如果矩陣X是m×n階，并且m＞n，則僅計算出Q的前n列?！ぃ跶，R，E］=qr(X，0)：此調(diào)用格式將生成一種“經(jīng)濟”的分解，其中E是一個置換向量，使得QR=A(:，E)，選擇列置換向量使得abs(diag(R))遞減。

【例6-3】qr分解。

X=［3-125；12-17；-21-24］％輸入矩陣X

X=

3

-1

2

5

1

2

-1

7

-2

1

-2

4［Q，R］＝qr(X)%進行qr分解

Q＝

-0.8018

0.1543

0.5754

-0.2673

-0.9567

-0.1155

0.5345

-0.2469

0.8083

R＝

-3.7417

0.8018

-2.4054

-3.7417

0

-2.3146

1.7591

-6.9128

0

0

-0.3464

5.3116［Q，R，E］=qr(X，0)

Q＝

-0.5270

0.6463

0.5518

-0.7379

0.0259

-0.6754

-0.4216

-0.7626

0.4905

R＝

-9.4868

-1.4754

-1.3703

0.5270

0

3.4383

-1.4607

2.8437

0

0

-1.4102

0.7971

E＝

4123

4.奇異值分解奇異值分解在矩陣分析中占有極為重要的位置。奇異值分解的定義為:[HJ2.3mm]對于任意的矩陣A∈Cm×n，存在酉矩陣(UnitaryMatrix)，U=［u1，u2，…，um］，V=［v1，v2，…，vn］使得:

UTAV=diag(σ1，σ2，…，σp)σ1≥σ2≥…≥σp≥0，p=min{m，n}其中，σi、ui、vi分別稱為矩陣A的第i個奇異值、左奇異向量和右奇異向量，而它們的組合就稱為奇異值分解三對組。這里的上標“T”表示共軛轉(zhuǎn)置。

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

·［U，S，V］=svd(X):產(chǎn)生一個與矩陣X具有相同維數(shù)的矩陣S，其對角線元素為遞減的非負值，同時得到酉矩陣U和V，使得X=U*S*V′。

·S=svd(X):得到矩陣X的奇異值組成的向量。

·［U，S，V］=svd(X，0):得到一個“經(jīng)濟大小”的分解結(jié)果，如果X是m×n矩陣且m＞n，則只計算U矩陣的前n行，且S矩陣是n×n階的。

【例6-4】奇異值分解。

X=［314；224；1-3-2；123］%輸入矩陣X

X=

3

1

4

2

2

4

1

-3

-2

1

2

3［U，S，V］=svd(X)%進行奇異值分解

U=

0.5873

-0.5075

-0.5954

0.2073

0.5951

-0.0691

0.4057

-0.6903

-0.3132

-0.8424

0.4331

0.06881

0.4503

0.1674

0.5417

0.68981

S=

8.2230

0

0

0

3.2221

0

0

0

0

0

0

0

V=

0.3757

-0.7249

0.5774

0.4400

0.6878

0.5774

0.8157

-0.0371

-0.5774［U，S，V］=svd(X，0)%進行奇異值分解

U=

0.5873

-0.5075

-0.5954

0.5951

-0.0691

0.4057

-0.3132

-0.8424

0.4331

0.4503

0.1674

0.5417

S=

8.2230

0

0

0

3.2221

0

0

0

0

V=

0.3757

-0.7249

0.5774

0.4400

0.6878

0.5774

0.8157

-0.0371

-0.5774奇異值分解也是矩陣求秩運算的基礎，對矩陣A進行奇異值分解S＝svd(A)，得到向量S的非零元素的個數(shù)就是矩陣A的秩。如:

X=［314；224；1-3-2；123］；

S=svd(X)

S=

8.2230

3.2221

0可見矩陣X的秩為2，用求秩運算rank(X)可以驗證這一結(jié)果。

5.Schur分解

schur分解的定義為:

A==USUT

其中，A是方陣，U是一個正交矩陣，S是一個塊上三角矩陣，由對角線上的1×1和2×2塊組成。特征值由S的對角線和塊給出，而U的列比一系列特征向量給出了更多的數(shù)值特性。schur分解可以對缺陷矩陣進行。矩陣的復schur形式矩陣是一個特征值在對角線上的上三角陣；而實schur形式是實特征值在對角線上，復特征值以2×2的塊矩陣排列在對角線上。用函數(shù)T=schur(X)求矩陣X的schur形式矩陣時，如果X是實矩陣，則函數(shù)返回實數(shù)schur形式矩陣；如果X是復矩陣，則函數(shù)返回復數(shù)schur形式矩陣。函數(shù)rsf2csf()可以將實數(shù)形式轉(zhuǎn)化為復數(shù)形式。矩陣X必須是方陣。函數(shù)schur的調(diào)用格式為:

·T=schur(X):僅僅返回schur形式矩陣T。

·［U，T］=schur(X):得到schur形式矩陣T和酉矩陣U，使得X＝U*T*U′和U′*U＝eye(size(X))。

【例6-5】schur分解。

A=［1

2

-1；3

4

0；1

2

3］％輸入方陣A

A=

1

2

-1

3

4

0

1

2

3［U，T］=schur(A)％進行schur分解

U=

0.8262

-0.2294

-0.5145

-0.5571

-0.4680

-0.6860

0.0835

-0.8534

0.5145

T=

-0.4495

0.8575

0.3760

-0.0000

4.4495

2.8501

0

0

4.0000［UU，TT］=rsf2csf(U，T)％對U和T進行轉(zhuǎn)換

UU=

0.8262

-0.2294

-0.5145

-0.5571

-0.4680

-0.6860

0.0835

-0.8534

0.5145

TT=

-0.44495

0.8575

0.3760

0

4.4495

2.8501

0

0

4.0000

6.矩陣的hessenberg分解

矩陣的hessenberg分解由函數(shù)hess完成，該函數(shù)命令的調(diào)用格式為:

·H=hess(A):可求矩陣的hessenberg形式矩陣H，H的第一子對角線以下的元素為0。如果矩陣A是對稱的或是hermitian矩陣，則H是對角三角陣。矩陣A必須是方陣。

·［P，H］=hess(A):得到一個酉矩陣P和一個hessenberg形式矩陣H，使得A＝P*H*P′和P′*P＝eye(slze(A))。

【例6-6】hessenberg分解。

A=［1

2

-1；3

4

0；1

2

3］％輸入方陣A

A=

1

2

-1

3

4

0

1

2

3［P，H］=hess(A)％進行hessenberg分解

P=

1.0000

0

0

0

-0.9487

-0.3162

0

-0.3162

0.9487

H=

1.0000

-1.5811

-1.5811

-3.1623

4.5000

0.5000

0

-1.5000

2.50006.1.2線性方程組的求解在工程教學和計算中，線性方程組的求解是一個很重要的問題。在矩陣表示方法中，線性方程組可以表示為給定兩個矩陣A和B，求X的唯一解，使得:

AX=B或XA=B在MATLAB中，求解線性方程組時，主要采用前面章節(jié)介紹的除法運算符“/”和“＼”求解。如:

·X＝A＼B表示求矩陣方程AX＝B的解；

·X＝B/A表示求矩陣方程XA=B的解。對方程X＝A＼B，要求矩陣A和B有相同的行數(shù)，X和B有相同的列數(shù)，它的行數(shù)等于矩陣A的列數(shù)。方程X＝B/A同理。如果矩陣A不是方陣，其維數(shù)是m×n，則有:

·當m=n時，為恰定方程，尋求精確解；

·當m>n時，為超定方程，尋求最小二乘解；

·當m<n時，為不定方程，尋求基本解，其中至多有m個非零元素。針對不同的情況，MATLAB將采用不同的算法來求解。6.1.3恰定方程組

恰定方程組由n個未知數(shù)的n個方程構(gòu)成，方程有唯一的一組解。其一般形式可用矩陣、向量寫成如下形式:

Ax=b其中，A是方陣，b是一個列向量。在線性代數(shù)教科書中，最常用的方程解法有:

·利用cramer公式求解法；

·利用矩陣求逆解法，即x=A-1b；

·利用gaussian消去法；

·利用lu法求解。一般來說，對于維數(shù)不高、條件數(shù)不大的矩陣，上面四種解法所得的結(jié)果差別不大。前三種解法的真正意義是在其理論上，而不是實際的數(shù)值計算。在MATLAB中，出于對算法穩(wěn)定性的考慮，行列式及其逆的計算大都在lu分解的基礎上進行。在MATLAB中，求解這類方程組的命令十分簡單，直接采用表達式:x=A＼b。

MATLAB的指令解釋器在確認變量A非奇異后，就對它進行l(wèi)u分解，并最終給出解x；若矩陣A的條件數(shù)很大，則MATLAB會提醒用戶注意所得解的可靠性。如果矩陣A是奇異的，則Ax=b的解不存在，或者存在但不唯一；如果矩陣A接近奇異，則MATLAB將給出警告信息；如果發(fā)現(xiàn)A是奇異的，則計算結(jié)果為inf，并且給出警告信息；如果矩陣A是病態(tài)矩陣，則也會給出警告信息。注意:在求解方程時，盡量不要用inv(A)*b命令，而應采用A＼b的解法。因為后者的計算速度比前者快、精度高，尤其當矩陣A的維數(shù)比較大時。另外，除法命令的適用性較強，對于非方陣的A，也能給出最小二乘解。

【例6-7】“左除”法與“求逆”法解恰定方程組時的性能比較。

rand(′state′，12)；％選定隨機種子，產(chǎn)生隨機矩陣A

A=rand(100，100)＋1.e8；％生成(100×100)均勻分布的隨機陣％加大數(shù)值是要使矩陣A的條件數(shù)增大

X=ones(100，1) ％生成100行1列的向量

b＝A*x；

flops(0) ％浮點運算計數(shù)器置0

tic ％啟動計時器StopwatchTime，開始計時

y=inv(A)*b； ％求逆法解方程運算次數(shù)

tic=toc ％關(guān)閉計時器，并顯示解方程所用的時間

ci=flops ％求逆法解方程運算次數(shù)

erri=norm(y-x) ％結(jié)果與精確解的范-2誤差

errif=norm(A*y-b) ％方程的范-2誤差％％％

flops(0)

tic

y=A＼b； ％用除法解方程

td=toc ％關(guān)閉計時器，并顯示解方程所用的時間

cd=flops ％除法解方程運算次數(shù)

errd=norm(y-x) ％結(jié)果與精確解的范-2誤差

errdf=norm(A*y-b) ％方程的范-2誤差運算結(jié)果:

ti=

0.0600 ％求逆法解方程所用的時間

ci=

2070322 ％求逆法解方程所用的運算次數(shù)

erri=

3.0708e-004

errif=

6.6280e+004

td= ％除法解方程所用的時間

0

cd= ％除法解方程所用的運算次數(shù)

741872

errd=

3.2243e-004

errdf=

2.0095e-005

從本例執(zhí)行的結(jié)果可知，求逆法解方程所需運算次數(shù)是除法解方程的2.5倍，時間上除法幾乎是“機器零”。6.1.4超定方程組

對于方程組Ax=b，A為n×m矩陣，如果A列滿秩，且n>m，則方程沒有精確解，此時稱方程組為超定方程組。線性超定方程經(jīng)常遇到的問題是數(shù)據(jù)的曲線擬合。對于超定方程，在MATLAB中，利用左除命令(x=A＼b)來尋求它的最小二乘解；還可以用廣義逆來求解，即x=pinv(A)，所得的解不一定滿足Ax=b，x只是最小二乘意義上的解。左除的方法建立在奇異值分解基礎之上，由此獲得的解最可靠；廣義逆法建立在對原超定方程直接進行householder變換的基礎上，其算法可靠性稍遜于奇異值分解，但速度較快?！纠?-8】求超定方程組的解

A=［2

-1

3；3

1

-5；4

-1

1；1

3-1

3］％輸入矩陣A

A=

2

-1

3

3

1

-5

4

-1

1

1

3-

13

b=［3

0

3

-6］′；

rank(A)

ans=

3

x1=A＼b％左除解方程

x2=pinv(A)b ％廣義逆求解

x1=

1.0000

2.0000

1.0000

x2=

1.0000

2.0000

1.0000

A*x1-b

ans=

1.0e-014*

-0.0888

-0.0888

-0.1332

0可見x1并不是方程Ax=b的精確解，用x2=pinv(A)*b所得的解與x1相同。6.1.5欠定方程組

欠定方程組未知量個數(shù)多于方程個數(shù)，但理論上有無窮個解。MATLAB將尋求一個基本解，其中最多只能有m個非零元素。特解由列主元qr分解求得?！纠?-9】解欠定方程

A=［1-211；1-21-1；1-215］ ％輸入矩陣A

A=

1

-2

1

1

1

-2

1

-1

1

-2

1

5

b=［1-15］＇；

x1=A＼b％左除法解方程

Warning:Rankdeficient，rank=2tol=4.6151e-015

x1=

0

-0.0000

0

1.0000

x2=pinv(A)*b ％用廣義逆解方程

x2=

0

-0.0000

0.0000

1.00006.1.6方程組的非負最小二乘解在某種情況下，所求的線性方程組的解出現(xiàn)負數(shù)是沒有意義的。雖然方程組可以得到精確解，但卻不能取負值解。在這種情況下，其非負最小二乘解比方程的精確解更有意義。在MATLAB中，求非負最小二乘解通常用函數(shù)nnls，其調(diào)用格式為:

·X=nnls(A，b):返回方程Ax=b的最小二乘解，方程的求解過程被限制在x≥0的條件下。

·X=nnls(A，b，TOL):指定誤差TOL來求解，TOL的默認值為TOL=max(size(A))*norm(A，1)*eps，矩陣的-1范數(shù)越大，求解的誤差越大。

·［X，W］=nnls(A，b):當x(i)=0時，w(i)<0；當x(i)>0時，w(i)≈0，同時返回一個雙向量w。

【例6-10】求方程組的非負最小二乘解。

A=［3.4336 -0.5238

0.6710

-0.5238 3.2833

-0.7302

0.6710 -0.7302

4.0261］；％輸入矩陣A

b=［-1.0000 1.5000

2.5000］；［X，W］=nnls(A，b)％求方程的非負最小二乘解

X=

0

0.6563

0.6998

W=

-3.6820

-0.0000

-0.0000

x1=A＼b ％用除法解方程

x1=

-0.3569

0.5744

0.7846

A*X-b ％驗證非負最小二乘解

ans=

1.1258

0.1437

-0.616

A*x1-b ％驗證除法解

ans=

1.0e-015*

-0.2220

0.4441

06.1.7方程解的精度人們在用計算機計算一個問題時，最關(guān)心的是該問題的數(shù)值解。而所得的解是否可靠，除非計算所用的是一些特殊的整數(shù)或有理數(shù)，否則計算中的圓整誤差、截斷誤差等都不可避免。在MATLAB中，數(shù)與數(shù)之間的最小分辨率用eps表示。表達式中任何數(shù)的相對誤差都不可能小于eps。此外，在程序執(zhí)行過程中，都不可避免地使這個最小誤差比放大。范數(shù)和條件數(shù)對求解過程中誤差放大現(xiàn)象有重要的定量描述，在這里僅重點說明方程的精度問題。在數(shù)值分析中，方程精度的估計可利用解的相對誤差來進行，其公式為其中，K是矩陣的條件數(shù)，在MATLAB中的命令為cond()。由于eps是機器精度，所以可以用K*eps的大小粗略判斷所得的方程解是否可靠。

【例6-11】對方程Ax=b的近似解和精確解進行比較，其中A是Hilbert矩陣。

N=［6，14］；％計算的矩陣階數(shù)

fori=1:length(N)

n=N(i) ％所要計算的矩陣階數(shù)放入n中

A=hilb(n)； ％產(chǎn)生n階Hilbert矩陣

Ai=invhilb(n)； ％產(chǎn)生完全準確的n階逆Hilbert矩陣

b=ones(n，1)； ％產(chǎn)生n階全是1的向量

x_app=A＼b； ％利用左除來求近似解

x_exa=Ai*b； ％利用準確求逆來求方程的準確解

fdb=norm(A*x_app-b)；

fb=norm(b)；

fdx=norm(x_app-x_exa)；

fx=norm(x_app)；

err_actual(i)=fdx/fx； ％實際的相對誤差

K=cond(A)；

err_app(i)=K*eps； ％最大可能的近似相對誤差

err_max(i)=K*fdb/fb； ％最大可能的相對誤差

end

disp(′Hilbert矩陣階數(shù)′)，disp(N)

formatshorte

disp(′實際的相對誤差′)，disp(err_actual)，disp(′′)

disp(′最大可能的近似相對誤差′)，disp(err_app)，disp(′′)

disp(′最大可能的相對誤差′)，disp(err_max)，disp(′′)運行結(jié)果如下:

Hilbert矩陣階數(shù)

6 14實際的相對誤差

5.0339e-011 3.9641e+000最大可能的近似相對誤差

3.3198e-009 9.4718e+001最大可能的相對誤差

6.0095e-007 6.6239e+009

從該例可以看出，14階的Hilbert矩陣是嚴重錯誤的。6.1.8用函數(shù)零點求方程的解

在MATLAB中，用函數(shù)零點法求方程或方程組解的函數(shù)命令有兩個，即fzero和fsolve。在此，首先需對方程和方程組進行轉(zhuǎn)化，比如將方程f(x)=g(x)轉(zhuǎn)化為F(x)=f(x)-g(x)=0，方程組也是如此；然后將函數(shù)F(x)寫成MATLAB的m函數(shù)，以便在fzero和fsolve命令中調(diào)用。本小節(jié)重點講述其使用方法。

1.一元方程轉(zhuǎn)化的函數(shù)零點求法在一元方程中，任意函數(shù)方程F(x)=0可能有零點，也可能沒有零點；可能有一個零點，也可能有多個零點，很難有一個通用的解法。一般來說，其求解的過程為:先猜測一個初始零點或者該零點大概所在的區(qū)間，然后通過計算使猜測值不斷精確化，或使猜測區(qū)間不斷收縮，直至達到預先指定的精度為止。在猜測一個初始的零點時，一般用MATLAB的作圖命令來獲取初始近似解。其具體步驟為:先確定一個零點可能存在的自變量區(qū)間，然后利用fplot命令繪出f(x)在該區(qū)間的圖形，用眼觀察F(x)與橫軸的交點坐標，或者更細些用zoom命令對交點局部放大來讀數(shù)，或借助ginput命令來獲得更精確的交點坐標。

1)利用fzero函數(shù)求一元函數(shù)零點命令fzero的調(diào)用格式為:

(1)x=fzero(fun，x0):求一元函數(shù)零點命令的最簡形式；

(2)［x，fval，exitflag］=fzero(fun，x0，options，P1，P2，…):從MATLAB5.3版本起，求一元函數(shù)零點命令最完整的格式。其中參數(shù)x0是初始猜測的零點，同時也是預定搜索零點的大致位置。它可以是標量或是二元向量，當x0是標量時，該命令將在它的兩側(cè)尋找一個與之最靠近的零點:當x0是一個二元向量［a，b］時，該命令將在區(qū)間［a，b］內(nèi)尋找零點。options是優(yōu)化迭代所采用的參數(shù)選項。它是MATLAB設計優(yōu)化程序fzero、fsolve、fminbnd、fminsearch和fminunc時都需要的一個“模塊”，采用結(jié)構(gòu)數(shù)組存放優(yōu)化參數(shù)。如果沒有優(yōu)化參數(shù)要設置，則可以在optons的位置上用“［］”作為占位符。在此命令中，options的缺省設置可以用命令options=optimset(′fzero′)獲得。它有兩個域:Display和TolX。其中:

options.Display:顯示設置，有三個選取值［off|iter|{final}］；

options.TolX:自變量計算的終止誤差，可選?。踥psitivescalar］，缺省值為2.2204e-016。若要修改缺省值，可用options.TolX=0.001或用options=optim-set(′TolX′，0.001)來改動。

P1和P2是向函數(shù)fun傳遞的附加參數(shù)。它的具體取名和函數(shù)fun中一致。

x是輸出參數(shù)，為所求的零點自變量值。

fval是輸出參數(shù)，為函數(shù)fun在x處的值。

exitflag用于描述函數(shù)fun的退出情況，表明程序的終止條件。若exitflag>0，則表示找到函數(shù)零點后退出；若exitflag<0，則表示沒有找到零點或在搜索過程中遇到了無窮大的函數(shù)值。

【例6-12】試用fzero命令求解函數(shù)f(x)=x4-4x-5的零點。在MATLAB內(nèi)，建立M文件C6L12。

(1)建立函數(shù)f(x)M文件。

funtiony=fun1(x)%例6-12用M函數(shù)文件

y=x.^4-4*x-5；

(2)建立水平橫軸的M文件。

functiony=fun2(x) %例6-12用M函數(shù)文件(繪出水平軸的函數(shù))

y=0；

(3)用作圖法估計函數(shù)的零點位置。

fplot(′fun1′，［-5，5］)

holdon

fplot(′fun2′，［-5，5］，′r′)程序的運行結(jié)果如圖6-1所示。圖6-1函數(shù)零點估計圖

2)用zoom和ginput命令獲得零點的初始近似值在MATLAB的繪圖窗口中，可實現(xiàn)zoom和ginput命令進行圖形和鼠標的交互操作，或直接在程序中輸入下列命令，便可得到如圖6-2所示的局部放大圖及鼠標操作線。

zoomon％局部放大命令［tt，yy］=ginput(2) ％用鼠標獲取2個零點猜測值

zoomoff ％恢復原來圖形大小顯示所得零點初始猜測值，結(jié)果為:

tt=

-1.0138

1.9124

yy=

0.5848

0.5848圖6-2局部放大圖和利用鼠標取值圖

3)用函數(shù)fzero求函數(shù)的精確零點［x，fval，exitflag］=fzero(′fun1′，tt(1)，［］)％靠近tt(1)點處的精確零點［x，fval，exitflag］=fzero(′fun1′，tt(2)，［］) ％靠近tt(2)處的精確零點

結(jié)果為:

Zerofoundintheinterval:［-0.98515，-1.0138］

x=

-1

fval=

0

fxitflag=

1

Zerofoundintheinterval:［1.8584，1.9124］

x=

1.8812

fval=

-6.2172e-015

exitflag=

1

2.多元函數(shù)的零點解非線性方程組

對于采用多元函數(shù)的零點解非線性方程組，MATLAB提供有fsolve命令來完成這類工作。一般來說，多元函數(shù)的零點問題更難解決。而一旦知道零點的大致位置，就可用數(shù)值方法搜索到精確零點。同一元函數(shù)用函數(shù)曲線與橫軸的交點來獲得零點的大致位置一樣，二元函數(shù)的零等位線有助于問題的解決。但更高維的初始零點則很難通過圖形獲得。有了初始零點后，即可借助函數(shù)fsolve命令來求零點的精確解。其調(diào)用格式為:

·x=fsolve(fun，x0):解非線性方程組最簡單的調(diào)用格式。該式中除兩個輸入?yún)?shù)外，其余輸入輸出參數(shù)都可以缺省。

·［x，fval，exitflag，output，jacob］=fsolve(fun，x0，options，P1，P2…):解非線性方程組最完整的調(diào)用格式。各參數(shù)的含義如下:

·x0：零點數(shù)是猜測值的向量。

·options:fsolve的優(yōu)化迭代所采用參數(shù)的結(jié)構(gòu)數(shù)組。它的缺省值可用options=optimset(′fsolve′)獲得。它含有5個域:Display、MaxFunEvals、Maxlter、TolFun和TolX。其關(guān)鍵參數(shù)為options.Display。

options.Display:顯示設置，有3個選取值［off|iter|{final}］；

options.MaxfunEvals:允許函數(shù)計算的最多次數(shù)，可選?。踦ositiveinteger］，缺省值是自變量的100倍；

options.Maxlter:允許的最多迭代次數(shù)，可選?。踦ositiveinteger］，缺省值為400；

options.TolFun:函數(shù)計算的終止誤差，可選取［positivescalar］，缺省值為1.0000e-006；

options.TolX:自變量計算的終止誤差，可選?。踦ositivescalar］，缺省值為1.0000e-006。

·P1和P2:向函數(shù)fun傳遞的參數(shù)。

·x和fval:輸出參數(shù)，分別是所求零點的自變量值和函數(shù)值?！xitflag:若exitflag>0，則表示找到零點后退出；若exitflag<0，則表示沒有找到零點或在搜索過程中遇到了無窮大的函數(shù)值。

·output:輸出本命令所用的計算方法、迭代次數(shù)等信息，它是結(jié)構(gòu)數(shù)組。

·jacob:函數(shù)在x處的jacobian。注意:上兩例中輸入?yún)?shù)fun有三種編寫方式:字符串方式、內(nèi)聯(lián)函數(shù)方式和M函數(shù)文件方式。

【例6-13】試求方程組的解。

(1)觀察量函數(shù)0等位線的交點情況。

x=-1:0.5:1；

y=x；［X，Y］=meshgrid(x，y)；％產(chǎn)生x-y平面上的網(wǎng)點坐標

fun1=sin(X)+Y； ％函數(shù)fun1

fun2=X+6*Y ； ％函數(shù)fun2

v=［-0.2，0，0.2］ ； ％指定三個等位線，是為了更可靠地判斷0等位線的位置

contour(X，Y，fun1，v) ％繪出fun1的三條等位線

holdon

contour(X，Y，fun2，v) ％繪出fun2的三條等位線

holdoff

(2)從圖6-3中用鼠標獲取零點的初始近似值。［x0，y0］=ginput(1)； ％用鼠標在圖形上取一個點的坐標(三線組中間那條線)

disp(［x0，y0］) ％顯示初始零點猜測值

0.0000 -0.0058圖6-3兩個二元函數(shù)0等位線的交點圖有了初始零點后，便可利用fsolve函數(shù)命令求精確解。

fun=′［sin(x(1))+x(2)，x(1)+6*x(2)］′%用字符串表達式形式命令。注意自變量必須寫成x(1)和x(2)

XXYY=fsolve(fun，［x0(1)，y0(1)］)%解此非線性方程組

disp(XXYY)

XXYY=

1.0e-016*

0.169 -0.0347

1.0e-016*

0.1691-0.0347檢驗:

FF1=sin(XXYY(1))+XXYY(2) ％檢驗方程1

FF2=XXYY(1)+6*XXYY(2) ％檢驗方程2

FF1=

1.3444e-017

FF2=

-3.9031e-018上面的FF1和FF2表明，所求的零點相當準確。說明:在使用函數(shù)fsolve命令時，fun也可是內(nèi)聯(lián)函數(shù)的形式，此時

fun=inline(′［sin(x(1))+x(2)，x(1)+6*x(2)］′，′x′)；%′x′項必須有

XXYY=fsolve(fun，［x0(1)，y0(1)］)

fun函數(shù)也可用M函數(shù)文件的形式，此時應先用fun.m表示被解函數(shù)(并在搜索路徑上):

functionyy=fun(x)

yy(1)=sin(x(1))+x(2)′

yy(2)=x(1)+6*x(2)；

XXYY=fsolve(′fun′，［x0(1)，y0(1)］)6.1.9符號方程及方程組的求解

在工程教學或工程實踐中，會遇到一些帶有符號的方程，或某些問題需求它的符號解，這類問題就是符號方程和方程組的求解問題。一般在實際中遇到的代數(shù)方程組包括線性方程組、非線性方程組和超越方程組。在MATLAB中，用于求解代數(shù)方程組的命令有l(wèi)insolve命令和solve命令。

1.線性方程組的符號解

矩陣計算是求解線性方程組最簡單有效的方法。從MATLAB5.x版本起，不管數(shù)據(jù)對象是數(shù)值還是符號，實現(xiàn)矩陣運算的指令幾乎相同。因而，求線性方程組符號解時，可以套用數(shù)值解的命令的編寫方法進行。而在MATLAB中，函數(shù)命令linsolve專門用來求解線性方程組，其使用等同于數(shù)值計算方法。對方程A*X=B，函數(shù)命令linsolve的調(diào)用格式為:

X=linsolve(A，B)等同于

X=sym(A)＼sym(B)矩陣A必須至少是行滿秩的。當A的列數(shù)大于行數(shù)時，將給出解不唯一的警告提示。使用該格式可得到方程組的特解X。方程組的通解XX為

XX==X+′k′*null(A)其中，k是任意常數(shù)；null命令將求出A的“零空間”的基。

【例6-14】求給定線性方程組的解。

A=sym(′［12/31/2；351；121］′) %輸入矩陣A

A=［1，2/3，1/2］［3，5，1］［1，2，1］

b=sym(′［12；1/33；11/7］′)；

X=linsolve(A，b)%計算線性方程的解

X=［23/39，283/91］［-8/13，-102/91］［64/39，-66/91］

【例6-15】求欠定方程組。

A=sym(′［12/31/21；3511；1211］′)％輸入矩陣A

A=［1，2/3，1/2，1］［3，5，1，1］［1，2，1，1］

b=sym(′［1；1/3；1］′)；

X=linsolve(A，b) %特解

X=［-1/3］［0］［0］［4/3］

XX=X+′k′null(A)%通解

XX=［-1/3-3/2*k］［k］［-8/3*k］［4/3+13/6*k］說明:欠定方程組的解不唯一，而MATLAB總是給出一個特解，并且解向量中含非零元素最少。

2.一般代數(shù)方程組的符號解

slove命令可以解一般代數(shù)方程，包括線性方程、非線性方程和超越方程。當方程不存在符號解，且又無其他自由參數(shù)時，函數(shù)solve將給出數(shù)值解。其命令調(diào)用格式為:

·solve(′eqn1′，′eqn2′，…，′eqnN′):對N個方程的默認變量求解。

·solve(′eqn1′，′eqn2′，…，′eqnN′，′var1，var2，…，varN′):對N個方程的var1，var2，…，varN變量求解，但該式要注意變量的英文字母順序，并且在變量前不可有空格。

·S=solve(′eqn1′，′eqn2′，…，′eqnN′，′var1′，′var2′，…，′varN′):對N個方程的′var1′，′var2′，…，′varN′變量求解。

·［x1，x2，…，xn］=solve(′eqn1′，′eqn2′，…，′eqnN′，′var1′，′var2′，…，′varN′):對變量

var1，var2，…，varN求解，并且求解的結(jié)果分別賦給x1，x2，…，xn。此式中，MATLAB是按照變量var1，var2，…，varN在英文字母中的順序賦值給x1，x2，…，xn的。提示:′eqn1′，′eqn2′，…，′eqnN′是字符串表達的方程，或是字符串表達式；′eqnN′，′var1′，′var2′，…，′varN′是字符串表達的求解變量名。第三個命令中，S是一個結(jié)構(gòu)數(shù)組，如果要顯示結(jié)果，則必須使用S.var1，S.var2，…，S.varN的引用形式。在得不到“封閉性解析解”時，將給出數(shù)值解。

【例6-16】求非線性方程組的解。［x，y，z］=solve(′x^2+sqrt(2)*x+2=0′，′x+3*z=4′，′y*z=-1′，′x′，′y′，′z′)

x=［(-1/2+1/2*i*3^(1/2))*2^(1/2)］［(-1/2-1/2*i*3^(1/2))*2^(1/2)］

y=［-51/73+3/73*i*3^(1/2)-27/146*(-1/2+1/2*i*3^(1/2))*2^(1/2)-3/146*2^(1/2)］［-51/73-3/73*i*3^(1/2)-27/146*(-1/2-1/2*i*3^(1/2))*2^(1/2)-3/146*2^(1/2)］

z=［-1/3*(-1/2+1/2*i*3^(1/2))*2^(1/2)+4/3］［-1/3*(-1/2-1/2*i*3^(1/2))*2^(1/2)+4/3］

S=solve(′x^2+sqrt(2)*x+2=0′，′x+3*z=4′，′y*z=-1′，′x′，′y′，′z′)

S=

x:［2x1sym］

y:［2x1sym］

z:［2x1sym］

disp(′S.x′)，disp(S.x)％顯示結(jié)構(gòu)數(shù)組S中x的內(nèi)容

S.x［(-1/2+1/2*i*3^(1/2))*2^(1/2)］［(-1/2-1/2*i*3^(1/2))*2^(1/2)］

disp(′S.y′)，disp(S.y)％顯示結(jié)構(gòu)數(shù)組S中y的內(nèi)容

S.y［-51/73+3/73*i*3^(1/2)-27/146*(-1/2+1/2*i*3^(1/2))*2^(1/2)-3/146*2^(1/2)］［-51/73-3/73*i*3^(1/2)-27/146*(-1/2-1/2*i*3^(1/2))*2^(1/2)-3/146*2^(1/2)］

disp(′S.z′)，disp(S.z)％顯示結(jié)構(gòu)數(shù)組S中z的內(nèi)容

S.z［-1/3*(-1/2+1/2*i*3^(1/2))*2^(1/2)+4/3］［-1/3*(-1/2-1/2*i*3^(1/2))*2^(1/2)+4/3］［x，y］=solve(′x^x=2′，′x/y=3′)

x=

log(2)/lambertw(log(2))

y=

1/3*log(2)/lambertw(log(2))上述符號解中的lambertw(ω)代表ω函數(shù)ω(x)eω(x)=x的解。有關(guān)詳細內(nèi)容，請用mhelp命令訪問MAPLE庫的lambertw條目。

【例6-17】解超越方程組的解。6.1.10矩陣的特征值和特征向量

特征值和特征向量是線性代數(shù)中非常重要的概念，在實際的工程應用和求解數(shù)學問題中占有非常重要的地位。比如工程中的振動問題、穩(wěn)定問題等，從數(shù)學關(guān)系上常常歸結(jié)為矩陣的特征值和特征向量的求解問題。在求解微分方程組以及簡化矩陣時都要用到特征理論。而在與振動有關(guān)的各學科中，特征值和特征向量的問題都有著廣泛的應用。矩陣A與向量x相乘，即表示矩陣對向量的變換(Transformation)。對于任何一個矩陣來說，總存在一些特殊的向量，在變換的作用下，向量的方向不變，僅是長短發(fā)生變化。這種向量就是所謂的特征向量(Eigenvector)。它滿足方程:Ax=λx其中，A是n×n的方陣，λ稱為特征值(Eigenvalue)，是個標量。x是對應λ的一個長度為n的列向量。該方程就稱為特征方程(EigenvalueEquation)。上式的廣義特征方程是:Ax=λBx。此式中，A和B都是n×n的矩陣，λ是標量。當B＝1時，就退化為Ax=λx。而廣義特征值問題就是方程Ax=λBx的非平凡解問題。一些有缺陷的矩陣不能對角化，不能進行特征值分解。由A的特征值構(gòu)成的對角矩陣以及由對應的特征向量構(gòu)成的矩陣V的各列，必須滿足:AX＝VD。如果V是非奇異的，則這就是矩陣A的特征值分解。

1.矩陣的數(shù)值特征值和數(shù)值特征向量在數(shù)學教科書上，最常見的求解特征方程Ax=λx的方法是:先根據(jù)|A-λxi|=0求特征值λxi(i＝1，2，…，n)，然后由Ax=λx求對應的特征向量xi。而在MATLAB中，其計算特征值和特征向量的算法取自EISPACK程序庫，相應的計算命令較為簡單，調(diào)用格式如下:

·D＝eig(A):僅計算A的特征值組成的列向量。

·［V，D］=eig(A):生成由矩陣A的特征值、矩陣D和特征向量構(gòu)成的矩陣V，使A*V=V*D。

·［V，D］=eig(A，′nobalance′):計算時不采用預先平衡。通常，預先平衡增加了特征值和特征向量的計算精度，然而，當一個矩陣包含有與截斷誤差數(shù)量級相差不遠的元素時，平衡過程有可能將它們放大，從而導致錯誤的特征值。該指令可使精度增加。

·D＝eig(A，B):如果A和B是方陣，則生成廣義特征值D?！ぃ踁，D］=eig(A，B):計算廣義特征值矩陣D和廣義特征向量矩陣V，使得A*V=B*V*D

【例6-18】計算矩陣A的特征值和特征向量。

A=［-110；-430；102］；

B=［123；139；402］；

D＝eig(A)％計算矩陣A的特征值

D＝

2

1

1［V，D］＝eig(A)

V＝

0

0.4082

-0.4082

0

0.8165

-0.8165

1.0000

-0.4082

0.4082

D＝

2

0

0

0

1

1

0

0

1［V，D］＝eig(A，B)

V＝

0.5940-0.0854i

-0.0356+0.5465i

-0.5983-0.0466i

-0.2068+0.3408

-0.0538+0.8255i

0.2285+0.3267i

-0.4748+0.5054i

0.0081-0.1251i

0.4410-0.5352i

D＝

0.2524-0.5295i

0

0

0

0.1531-0.0000i

0

0

0

0.2524+0.5292i

A*V

ans＝

-0.8008+0.4262i

-0.0182+0.2790i

0.8296+0.3733i

-2.9965+1.3640i

-0.0189+0.2905i

3.0789+1.1663i

-0.3556+1.0962i

-0.0193+0.2964i

0.2836-1.1170i

B*V*D

ans＝

-0.8008+0.4262i

-0.0182+0.2790i

0.8269+0.3733i

-2.9965+1.3640i

-0.0189+0.2905i

3.0789+1.1663i

-0.3556+1.0962i

-0.0193+0.2964i

0.2836-1.1170i

2.符號特征值和特征向量

在MATLAB中，也是使用函數(shù)eig計算方陣A的符號特征值和特征向量。其調(diào)用格式為:

Lambda＝eig(A)［V，D］=eig(A)計算任意精度的矩陣特征值和特征向量的調(diào)用格式為:

Lambds＝eig(vpa(A))［V，D］=eig(vpa(A))上式中，各參數(shù)的具體含義與數(shù)值特征值和特征向量相同，在此不再贅述。

【例6-19】計算符號特征值和特征向量。

A＝［8/91/21/3；1/21/31/4；1/31/41/5］；

A＝sym(A)%轉(zhuǎn)化為符號矩陣

A＝［8/9，1/2，1/3］［1/2，1/3，1/4］［1/3，1/4，1/5］［V，D］＝eig(A)

V＝［28/153+2/153*12589^(1/2)，28/153-2/153*12589^(1/2)，1］［1，1，-4］［292/255-1/255*12589^(1/2)，299/255+1/255*12589^(1/2)，10/3］

D＝［32/45+1/180*12589^(1/2)，0，0］［0，32/45-1/180*12589^(1/2)，0］［0，0，0］

eig(vpa(A))

ans＝［-.31796273547608369503007e-31］［.87773816609932252988056841486003e-1］［1.3344484056122899692341653807363］［V，D］=eig(vpa(A))

V＝［.18860838403857944292978827467058，.80318501188318836927227766194002，-.56508469644519509785929796275283］［.75443353615431777171915309868189，.48687247435830155621160586316434，.44021044199100581223229247499660］［.62869461346193147643262758223468，.34329146566500535667224373721531，.69777793932276550403654035882202］

D=［-.33889316547608368542146e-31，0，0］［0，1.3344484056122899692341653907363，0］［0，0，.087773816609932252988056841485995e-1］6.1.11矩陣的對角化和其他矩陣函數(shù)矩陣運算經(jīng)過矩陣對角化后變得十分簡單，特別對于維數(shù)較多的矩陣更是如此。矩陣對角化在實際的工程教學和實踐中有著非常廣泛的用途。幾乎所有的對角化都基于特征值和特征向量的求解，特征值和特征向量的求解目的也是為了對角化。

1.矩陣的PAP對角化由線性代數(shù)可知，對于任意可對角化的矩陣A，都存在一個可逆矩陣P，使得P-1AP為對角陣，并且對角陣的對角線元素為矩陣A的特征值。在MATLAB中，求出矩陣特征值和特征向量D和V后，即可滿足上述條件(P＝V)。

【例6-20】矩陣的PAP對角化。

A＝［32-1；-2-22；361］；［P，D］＝eig(A)

P＝

0.8890

0.2673

-0.0531

-0.2540

-0.5345

0.4677

0.3810

0.8081

0.8823

D＝

2.0000

0

0

0

-4.0000

0

0

02.000

0

Inv(P)*A*P

ans＝

2.0000

0.0000

-0.0000

-0.0000

-4.0000

0.0000

-0.0000

-0.0000

2.0000

2.實對稱矩陣的QRQ對角化

實對稱矩陣A都可對角化，并且都存在正交矩陣Q，使得Q-1AQ即(Q′AQ)為對角陣。對角陣的對角線元素均為矩陣A的特征值。由線性代數(shù)的知識可知，對于實對稱矩陣A，特征值分解函數(shù)eig(A)返回的特征向量矩陣就是正交矩陣。

【例6-21】實矩陣的QRQ對角化。

A=［2

2

-2；2

5

-；-2

-4

5］；［Q，D］＝eig(A)

Q＝

0.8944

0.3333

-0.2981

-0.4472

0.6667

-0.5963

0-0.6667-0.7454

D＝

1.0000

0

0

0

10.0000

0

0

0

1.0000

Q*A*Q

ans＝

1.0000

0

0

0

10.0000

0

-0.0000

-0.0000

1.0000

inv(Q)*A*Q

ans＝

1.0000

-0.00000.0000

-0.0000

10.0000-0.0000

0.0000

-0.0000

1.0000

3.約旦(Jordan)標準型矩陣當用相似變換對角化矩陣時，就會產(chǎn)生約旦標準型矩陣。即給定矩陣A，尋找非奇異矩陣V使inv(V)*A*V(更簡潔地說，是使J=V＼A*V)盡可能地接近對角陣。約旦標準型是特征值矩陣，它的變換矩陣的每一列就是特征向量。對于有重特征值的非對稱矩陣，不能進行對角化。約旦標準型對角線上的元素是特征值，但一些對角線以上的元素是1，而不是0。在MATLB中，計算約旦標準型的函數(shù)是jordan，其調(diào)用格式為:

J=jordan［V，J］=Jordan(A)式中:J是約旦標準型，V是相似變換矩陣，使得V＼A*V=J。V的列是A的一般化特征向量。

【例6-22】求矩陣的約旦標準型。

A=［111；021；131］；［V，J］=Jordan(A)

V＝

0.7500

0.0947

0.1553

-0.2500

0.0947

0.1553

0.2500

0.1479

-0.3979

J＝

1.0000

0

0

0

3.5616

0

0

0

-0.5616

d＝eig(A)