MATLAB的其他函數(shù)庫

MATLAB的其他函數(shù)庫4.1數(shù)據(jù)分析和傅里葉變換函數(shù)庫4.2矩陣的分解與變換函數(shù)庫4.3多項式函數(shù)庫4.4函數(shù)功能和數(shù)值分析函數(shù)庫4.5字符串函數(shù)庫4.6符號數(shù)學(xué)函數(shù)庫4.7系統(tǒng)仿真函數(shù)庫4.1數(shù)據(jù)分析和傅里葉變換函數(shù)庫4.1.1基本的數(shù)據(jù)分析

MATLAB的基本數(shù)據(jù)處理功能是按列向進(jìn)行的，因此要求待處理的數(shù)據(jù)矩陣按列向分類,而不同的行向則表示數(shù)據(jù)的不同樣本。例如，10個學(xué)生的身高及三門課程分?jǐn)?shù)列表如下:data=15449836715899817515510068861456375961456375961415565751555664851478987771479654100145607667進(jìn)行簡單數(shù)據(jù)處理的命令見表4-1。其中大部分命令的意義很明確，不需解釋。

std標(biāo)準(zhǔn)差是指列中N個元素與該列平均值之差的平方和按N-1點取平均值的開方，即

數(shù)列求和命令sum相當(dāng)于矩形法求和，可用來近似求積分，用梯形法求和命令trapz求積分更為精確。梯形法是把相鄰兩點數(shù)據(jù)的平均值作為數(shù)據(jù)點，十個數(shù)據(jù)只能產(chǎn)生9個數(shù)據(jù)點。如果數(shù)組長度為N，則sum和trapz的關(guān)系為

trapz(data)=sum(data)-0.5(data(1)+data(N))其差額為半個首點和半個末點的數(shù)據(jù)和，將相加以后的結(jié)果乘以步長，才近似表示了這些數(shù)據(jù)包絡(luò)下的面積。有些數(shù)據(jù)處理命令的結(jié)果不是一個標(biāo)量而是一個列向量，為了節(jié)省篇幅，我們只取數(shù)據(jù)中的前三行，其結(jié)果見表4-2。注意其結(jié)果一般與原數(shù)據(jù)具有同樣的行數(shù)，只有求差分(diff)才會減少一行，因為它是求相鄰行之間的差。另外，cumtrapz函數(shù)是用梯形法累計求面積,和trapz相仿，它也會使數(shù)據(jù)長度減少一。4.1.2用于場論的數(shù)據(jù)分析函數(shù)用于場論的命令有以下幾個：

·gradient：用來求二維和三維場的近似梯度，例如根據(jù)電位分布求電場就可用這個函數(shù)。

·del2：是二維和三維場的拉普拉斯算子。

·cross：為兩個向量的矢量積。

·dot：為兩個向量的數(shù)量積。設(shè)i，j，k為沿x，y，z方向的單位向量，則對于兩向量a=axi+ayj+azk和b=bxi+byj+bzk而言：向量的矢量積為(叉乘)a×b=(aybz-azby)i+(azbx-axbz)j+(axby-aybx)k向量的數(shù)量積為(點乘)a·b=axbx+ayby+azbz

在MATLAB中這兩個向量可表為: a=［ax，ay，az］;b=［bx，by，bz］; cross(a，b)=[ay*bz-az*by，az*bx-ax*bz，ax*by-ay*bx]; dot(a，b)=a*b′;4.1.3用于隨機(jī)數(shù)據(jù)分析的函數(shù)

MATLAB有兩個產(chǎn)生隨機(jī)數(shù)的命令:一個是rand(m，n)，它產(chǎn)生在0與1之間均勻分布的m行n列隨機(jī)數(shù)矩陣，其均值為0.5，標(biāo)準(zhǔn)差(或均方根差)為0.2887;另一個是randn(m，n)，它產(chǎn)生正態(tài)分布的m行n列隨機(jī)數(shù)矩陣，其均值為0，標(biāo)準(zhǔn)差為1。其分布情況可用直方圖命令hist(x，N)來顯示，其中N表示直方圖橫坐標(biāo)的分割數(shù)，默認(rèn)值為10。例如

x=rand(1，1000);hist(x) y=randn(1，1000);hist(y，50)得出的兩組圖形分別如圖4-1(a)、(b)所示。hist(x)是把1000個x中處于0~0.1，0.1~0.2，…，0.9~1.0各個區(qū)域中的數(shù)目分別清點出來，畫成直方圖。如果x真是均勻分布的，那么這個圖應(yīng)該是水平直線。實際上，隨機(jī)數(shù)規(guī)律是按統(tǒng)計方法確定的，所以各區(qū)域的數(shù)量仍參差不齊，只有數(shù)據(jù)量無限增加時，此規(guī)律才越來越明顯地表現(xiàn)出來。hist(y，50)則把y的最小值和最大值之間分成50份進(jìn)行統(tǒng)計，得到一個鐘形的，即正態(tài)分布的曲線。圖4-1均勻分布與正態(tài)分布隨機(jī)數(shù)直方圖4.1.4用于相關(guān)分析和傅里葉分析的函數(shù)相關(guān)分析(包括卷積)和傅里葉分析分別用于信號的時域和頻域處理。這里雖然只給出了十幾個函數(shù)，實際上它們是整個信號處理計算的基礎(chǔ)。

(1)corrcoef給出兩個同長信號的相關(guān)系數(shù)，例如對前面兩個隨機(jī)序列，鍵入

R=corrcoef(x，y)得

R=［1-0.0508 -0.05081］主對角線上是x和y的自相關(guān)系數(shù)，它必定為1；此處可看出x和y的互相關(guān)很弱。

(2)cov(x，y)給出x、y的協(xié)方差矩陣，對上述x、y，有

cov(x，y)=［0.0785-0.0148 -0.01481.0782］其主對角線上的值分別為x和y的均方差，即標(biāo)準(zhǔn)差的平方(因為是隨機(jī)數(shù)，所以它不會嚴(yán)格等于理論值)。

(3)conv(x，y)給出x，y的卷積。如果x是輸入信號，y是線性系統(tǒng)的脈沖過渡函數(shù)，則x，y的卷積就給出系統(tǒng)的輸出信號。卷積函數(shù)也用于多項式相乘，見4.3.1小節(jié)。

