《MATLAB程序設(shè)計(jì)基礎(chǔ)教程》課件第7章

第7章數(shù)值計(jì)算與分析

7.1MATLAB多項(xiàng)式7.2有理多項(xiàng)式的運(yùn)算7.3多項(xiàng)式估值與擬合7.4數(shù)據(jù)插值7.5數(shù)值分析7.6代數(shù)方程組求解7.7微分方程的數(shù)值解

7.1MATLAB多項(xiàng)式

7.1.1概述

多項(xiàng)式在數(shù)學(xué)中有著極為重要的作用，同時(shí)多項(xiàng)式的運(yùn)算也是工程和應(yīng)用中經(jīng)常遇到的問(wèn)題。曲線擬合、插值運(yùn)算在具有數(shù)據(jù)統(tǒng)計(jì)、信號(hào)處理和圖像處理等領(lǐng)域應(yīng)用十分廣泛，是工程中經(jīng)常要用到的技術(shù)之一。每當(dāng)難以對(duì)一個(gè)函數(shù)進(jìn)行積分、微分或者解析上確定一些特殊的值時(shí)，就可以借助計(jì)算機(jī)進(jìn)行數(shù)值分析。

MATLAB提供了一些專門用于處理多項(xiàng)式的函數(shù)，包括多項(xiàng)式求根、多項(xiàng)式的四則運(yùn)算及多項(xiàng)式的微積分等，如表7-1所示。

1．多項(xiàng)式與行向量

在MATLAB中，多項(xiàng)式可用一個(gè)行向量表示，向量中的元素為該多項(xiàng)式的系數(shù)，按照降序排列。如多項(xiàng)式p(x)?=?9x3?+?7x2?+?4x?+?3可以表示為向量p=[9743]。

用戶可以先用創(chuàng)建向量的方式創(chuàng)建一個(gè)多項(xiàng)式，然后按照降序排列其表示為多項(xiàng)式，也可用ployzstr(p)函數(shù)將向量p轉(zhuǎn)換為多項(xiàng)式。

2．多項(xiàng)式的運(yùn)算

由于多項(xiàng)式是利用向量來(lái)表示的，因此多項(xiàng)式的四則運(yùn)算可以轉(zhuǎn)化為向量的運(yùn)算。

因?yàn)槎囗?xiàng)式的加減為對(duì)應(yīng)項(xiàng)系數(shù)的加減，所以可以通過(guò)向量的加減來(lái)實(shí)現(xiàn)。但是在向量的加減中兩個(gè)向量需要有相同的長(zhǎng)度，因此在進(jìn)行多項(xiàng)式加減時(shí)，需要將短的向量前面補(bǔ)0。

多項(xiàng)式的乘法實(shí)際上是多項(xiàng)式系數(shù)向量之間的卷積運(yùn)算，可以通過(guò)MATLAB中的卷積函數(shù)conv()來(lái)完成。

多項(xiàng)式的除法為乘法的逆運(yùn)算，可以通過(guò)反卷積函數(shù)deconv()來(lái)實(shí)現(xiàn)。除多項(xiàng)式的四則運(yùn)算外，MATLAB還提供了多項(xiàng)式的一些其他運(yùn)算。7.1.2多項(xiàng)式與根

1．用roots()函數(shù)求多項(xiàng)式的根

找出多項(xiàng)式(polynomial)的根，即多項(xiàng)式為零的值，這可能是許多學(xué)科共同的問(wèn)題。MATLAB求解這個(gè)問(wèn)題，并提供其他的多項(xiàng)式操作工具。

在MATLAB中，首先用一個(gè)行向量表示多項(xiàng)式的系數(shù)，該系數(shù)向量按降序排列。

用函數(shù)roots()找出多項(xiàng)式的根的函數(shù)用法如下：

r=roots(c)：返回一個(gè)列向量，其元素就是多項(xiàng)式c的根。行向量c包含多項(xiàng)式按降序指數(shù)排列的系數(shù)，如果c有n+1個(gè)元素，則表示多項(xiàng)式為(7.1.1)

例7-1-1

求多項(xiàng)式x4-12x3?+?25x?+?116的根。

解輸入多項(xiàng)式系數(shù)向量：

>>C=[1-12025116];注意，向量中必須包括系數(shù)為零的項(xiàng)。

>>r=roots(C)

r=

11.7473

2.7028

-1.2251+1.4672i

-1.2251-1.4672i

2．用poly()函數(shù)根重組多項(xiàng)式

poly()函數(shù)與roots()函數(shù)是一對(duì)逆運(yùn)算函數(shù)。使用poly()函數(shù)，可用根重組多項(xiàng)式。函數(shù)用法如下：

(1)?p=poly(r)：r是一個(gè)多項(xiàng)式的根向量，該函數(shù)返回多項(xiàng)式系數(shù)的行向量。例如：

>>pp=poly(r)

pp=

1.0000-12.0000-0.000025.0000116.0000

(2)

p=poly(A)：A是一個(gè)n

×

n矩陣，返回一個(gè)n+1個(gè)元素的行向量，其元素是多項(xiàng)式的系數(shù)。例如：

>>A=[123

456

780]

>>p=poly(A)

p=

1.0000-6.0000-72.0000-27.0000該語(yǔ)句找出多項(xiàng)式的根，在運(yùn)算中實(shí)際是求A的本征值(eigenvalues)，因此也可以使用eig()函數(shù)求A的本征值。

因?yàn)镸ATLAB無(wú)隙地處理復(fù)數(shù)，當(dāng)用根重組多項(xiàng)式時(shí)，如果一些根有虛部，由于截?cái)嗨a(chǎn)生的誤差，則poly()的結(jié)果有一些小的虛部，這是很普通的。要消除虛假的虛部，只要使用函數(shù)real抽取實(shí)部即可。7.1.3卷積運(yùn)算與多項(xiàng)式乘法

conv()函數(shù)用于卷積與多項(xiàng)式乘法。函數(shù)的定義和用法如下：

w=conv(u,v)：計(jì)算向量u和v的卷積。在代數(shù)計(jì)算中，多項(xiàng)式的乘法與此相同，多項(xiàng)式的系數(shù)就是u和v的元素。

定義：若m=length(u)、n=length(v)，則w是長(zhǎng)度為m?+?n-1的向量。第k個(gè)元素是：總和w(k)是所有u(j)和v(k?+?1-j)的乘積的和，j的取值范圍是j=max(1,k+1-n):min(k,m)。當(dāng)m=n時(shí)，即兩向量長(zhǎng)度相同，則w是長(zhǎng)度為2n?-1的向量，結(jié)果是：

w(1)=u(1)*v(1)

w(2)=u(1)*v(2)+u(2)*v(1)

w(3)=u(1)*v(3)+u(2)*v(2)+u(3)*v(1)

…

w(n)=u(1)*v(n)?+?u(2)*v(n?-1)?+?…?+?u(n)*v(1)

…

w(2*n-1)=u(n)*v(n)

conv()函數(shù)支持多項(xiàng)式乘法，即執(zhí)行兩個(gè)數(shù)組的卷積，兩個(gè)以上的多項(xiàng)式的乘法需要重復(fù)使用conv。7.1.4反卷積運(yùn)算與多項(xiàng)式除法

函數(shù)deconv()用于反卷積與多項(xiàng)式除法。函數(shù)的定義和用法如下：

[q,r]=deconv(v,u)：使用于長(zhǎng)除法，輸出向量u和v的反卷積。返回的商是向量q，余數(shù)是向量r。例如：

>>u=[1234]

v=[102030]

卷積結(jié)果：

>>x=conv(u,v)

x=1040100160170120反卷積結(jié)果：

>>[q,r]=deconv(x,u)

q=102030

r=000000

>>[q,r]=deconv(x,v)

q=1234

r=000000在一些特殊情況，一個(gè)多項(xiàng)式需要除以另一個(gè)多項(xiàng)式。如果u和v是多項(xiàng)式系數(shù)向量，卷積等于多項(xiàng)式乘，則反卷積等于多項(xiàng)式除，v被u除，商是q，余數(shù)是r。

在MATLAB中，多項(xiàng)式除由函數(shù)deconv完成。由上面的多項(xiàng)式a和c，得出b：

