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

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

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

7.1MATLAB多項式

7.1.1概述

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

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

1．多項式與行向量

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

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

2．多項式的運算

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

因為多項式的加減為對應(yīng)項系數(shù)的加減，所以可以通過向量的加減來實現(xiàn)。但是在向量的加減中兩個向量需要有相同的長度，因此在進行多項式加減時，需要將短的向量前面補0。

多項式的乘法實際上是多項式系數(shù)向量之間的卷積運算，可以通過MATLAB中的卷積函數(shù)conv()來完成。

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

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

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

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

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

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

例7-1-1

求多項式x4-12x3?+?25x?+?116的根。

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

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

>>r=roots(C)

r=

11.7473

2.7028

-1.2251+1.4672i

-1.2251-1.4672i

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

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

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

>>pp=poly(r)

pp=

1.0000-12.0000-0.000025.0000116.0000

(2)

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

×

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

>>A=[123

456

780]

>>p=poly(A)

p=

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

因為MATLAB無隙地處理復(fù)數(shù)，當(dāng)用根重組多項式時，如果一些根有虛部，由于截斷所產(chǎn)生的誤差，則poly()的結(jié)果有一些小的虛部，這是很普通的。要消除虛假的虛部，只要使用函數(shù)real抽取實部即可。7.1.3卷積運算與多項式乘法

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

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

定義：若m=length(u)、n=length(v)，則w是長度為m?+?n-1的向量。第k個元素是：總和w(k)是所有u(j)和v(k?+?1-j)的乘積的和，j的取值范圍是j=max(1,k+1-n):min(k,m)。當(dāng)m=n時，即兩向量長度相同，則w是長度為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ù)支持多項式乘法，即執(zhí)行兩個數(shù)組的卷積，兩個以上的多項式的乘法需要重復(fù)使用conv。7.1.4反卷積運算與多項式除法

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

[q,r]=deconv(v,u)：使用于長除法，輸出向量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在一些特殊情況，一個多項式需要除以另一個多項式。如果u和v是多項式系數(shù)向量，卷積等于多項式乘，則反卷積等于多項式除，v被u除，商是q，余數(shù)是r。

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

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

b=149-15

r=0000000

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

對于多項式的加法，MATLAB沒有提供一個直接的函數(shù)。如果兩個多項式向量大小相同，則可使用標(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)兩個多項式階次不同時，低階的多項式必須用首零填補，使其與高階多項式有同樣的階次。這里要求是首零而不是尾零，是因為相關(guān)的系數(shù)像x冪次一樣，必須整齊。

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

>>e=c+[000d]

e=

162021193-71

結(jié)果是：

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

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

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

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

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

例7-1-3

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

解(1)輸入多項式系數(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)換為多項式。

解程序如下：

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有理多項式的運算

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

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

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

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

b=[b1b2…bm]

a=[a1a2…an]

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

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