(4)filter(b，a，x)也是根據(jù)輸入信號x和線性系統(tǒng)特性求輸出信號的函數(shù)。它與conv的不同在于系統(tǒng)的特性是以傳遞函數(shù)的分子多項式系數(shù)向量b和分母多項式系數(shù)向量a給出，而不是以脈沖過渡函數(shù)的形式給出的。

(5)X=fft(x，N)求出時域信號x的離散傅里葉變換X。N為規(guī)定的點數(shù)。N的默認(rèn)值為所給x的長度。當(dāng)N取2的整數(shù)冪時變換的速度最快。通常取大于又最靠近x的整數(shù)冪，即令N＝2^nextpow2(length(x))。例如x的長度為12，nextpow2(12)=4，N=2^4=16，多出的各點補以零。一般情況下，fft求出的函數(shù)為復(fù)數(shù)，可用abs及angle分別求其幅度和相位。在畫頻譜圖時往往最關(guān)心其幅頻特性。

【例4-1-1】給出一個信號t=0:.001:3;u=sin(300*t)+2*cos(200*t)它的幅頻特性可用下列語句求得:U=fft(u);plot(abs(U))得出的頻譜曲線如圖4-2(a)所示，它對采樣頻率呈對稱形式。為了看得更清楚，把坐標(biāo)間隔縮小，鍵入axis(［0，300，0，3000］)，得出的頻譜曲線如圖4-2(b)所示。圖4-2例4-1-1中信號的幅頻曲線

(6)x=ifft(X)為傅里葉反變換函數(shù)，其用法與fft相仿。

(7)sound(u，s)會在音箱中產(chǎn)生u所對應(yīng)的聲音。s規(guī)定重放的速度，其缺省值為8192(b/s)。數(shù)據(jù)分析和傅里葉變換函數(shù)庫見表4-3。4.2矩陣的分解與變換函數(shù)庫4.2.1線性方程組的系數(shù)矩陣在2.2節(jié)中我們提到了可以用矩陣除法來解線性方程組，本節(jié)將討論有關(guān)解線性方程組的一些深入的問題及其工具函數(shù)。這些函數(shù)見表4-4中的矩陣分析和線性方程部分。

det(a)用以求方陣a的行列式。若det(a)不等于零，則a的逆陣inv(a)存在。線性方程組的系數(shù)矩陣只有滿足這個條件，它的解才存在。rank(a)用以求任意矩陣a的秩，也就是它所能劃分出的行列式不為零的最大方陣的邊長。trace(a)求出矩陣主對角線上元素的和。如果det(a)雖不等于零，但數(shù)值很小，近似于零，則這樣的線性方程組稱為病態(tài)的線性方程組，其解的精度比較低。為了評價線性方程組系數(shù)矩陣的病態(tài)程度，用了條件數(shù)(ConditionNumber)的概念。條件數(shù)愈大，方程病態(tài)愈重，解的精度愈低。當(dāng)系數(shù)矩陣的條件數(shù)很大(達(dá)到1016)而又拿它作除數(shù)(即解此線性方程組)時，MATLAB會提出警告:“條件數(shù)太大，結(jié)果可能不準(zhǔn)確”，求條件數(shù)的函數(shù)為cond(a)，a可以不是方陣。在線性方程組A*x=B中，系數(shù)矩陣A的行數(shù)m表示方程的數(shù)目，其列數(shù)n表示未知數(shù)的數(shù)目。在正常情況下，方程數(shù)等于未知數(shù)數(shù)，即n=m，A為方陣，A＼B意味著inv(A)*B。實際上，對于方程數(shù)大于未知數(shù)數(shù)(m>n)的超定方程組，以及方程數(shù)小于未知數(shù)數(shù)(m<n)的不定方程組，MATLAB中A＼B的算式都仍然合法。前者是最小二乘解;而后者則是令x中n-m個元素為零的一個特殊解。對前一種情況下，因為A不是方陣，其逆inv(A)不存在，解的MATLAB算式均為x=inv(A′*A)*(A′*B)。把pinv(A)=inv(A′*A)*A′定義為廣義逆函數(shù)，則A＼B就等于pinv(A)*B。因此，MATLAB中引入了廣義逆的概念，除數(shù)A就可以不是方陣。

【例4-2-1】求下列矩陣的行列式及逆陣等特性。

a=2900041475517874得

det(a)=-275rank(a)=4inv(a)=-0.07270.42550.7855-0.62180.1273-0.0945-0.17450.13820.0000-0.6000-0.80000.8000-0.12730.49450.3745-0.3382trace(a)=15cond(a)=33.47634.2.2矩陣的分解矩陣可以分解為幾個具有特殊構(gòu)造性質(zhì)的矩陣的乘積，這是分析矩陣的一種重要手段，而這種分解在數(shù)學(xué)計算上通常又是非常繁瑣的工作。MATLAB提供了一些現(xiàn)成的函數(shù)可供調(diào)用，主要有以下三種。

(1)三角分解(lu分解)。它把一個任意矩陣分解為一個準(zhǔn)下三角矩陣和一個上三角矩陣的乘積。由于此函數(shù)有兩個輸出矩陣，因此其左端應(yīng)有兩個變量l和u，鍵入 ［l，u］=lu(a)得

l=0.28571.00000000.52830.68381.00001.00000001.00000.39621.00000u=7.00005.00005.00001.000007.5714-1.4286-0.2857002.56603.11320002.0221

l被稱為準(zhǔn)下三角陣，因為它必須交換幾行才能成為真的下三角陣，它的行列式絕對值等于1(可正可負(fù));上三角陣u的行列式等于a的行列式。

(2)正交分解(qr分解)。qr(a)把任意矩陣a分解為一個正交方陣q和一個與a有同樣階數(shù)的上三角矩陣r的乘積。該方陣q的邊長為m，且其行列式的值為1。

【例4-2-2】求下列3×5矩陣b的qr分解。

b=0.21900.67930.51940.05350.00770.04700.93470.83100.52970.38340.67890.38350.03460.67110.0668鍵入 ［q，r］=qr(b)得