>>[b,r]=deconv(c,a)

b=149-15

r=0000000

這個(gè)結(jié)果是c被a除，給出商多項(xiàng)式b和余數(shù)r，在現(xiàn)在情況下r是零，因?yàn)閎和a的乘積恰好是c。7.1.5多項(xiàng)式加法

對(duì)于多項(xiàng)式的加法，MATLAB沒(méi)有提供一個(gè)直接的函數(shù)。如果兩個(gè)多項(xiàng)式向量大小相同，則可使用標(biāo)準(zhǔn)的數(shù)組加法。例如，將例7-1-2中a(x)和b(x)相加，由于a?=?[1234]，

b?=?[149-15]，故

>>d=a+b

d=

2612-11

結(jié)果是

d(x)?=2x3?+?6x2?+?12x-11當(dāng)兩個(gè)多項(xiàng)式階次不同時(shí)，低階的多項(xiàng)式必須用首零填補(bǔ)，使其與高階多項(xiàng)式有同樣的階次。這里要求是首零而不是尾零，是因?yàn)橄嚓P(guān)的系數(shù)像x冪次一樣，必須整齊。

例如，考慮將例7-1-2的多項(xiàng)式c和d相加：

>>e=c+[000d]

e=

162021193-71

結(jié)果是：

e(x)?=x6?+?6x5?+?20x4?+?21x3?+?19x2?+?3x-717.1.6多項(xiàng)式求導(dǎo)數(shù)

函數(shù)polyder()用于多項(xiàng)式的求導(dǎo)。該函數(shù)可以用于求解一個(gè)多項(xiàng)式的導(dǎo)數(shù)、兩個(gè)多項(xiàng)式乘積的導(dǎo)數(shù)和兩個(gè)多項(xiàng)式商的導(dǎo)數(shù)。該函數(shù)的用法如下：

(1)?q=polyder(p)：該命令計(jì)算多項(xiàng)式p的導(dǎo)數(shù)。

(2)?c=polyder(a,b)：該命令實(shí)現(xiàn)多項(xiàng)式a、b的積的導(dǎo)數(shù)。

(3)?[q,d]=polyder(a,b)：該命令實(shí)現(xiàn)多項(xiàng)式a、b的商的導(dǎo)數(shù)，q、d為最后的結(jié)果。

例7-1-3

求多項(xiàng)式g(x)=x6+6x5+20x4+21x3+9x2+3x+11的導(dǎo)數(shù)。

解(1)輸入多項(xiàng)式系數(shù)向量：

>>g=[1620219311];

(2)求導(dǎo)：

>>h=polyder(g)

h=6308063183

結(jié)果是

h(x)?=6x5?+?30x4?+?80x3?+?63x2?+?18x?+?3例7-1-4用polyzstr()函數(shù)將向量轉(zhuǎn)換為多項(xiàng)式。

解程序如下：

a=[12345];

>>a=[12345];

>>poly2str(a,'x')

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

>>b=polyder(a)

b=4664

>>poly2str(b,'x')

ans=4x^3+6x^2+6x+4 7.2有理多項(xiàng)式的運(yùn)算

7.2.1使用residue()函數(shù)展開(kāi)部分分式

residue()函數(shù)是執(zhí)行部分分式展開(kāi)和多項(xiàng)式系數(shù)之間的轉(zhuǎn)換。語(yǔ)法如下：

[r,p,k]=residue(b,a)：返回部分分式展開(kāi)式中的極點(diǎn)p(poles)、相應(yīng)極點(diǎn)對(duì)應(yīng)的留數(shù)r(residues)和余項(xiàng)k(directterm)，a、b為多項(xiàng)式的系數(shù)。定義如下：設(shè)s的有理分式為

(7.2.1)式中，b和a分別表示F(s)分子和分母的系數(shù)，即

b=[b1b2…bm]

a=[a1a2…an]

(1)如果沒(méi)有多重根，將按下式給出F(s)部分分式展開(kāi)式中的留數(shù)、極點(diǎn)和余項(xiàng)：(7.2.2)式中，n=length(r)=length(p)=length(a)-1。①k(s)是余項(xiàng)，如果length(b)<length(a)，則k返回一個(gè)空向量，余項(xiàng)k為零。

例7-2-1試將下列函數(shù)展開(kāi)成部分分式

解對(duì)于該函數(shù)，有

>>num=[0

1

4

6];

>>den=[1

3

3

1];

>>[r,p,k]=residue(num,den)

將得到如下結(jié)果：

r=

1.0000

2.0000

3.0000

p=

-1.0000

-1.0000

-1.0000

k=

[

]０所以可得注意，本例的余項(xiàng)k為零。②否則，余項(xiàng)k不為零，其長(zhǎng)度為length(k)?=length(b)?-length(a)?+?1。例7-2-2試求下列函數(shù)的部分分式展開(kāi)式

解由此函數(shù)得

>>num=[1

11

39

52

26];

>>den=[1

10

35

50

24];

>>[r,p,k]=residue(num,den)

于是，得到下列結(jié)果：

r=

1.0000

2.5000

-3.0000

0.5000

p=

-4.0000

-3.0000

-2.0000

-1.0000

k=1則得

(2)如果F(s)中含多重極點(diǎn)，即p(j)?=?…?=?p(j?+?m-1)，則部分分式展開(kāi)式將包括下列各項(xiàng)(7.2.3)式中，pj為一個(gè)m重極點(diǎn)。7.2.2residue()函數(shù)的逆運(yùn)算

residue()函數(shù)也可執(zhí)行逆運(yùn)算[b,a]=residue(r,p,k)。

例7-2-3

試求下列函數(shù)的部分分式展開(kāi)式

解

>>num=10*[12];

>>den=poly([-1;-3;-4]);

>>[res,poles,k]=residue(num,den)

res=

-6.6667

5.0000

1.6667

poles=

-4.0000

-3.0000

-1.0000

k=

[]根據(jù)上述余數(shù)、極點(diǎn)和留數(shù)，該部分分式展開(kāi)的常數(shù)項(xiàng)結(jié)果是

現(xiàn)在根據(jù)上述余數(shù)、極點(diǎn)和留數(shù)執(zhí)行逆運(yùn)算。

>>[b,a]=residue(res,poles,k)

b=-0.000010.000020.0000

a=

1.00008.000019.000012.0000

>>roots(a)

ans=

-4.0000

-3.0000