解對于該函數(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=

[

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

解由此函數(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)中含多重極點，即p(j)?=?…?=?p(j?+?m-1)，則部分分式展開式將包括下列各項(7.2.3)式中，pj為一個m重極點。7.2.2residue()函數(shù)的逆運算

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

例7-2-3

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

解

>>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ù)、極點和留數(shù)，該部分分式展開的常數(shù)項結(jié)果是

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

>>[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在截斷誤差內(nèi)，這與我們開始時的分子和分母多項式一致。即7.2.3polyder()函數(shù)對有理多項式的求導(dǎo)

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

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

例7-2-4

對上述多項式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多項式估值與擬合

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

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

對于給定的多項式，利用函數(shù)polyval()可以計算該多項式在任意點的值，其語法如下：

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

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

>>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)生誤差估計值y±delta。

例7-3-1

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

解

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

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

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

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

(3)計算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多項式估值曲線7.3.2曲線擬合與polyfit()函數(shù)

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

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

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

不要求曲線經(jīng)過所有的點；

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

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

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

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

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

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

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

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

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

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

(3)?[p,S,mu]=polyfit(x,y,n)：以下列值取代x，找出多項式系數(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)輸入多項式系數(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ù)的多項式的階次或度。如果我們選n?=?2作為階次，可得到一個2階多項式。

>>n=2;

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

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

p=-?8.040419.5507-0.7627

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

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

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

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

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

為了計算在xi數(shù)據(jù)點的多項式值，調(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ù)點，在數(shù)據(jù)點之間，再用直線重畫原始數(shù)據(jù)，然后用點“.”線畫出多項式數(shù)據(jù)xi和2階擬合曲線z。

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

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

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

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

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

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

多項式插值。

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

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

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

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

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

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

2．一維線性插值

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

>>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('線性插值')

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

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

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

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

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

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

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

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

extrap表示外插運算；

extrapval表示外插運算返回值，常用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小時內(nèi)，每一小時測量一次室外溫度，數(shù)據(jù)存儲在兩個MATLAB變量hours和temps中，請繪制出數(shù)據(jù)點線性插值的曲線。

解使用一維插值函數(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小時溫度曲線')

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

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

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

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

t=22.3000

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

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

t=11.000014.050015.490027.1400

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

除了采用直線連接數(shù)據(jù)點外，還可采用某些更光滑的曲線來擬合數(shù)據(jù)點。常用的方法是用一個3階多項式，即3次多項式，對相繼數(shù)據(jù)點之間的各段建模，每個3次多項式的頭兩個導(dǎo)數(shù)與該數(shù)據(jù)點相一致。這種類型的插值被稱為3次樣條或簡稱為樣條。函數(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小時溫度樣條曲線

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

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

t=22.4187

求出各時間的溫度：

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

t=11.020214.146315.534226.8163

樣條插值得到的結(jié)果與上面所示的線性插值的結(jié)果不同。因為插值是一個估計或猜測的過程，應(yīng)用不同的估計規(guī)則會導(dǎo)致不同的結(jié)果。使用樣條插值可獲得一個更平滑的插值曲線，但不一定是更精確的溫度估計。

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

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

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

二維插值比一維插值復(fù)雜，因為有更多的量要保持跟蹤。二維插值使用interp2()函數(shù)，已知點集在三維空間中，這些點的插值就二維插值問題，這在圖像處理中有廣泛的應(yīng)用，其使用語法如下：

(1)

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

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

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

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

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

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

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

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

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

解

[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)：二維插值的默認形式，即雙線性插值。假設(shè)x=1:n和y=1:m，[m,n]=size(z)。

(4)?zi=interp2(i,ntimes)：在每兩個元素之間插入插值，按ntimes遞歸擴展z，interp2(Z)等同于interp2(Z,1)。例7-4-4

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

解采集的數(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處的平板溫度同時，我們可以在三維空間立體觀察溫度分布情況。首先在三維坐標(biāo)畫出原始數(shù)據(jù)，觀察該原始數(shù)據(jù)表示的分布情況。

>>mesh(width,depth,temps)

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

>>title(['深度depth='num2str(d)'處的溫度']),axis('ij')圖7-9平板溫度的立體分布然后在兩個方向上插值，以平滑數(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ù)用一個整數(shù)系數(shù)(插值)增加采樣率，插值的方法是在原序列中插入一些0值，以提高采樣率。語法如下：

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

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

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

例7-4-5

抽樣信號的插值。

解

t=0:0.001:1;%時間向量

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

y=interp(x,4);

stem(x(1:30));

title('原抽樣信號');

figure

stem(y(1:120));

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

1．基本特征

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

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

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

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

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

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

2．分段多項式

在最簡單的用法中，spline()獲取數(shù)據(jù)x和y以及期望值xi，尋找擬合x和y的三次樣條內(nèi)插多項式，然后計算這些多項式，對每個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)插值的情況。不過，如果用戶需要多組內(nèi)插值數(shù)據(jù)時，例如需要從相同數(shù)據(jù)集里獲取另一組內(nèi)插值，再次計算三次樣條系數(shù)是沒有意義的。在這種情況下，可以調(diào)用僅帶前兩個參量的spline()：

>>pp=spline(x,y)

pp=

form:'pp'

breaks:[0123456789101112]

coefs:[12x4double]

pieces:12

order:4

dim:1上述給定的三次樣條pp形式存儲了斷點和多項式系數(shù)，以及關(guān)于三次樣條表示的其他信息。因為所有信息都被存儲在單個向量里，所以這種形式在MATLAB中是一種方便的數(shù)據(jù)結(jié)構(gòu)。

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

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

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

>>xi2=linspace(10,12);

>>yi2=ppval(pp,xi2);

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

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

例如，在限定的區(qū)間[10,15]內(nèi)，再次計算該三次樣條，注意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重新計算區(qū)域外樣條插值擬合結(jié)果

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

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

[breaks,coefs,l,k,d]=unmkpp(pp)：從分段多項式pp中提取斷點breaks、系數(shù)coefs、分段數(shù)l、階次k和維數(shù)d。運用上述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是斷點，coefs是矩陣，它的第i行是第i個三次多項式，npolys是多項式的分段數(shù)目，ncoefs是每個多項式系數(shù)的數(shù)目(階次)。注意，這種形式非常一般，樣條多項式不必是三次，這對于樣條的積分和微分是很有用的。

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

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

(2)?pp=mkpp(breaks,coefs,d)：指示的分段多項式pp是d值向量，即其系數(shù)的每個值是長度為d的向量，breaks是長度L+1的增向量。coefs是d?×?L?×?k?數(shù)組，coefs(r,i,:)包含分段多項式的第r個元素、第i個多項式段的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)解，也就是在某個區(qū)間內(nèi)有條件約束或者無條件約束地找到函數(shù)的最大值或者最小值。

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

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

MATLAB沒有求最大值的命令，因為f(x)的最大值等于-f(x)的最小值，所以利用相反函數(shù)h(x)=-f(x)的最小值可以求得最大值，故fminbnd()可用來求最小值和最大值。

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為可選項，用于指定程序的其他參數(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附近有一個最大值，并且在tmin=4附近有一個最小值。而這些點的解析值為：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

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

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

tmin=3.9270

ymin=-0.0279

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

MATLAB提供的函數(shù)fminsearch()用于計算多元函數(shù)的極小值。fminsearch()函數(shù)內(nèi)部應(yīng)用了Nelder-Mead單一搜索算法，通過調(diào)整x的各個元素的值來尋找f(x)的極小值。該算法雖然對于平滑函數(shù)搜索效率沒有其他算法高，但它不需要梯度信息，從而擴展了其應(yīng)用范圍。因此，該算法特別適用于不太平滑、難以計算梯度信息或梯度信息價值不大的函數(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ù)是一個經(jīng)典的例子，用于測試fminsearch()函數(shù)計算多元函數(shù)的極小值。banana()函數(shù)的定義為可知函數(shù)的最小值在[11]點，最小值為0，傳統(tǒng)的起始點選在[-1.21]。定義一個匿名函數(shù)banana如下：

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

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

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

x=

1.00001.0000

fval=

8.1777e-010

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

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

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

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

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

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

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

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

例7-5-2

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

解

(1)定義一個匿名函數(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ù)的零點。為了說明該函數(shù)的使用，讓我們再運用humps例子。

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

xzero=1.2995

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

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

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

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

1．trapz()函數(shù)

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

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

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

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

例7-5-3

使用梯形法計算積分。

解函數(shù)。

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

>>Y=sin(X);

>>Z=trapz(X,Y)

Z=1.9998

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

MATLAB中一元函數(shù)的積分可以用quad()和quadl()兩個函數(shù)來實現(xiàn)。函數(shù)quad()和quad8()是基于數(shù)學(xué)上的正方形概念來計算函數(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)方法計算函數(shù)fun()在區(qū)間[a,b]的定積分，其精確度為1e-6，fun可以是M文件，也可以是匿名函數(shù)，即定義為(7.5.1)例如，上例對正弦函數(shù)的積分：

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

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

F=2.0000

(2)

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

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

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

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

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

area=29.3262

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

area=29.3262

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

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

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

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

4．二重積分和三重積分

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

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

dblquad()函數(shù)通過調(diào)用quad()函數(shù)來計算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，默認值是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