q=-0.3063-0.4667-0.8297-0.0658-0.85910.5076-0.94970.21010.2324r=-0.7149-0.6338-0.2466-0.6886-0.09110-1.0395-0.9490-0.3390-0.318900-0.00110.38050.2038

(3)奇異值分解(svd分解)。svd(a)把任意m×n矩陣a分解為三個矩陣的乘積，即a=u*s*v。其中u，v分別為m×m和n×n的正交方陣，s則為m×n的矩陣。其左上方為一對角方陣，對角線上的元素就是矩陣a的奇異值，其余元素全為零。例如對上述b作奇異值分解，鍵入 ［u，s，v］=svd(b)得

u=0.46230.22730.8571 0.78220.3507-0.5149 0.4176-0.90850.0157s=1.7539000000.7995000000.353400v=0.2403-0.68850.4927-0.474800.68720.16740.30270.41930.48190.51580.47290.0506-0.3723-0.60760.4102-0.5151-0.61220.3194-0.29940.18900.0944-0.5369-0.59850.5558矩陣最大奇異值和最小奇異值之比就是它的條件數(shù)，即

cond(b)=max(diag(s))/min(diag(s))4.2.3矩陣的特征值分析

eig(a)用來求方陣a的特征根和特征向量，其輸出有兩個，即特征向量e和特征根r。鍵入 ［e，r］=eig(a)得

e=-0.2568-0.3834+0.4681i-0.3834-0.4681i0.6167-0.3481-0.2177-0.2869i-0.2177+0.2869i-0.1850-0.46820.5152+0.2228i0.5152-0.2228i-0.6624-0.77050.4217-0.1060i0.4217+0.1060i0.3829r=14.2004 0 0 00 0.7495+5.2088i 0 0 0 0 0.7495-5.2088i0 0 0 0-0.6993特征根是特征方程的根，矩陣的特征方程系數(shù)可用poly函數(shù)求出，用roots命令可求出其特征根。例如，鍵入

p=poly(a)得

p=1.0000-15.000038.0000-359.0000-275.0000而

roots(p)=14.20040.7495+5.2088i0.7495-5.2088i-0.6993結(jié)果與eig函數(shù)求出的特征根相同，但roots函數(shù)不能求特征向量。4.2.4特殊矩陣庫有一些特殊構(gòu)造的矩陣在矩陣變換中很有用處，它們組成一個專門的函數(shù)庫。讀者在應(yīng)用中遇到有關(guān)的矩陣，即可在此調(diào)用。其調(diào)用的參數(shù)和調(diào)用方法均可從help文本中查得，特殊函數(shù)庫見表4-5。4.3多項式函數(shù)庫一元高次代數(shù)多項式a(x)=a1xn+a2xn-1+…+anx+an+1，在MATLAB中可以用它的系數(shù)向量

a=［a(1)，a(2)，…，a(n)，a(n+1)］來表示。其方次已隱含在系數(shù)元素離向量右端元素的下標(biāo)差中。要注意，如果x的某次冪的系數(shù)為零，這個零必須列入系數(shù)向量。在微分方程中，通過運算微積可把線性系統(tǒng)的特性表示為兩個算子s的多項式之比。因此多項式在近代信息和控制理論中有著十分重要的地位。MATLAB的多項式和插值函數(shù)庫見表4-6。4.3.1多項式的四則運算

【例4-3-1】設(shè)有兩個多項式a(x)=2x3+4x2+6x+8及b(x)=3x2+6x+9，要求對此兩個多項式作如下運算:

(1)多項式相乘(conv):conv函數(shù)本來是卷積(Convolution)的意思，但它也符合多項式相乘的運算規(guī)則。可以想象把系數(shù)向量a正常排列，而把b反轉(zhuǎn)，先將a(1)與b(1)

對齊: a(1)a(2)a(3)a(4)b(3)b(2)b(1)把上下對應(yīng)的項相乘，a(1)*b(1)得出多項式乘積c的最高次項系數(shù)c(1);把b右移一位，把上下對應(yīng)的項相乘并求和，a(1)*b(2)+a(2)*b(1)為次高次項系數(shù)c(2);依此類推，可得到乘積c的全部系數(shù)。這種運算和卷積運算的規(guī)則完全相同，故也用conv。鍵入

a=［2，4，6，8］，b=［3，6，9］，c=conv(a，b)得

a=2468 b=369 c=624609610272

(2)多項式相加:MATLAB規(guī)定，只有長度相同的向量才能相加，因此必須在短的向量前面補以若干個零元素，才能用MATLAB的矩陣加法運算符。鍵入

d=a+［0，b］得

d=271217這種手工數(shù)兩個多項式的長度再補零的方法是不可取的，最好讓計算機(jī)自動完成。為此可編一個子程序polyadd.m，其內(nèi)容為:functiony=polyadd(x1，x2)n1=length(x1);n2=length(x2);

ifn1>n2x2=［zeros(1，n1-n2)，x2］;elseifn1<n2x1=［zeros(1，n2-n1)，x1］;end，y=x1+x2;這樣，多項式相加就可寫成:c=polyadd(a，b)，相減可另編一個子程序，或在polyadd的輸入變元中加負(fù)號來實現(xiàn)。

(3)多項式相除:相除是相乘的逆運算，用deconv實現(xiàn)。但除法不一定除得盡，會有余子式，因此鍵入

［q，r］=deconv(c，a) q=369 r=000000其中，q是商式，r是余子式。因為用的是相乘的數(shù)據(jù)a和c，恰好除凈。如令a1=a+1，則有: a1=3579

［q1，r1］=deconv(c，a1) q1=2.00004.66677.5556 r1=0007.55567.11114.0000余式為7.5556x2+7.1111x+4，可以用商式與除式相乘，再加上余式的方法來檢驗:

c1=conv(q1，a1)+r1得

c1=624609610272與c相同。4.3.2多項式求導(dǎo)、求根和求值

(1)多項式求導(dǎo)數(shù)(polyder):鍵入

e=polyder(c)得

e=3096180192102

(2)多項式求根(roots和poly函數(shù)):鍵入

ra=roots(a);rb=roots(b);rc=roots(c);ra，rb，rc得