-1.0000在截?cái)嗾`差內(nèi)，這與我們開(kāi)始時(shí)的分子和分母多項(xiàng)式一致。即7.2.3polyder()函數(shù)對(duì)有理多項(xiàng)式的求導(dǎo)

函數(shù)polyder()可對(duì)有理多項(xiàng)式求導(dǎo)，其命令格式如下：

[b,a]=polyder(num,den)：如果給出兩個(gè)輸入，則它可對(duì)有理多項(xiàng)式求導(dǎo)，即求多項(xiàng)式num、den商的微分。其中，num為分子，den為分母。

例7-2-4

對(duì)上述多項(xiàng)式F(s)求導(dǎo)。

解

>>num=[1020];den=[181912]

>>[b,a]=polyder(num,den)

b=-20-140-320-260

a=116102328553456144

即求導(dǎo)結(jié)果為 7.3多項(xiàng)式估值與擬合

7.3.1多項(xiàng)式擬合的估值與polyval()函數(shù)

在MATLAB中，根據(jù)多項(xiàng)式系數(shù)的行向量，既可對(duì)多項(xiàng)式進(jìn)行加、減、乘、除和求導(dǎo)，也可以對(duì)它們進(jìn)行估值，這由polyval()或polyvalm()函數(shù)來(lái)完成。

對(duì)于給定的多項(xiàng)式，利用函數(shù)polyval()可以計(jì)算該多項(xiàng)式在任意點(diǎn)的值，其語(yǔ)法如下：

(1)?y=polyval(p,x)：表示返回x在指數(shù)為n的多項(xiàng)式值。輸入變量p是長(zhǎng)度為length=n+1的向量，其元素是被估值的降冪多項(xiàng)式的系數(shù)，即

x可以是矩陣或向量，polyval()函數(shù)對(duì)x的每一個(gè)元素進(jìn)行估值。當(dāng)x是矩陣時(shí)，可使用polyvalm(p,x)函數(shù)。例如：多項(xiàng)式y(tǒng)(x)?=?3x2?+?2x?+?1，估算x在5、7、9時(shí)的值：

>>y=[321];

>>polyval(y,[579])

ans=

86162262

(3)?[y,delta]=polyval(p,x,S)?和[y,delta]=polyval(p,x,S,mu)：使用可選的輸出結(jié)構(gòu)S(由polyfit()函數(shù)生成)產(chǎn)生誤差估計(jì)值y±delta。

例7-3-1

求多項(xiàng)式p(x)=x3?+?4x2-7x-10的估值曲線。

解

(1)輸入多項(xiàng)式系數(shù)向量：

>>p=[14-7-10];

(2)選擇生成數(shù)據(jù)點(diǎn)：

>>x=linspace(-1,3);

(3)計(jì)算x值上的p(x)，把結(jié)果存在v里：

>>v=polyval(p,x);

(4)然后用函數(shù)plot()繪出結(jié)果：

>>plot(x,v),title('x^3+4x^2-7x-10'),xlabel('x')圖7-1多項(xiàng)式估值曲線7.3.2曲線擬合與polyfit()函數(shù)

在許多應(yīng)用領(lǐng)域中，有些數(shù)據(jù)通常是各種物理問(wèn)題和統(tǒng)計(jì)問(wèn)題有關(guān)量的多次觀測(cè)值或由實(shí)驗(yàn)總結(jié)出來(lái)的離散值。離散點(diǎn)組或數(shù)據(jù)往往是零散的，這不僅不便于處理，而且通常不能確切和充分地體現(xiàn)出其固有的規(guī)律。

這種問(wèn)題可由適當(dāng)?shù)慕馕霰磉_(dá)式來(lái)解決，即需要用一個(gè)解析函數(shù)來(lái)描述數(shù)據(jù)。對(duì)這個(gè)問(wèn)題可以有兩種方法解決：曲線擬合和插值法。在曲線擬合(curvefitting)或回歸法中，設(shè)法找出某條光滑曲線，它可以最佳地?cái)M合數(shù)據(jù)，但不必經(jīng)過(guò)任何數(shù)據(jù)點(diǎn)。在插值法里，數(shù)據(jù)假定是正確的，要求以某種方法描述數(shù)據(jù)點(diǎn)之間所發(fā)生的情況。

最佳擬合可用許多不同的方法定義，并存在無(wú)窮數(shù)目的曲線。因此，我們需要一種新的逼近原函數(shù)的手段：

不要求曲線經(jīng)過(guò)所有的點(diǎn)；

盡可能表現(xiàn)數(shù)據(jù)的趨勢(shì)，盡量靠近這些數(shù)據(jù)點(diǎn)。曲線擬合需要回答兩個(gè)基本問(wèn)題：最佳擬合的意義是什么？應(yīng)該用什么樣的曲線？

在數(shù)值分析中，曲線擬合就是用解析表達(dá)式逼近離散數(shù)據(jù)，即離散數(shù)據(jù)的公式化，“最佳擬合”被定義為在數(shù)據(jù)點(diǎn)的最小誤差平方和。

“最小二乘”法(又稱“最小平方”法)是一種數(shù)學(xué)優(yōu)化技術(shù)。它通過(guò)最小化誤差的平方和尋找數(shù)據(jù)的最佳函數(shù)匹配。利用最小二乘法可以簡(jiǎn)便地求得未知的數(shù)據(jù)，并使得這些求得的數(shù)據(jù)與實(shí)際數(shù)據(jù)之間誤差的平方和為最小。最小二乘法可用于曲線擬合，其他一些優(yōu)化問(wèn)題也通過(guò)最小化能量或最大化熵用最小二乘法來(lái)表示。由于“最小二乘”就是使誤差平方和最小，當(dāng)所用的曲線限定為多項(xiàng)式時(shí)，曲線擬合是相當(dāng)簡(jiǎn)捷的。數(shù)學(xué)上稱之為多項(xiàng)式的“最小二乘”曲線擬合。

擬合曲線(虛線)和標(biāo)志的數(shù)據(jù)點(diǎn)之間的垂直距離是在該點(diǎn)的誤差。對(duì)各數(shù)據(jù)點(diǎn)距離求平方，就是誤差平方，然后把誤差平方全加起來(lái)，就是誤差平方和。這條虛線是使誤差平方和盡可能小的曲線，即是最佳擬合。

曲線擬合是工程中經(jīng)常要用到的技術(shù)之一。MATLAB提供了曲線擬合工具箱滿足用戶要求，并可以使用polyfit()函數(shù)求解最小二乘曲線擬合問(wèn)題。

polyfit()函數(shù)給出在最小二乘意義下的最佳擬合系數(shù)。該函數(shù)的用法如下：

(1)?p=polyfit(x,y,n)：找出擬合于數(shù)據(jù)的n階(或度)多項(xiàng)式p(x)的系數(shù)。以最小二乘，從p(x(i))到y(tǒng)(i)，結(jié)果p是長(zhǎng)度為length=n+1的行向量，其元素是降冪的多項(xiàng)式系數(shù)：p(x)?=?p1xn?+?p2xn-1?+?…?+?pnx?+?pn+1

其中，x、y分別為待擬合數(shù)據(jù)的x坐標(biāo)和y坐標(biāo)，n用于指定返回多項(xiàng)式的次數(shù)。當(dāng)n=1時(shí)，得到最簡(jiǎn)單的線性近似，通常稱為一元“線性回歸”。

注意：如果在回歸分析中，只包括一個(gè)自變量一個(gè)因變量，且二者的關(guān)系可用一條直線近似表示，這種回歸分析稱為一元線性回歸分析。

(2)?[p,S]=polyfit(x,y,n)：返回多項(xiàng)式系數(shù)p和一個(gè)結(jié)構(gòu)S，用于誤差估計(jì)或預(yù)測(cè)。

(3)?[p,S,mu]=polyfit(x,y,n)：以下列值取代x，找出多項(xiàng)式系數(shù)：(7.3.1)例7-3-2

已知待擬合數(shù)據(jù)的x坐標(biāo)和y坐標(biāo)如下：

x=[00.10.20.30.40.50.60.70.80.91];

y=[-0.4050.9681.286.185.087.447.6010.538.469.5511.5];

求其最小二乘曲線擬合。

解：(1)輸入多項(xiàng)式系數(shù)向量：

>>x=[00.10.20.30.40.50.60.70.80.91];

>>y=[-0.4050.9681.286.185.087.447.6010.538.469.5511.5];

(2)為了使用polyfit()函數(shù)，我們必須給函數(shù)賦予上面的數(shù)據(jù)和我們希望最佳擬合數(shù)據(jù)的多項(xiàng)式的階次或度。如果我們選n?=?2作為階次，可得到一個(gè)2階多項(xiàng)式。

>>n=2;

(3)求最小二乘曲線擬合。

>>p=polyfit(x,y,n)

p=-?8.040419.5507-0.7627

polyfit()函數(shù)的輸出是一個(gè)多項(xiàng)式系數(shù)的行向量，其解是：

y=-8.0404x2?+?19.5507x-?0.7627

(4)為了將曲線擬合解與數(shù)據(jù)點(diǎn)比較，讓我們把二者都繪成圖。

選擇生成數(shù)據(jù)點(diǎn)。

>>xi=linspace(0,1,100);

為了計(jì)算在xi數(shù)據(jù)點(diǎn)的多項(xiàng)式值，調(diào)用MATLAB的函數(shù)polyval()。

>>z=polyval(p,xi);

繪制曲線。

>>plot(x,y,'o',x,y,xi,z,'.')

>>xlabel('x'),ylabel('y=f(x)'),title('2階曲線擬合')在此畫出了原始數(shù)據(jù)x和y，用“o”標(biāo)出該數(shù)據(jù)點(diǎn)，在數(shù)據(jù)點(diǎn)之間，再用直線重畫原始數(shù)據(jù)，然后用點(diǎn)“.”線畫出多項(xiàng)式數(shù)據(jù)xi和2階擬合曲線z。

根據(jù)上述步驟繪制的結(jié)果顯示于圖7-2中。圖7-22階曲線擬合 7.4數(shù)據(jù)插值

7.4.1一維插值與interp1()函數(shù)

插值定義為對(duì)數(shù)據(jù)點(diǎn)之間函數(shù)的估值方法，這些數(shù)據(jù)點(diǎn)是由某些集合給定的。當(dāng)人們不能很快地求出所需中間點(diǎn)的函數(shù)值時(shí)，插值是一個(gè)有價(jià)值的工具。

具體地說(shuō)，插值運(yùn)算就是根據(jù)數(shù)據(jù)的分布規(guī)律，找到一個(gè)函數(shù)表達(dá)式可以連接已知的各點(diǎn)，并用此函數(shù)表達(dá)式預(yù)測(cè)兩點(diǎn)之間任意位置上的函數(shù)值。例如，當(dāng)數(shù)據(jù)點(diǎn)是某些實(shí)驗(yàn)測(cè)量的結(jié)果時(shí)，就有需要使用插值運(yùn)算的情況。

MATLAB提供了對(duì)數(shù)組的任意一維進(jìn)行插值的工具，這些工具大多需要用到多維數(shù)

組的操作，一維插值在曲線擬合和數(shù)據(jù)分析中具有重要的地位。MATLAB中的一維插值主要有：

多項(xiàng)式插值。

快速傅里葉變換(FFT)插值。

本節(jié)介紹對(duì)數(shù)據(jù)的一維多項(xiàng)式插值。

1．內(nèi)插運(yùn)算與外插運(yùn)算

(1)只對(duì)已知數(shù)據(jù)點(diǎn)集內(nèi)部的點(diǎn)進(jìn)行的插值運(yùn)算稱為內(nèi)插，可比較準(zhǔn)確地估測(cè)插值點(diǎn)上的函數(shù)值。

(2)當(dāng)插值點(diǎn)落在已知數(shù)據(jù)集的外部時(shí)的插值稱為外插，要估計(jì)外插函數(shù)值比較困難。MATLAB的外插結(jié)果偏差較大。

MATLAB對(duì)已知數(shù)據(jù)集外部點(diǎn)上函數(shù)值的預(yù)測(cè)都返回NaN，但可通過(guò)為interp1()函數(shù)添加?'extrap'?參數(shù)指明該函數(shù)也用于外插。

2．一維線性插值

最簡(jiǎn)單插值的例子是MATLAB的作圖，按缺省情況，MATLAB作圖用直線連接所用的數(shù)據(jù)點(diǎn)，這個(gè)線性插值猜測(cè)中間值落在數(shù)據(jù)點(diǎn)之間的直線上。例如：

>>x1=linspace(0,2*pi,60);

>>x2=linspace(0,2*pi,6);

>>plot(x1,sin(x1),x2,sin(x2),'-')

>>xlabel('x'),ylabel('sin(x)'),title('線性插值')

運(yùn)算結(jié)果如圖7-3所示，在數(shù)據(jù)點(diǎn)之間用60個(gè)點(diǎn)的曲線，比只用6個(gè)點(diǎn)的曲線更光滑和精確。線性插值落在數(shù)據(jù)點(diǎn)之間的直線上。當(dāng)數(shù)據(jù)點(diǎn)個(gè)數(shù)增加和它們之間距離的減小時(shí)，線性插值就更精確。圖7-3線性插值

3．一維線性插值與interp1()函數(shù)

如同曲線擬合一樣，插值要作決策，根據(jù)所作的假設(shè)，有多種插值可供選擇，而且可以在一維以上空間中進(jìn)行插值。如果有反映兩個(gè)變量函數(shù)的插值，即已知的數(shù)據(jù)集是平面上的一組離散點(diǎn)集z=f(x,y)，那么就可以在x值之間和在y值之間找出z的中間值進(jìn)行插值，其相應(yīng)的插值就是一維插值，MATLAB中一維插值函數(shù)是interp1()。該函數(shù)的調(diào)用格式為

yi=interp1([x,]y,xi,[method],['extrap'],[extrapval])其中：

[]表示可選，缺省時(shí)為線性插值；

x、y表示采用數(shù)據(jù)的x坐標(biāo)(獨(dú)立變量)和y坐標(biāo)(因變量)；

xi表示待插值位置的數(shù)值數(shù)組；

method表示采用的插值方法，該語(yǔ)句中的method參數(shù)可以選擇的內(nèi)容如表7-2所示；

extrap表示外插運(yùn)算；

extrapval表示外插運(yùn)算返回值，常用0或NaN。

例7-4-1

上述線性插值的正弦函數(shù)使用一維插值函數(shù)interp1()的程序如下：

解

x=linspace(0,2*pi,6);

y=sin(x);

xi=0:.25:6;

yi=interp1(x,y,xi);

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

title('y=sin(x)'),xlabel('x')圖7-4一維插值函數(shù)interp1()的線性插值

例7-4-2

在每天12小時(shí)內(nèi)，每一小時(shí)測(cè)量一次室外溫度，數(shù)據(jù)存儲(chǔ)在兩個(gè)MATLAB變量hours和temps中，請(qǐng)繪制出數(shù)據(jù)點(diǎn)線性插值的曲線。

解使用一維插值函數(shù)interp1()的程序如下：

%ex7_4_2.m

hours=1:12;

temps=[11109.28.510.513.015.116.419.225.426.828.5];

hi=[1:0.5:12];

t=interp1(hours,temps,hi);

plot(hours,temps,'o',hi,t,'-');

title('12小時(shí)溫度曲線')

xlabel('Time(h)'),ylabel('溫度

(C)')圖7-5在12小時(shí)內(nèi)的室外溫度線性插值曲線我們可用上圖對(duì)在1點(diǎn)鐘到12點(diǎn)鐘的12小時(shí)內(nèi)的室外溫度作出溫度分布分析、對(duì)溫度變化趨勢(shì)進(jìn)行解釋、可用interp1()函數(shù)計(jì)算在該時(shí)間段內(nèi)任意給定時(shí)間的溫度。

例如：(1)求早上9:30的室外溫度：

>>t=interp1(hours,temps,9.5)

t=22.3000

(2)求其他各指定時(shí)間的室外溫度：

>>t=interp1(hours,temps,[5.26.57.311.2])

t=11.000014.050015.490027.1400

4．一維樣條插值與interp1()函數(shù)

除了采用直線連接數(shù)據(jù)點(diǎn)外，還可采用某些更光滑的曲線來(lái)擬合數(shù)據(jù)點(diǎn)。常用的方法是用一個(gè)3階多項(xiàng)式，即3次多項(xiàng)式，對(duì)相繼數(shù)據(jù)點(diǎn)之間的各段建模，每個(gè)3次多項(xiàng)式的頭兩個(gè)導(dǎo)數(shù)與該數(shù)據(jù)點(diǎn)相一致。這種類型的插值被稱為3次樣條或簡(jiǎn)稱為樣條。函數(shù)interp1()也能執(zhí)行3次樣條插值。例如：

>>hours=1:12;

>>temps=[11109.28.510.513.015.116.419.225.426.828.5];

>>hi=[1:0.5:12];

>>t=interp1(hours,temps,hi,'spline');

>>plot(hours,temps,'o',hi,t,'-');

>>title('12小時(shí)溫度樣條曲線

');xlabel('Time(h)'),ylabel('溫度(C)')圖7-6在12小時(shí)內(nèi)的室外溫度樣條插值曲線求出早上9:30的溫度：

>>t=interp1(hours,temps,9.5,'spline')

t=22.4187

求出各時(shí)間的溫度：

>>t=interp1(hours,temps,[5.26.57.311.2],'spline')

t=11.020214.146315.534226.8163

樣條插值得到的結(jié)果與上面所示的線性插值的結(jié)果不同。因?yàn)椴逯凳且粋€(gè)估計(jì)或猜測(cè)的過(guò)程，應(yīng)用不同的估計(jì)規(guī)則會(huì)導(dǎo)致不同的結(jié)果。使用樣條插值可獲得一個(gè)更平滑的插值曲線，但不一定是更精確的溫度估計(jì)。

獨(dú)立變量必須是單調(diào)的，即獨(dú)立變量在值上必須總是增加或總是減小的。在上述例子里，自變量hours從1點(diǎn)到12點(diǎn)是單調(diào)的。然而，如果我們已經(jīng)定義獨(dú)立變量為12小時(shí)制一天的實(shí)際時(shí)間，則獨(dú)立變量將不是單調(diào)的，因?yàn)闀r(shí)間從1增加到12后跌到1，再增加到12。如果用24小時(shí)制，也是單調(diào)的。否則，自變量不單調(diào)將會(huì)返回一個(gè)錯(cuò)誤。

不能要求有獨(dú)立變量范圍以外的結(jié)果，例如，interp1(hours,temps,13.5)導(dǎo)致一個(gè)錯(cuò)誤，因?yàn)閔ours在1～12之間變化。7.4.2二維插值與interp2()函數(shù)

二維插值(Two-dimensionaldatainterpolation)是對(duì)兩自變量的函數(shù)z=f(x,y)進(jìn)行插值，其思路與一維插值基本相同。

二維插值比一維插值復(fù)雜，因?yàn)橛懈嗟牧恳３指?。二維插值使用interp2()函數(shù)，已知點(diǎn)集在三維空間中，這些點(diǎn)的插值就二維插值問(wèn)題，這在圖像處理中有廣泛的應(yīng)用，其使用語(yǔ)法如下：

(1)

zi=interp2(x,y,z,xi,yi,method)：二維插值的基本形式。其中，x和y是兩個(gè)獨(dú)立變量，z是一個(gè)應(yīng)變量矩陣，x、y指定z點(diǎn)的數(shù)據(jù)。x和y對(duì)z的關(guān)系是：z(i,:)=f(x,y(i))和z(:,j)=f(x(j),y)。也就是說(shuō)，當(dāng)x變化時(shí)，z的第i行與y的第i個(gè)元素y(i)相關(guān)；當(dāng)y變化時(shí)，z的第j列與x的第j個(gè)元素x(j)相關(guān)。

在所有的情況下，假定獨(dú)立變量x和y是線性間隔和單調(diào)的，并有相同的格式(柵格)。

xi是沿x軸插值的一個(gè)數(shù)值數(shù)組；yi是沿y軸插值的一個(gè)數(shù)值數(shù)組。

method是可選的參數(shù)?？梢允?linear'、'cubic'、'spline'或'nearest'。代表的意義如下：

linear方法是默認(rèn)值，代表雙線性插值(雙線性內(nèi)插)，僅用作連接圖上數(shù)據(jù)點(diǎn)。

nearest方法只選擇最接近各估計(jì)點(diǎn)的粗略數(shù)據(jù)點(diǎn)。

'spline'：3次樣條插值。'cubic'：雙3次樣條插值。在這種情況下，cubic不意味著3次樣條，而是使用3次多項(xiàng)式的另一種算法。

(2)?zi?=?interp2(x,y,z,xi,yi)：二維插值的默認(rèn)形式，即雙線性插值。

例7-4-3peaks()函數(shù)經(jīng)二維插值后，獲得更細(xì)密的數(shù)據(jù)點(diǎn)。程序如下：

解

[X,Y]=meshgrid(-3:.25:3);

Z=peaks(X,Y);

[XI,YI]=meshgrid(-3:.125:3);

ZI=interp2(X,Y,Z,XI,YI);

mesh(X,Y,Z),hold,mesh(XI,YI,ZI+15)

holdoff

axis([-33-33-520])圖7-7二維插值后的peaks()函數(shù)

(3)?zi?=?interp2(z,xi,yi)：二維插值的默認(rèn)形式，即雙線性插值。假設(shè)x=1:n和y=1:m，[m,n]=size(z)。

(4)?zi=interp2(i,ntimes)：在每?jī)蓚€(gè)元素之間插入插值，按ntimes遞歸擴(kuò)展z，interp2(Z)等同于interp2(Z,1)。例7-4-4

要對(duì)平板上的溫度分布進(jìn)行估計(jì)，給定的溫度值取自平板表面均勻分布的格柵。格柵的寬度用自變量width表示，在數(shù)據(jù)矩陣中代表矩陣的行向量，在二維插值中是x維；格柵的高度(深度)用自變量depth表示，在數(shù)據(jù)矩陣中代表矩陣的列向量，在二維插值中是y維。因變量矩陣temps表示整個(gè)平板的溫度分布，為了估計(jì)在數(shù)據(jù)點(diǎn)之間的點(diǎn)的溫度，必須對(duì)它們進(jìn)行插值。

解采集的數(shù)據(jù)和插值程序如下：

width=1:5;depth=1:3;

temps=[8081808284;7963616581;8484828586];

wi=1:0.2:5;d=2;

zlinear=interp2(width,depth,temps,wi,d);%線性插值

zcubic=interp2(width,depth,temps,wi,d,'cubic');%雙3次樣條插值

plot(wi,zlinear,'-',wi,zcubic)%plotresults

xlabel('板的寬度(Width)'),ylabel('板的溫度(DegreesCelsius)')

title(['深度depth='num2str(d)'處的溫度'])圖7-8在深度d=2處的平板溫度同時(shí)，我們可以在三維空間立體觀察溫度分布情況。首先在三維坐標(biāo)畫出原始數(shù)據(jù)，觀察該原始數(shù)據(jù)表示的分布情況。

>>mesh(width,depth,temps)

>>xlabel('板的寬度(Width)'),ylabel('板的深度(Depth)'),zlabel('板的溫度(DegreesCelsius)')

>>title(['深度depth='num2str(d)'處的溫度']),axis('ij')圖7-9平板溫度的立體分布然后在兩個(gè)方向上插值，以平滑數(shù)據(jù)。

>>[wi,di]=meshgrid(1:0.2:5,1:0.2:3,temps);

>>zcubic=interp2(width,depth,temps,wi,di,'cubic');

>>mesh(wi,di,zcubic)

>>xlabel('板的寬度(Width)'),ylabel('板的深度(Depth)'),zlabel('板的溫度(DegreesCelsius)')

>>title(['深度depth='num2str(d)'處的溫度']),axis('ij')

二維插值后的平板溫度的立體分布情況如圖7-10所示。圖7-10二維插值后的平板溫度的立體分布7.4.3抽樣插值與interp()函數(shù)

interp()函數(shù)用一個(gè)整數(shù)系數(shù)(插值)增加采樣率，插值的方法是在原序列中插入一些0值，以提高采樣率。語(yǔ)法如下：

(1)?y=interp(x,r)：用一個(gè)整數(shù)系數(shù)(插值)?r增加x的采樣率，插值后的向量y的長(zhǎng)度數(shù)倍于原始輸入信號(hào)x。

(2)?y=interp(x,r,length,alpha)：指定濾波器的長(zhǎng)度length，截止頻率alpha，默認(rèn)值是length=4，alpha=0.5。

(3)?[y,b]=interp(x,r,l,alpha)：返回插值y和濾波器系數(shù)b。

例7-4-5

抽樣信號(hào)的插值。

解

t=0:0.001:1;%時(shí)間向量

x=sin(2*pi*30*t)+sin(2*pi*60*t);

y=interp(x,4);

stem(x(1:30));

title('原抽樣信號(hào)');

figure

stem(y(1:120));

title('插值抽樣信號(hào)');圖7-11原抽樣信號(hào)圖7-12插值抽樣信號(hào)7.4.4三次樣條與spline()函數(shù)

1．基本特征

在三次樣條中，要找出一個(gè)三次多項(xiàng)式，以逼近每對(duì)數(shù)據(jù)點(diǎn)間的曲線。在樣條術(shù)語(yǔ)中，這些數(shù)據(jù)點(diǎn)稱之為斷點(diǎn)。因?yàn)閮牲c(diǎn)只能決定一條直線，而在兩點(diǎn)間的曲線可用無(wú)限多的三次多項(xiàng)式近似。因此，為使結(jié)果具有唯一性，在三次樣條中，增加了三次多項(xiàng)式的約束條件。通過(guò)限定每個(gè)三次多項(xiàng)式的一階和二階導(dǎo)數(shù)，使其在斷點(diǎn)處相等，就可以較好地確定所有內(nèi)部三次多項(xiàng)式。此外，近似多項(xiàng)式通過(guò)這些斷點(diǎn)的斜率和曲率是連續(xù)的。然而，第一個(gè)和最后一個(gè)三次多項(xiàng)式在第一個(gè)和最后一個(gè)斷點(diǎn)以外，沒(méi)有伴隨多項(xiàng)式。因此必須通過(guò)其他方法確定其余的約束。

MATLAB提供了spline()函數(shù)逼近樣條。函數(shù)spline()所采用的方法，也是最常用的方法，就是采用非結(jié)(not-a-knot)條件。這個(gè)條件強(qiáng)迫第一個(gè)和第二個(gè)三次多項(xiàng)式的三階導(dǎo)數(shù)相等。對(duì)最后一個(gè)和倒數(shù)第二個(gè)三次多項(xiàng)式也做同樣地處理。

基于上述描述，人們可能猜想到，尋找三次樣條多項(xiàng)式需要求解大量的線性方程。實(shí)際上，給定N個(gè)斷點(diǎn)，就要尋找N-1個(gè)三次多項(xiàng)式，每個(gè)多項(xiàng)式有4個(gè)未知系數(shù)。這樣，所求解的方程組包含有4?×?(N-1)個(gè)未知數(shù)。把每個(gè)三次多項(xiàng)式列成特殊形式，并且運(yùn)用各種約束，通過(guò)求解N個(gè)具有N個(gè)未知系數(shù)的方程組，就能確定三次多項(xiàng)式。這樣，如果有50個(gè)斷點(diǎn)，就有50個(gè)具有50個(gè)未知系數(shù)的方程組。所幸的是，若使用稀疏矩陣，則這些方程式就能夠簡(jiǎn)明地列出并求解，這就是函數(shù)spline()所使用的計(jì)算未知系數(shù)的方法。spline()函數(shù)的使用方式如下：

(1)?cs=spline(x,y)：以已知的斷點(diǎn)序列x，在所有i并滿足“not-a-knot”(非節(jié)點(diǎn))條件的x(i)點(diǎn)上的y(i)值，返回三次樣條的ppform(pp形式或分段多項(xiàng)式形式)。給出的結(jié)果與在樣條工具箱中使用命令cs=csapi(x,y)相同。

(2)?yi=spline(x,y,xi)：使用三次樣條插值來(lái)尋找y在矢量xi指定的點(diǎn)上的函數(shù)值yi。y是x的函數(shù)，即y?=?f(x)，x和y是一一對(duì)應(yīng)關(guān)系，均是有多個(gè)元素的向量，如果y是一個(gè)矩陣，數(shù)據(jù)結(jié)果是一個(gè)向量，插值按y的每一行算出。這種情況下，length(x)必須等于size(y,2)，yy是size(y,1)*length(xx)向量。

該語(yǔ)法與一維插值函數(shù)yi=interp1(x,y,xi,'spline')結(jié)果相同。

2．分段多項(xiàng)式

在最簡(jiǎn)單的用法中，spline()獲取數(shù)據(jù)x和y以及期望值xi，尋找擬合x(chóng)和y的三次樣條內(nèi)插多項(xiàng)式，然后計(jì)算這些多項(xiàng)式，對(duì)每個(gè)xi的值，尋找相應(yīng)的yi。例如：

>>x=0:12;

>>y=tan(pi*x/25);

>>xi=linspace(0,12);

>>yi=spline(x,y,xi);

>>plot(x,y,'o',xi,yi),title('Spline三次樣條內(nèi)插')圖7-13樣條擬合結(jié)果這種方法適合于只需要一組內(nèi)插值的情況。不過(guò)，如果用戶需要多組內(nèi)插值數(shù)據(jù)時(shí)，例如需要從相同數(shù)據(jù)集里獲取另一組內(nèi)插值，再次計(jì)算三次樣條系數(shù)是沒(méi)有意義的。在這種情況下，可以調(diào)用僅帶前兩個(gè)參量的spline()：

>>pp=spline(x,y)

pp=

form:'pp'

breaks:[0123456789101112]

coefs:[12x4double]

pieces:12

order:4

dim:1上述給定的三次樣條pp形式存儲(chǔ)了斷點(diǎn)和多項(xiàng)式系數(shù)，以及關(guān)于三次樣條表示的其他信息。因?yàn)樗行畔⒍急淮鎯?chǔ)在單個(gè)向量里，所以這種形式在MATLAB中是一種方便的數(shù)據(jù)結(jié)構(gòu)。

當(dāng)采用這種方式調(diào)用時(shí)，spline返回一個(gè)稱之為三次樣條的pp形式或分段多項(xiàng)式形式的數(shù)組。這個(gè)數(shù)組pp包含了對(duì)于任意一組所期望的內(nèi)插值和計(jì)算三次樣條所必須的全部信息。給定pp形式，可以在樣條工具箱中使用，函數(shù)ppval可以使用pp反復(fù)計(jì)算該三次樣條。例如：

>>yi=ppval(pp,xi);計(jì)算前面計(jì)算過(guò)的yi，與圖7-12所示的結(jié)果相同。

利用分段多項(xiàng)式形式可以很方便地在一個(gè)更精確的區(qū)域[10,12]內(nèi)重新計(jì)算該三次樣條插值。例如：

>>xi2=linspace(10,12);

>>yi2=ppval(pp,xi2);

>>plot(x,y,'o',xi2,yi2),title('Spline三次樣條內(nèi)插')

在一個(gè)更精確的區(qū)域[10,12]內(nèi)，重新計(jì)算該三次樣條插值，樣條擬合結(jié)果如圖7-14所示。圖7-14重新計(jì)算樣條插值擬合結(jié)果運(yùn)用pp形式，還可以計(jì)算區(qū)域之外的插值。當(dāng)數(shù)據(jù)出現(xiàn)在最后一個(gè)斷點(diǎn)之后或第一個(gè)斷點(diǎn)之前時(shí)，則分別運(yùn)用最后一個(gè)或第一個(gè)三次多項(xiàng)式來(lái)尋找內(nèi)插值。

例如，在限定的區(qū)間[10,15]內(nèi)，再次計(jì)算該三次樣條，注意pp的有效區(qū)域是[0,12]。

>>xi3=10:15;

>>yi3=ppval(pp,xi3);

>>plot(x,y,'o',xi3,yi3),title('Spline三次樣條內(nèi)插')

樣條擬合結(jié)果如圖7-15所示。圖7-15重新計(jì)算區(qū)域外樣條插值擬合結(jié)果

3．pp形式的分解與重構(gòu)

當(dāng)要計(jì)算三次樣條表示時(shí)，必須把pp形式分解成它的各個(gè)表示段。在MATLAB中，通過(guò)函數(shù)unmkpp()完成這一過(guò)程，使用方法如下：

[breaks,coefs,l,k,d]=unmkpp(pp)：從分段多項(xiàng)式pp中提取斷點(diǎn)breaks、系數(shù)coefs、分段數(shù)l、階次k和維數(shù)d。運(yùn)用上述pp形式，該函數(shù)給出如下結(jié)果：

>>[breaks,coefs,npolys,ncoefs,d]=unmkpp(pp)

breaks=0123456789101112

coefs=

0.0007-0.00010.12570

0.0007?0.00200.12760.1263

0.0010?0.00420.13390.2568

0.0012?0.00720.14540.3959

0.0024?0.01090.16350.5498

0.0019?0.01810.19250.7265

0.0116?0.02370.23440.9391

-0.0083?0.05860.31671.2088

0.1068?0.03360.40891.5757

-0.1982?0.35420.79672.1251

1.4948-0.24060.91023.0777

1.4948?4.24394.91365.2422

npolys=12

ncoefs=4

d=1這里，breaks是斷點(diǎn)，coefs是矩陣，它的第i行是第i個(gè)三次多項(xiàng)式，npolys是多項(xiàng)式的分段數(shù)目，ncoefs是每個(gè)多項(xiàng)式系數(shù)的數(shù)目(階次)。注意，這種形式非常一般，樣條多項(xiàng)式不必是三次，這對(duì)于樣條的積分和微分是很有用的。

如果給定上述分散形式，函數(shù)mkpp()可恢復(fù)pp形式。

(1)?pp=mkpp(breaks,coefs)：以其斷點(diǎn)breaks和系數(shù)coefs生成分段多項(xiàng)式的pp。breaks是一個(gè)長(zhǎng)度為L(zhǎng)+1的向量，從開(kāi)始和結(jié)束的每個(gè)元素都嚴(yán)格按L間隔增量。coefs是L?×?k的矩陣，每一行coefs(i,:)都包含k階多項(xiàng)式有效區(qū)間[breaks(i),breaks(i+1)]中的系數(shù)，從最高到最低的冪指數(shù)排列。

(2)?pp=mkpp(breaks,coefs,d)：指示的分段多項(xiàng)式pp是d值向量，即其系數(shù)的每個(gè)值是長(zhǎng)度為d的向量，breaks是長(zhǎng)度L+1的增向量。coefs是d?×?L?×?k?數(shù)組，coefs(r,i,:)包含分段多項(xiàng)式的第r個(gè)元素、第i個(gè)多項(xiàng)式段的k系數(shù)。

例如：

>>pp=mkpp(breaks,coefs)

pp=

form:'pp'

breaks:[0123456789101112]

coefs:[12x4double]

pieces:12

order:4

dim:1

7.5數(shù)值分析

7.5.1求極值

最優(yōu)化是求最優(yōu)解，也就是在某個(gè)區(qū)間內(nèi)有條件約束或者無(wú)條件約束地找到函數(shù)的最大值或者最小值。

在許多應(yīng)用中，需要確定函數(shù)的極值，即最大值(峰值)和最小值(谷值)。數(shù)學(xué)上，可通過(guò)確定函數(shù)導(dǎo)數(shù)為零的點(diǎn)，解析出這些極值點(diǎn)。然而，很容易求導(dǎo)的函數(shù)，常常很難找到導(dǎo)數(shù)為零的點(diǎn)。在不可能從解析上求得導(dǎo)數(shù)的情況下，必須從數(shù)值上尋找函數(shù)的極值點(diǎn)。

MATLAB使用數(shù)字方法，即迭代算法求函數(shù)的最小值，也就是有些步驟要重復(fù)許多次。MATLAB提供了兩個(gè)完成此功能的函數(shù)fminbnd()和fmins()，這兩個(gè)函數(shù)分別尋找一維或n維函數(shù)的最小值。但在MATLABR2010版本中不再包含下列函數(shù)：fmin、fmins、icubic、interp4、interp5、interp6、meshdom、nnls和saxis。

MATLAB沒(méi)有求最大值的命令，因?yàn)閒(x)的最大值等于-f(x)的最小值，所以利用相反函數(shù)h(x)=-f(x)的最小值可以求得最大值，故fminbnd()可用來(lái)求最小值和最大值。

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

(1)?fminbnd()函數(shù)求得在給定區(qū)間內(nèi)的函數(shù)極小值的自變量值。該函數(shù)的調(diào)用格式為

x=fminbnd(fun,x1,x2)

x=fminbnd(fun,x1,x2,options)

其中：fun為函數(shù)句柄，可用于M文件函數(shù)或匿名函數(shù)；x1和x2可分別用于指定區(qū)間的左右邊界；options為可選項(xiàng)，用于指定程序的其他參數(shù)，其取值如表7-4所示。

(2)?feval()函數(shù)求得在給定區(qū)間內(nèi)的函數(shù)極小值或極大值。該函數(shù)的調(diào)用格式為

y=feval(fun,x)

其中，fun為函數(shù)句柄，可用于M文件函數(shù)或匿名函數(shù)，x是極小值或極大值的自變量值。

(3)?[x,fval]=fminbnd(...)：直接求出函數(shù)極值的自變量值和函數(shù)的極值。

例7-5-1

已知函數(shù)y

=

2.0e-tsin(t)，求其極大值和極小值。

解

(1)繪制函數(shù)圖形：

>>y=@(t)2.0*exp(-t).*sin(t);

>>fplot(y,[08]);

>>title('函數(shù)極值'),xlabel('(t)'),ylabel('y=2.0*exp(-t).*sin(t)');

在區(qū)間0≤t≤8內(nèi)繪出上述函數(shù)，得到圖7-16所示的圖形。

由圖可知，在tmax=0.7附近有一個(gè)最大值，并且在tmin=4附近有一個(gè)最小值。而這些點(diǎn)的解析值為：tmax=p/4≈0.785和tmin=5p/4≈3.93。圖7-16繪制函數(shù)曲線圖形

(2)使用fminbnd()函數(shù)求函數(shù)在給定區(qū)間內(nèi)的極小值的自變量值：

>>tmin=fminbnd(y,0,8)

tmin=3.9270

其誤差為

>>emin=5*pi/4-tmin

emin=1.6528e-005

(3)使用fminbnd()函數(shù)求反函數(shù)在給定區(qū)間內(nèi)的極小值的自變量值，作為該函數(shù)的極大值的自變量值：

>>ym=@(t)(-2.0*exp(-t).*sin(t));%定義反函數(shù)

>>tmax=fminbnd(ym,0,8)

tmax=0.7854

其誤差為

>>emax=pi/4-tmax

emax=-1.0079e-005

(4)根據(jù)極小值的自變量值tmin和極大值的自變量值tmax，使用feval()函數(shù)，求出函數(shù)的極小值和極大值。

>>ymin=feval(y,tmin)

ymin=-0.0279

>>ymax=feval(y,tmax)

ymax=0.6448

上述過(guò)程可以直接用下列方式求出極值的自變量值和函數(shù)的極值：

>>[tminymin]=fminbnd(y,0,8)

tmin=3.9270

ymin=-0.0279

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

MATLAB提供的函數(shù)fminsearch()用于計(jì)算多元函數(shù)的極小值。fminsearch()函數(shù)內(nèi)部應(yīng)用了Nelder-Mead單一搜索算法，通過(guò)調(diào)整x的各個(gè)元素的值來(lái)尋找f(x)的極小值。該算法雖然對(duì)于平滑函數(shù)搜索效率沒(méi)有其他算法高，但它不需要梯度信息，從而擴(kuò)展了其應(yīng)用范圍。因此，該算法特別適用于不太平滑、難以計(jì)算梯度信息或梯度信息價(jià)值不大的函數(shù)。

fminsearch()函數(shù)與fminbnd()函數(shù)的用法基本相同，不同之處在于：fminbnd()函數(shù)的輸入?yún)?shù)為尋找最小值的區(qū)間，并且只能用于求解一元函數(shù)的極值，fminsearch()函數(shù)的輸入?yún)?shù)為初始值。該函數(shù)的調(diào)用格式為

x=fminsearch(fun,x0)

x=fminsearch(fun,x0,options)

例：Rosenbrock的banana()函數(shù)是一個(gè)經(jīng)典的例子，用于測(cè)試fminsearch()函數(shù)計(jì)算多元函數(shù)的極小值。banana()函數(shù)的定義為可知函數(shù)的最小值在[11]點(diǎn)，最小值為0，傳統(tǒng)的起始點(diǎn)選在[-1.21]。定義一個(gè)匿名函數(shù)banana如下：

>>banana=@(x)100*(x(2)-x(1)^2)^2+(1-x(1))^2;

使用fminsearch()函數(shù)計(jì)算多元函數(shù)的極小值：

>>[x,fval]=fminsearch(banana,[-1.2,1])

x=

1.00001.0000

fval=

8.1777e-010

函數(shù)的最小值在[11]點(diǎn)，函數(shù)的極小值近似于0，結(jié)果與之相符。7.5.2求零點(diǎn)

求函數(shù)f(x)的零值就等于求方程f(x)?=0的解，對(duì)于多項(xiàng)式可以用roots()函數(shù)來(lái)求解。單變量函數(shù)的零值可以用fzero()函數(shù)來(lái)求解。fzero()用迭代法來(lái)求解，使得初始的估計(jì)值接近理想的零值。

尋找一元函數(shù)零點(diǎn)時(shí)，可以指定一個(gè)初始點(diǎn)，或者指定一個(gè)區(qū)間，然后使用fzero()函數(shù)來(lái)求一元函數(shù)的零點(diǎn)。當(dāng)指定一個(gè)初始點(diǎn)時(shí)，此函數(shù)在初始點(diǎn)附近尋找一個(gè)使函數(shù)值變號(hào)的區(qū)間，如果沒(méi)有找到這樣的區(qū)間，則函數(shù)返回NaN。該函數(shù)的調(diào)用格式為

x=fzero(fun,x0)，x=fzero(fun,[x1,x2])：尋找x0附近或者區(qū)間[x1,x2]內(nèi)fun的零點(diǎn)，返回該點(diǎn)的x坐標(biāo)；

x=fzero(fun,x0,options)，x=fzero(fun,[x1,x2],options)：通過(guò)options設(shè)置參數(shù)；

[x,fval]=fzero(...)：返回零點(diǎn)的同時(shí)返回該點(diǎn)的函數(shù)值；

[x,fval,exitflag]=fzero(...)：返回零點(diǎn)、該點(diǎn)的函數(shù)值及程序退出的標(biāo)志；

[x,fval,exitflag,output]=fzero(...)：返回零點(diǎn)、該點(diǎn)的函數(shù)值、程序退出的標(biāo)志及選定的輸出結(jié)果。

例7-5-2

在MATLAB中提供一個(gè)名為humps()的函數(shù)，定義為求該函數(shù)的零點(diǎn)。

解

(1)定義一個(gè)匿名函數(shù)：

>>y=@(x)(1./((x-.3).^2+.01)+1./((x-.9).^2+.04)-6);

(2)繪制出該函數(shù)圖形：

>>fplot(y,[02])

>>line([02],[00])

>>title('humps函數(shù)')圖7-17humps()函數(shù)的圖形

(3)使用fzero()函數(shù)尋找該一維函數(shù)的零點(diǎn)。為了說(shuō)明該函數(shù)的使用，讓我們?cè)龠\(yùn)用humps例子。

>>xzero=fzero('humps',[02])

xzero=1.2995

所以，fzero()函數(shù)找出humps函數(shù)的零點(diǎn)接近于1.3，即1.2995。

fzero()函數(shù)不僅能尋找零點(diǎn)，還可以尋找函數(shù)等于任何常數(shù)值的點(diǎn)。例如，為了尋找f(x)=c的點(diǎn)，定義函數(shù)g(x)=f(x)-c，然后，在fzero()函數(shù)中使用g(x)，找出g(x)為零的x值，即f(x)=c。7.5.3數(shù)值積分

一個(gè)函數(shù)的積分或面積也是它的另一個(gè)有用的屬性。MATLAB中提供了用于積分的函數(shù)，包括一元函數(shù)的自適應(yīng)數(shù)值積分、一元函數(shù)的矢量積分、二重積分和三重積分。

這些函數(shù)如表7-5所示。MATLAT提供了在有限區(qū)間內(nèi)，數(shù)值計(jì)算某函數(shù)下的面積的三種函數(shù)：trapz()、quad()和quad8()。

1．trapz()函數(shù)

函數(shù)trapz()通過(guò)計(jì)算若干梯形面積的和來(lái)近似某函數(shù)的積分，這些梯形是通過(guò)使用函數(shù)humps的數(shù)據(jù)點(diǎn)形成。調(diào)用格式如下：

Z=trapz(Y)：使用梯形法(以單位間距)計(jì)算Y的積分的近似值。若要計(jì)算多于1個(gè)單位間距的積分，就將Z乘以間距增量。如果Y是向量，trapz(Y)是Y的積分。如果Y是一個(gè)矩陣，trapz(Y)是一個(gè)行向量，每個(gè)元素是Y的每一列的積分值。

Z=trapz(X,Y)：使用梯形積分計(jì)算Y關(guān)于X的積分。X表示橫坐標(biāo)向量，Y為對(duì)應(yīng)的縱坐標(biāo)向量。要求X與Y的長(zhǎng)度必須相等。

Z=trapz(...,dim)：使用梯形積分計(jì)算Y跨維積分，維度由標(biāo)量dim指定。

例7-5-3

使用梯形法計(jì)算積分。

解函數(shù)。

>>X=0:pi/100:pi;

>>Y=sin(X);

>>Z=trapz(X,Y)

Z=1.9998

2．一元函數(shù)的積分

MATLAB中一元函數(shù)的積分可以用quad()和quadl()兩個(gè)函數(shù)來(lái)實(shí)現(xiàn)。函數(shù)quad()和quad8()是基于數(shù)學(xué)上的正方形概念來(lái)計(jì)算函數(shù)的面積，函數(shù)quad()采用低階的自適應(yīng)遞歸Simpson方法，quad8()函數(shù)采用高階的自適應(yīng)遞歸Simpson方法，而且quad8()比quad()更精確。quadl()采用高階自適應(yīng)Lobatto方法，該函數(shù)是quad8()函數(shù)的替代。

quad()函數(shù)的調(diào)用格式如下：

(1)?q=quad(fun,a,b)：采用遞歸自適應(yīng)方法計(jì)算函數(shù)fun()在區(qū)間[a,b]的定積分，其精確度為1e-6，fun可以是M文件，也可以是匿名函數(shù)，即定義為(7.5.1)例如，上例對(duì)正弦函數(shù)的積分：

>>f=@(X)sin(X);

>>F=quad(f,0,pi)

F=2.0000

(2)

q=quad(fun,a,b,tol)：指定允許誤差，指定的誤差tol需大于1e-6。該命令運(yùn)行更快，但是得到的結(jié)果精確度降低。

(3)?q=quad(fun,a,b,tol,trace)：跟蹤迭代過(guò)程，輸出[fcntab-aQ]的值，分別為計(jì)算函數(shù)值的次數(shù)fcnt、當(dāng)前積分區(qū)間的左邊界a、步長(zhǎng)b-a和該區(qū)間內(nèi)的積分值Q。

(4)?[q,fcnt]=quadl(fun,a,b,...)：輸出函數(shù)值的同時(shí)輸出計(jì)算函數(shù)值的次數(shù)。

例如，估計(jì)humps()函數(shù)在[0,2]區(qū)間的面積：

>>area=quad('humps',0,2)

area=29.3262

>>area=quad8('humps',0,2)

area=29.3262

這兩個(gè)函數(shù)返回完全相同的估計(jì)面積。

3．一元函數(shù)的矢量積分

矢量積分相當(dāng)于多個(gè)一元函數(shù)積分。當(dāng)被積函數(shù)中含有參數(shù)，需要對(duì)該參數(shù)的不同值計(jì)算該函數(shù)的積分時(shí)，可以使用一元函數(shù)的矢量積分。

矢量積分返回一個(gè)向量，每個(gè)元素的值為一個(gè)一元函數(shù)的積分值。quadv()函數(shù)與quad()和quadl()函數(shù)相似，可以設(shè)置積分參數(shù)和結(jié)果輸出。

4．二重積分和三重積分

在MATLAB中，二重積分和三重積分分別由函數(shù)dblquad()和函數(shù)triplequad()來(lái)實(shí)現(xiàn)。首先介紹函數(shù)dblquad()，該函數(shù)的基本格式如下：

(1)?q=dblquad(fun,xmin,xmax,ymin,ymax)，函數(shù)的參數(shù)分別為函數(shù)句柄、兩個(gè)自變量的積分上下限，返回二重積分結(jié)果。二重積分的定義為(7.5.2)

dblquad()函數(shù)通過(guò)調(diào)用quad()函數(shù)來(lái)計(jì)算fun(x,y)在矩形區(qū)域xmin≤x≤xmax、ymin≤y≤ymax的雙重積分，fun()可以是函數(shù)句柄、M文件函數(shù)或匿名函數(shù)。

(2)?q=dblquad(fun,xmin,xmax,ymin,ymax,tol)：指定積分結(jié)果的精度tol，默認(rèn)值是tol=1.0e-6。

(3)?q=dblquad(fun,xmin,xmax,ymin,ymax,tol,method)：指定結(jié)果精度tol和積分方法method，method的取值可以是@quadl，也可以是用戶自定義的積分函數(shù)句柄，該函數(shù)的調(diào)用格式必須與quad()的調(diào)用格式相同。

triplequad()函數(shù)的調(diào)用格式和dblq