ra=-1.6506-0.1747+1.5469i-0.1747-1.5469i rb=-1.0000+1.4142i-1.0000-1.4142irc是ra和rb的并集，這是完全可以預(yù)計到的，為節(jié)省篇幅，不再列出。由根求多項式系數(shù)是roots的逆運算，其函數(shù)名也是poly，有

a=poly(ra);b=poly(rb)從poly函數(shù)的用法可以看出MATLAB的智能特點:當(dāng)a是向量時，poly把它看做根來組成多項式;當(dāng)a是方陣時，poly用它組成方陣的特征多項式。

(3)多項式求值(polyval):將多項式a中的自變量x賦予值xv時，該多項式的值可用

F=polyval(a，xv)求得，其中xv可以是復(fù)數(shù)，而且可以是矩陣或數(shù)組，此時polyval對輸入變元作元素群運算，這對于求線性系統(tǒng)的頻率特性特別方便。polyvalm則對輸入的變元陣方陣作矩陣多項式運算。

【例4-3-2】設(shè)a為系統(tǒng)分母系數(shù)向量，b為系統(tǒng)分子系數(shù)向量，求此系統(tǒng)的頻率響應(yīng)并畫出頻率特性。先令頻率數(shù)組w取線性間隔: w=linspace(0,10); %在w等于0~10之間按線性間隔取100點(默認(rèn)值) A=polyval(a，j*w); %分別求分母分子多項式的值(為復(fù)數(shù)數(shù)組) B=polyval(b，j*w); subplot(2,1,1);plot(w，abs(B./A)) %畫兩者元素群相除所得的幅頻特性

subplot(2,1,2);plot(w，angle(B./A)) %畫相頻特性頻率特性通常在對數(shù)坐標(biāo)中繪制。因此輸入頻率數(shù)組取對數(shù)等間隔: w1=logspace(-1,1) %在w1從10-1~10之間，按對數(shù)分割為50點(默認(rèn)值) F=polyval(b,j*w1)./polyval(a，j*w1); %求出這些點上的頻率響應(yīng)(復(fù)數(shù)) subplot(2,1,1)，loglog(w1，abs(F)) %在雙對數(shù)坐標(biāo)中畫出幅頻特性

subplot(2,1,2);semilogx(w1，angle(F)) %在雙對數(shù)坐標(biāo)(x)中畫出相頻特性所得曲線如圖4-3所示。圖4-3線性坐標(biāo)和對數(shù)坐標(biāo)中的頻率特性(a)線性坐標(biāo)中的頻率特性；(b)對數(shù)坐標(biāo)中的頻率特性4.3.3多項式擬合

p=polyfit(x，y，n)用于多項式曲線擬合。其中x，y是已知的N個數(shù)據(jù)點坐標(biāo)向量，當(dāng)然其長度均為N。n是用來擬合的多項式次數(shù)，p是求出的多項式的系數(shù)，n次多項式應(yīng)該有n+1個系數(shù)，故p的長度為n+1。擬合的準(zhǔn)則是最小二乘法。

【例4-3-3】設(shè)原始數(shù)據(jù)為x在11個點上測得的y值: x=0:0.1:1; y=［-0.447，1.978，3.28，6.16，7.08，7.34，7.66，9.56，9.48，9.30，11.2］;線性擬合:a1=polyfit(x，y，1);求出a1后，可求出xi=linspace(0，1)及100個點上的yi1值并繪圖: yi1=polyval(a1，xi);plot(x，y，‘o’，xi，yi1，‘b’)，pause其中，原始數(shù)據(jù)用圓圈標(biāo)出，而擬合曲線為藍(lán)色。依此類推，有:二次擬合:a2=polyfit(x,y,2);yi2=polyval(a2,xi); plot(x,y,′o′,xi,yi2,′m′)三次擬合:a3=polyfit(x,y,3);yi3=polyval(a3,xi); plot(x,y,′o′,xi,yi3,′r′)九次擬合:a9=polyfit(x,y,9);yi9=polyval(a9,xi); plot(x,y,′o′,xi,yi9,′c′)十次擬合:a10=polyfit(x,y,10);yi10=polyval(a10,xi); plot(x,y,′o′,xi,yi10,′g′)所得的曲線如圖4-4所示。給定11點的最大擬合階次為10，此時擬合曲線將通過全部給定點?？梢钥闯?，擬合曲線的階次太高會造成曲線振蕩，反而看不出函數(shù)關(guān)系的基本規(guī)律，并不一定好。圖4-4不同的逼近次數(shù)產(chǎn)生的不同曲線4.3.4多項式插值插值和擬合的不同在于:(1)插值函數(shù)通常是分段的，因而人們關(guān)心的不是函數(shù)的表達(dá)式，而是插值得的數(shù)據(jù)點;(2)插值函數(shù)應(yīng)通過給定的數(shù)據(jù)點x、y。插值函數(shù)一般可表示為

yi=interp1(x，y，xi，′method′)其中，xi為插值范圍內(nèi)的任意點集x坐標(biāo)，yi是插值后的對應(yīng)數(shù)據(jù)點集的y坐標(biāo)?！鋗ethod′為插值函數(shù)的類型選項，有′linear′(線性，默認(rèn)項)、′cubic′(三次)和′cubicspline′(三次樣條)等三種。

(1)一維插值函數(shù)interp1。

【例4-3-4】仍取例4-3-3中的x、y、xi，求其線性和三次插值曲線。線性插值:yi1=interp1(x，y，xi);plot(x，y，′o′，xi，yi1)三次插值:yi2=interp1(x，y，xi，′spline′);plot(x，y，′o′，xi，yi2，′g′)所得曲線如圖4-5及圖4-6所示，三次插值的結(jié)果比較光滑。圖4-5線性插值曲線圖4-6三次插值曲線

(2)二維插值函數(shù)zi=interp2(x，y，z，xi，yi，′method′)，其變元的意義可以類推。

【例4-3-5】已知某矩形溫箱中3×5個測試點上的溫度，求全箱的溫度分布。

解：給定width=1:5;depth=1:3; temps=［8281808284;7963616581;8484828586］;計算沿寬度和深度細(xì)分網(wǎng)格di=1:0.2:3;wi=1:0.2:5;交點上的各點溫度。程序如下: tc=interp2(width，depth，temps，wi，di′，′cubic′); %求各點溫度

mesh(wi，di，tc) %畫三維曲面所得溫度分布圖形如圖4-7所示。注意interp2中所用的wi和di′是寬度和深度方向的細(xì)分坐標(biāo)向量，di必須變換為列向量，插值函數(shù)運算時會自動將它們轉(zhuǎn)變?yōu)閷挾瘸松疃绕矫嫔系木W(wǎng)格，并計算網(wǎng)格點上的溫度。圖4-7二維插值的曲面4.3.5線性微分方程的解線性常微分方程的解可用拉普拉斯算子s表示為Y(s)=B(s)/A(s)其中，B(s)和A(s)都是s的多項式，分母多項式的次數(shù)n通常高于分子多項式的次數(shù)m。在時間域的解y(t)是Y(s)的拉普拉斯反變換。求反變換的重要方法之一是部分分式法，即將上述多項式分解為多個s的一次分式之和。表4-6中的留數(shù)函數(shù)residue可以完成這一任務(wù)。步驟為:

(1)用［r，p，k］=residue(b，a)求出Y(s)的極點數(shù)組p和留數(shù)數(shù)組r，因而Y(s)可表示為

(2)此時可以很簡單地求出它的反變換

y(t)=r(1)*exp(p(1)*t)+r(2)*exp(p(2)*t)+r(3)*exp(p(3)*t) +r(4)*exp(p(4)*t)+…

【例4-3-6】求解線性常微分方程+5y″+4y′+7y=3u″+0.5u′+4u在輸入u(t)為單位脈沖及單位階躍信號時的解析解。

解:

用Laplace變換(脈沖輸入u(s)=1，階躍輸入u(s)=1/s)求輸出y(s)，得再對y(s)求反變換，得y(t)。

(1)在脈沖輸入時的響應(yīng):a=［1，5，4，7］;b=［3，0.5，4］;［r，p，k］=residue(b，a)得

r=3.2288-0.1144+0.0730i-0.1144-0.0730ip=-4.4548-0.2726+1.2235i-0.2726-1.2235ik=［］求時域解，先設(shè)定時間數(shù)組t=0:0.2:10;然后列出

yi=r(1)*exp(p(1)*t)+r(2)*exp(p(2)*t)+r(3)*exp(p(3)*t);plot(t,yi)所得曲線如圖4-8(a)所示。

(2)在階躍輸入時的響應(yīng):此時由于分母乘了個s，因此a將提高一階，右端多加一個零。

a=［1,5,4,7,0］;b=［3,0.5,4］;［r,p,k］=residue(b，a)

r=-0.72480.0767+0.0764i0.0767-0.0764i0.5714+0.0000ip=-4.4548-0.2726+1.2235i-0.2726-1.2235i0k=［］

ys=r(1)*exp(p(1)*t)+r(2)*exp(p(2)*t)+r(3)*exp(p(3)*t)+r(4);plot(t，ys)所得曲線如圖4-8(b)所示。圖4-8脈沖過渡響應(yīng)和階躍過渡響應(yīng)4.4函數(shù)功能和數(shù)值分析函數(shù)庫4.4.1本函數(shù)庫的主要子程序本函數(shù)庫的主要子程序見表4-7，我們將它分為兩類來探討。第一類是對任意非線性函數(shù)的分析，包括求極值、過零點等;第二類是求任意函數(shù)的數(shù)值積分，包括定積分和微分方程的數(shù)值解等。它們的共同特點是必須會自行定義函數(shù)。這在本書2.6節(jié)中已經(jīng)學(xué)過。就取其中的函數(shù)humps.m來說明這些子程序的用法，它定義了下列非線性函數(shù)，即4.4.2非線性函數(shù)的分析

1)繪制函數(shù)曲線其格式為:fplot(′函數(shù)名′，［初值x0，終值xf］)，例如，要畫出humps函數(shù)在x=0~2之間的曲線，鍵入

fplot(′humps′，［0，2］)，grid得出圖4-9所示的曲線。

fplot函數(shù)對于快速了解一些復(fù)雜特殊函數(shù)的波形很有用處。例如求第一類bessel函數(shù)(見表4-5)，可用

fplot(′besselj(alpha，x)′，［0，10］)設(shè)alpha為1，2，5時所得曲線如圖4-10所示。圖4-9humps函數(shù)的曲線圖4-10第一類bessel函數(shù)的曲線

2)求函數(shù)極值其格式為:fminbnd(′函數(shù)名′，初值x0，終值xf)，例如，求humps函數(shù)在x=0～1.5之間的極小值，鍵入

m=fminbnd(′humps′，0，1.5)得

m=0.6370

3)求函數(shù)零點其格式為:fzero(′函數(shù)名′，初猜值x0)，例如，求humps函數(shù)在x=1附近的過零點，鍵入

z=fzero(′humps′，1)得

z=1.2995以上給出的是這些函數(shù)調(diào)用的典型格式，還有其他選項可作為變元，例如

fplot(′tan′，［-2*pi2*pi-2*pi2*pi］，′*′)，grid在第二項變元中增加了y軸的上下限，第三項變元是線型。所得圖形如圖4-11(a)所示，讀者可從helpfplot中得到進(jìn)一步的信息。圖4-11由fplot和ezplot畫出的曲線還有一個簡便畫函數(shù)圖的命令ezplot(讀作easyplot)，它連自變量范圍都無需規(guī)定，其默認(rèn)的自變量范圍為［-2π，2π］，因此只要鍵入

ezplottan(x)，grid也可得到類似于圖4-11(a)的曲線，只是*號變?yōu)榱藢嵕€。若鍵入

ezplottan(sin(x))-sin(tan(x))所得圖形如圖4-11(b)所示?？梢钥闯?，圖上還自動作出了標(biāo)注。4.4.3函數(shù)和微分方程的數(shù)值積分

1)定積分子程序(quad、quad8、quadl等)其格式為:quad(′函數(shù)名′，初值x0，終值xf)，例如，求humps函數(shù)在x=1～2之間的定積分，鍵入

s=quad(′humps′，1，2)得

s=-0.5321不難用定積分函數(shù)來求不定積分的數(shù)值解。只要固定積分下限，用for循環(huán)，把積分上限逐步增加即可。例如要求humps函數(shù)以x=0為下限的不定積分，可編寫下列程序:fori=1:20x(i)=0.1*i;y(i)=quad(′humps′，0，x(i));end，plot(x，y)得出的曲線如圖4-12所示。可以與圖4-9對照確認(rèn)它是humps曲線的積分。

quad、quad8、quadl等調(diào)用格式相同，只是內(nèi)部算法不同。圖4-12humps函數(shù)的積分曲線

2)微分方程數(shù)字解(ode23、ode45等)如果微分方程可化為一階微分方程組的形式，即

dy/dx=f(x，y)其中，x是標(biāo)量，y可以是一個列向量。f(x，y)是以x，y為變元的函數(shù)，用MATLAB函數(shù)文件表述。設(shè)文件名為yprime.m，則求此微分方程的數(shù)值解的子程序調(diào)用格式為

[x,y]=ode23(′yprime′,自變量區(qū)間xspan,因變量初值y0)對于humps函數(shù)，不能直接用ode23作數(shù)值積分，其原因在于humps只有一個輸入變元x，微分方程數(shù)值解的函數(shù)ode23等要求被調(diào)用的函數(shù)有兩個輸入變元。如果我們把humps函數(shù)文件加一個虛的變元y，即把它的第一句換成functionyp=humps1(x，y)，并將此函數(shù)另存為一個humps1.m文件，則［x，y］=ode23(′humps1′，［0，2］，1);plot(x，y)表示在初值y0=1，xspan為［0，2］的條件下，求微分方程的數(shù)值解，則可以得到與圖4-12相仿的曲線，只是向上平移了一個單位，因為這里設(shè)y0=1。在這個例子中，函數(shù)humps1中的y只是一個虛的變元，比較簡單。

【例4-4-1】求下列微分方程(范德堡方程)的數(shù)值解，即y″+r(y2-1)y′+y=0它可寫成導(dǎo)數(shù)在左端的兩個一階微分方程構(gòu)成的方程組，也可寫成矩陣方程：其中，。先要建立反映此微分方程組右端的函數(shù)文件vdpl.m，存入子目錄user中，其內(nèi)容為：

functionyprime=vdpl(x，y) globalr%r值由主程序通過全局變量傳送

yprime=［y(2);r*(1-y(1).^2).*y(2)-y(1)］; %兩行單列向量主程序如下: globalr，r=input(′輸入r，在0<r<10之間選擇′) xspan=input(‘［x0，xf］=’);y0=input(‘y0=［y10;y20］=’);

［x，y］=ode45(′vdpl′，xspan，y0);plot(x，y)在r=2，xspan=［0，30］，y0=［1;2］條件下得出的曲線如圖4-13所示。圖4-13范德堡方程積分的曲線

ode45是高階的數(shù)值積分函數(shù)，其步長可以取得較大，并能保證較高的精度，調(diào)用方法與ode23相仿。這些函數(shù)都有自動選擇步長的功能，以保證把誤差tol控制在0.001以下。如果要改變?nèi)菰S誤差tol，可在輸入變元中增加其他選項，詳情可從help命令中獲得。此外，還有ode23s，ode15s，ode113等內(nèi)部算法各異的數(shù)值積分函數(shù)。在圖形輸出方面，則有繪制相平面曲線的命令odephas2，odephas3等。4.5字符串函數(shù)庫

MATLAB的程序和標(biāo)識符都是用字符串來表示的。每個字符有它對應(yīng)的ASCII碼。4.4節(jié)中的函數(shù)，都把字符串當(dāng)做變元來看待，從而使程序更為簡潔。有時也需要把字符串當(dāng)做數(shù)碼來處理。字符串函數(shù)庫中的命令，都是為了增強(qiáng)這一功能而建立的。對初學(xué)者，這部分并不重要，但若要編寫出人機(jī)界面優(yōu)良、能調(diào)用各種函數(shù)和文件的高級程序，字符串函數(shù)庫是必不可少的。MATLAB的字符串函數(shù)庫見表4-8。4.5.1字符串的賦值語句s=′abyzABYZ0189′把字符串賦值給變量s，其結(jié)果是

s=abyzABYZ0189而

size(s)=112說明它是以行向量的形式存儲的。當(dāng)然它內(nèi)部帶有字符串的標(biāo)志，故在屏幕上顯示出字符。要找到s所對應(yīng)的ASCII碼，可用abs命令: abs(s)=97981211226566899048 495657從中可以知道英文大小寫字母和數(shù)字的十進(jìn)制ASCII碼值，還可用setstr命令作逆向變換: setstr(abs(s))=abyzABYZ0189要求出字母和數(shù)字的十六進(jìn)制ASCII碼值，可用dec2hex命令: dec2hex(abs(s))=6162797A4142595A30313839

MATLAB顯示時并沒有各碼之間的空格，這里加上空格是為了便于讀者閱讀?？梢园褞讉€字符串沿行向串接，構(gòu)成更長的字符串，如鍵入

s1=［′welcome′，s］得

s1=welcomeabyzABYZ0189可以把幾個長度相同的字符串沿列向并列，組成一個字符串矩陣，如

s2=［′a=5′;′b=2′;′c=a+b*b′］這時必須在前兩個字符串中增添若干個空格，保證三個字符串長度均為7，否則賦值無效。4.5.2字符串語句的執(zhí)行如果字符串的內(nèi)容是MATLAB語句，如上述的s2，則可以用eval命令來執(zhí)行它。如鍵入

fork=1:3eval(s2(k，:))，end結(jié)果為a=5，b=2，c=9。在編寫MATLAB的演示程序時，往往要通過人機(jī)交互，讓用戶輸入某種表達(dá)式(而不只是數(shù)據(jù))，然后按此表達(dá)式執(zhí)行，這種程序的格式為

st=input(′s=表達(dá)式′，′s′);eval(st)

input語句中的′s′表示把輸入當(dāng)作字符串來接受，因此用戶鍵入的字符串就不必加引號了。下面的例子說明了如何將eval函數(shù)和load函數(shù)一起使用，讀出10個具有連續(xù)文件名mydata1,mydata2，…，mydata10的數(shù)據(jù)文件: fori=1:10fname=′mydata′;eval(［′load′，fname，int2str(i)］)，end特別注意的是int2str(i)把數(shù)1:10轉(zhuǎn)換為字符1:10。在顯示屏上看不出兩者有什么差別,但要記住字符0:9的ASCII碼是48:57。數(shù)10只占一個MATLAB雙精度存儲單元，而字符串10卻占兩個存儲單元，并構(gòu)成一個單行兩列的矩陣。在eval后的輸入變元必須是一個構(gòu)成MATLAB語句的字符串，因此必須把三個字符串接起來，并且不要忘掉在load后面加一個空格。eval命令是在較高級的程序中常常遇到的，要學(xué)會使用。4.5.3字符串輸入輸出前面我們一直用disp函數(shù)來進(jìn)行字符串和數(shù)據(jù)的輸出。disp(′pi=′)將顯示引號內(nèi)的字符串，此處為pi=。disp(pi)將顯示變量pi的值3.1416或其他的八種顯示格式之一，由format命令確定。想把字符串pi=和變量pi的值顯示在一行上，試用disp(′pi=′，pi)，回答這是非法的。這時應(yīng)該用sprintf函數(shù)，它可把數(shù)據(jù)按要求的格式轉(zhuǎn)換為字符串，再把它與需要顯示的字符串組裝成一個長字符串，使顯示格式非常靈活，人機(jī)界面更為友好。如鍵入 st=sprintf(′圓周率pi=%8.5f′，pi);disp(st)結(jié)果為 圓周率pi=3.14159其中，%為數(shù)據(jù)格式符，f表示十進(jìn)制浮點，8.5表示數(shù)字的長度為8位，小數(shù)點后5位。從%到f之間的字符都是不顯示的，它只指出顯示數(shù)據(jù)pi的格式。

【例4-5-1】再舉一個用sprintf的例子，它為y的兩列規(guī)定了不同的顯示格式。

x=0:10:90;y=［x;sin(x*pi/180)］;disp(sprintf (′%10.2f%12.8f＼n′，y)) 0.000.00000000 10.000.17364818

80.000.98480775 90.001.00000000

sprintf命令是從C語言中的同名命令演化來的，sscanf則是它的逆命令，相仿的還有fprintf和fscanf(見表3-3)。4.6符號數(shù)學(xué)函數(shù)庫4.6.1符號數(shù)學(xué)函數(shù)庫的主要功能符號數(shù)學(xué)函數(shù)庫是利用MATLAB的界面來調(diào)用Maple軟件的工具的。Maple是加拿大Waterloo大學(xué)開發(fā)的推理軟件，在國際上使用最早，也最為成功。從90年代初，Mathworks公司就與Maple合作，將MATLAB的數(shù)值計算與符號數(shù)學(xué)功能相結(jié)合，構(gòu)成一個Symbolic工具箱，使得在MATLAB的界面下，可兼做符號數(shù)學(xué)。符號數(shù)學(xué)，顧名思義，就是以符號(如a，b，c，x，y，z)為對象的數(shù)學(xué)。它以推理為特點,區(qū)別于以數(shù)字為對象的MATLAB基本部分。符號數(shù)學(xué)是每門課都要用到的，因此Mathworks公司提供給大學(xué)生的版本(StudentEditionofMATLAB)中就包括了這個工具箱。國外的高等數(shù)學(xué)教科書中，都有部分例題和習(xí)題會用到這個工具箱。本書雖然要對它作簡要的介紹，但在“數(shù)學(xué)篇”中并沒有把這個工具箱作為重點，這是基于以下幾方面原因。

(1)學(xué)生從小學(xué)用筆算加減乘除，到中學(xué)用計算器，做了大量的四則運算和函數(shù)運算，確實不應(yīng)該在大學(xué)教學(xué)中再重復(fù)使用低級的老方法。把數(shù)值計算交給計算機(jī)做，不僅絕大部分大學(xué)老師能夠接受，而且也是他們應(yīng)負(fù)的責(zé)任。推理是一項基本功，學(xué)生在中小學(xué)中并未得到太多的訓(xùn)練，所以應(yīng)該在大學(xué)階段來學(xué)好這個本領(lǐng)。其實在大多數(shù)大學(xué)課程中，也沒有太復(fù)雜的推理。如果把公式推理也交給計算機(jī)，對教學(xué)是否有利，是一個有爭議的問題。

(2)作者十幾年來一直努力推動在大學(xué)課程中用計算機(jī)做工具的工作，但強(qiáng)調(diào)的一直是數(shù)值計算。因為數(shù)值計算的道理，從計算器到計算機(jī)，是一種從簡到繁的自然升級，老師和學(xué)生都比較清楚，沒有跳躍的感覺。而公式推理在計算機(jī)中是如何完成的，大部分人都不清楚。所以作者目前的看法是，凡是能夠用手工簡便推理的地方，一般還是少用計算機(jī)來推導(dǎo)。大量的使用符號數(shù)學(xué)工具還是放在科研中比較恰當(dāng)。

(3)符號數(shù)學(xué)工具箱中也包括一些數(shù)值計算和繪圖的功能，其中有些比MATLAB的基本部分優(yōu)越。比如：它的精度就可以超越16位十進(jìn)制，在解方程時允許未知量出現(xiàn)在方程右端，可以根據(jù)輸入函數(shù)的表達(dá)式直接繪圖，等等。這些把數(shù)值計算能力進(jìn)一步提高的內(nèi)容可以在大學(xué)本科中學(xué)習(xí);有些內(nèi)容過去是靠查表的，如積分表、傅里葉變換和拉普拉斯變換表等，現(xiàn)在也不妨用軟件來代替;在帶有研究性的任務(wù)(如數(shù)學(xué)建模等)中也可以較多地使用。即便如此，仍然要以懂得MATLAB的數(shù)值計算作為基礎(chǔ)。從這個思路出發(fā)，我們只對符號數(shù)學(xué)工具箱的部分功能做簡要的介紹。Symbolic工具箱的主要功能有以下8項，括號內(nèi)的英文字符串是Symbolic工具箱中的函數(shù)名。

(1)用符號定義各種數(shù)學(xué)運算和函數(shù)(syms，symop)等;

(2)對這些函數(shù)式進(jìn)行代數(shù)和三角運算，包括因式分解(factor)、展開(expand)、變量置換(subs)、復(fù)合函數(shù)(compose)等;

(3)微分和積分運算(diff，int)等;

(4)函數(shù)的整理和化簡(combine，simplify，simple)等;

(5)可變精度的運算，如可以設(shè)置任意多個有效計算位數(shù)進(jìn)行計算(vpa，digits)等;

(6)解方程，包括單變量的代數(shù)方程、多變量非線性的聯(lián)立代數(shù)方程(solve)、單變量微分方程、多變量聯(lián)立微分方程(dsolve)等;

(7)線性代數(shù)和矩陣運算(determ，linsolve)，包含在matfun函數(shù)庫中的絕大部分矩陣計算和變換的命令都適用于符號變量等;

(8)變換，包括拉普拉斯變換(laplace)、傅里葉變換(fourier)和Z變換(ztrans)等。表4-9列出了本科數(shù)學(xué)可能用到的符號運算工具箱中的主要函數(shù)。4.6.2符號數(shù)學(xué)式的基本表示方法符號數(shù)學(xué)是對字符串進(jìn)行運算的，在命令窗中，如果鍵入

f=3*x^2+5*x+2或

y=sin(x)系統(tǒng)會指出變量x無定義，因為它要求x必須是一個數(shù);MATLAB也可以接受形為f=′3*x^2+5*x+2′或y=′sin(x)′的語句，這時f和y都是一個字符串，但它沒有任何含義，因為MATLAB對字符串中的內(nèi)容不作任何分析。

Symbolic工具箱必須要能分析字符串的含義，為此，首先要對符號變量作出定義，用語句

x=sym(′x′);就定義了x是一個字符(串)變量，此后鍵入的算式

f=3*x^2+5*x+2或

y=sin(x)就具有了符號函數(shù)的意義，連表明字符串的引號都可省略。f和y已被定義(賦值)為符號自變量的函數(shù)，所以也自然成為字符(串)變量?？梢姡恍鑼λ闶降挠疫叿栕兞?即自變量)作定義，系統(tǒng)就會自動把因變量定義為符號變量。如果另外用syms命令定義因變量，那系統(tǒng)將把它看成自變量，而不理會因變量與自變量的函數(shù)關(guān)系了。如果一個數(shù)學(xué)符號表示式中有多個符號(包括自變量和系數(shù))，如

z=a*t^2+b*t+c則可以將簡化的多個符號變量定義語句放在此表示式的前面，如

symsabct在做符號運算時，比如求上述代數(shù)方程的根，就得知道哪個(些)是變量，其余的則是系數(shù)，這時可在語句中指定，比如r=solve('z=0'，t)表示把t作為求解的變量。如果不加說明，則系統(tǒng)將自動把離x最近的那個(些)符號當(dāng)作變量來求解。在本例中，鍵入r=solve(′z=0′)，自然會選到t作為變量。表4-10以表格方式列出了一些算例，從中可對符號數(shù)學(xué)工具箱的功能有一個印象。表中為了節(jié)省篇幅，盡量選了一些簡單的推導(dǎo)式，實際上還可以推導(dǎo)很繁的式子。在公式推導(dǎo)意義下使用MATLAB是很方便的，只是不給自變量賦以數(shù)值，而代之以

syms自變量1自變量2自變量3…以后的編程和普通MATLAB程序完全相同。其執(zhí)行的結(jié)果自然是表達(dá)式而不是數(shù)值解。如果要做進(jìn)一步的工作，例如化簡、代換、代入數(shù)值、解聯(lián)立方程等，那就需要對這個工具箱有較完整的了解。表4-10中列出的是本科數(shù)學(xué)中有可能用到的符號數(shù)學(xué)工具箱中的函數(shù)，供讀者參考，它的一些用法將在第5章中舉例說明。不管在課程中是否使用這個工具箱，用到了多少，應(yīng)該看到的是，計算機(jī)科學(xué)已經(jīng)成功地進(jìn)入了推理領(lǐng)域，而且已經(jīng)可以被普通的科技人員在微機(jī)上實現(xiàn)。我國數(shù)學(xué)家吳文俊教授在計算機(jī)推理領(lǐng)域做出了國際領(lǐng)先的成果，并因此榮獲國家最高科技獎，這也代表了國際科技發(fā)展的方向。沒有任何理由和方法能夠阻止人們?nèi)ソ佑|、學(xué)習(xí)、掌握和應(yīng)用這些推理軟件,教師要考慮到這一點，并探索它對大學(xué)教育改革可能產(chǎn)生的影響。教師首先要懂得，要會用，然后應(yīng)該讓學(xué)生了解這些軟件的功能，便于他們根據(jù)自己的需要來選擇。4.7系統(tǒng)仿真函數(shù)庫框圖是微分方程的等價物。即使是普通的微分方程，也可以用框圖來描述，在工程設(shè)計計算中，人們往往喜歡用框圖來表示復(fù)雜的系統(tǒng)，因為它形象地把各個環(huán)節(jié)區(qū)分開來，并且便于研究各個環(huán)節(jié)的參數(shù)對系統(tǒng)的影響。Simulink工具箱就提供了這樣的手段。以一個最簡單的彈簧-質(zhì)量系統(tǒng)為例，其結(jié)構(gòu)如圖4-14左下角所示。彈簧的左端固定在x1位置，彈簧的右端是一質(zhì)量為M的滑塊，彈簧剛度系數(shù)為K，滑塊M相對于初始平衡位置的位移為x2，它與地面間的阻尼力與滑塊的速度成正比，方向則相反，阻尼系數(shù)為C，于是滑塊受到彈簧彈力和阻尼力兩個力，其運動方程可寫成圖4-14Simulink工具箱所提供的多窗口系統(tǒng)仿真綜合圖形界面這個方程可以用圖4-14中帶兩個反饋的框圖部分來描述。它們分別表示由于速度和位置x2的變化，使作用于滑塊上的力