極值積分微分

13.1繪圖說到繪圖，只要計(jì)算函數(shù)在某一區(qū)間的值，并且畫出結(jié)果向量，這樣就得到了函數(shù)的圖形。在大多數(shù)情況下，這就足夠了。然而，有時一個函數(shù)在某一區(qū)間是平坦的并且無激勵，而在其它區(qū)間卻失控。在這種情況下，運(yùn)用傳統(tǒng)的繪圖方法會導(dǎo)致圖形與函數(shù)真正的特性相去甚遠(yuǎn)。MATLAB提供了一個稱為fplot的巧妙的繪圖函數(shù)。該函數(shù)細(xì)致地計(jì)算要繪圖的函數(shù)，并且確保在輸出的圖形中表示出所有的奇異點(diǎn)。該函數(shù)的輸入需要知道以字符串表示的被畫函數(shù)的名稱以及2元素數(shù)組表示的繪圖區(qū)間。例如：〉〉fplot（‘humps‘,[02]）〉〉title（‘FPLOTOFHUMPS‘）在0和2之間計(jì)算函數(shù)humps，并顯示該函數(shù)的圖形。（見圖13.1）。13.2極小化作圖除了提供視覺信息外，還常常需要確定一個函數(shù)的其它更多的特殊屬性。在許多應(yīng)用中，特別感興趣的是確定函數(shù)的極值，即最大值（峰值）和最小值（谷值）。數(shù)學(xué)上，可通過確定函數(shù)導(dǎo)數(shù)（斜率）為零的點(diǎn)，解析上求出這些極值點(diǎn)。檢驗(yàn)humps的圖形在峰值和谷值點(diǎn)上的斜率就很容易理解這個事實(shí)。顯然，如果定義的函數(shù)簡單，則這種方法常常奏效。然而，即使很多容易求導(dǎo)的函數(shù)，也常常很難找到導(dǎo)數(shù)為零的點(diǎn)。在這種情況下，以及很難或不可能解析上求得導(dǎo)數(shù)的情況下，必須數(shù)值上尋找函數(shù)的極值點(diǎn)。MATLAB提供了兩個完成

此功能的函數(shù)fmin和fmins。這兩個函數(shù)分別尋找一維或n維函數(shù)的最小值。這里僅討論fmin。有關(guān)fmins的詳細(xì)信息，參閱《MATLAB參考指南》。因?yàn)閒(x)的最大值等于-f(x)的最小值所以，上述fmin和fmins可用來求最大值和最小值。如果還不清楚，把上述圖形倒過來看，在這個狀態(tài)下，峰值變成了谷值，而谷值則變成了峰值。為了解釋求解一維函數(shù)的最小值和最大值，再考慮上述例子。從圖13.2可知，在xmax=0.7附近有一個最大值，并且在xmin=4附近有一個最小值。而這些點(diǎn)的解析值為：和。為了方便，用文本編輯器編寫一個腳本M文件，并用fmin尋出數(shù)值上極值點(diǎn)，給出函數(shù)主體如下：%exfmin.mfn二fn二‘2*exp(—x)*sin(x)‘;%definefunctionforminxmin=fmin(fn,2,5)%xmin=fmin(fn,2,5)%searchoverrange2<x<5emin=5*pi/emin=5*pi/4-xmin%finderrorx=xmin;%needxsincefnhasxasitsvariableyminx=xmin;%needxsincefnhasxasitsvariableymin二eval(fn)%evaluateatxminfx二fx二‘-2*exp(-x)*sin(x)‘;%defineformax:noteminussignxmax=fmin(fx,0,3)%searchoverrange0<x<3xmax=fmin(fx,0,3)x=xmax; %needxsincefnhasxasitsvariableymax二eval(fn) %evaluateatxmax下面是M文件的運(yùn)行結(jié)果：>>ex-fminxmin=3.9270emin=1.4523e-006ymin=-0.0279xmax=0.7854emax=-1.3781e-005ymax=這些結(jié)果與上述圖形非常吻合。注意，fmin的工作方式很像fplot。要計(jì)算的函數(shù)可用一個函數(shù)M文件表達(dá)或者只給出一個x為自變量的字符串。上述例子就是使用后一種方法。這個例子也使用了函數(shù)eval，它獲取一個字符串，并解釋它，如同在MATLAB提示符下輸入該字符串。由于要計(jì)算的函數(shù)以x為自變量的字符串形式給出，那么設(shè)置x等于xmin和xmax允許eval計(jì)算該函數(shù)找到y(tǒng)min和ymaxo最后，特別注意，求數(shù)值上的最小值包含一個搜索過程，fmin不斷計(jì)算函數(shù)值，尋求其最小值。如果計(jì)算的函數(shù)需要很大的計(jì)算量，或者該函數(shù)在搜索區(qū)間不止一個最小值則該搜索過程所花的時間比較長。在有些情況下，搜索過程根本找不到結(jié)果。當(dāng)fmin找不到最小值時，它會停止運(yùn)行并提供解釋。與函數(shù)fmin一樣，函數(shù)fmins搜索最小值。不過，fmins搜索向量的標(biāo)量函數(shù)的最小值。即fmins尋找這里x是函數(shù)f(.)的向量參數(shù)，函數(shù)f(.)返回標(biāo)量值。函數(shù)fmins利用單純形法求最小值，它不需要精確的梯度計(jì)算。任何一種優(yōu)化工具箱中具有更多擴(kuò)展的優(yōu)化算法13.3求零點(diǎn)正如人們對尋找函數(shù)的極點(diǎn)感興趣一樣有時尋找函數(shù)過零或等于其它常數(shù)的點(diǎn)也非常重要。一般試圖用解析的方法尋找這類點(diǎn)非常困難，而且很多時候是不可能的。在上述函數(shù)humps的圖中(如圖13.3所示)，該函數(shù)在x=1.2附近過零。圖13.3humps函數(shù)的圖形MATLAB再一次提供了該問題的數(shù)值解法。函數(shù)fzero尋找一維函數(shù)的零點(diǎn)。為了說明該函數(shù)的使用，讓我們再運(yùn)用humps例子。〉〉xzero二fzero(‘humps‘,1.2) %lookforazeronear1.2xzero二1.2995〉〉yzero二humps(xzero,1.2) %evaluateatxzeroyzero二3.5527e-15所以，humps的零點(diǎn)接近于1.3。如前所述，尋找零點(diǎn)的過程可能失敗。如果fzero沒有找到零點(diǎn)，它將停止運(yùn)行并提供解釋。當(dāng)調(diào)用函數(shù)fzero時，必須給出該函數(shù)的名稱。但由于某種原因，它不能接受以x為自變量的字符串來描述的函數(shù)。這樣，即使在fplot和fmin中都具有的這個特性，fzero將不工作。fzero不僅能尋找零點(diǎn)，它還可以尋找函數(shù)等于任何常數(shù)值的點(diǎn)。僅僅要求一個簡單的再定義。例如，為了尋找f(x)=c的點(diǎn)，定義函數(shù)g(x)=f(x)-c，然后，在fzero中使用g(x)，就會找出g(x)為零的x值，它發(fā)生在f(x)=c時。

一個函數(shù)的積分或面積也是它的另一個有用的屬性。MATLAT提供了在有限區(qū)間內(nèi)，數(shù)值計(jì)算某函數(shù)下的面積的三種函數(shù)：trap2,quad和quad8。函數(shù)trapz通過計(jì)算若干梯形面積的和來近似某函數(shù)的積分，這些梯形如圖13.4所示，是通過使用函數(shù)humps的數(shù)據(jù)點(diǎn)形成。圖13.4粗略的梯形逼近曲線下的面積示意圖從圖中可明顯地看出，單個梯形的面積在某一段欠估計(jì)了函數(shù)真正的面積，而在其它段又過估計(jì)了函數(shù)的真正面積。如同線性插值，當(dāng)梯形數(shù)目越多時，函數(shù)的近似面積越準(zhǔn)確。例如，在圖13.4中，如果我們大致增加一倍數(shù)目的梯形，我們得到如下頁(如圖13.5)所示的更好的近似結(jié)果。圖13.5較好的梯形逼近曲線下的面積示意圖對如上所示的兩個曲線，用trapz在區(qū)間T〈x〈2上計(jì)算y二humps(x)下面的面積：%roughapproximation%calltrapzjustlikethe〉〉x=T:%roughapproximation%calltrapzjustlikethe〉〉y二humps(x);〉〉area二trapz(x,y)plotcommandarea=〉〉x=T:〉〉x=T:0.07:2;%betterapproximation〉〉y二humps(x);〉〉area二trapz(x,y)area=26.6243自然地，上述兩個結(jié)果不同?；趯D形的觀察，粗略近似可能低估了實(shí)際面積。除非特別精確，沒有準(zhǔn)則說明哪種近似效果更好。很明顯，如果人們能夠以某種方式改變單個梯形的寬度，以適應(yīng)函數(shù)的特性，即當(dāng)函數(shù)變化快時，使得梯形的寬度變窄，這樣就能夠得到更精確的結(jié)果。MATLAB的函數(shù)quad和quad8是基于數(shù)學(xué)上的正方形概念來計(jì)算函數(shù)的面積，這些積分函數(shù)的操作方式一樣。為獲得更準(zhǔn)確的結(jié)果，兩個函數(shù)在所需的區(qū)間都要計(jì)算被積函數(shù)。此外，與簡單的梯形比較，這兩個函數(shù)進(jìn)行更高階的近似，而且quad8比quad更精確。這兩個函數(shù)的調(diào)用方法與fzero相同，即〉〉area二quad(‘humps‘,T,2)%findareabetweenTand2area=26.3450〉〉area二quad8(‘humps‘,T,2)area=26.3450注意，這兩個函數(shù)返回完全相同的估計(jì)面積，而且這個估計(jì)值在兩個trapz面積的估計(jì)值之間。有關(guān)MATLAB的積分函數(shù)的其它信息，參閱《MATLAB參考指南》或在線幫助。13.5微分與積分相反，數(shù)值微分非常困難。積分描述了一個函數(shù)的整體或宏觀性質(zhì)，而微分則描述一個函數(shù)在一點(diǎn)處的斜率，這是函數(shù)的微觀性質(zhì)。因此積分對函數(shù)的形狀在小范圍內(nèi)的改變不敏感。而微分卻很敏感。一個函數(shù)小的變化，容易產(chǎn)生相鄰點(diǎn)的斜率的大的改變。由于微分這個固有的困難，所以盡可能避免數(shù)值微分，特別是對實(shí)驗(yàn)獲得的數(shù)據(jù)進(jìn)行微分。在這種情況下，最好用最小二乘曲線擬合這種數(shù)據(jù)，然后對所得到的多項(xiàng)式進(jìn)行微分。或用另一種方法，對該數(shù)據(jù)進(jìn)行三次樣條擬合，然后尋找如第11章所討論的樣條微分。例如，再次考慮第11章曲線擬合的例子?！怠祒=[0.1.2.3.4.5.6.7.8.91]〉〉y二[-.4471.9783.286.167.087.347.669.569.489.3011.2];%data〉〉n=2;%orderoffit>>p=polyfit(x,y,n)%findpolynomialcoefficients-9.810820.1293-0.0317〉〉xi=linspace(0,1,100);〉〉z二polyval(p,xi); %evaluatepolynomial〉〉plot(x,y,‘o',x,y,xi,z,':')〉〉xlabel(‘x‘),ylabel(‘y=f(x)‘),title(‘SecondOrderCurveFitting‘)在這種情況下，運(yùn)用多項(xiàng)式微分函數(shù)polyder求得微分?！怠祊d二polyder(p)pd=-19.621720.1293圖13.6二次曲線擬合〈！［endif］〉的微分是dy/dx=-19.6217x+20.1293。由于一個多項(xiàng)式的微分是另一個低一階的多項(xiàng)式所以還可以計(jì)算并畫出該函數(shù)的微分?！怠祕二polyval(pd,xi);%evaluatederivative〉〉plot(xi,z)〉〉xlabel(‘x‘),ylabel(‘dy/dx‘),title(‘DerivativeofacurveFitPolynimial‘)(微分曲線如圖13.7所示)圖13.7曲線擬合多項(xiàng)式微分在這種情況下，擬合的多項(xiàng)式為二階，使其微分為一階多項(xiàng)式。這樣，微分為一條直線，它意味該微分與x成線性變化。給定一些描述某函數(shù)的數(shù)據(jù)，MATLAB提供了一個計(jì)算其非常粗略的微分的函數(shù)。這個函數(shù)命名為diff，它計(jì)算數(shù)組中元素間的差分。因?yàn)槲⒎侄x為：則y=f(x)的微分可近似為：這里h>0它是y的有限差分除以x的有限差分。因?yàn)閐iff計(jì)算數(shù)組元素間的差分，所以在MATLAB中，可近似求得函數(shù)的微分。繼續(xù)前一個例子:>>dy=diff(y)./diff(x);%computedifferencesandusearraydivision>>xd=x(1:length(x)T); %createnewxaxissincedyisshorterthany>>plot(xd,dy);〉〉title(‘ApproximateDerivativeUsingDIFF‘)〉〉ylabel(‘dy/dx‘),xlabel(‘x‘)圖13.8用diff得到的近似微分由于diff計(jì)算數(shù)組元素間的差分，所以，其所得輸出比原數(shù)組少了一個元素。這樣，畫微分曲線時，必須舍棄X數(shù)組中的一個元素。當(dāng)舍棄X的第一個元素時，上述過程給出向后差分近似，而舍棄X的最后一個元素，則給出向前差分近似。比較上述兩條曲線，顯而易見，用有限差分近似微分會導(dǎo)致很差的結(jié)果，特別是被噪聲污染了的數(shù)據(jù)。13.6微分方程一般微分方程式描述系統(tǒng)內(nèi)部變量的變化率如何受系統(tǒng)內(nèi)部變量和外部激勵，如輸入，的影響。當(dāng)常微分方程式能夠解析求解時，可用MATLAB的符號工具箱中的功能找到精確解。在本書的后面將介紹該工具箱的一些特點(diǎn)。在微分方程難以獲得解析解的情況下，可以方便地在數(shù)值上求解。為了說明起見，考慮描述振蕩器的經(jīng)典的范得波(VarderPol)微分方程。與所有的數(shù)值求解微分方程組的方法一樣高階微分方程式必須等價地變換成一階微分方程組。對于上述微分方程，通過重新定義兩個新的變量，來實(shí)現(xiàn)這種變換。令y1=x且y2=dy/dx則dy1/dt二y2根據(jù)這個微分方程組，可用MATLAB的函數(shù)ode23和ode45求出系統(tǒng)隨時間變化的運(yùn)動情況。調(diào)用這些函數(shù)時，需要編寫一個函數(shù)M文件，給定當(dāng)前時間及y1和y2的當(dāng)前值，該函數(shù)返回上述導(dǎo)數(shù)值。MATLAB中，這些導(dǎo)數(shù)由一個列向量給出。在本例中，這個列向量為yprime。同樣，y1和y2合并寫成列向量y。所得函數(shù)M文件是：functionyprime二vdpol(t,y);%VDP0L(t,y)returnsderivativesoftheVanderPolequation:%%x‘‘-mu*(l—x八2)*x‘+x=0 (‘=d/dx, ‘‘=d八2/dx八2)%%lety(l)=xandy(2)=x‘%%theny(1)‘=y(2)% y(2)‘=MU*(l-y(l廠2)*y(2)-y(l)globalMU%choose0〈MU〈10inCommandworkspaceyprime=[y(2)MU*(l-y(l廠2)*y(2)-y(l)];%outputmustbeacolumn給定這個完整地描述微分方程的函數(shù)，計(jì)算結(jié)果如下:〉〉globalMU%defineMUasaglobalvariableintheCommandWorkspace〉〉MU=2;%setglobalparametertodesiredvalue>>[t,y]=ode23(vdpol‘,0,30,[1;0])%to=0,tf=30,yo=[1;0]〉〉yl=y(:,1);%firstcolumnisy(1)versustimepointsint〉〉y2=y(:,2);%secondcolumnisy(2)〉〉plot(t,y1,t,y2,‘一‘)〉〉xlabel(‘Time,Second‘),ylabel(‘Y(l)andY(2)‘)〉〉title(‘VanderPolSolutionformu=2‘)所得的圖見圖13.9。圖13.9當(dāng)mu=2時的范得波方程的運(yùn)動曲線在圖13.9中，y2(虛線)是y1(實(shí)線)的導(dǎo)數(shù)。傳遞給ode23的參數(shù)由ode23(f_name,to,tf,yo,tol)描述。這里f_name是計(jì)算導(dǎo)數(shù)的M文件函數(shù)的字符串名，to是初始時間，tf是終止時間，yo是初始條件向量?？蛇x擇的參數(shù)tol(缺省值to1=1e-3)是所需的相對精度。在上例中，起始時間是第0秒，終止時間是第30秒，初始條件為y二[1;0]。兩個輸出參數(shù)是列向量t和矩陣y，其中向量t包含了估計(jì)響應(yīng)的時間點(diǎn)，而矩陣y的列數(shù)等于微分方程組的個數(shù)(本例為2),且其行數(shù)與t相同。t中的時間點(diǎn)不是等間隔的，因?yàn)闉榱吮３炙璧南鄬?，積分算法改變了步長。函數(shù)ode45的使用與ode23完全一樣。兩個函數(shù)的差別在于必須與所用的內(nèi)部算法相關(guān)。兩個函數(shù)都運(yùn)用了基本的龍格-庫塔(Runge-Kutta)數(shù)值積分法的變形。ode23運(yùn)用一個組合得2/3階龍格-庫塔-芬爾格(Runge-Kutta-Fehlerg)算法，而ode45運(yùn)用組合的4/5階龍格-庫塔-芬爾格算法。一般地，ode45可取較多的時間步。所以，要保持與ode23相同誤差時，在to和tf之間可取較少的時間步。然而，在同一時間，ode23每時間步至少調(diào)用f_name3次，而ode45每時間步至少調(diào)用f_name6次。正如使用高階多項(xiàng)式內(nèi)插常常得不到最好的結(jié)果一樣，ode45也不總是比ode23好。如果ode45產(chǎn)生的結(jié)果，對作圖間隔太大，則必須在更細(xì)的時間區(qū)間，對數(shù)據(jù)進(jìn)行內(nèi)插，比如用函數(shù)interpl。這個附加時間點(diǎn)會使ode23更有效。作為一條普遍規(guī)則，在所計(jì)算的導(dǎo)數(shù)中，如有重復(fù)的不連續(xù)點(diǎn)，為保持精度致使高階算法減少時間步長，這時低階算法更有效。正是由于這個原因，電子電路分析按缺省，就用一階算法編程，并且最多提供二階算法來解決暫態(tài)時間響應(yīng)問題。此外，通過對tol設(shè)置更小的值，要達(dá)到更高的精度，沒有必要使絕對誤差更小。tol設(shè)置每時間步的相對精度不一定引起絕對誤差減少?？傊?，不要盲目使用數(shù)值方法。對于給定的問題，在決定最好的方法之前，要試一試各種可能的方法。有關(guān)微分方程數(shù)值解法的更進(jìn)一步信息，請參考數(shù)值分析方面的書籍。有些參考書還提供了一些關(guān)于算法選擇和如何處理那些時間常數(shù)變化范圍大的病態(tài)方程的非M文件舉例：這里所介紹的《精通MATLAB工具箱》中的M文件可近似求解由采樣值給出的函數(shù)的積分和微分。這里假定這些函數(shù)本身不存在，且獨(dú)立變量也許不是線性間隔。例如，已裝載到MATLAB中要分析的數(shù)據(jù)來源于實(shí)驗(yàn)測試。對于所包含的數(shù)據(jù)缺乏函數(shù)描述，有許多種積分和微分的方法。如前所述，人們可以用最小二乘多項(xiàng)式擬合數(shù)據(jù)，然后在多項(xiàng)式的描述上進(jìn)行操作。另一種方法是尋找數(shù)據(jù)的三次樣條表示，然后運(yùn)用《精通MATLAB工具箱》中的函數(shù)spintgrl和spderiv來分別尋找積分和微分的樣條表示。這里所介紹的方法提供了另一種更簡單的方法。積分用梯形規(guī)則計(jì)算。用加權(quán)中心差分計(jì)算微分。此外，將函數(shù)設(shè)計(jì)成在矩陣形式下工作，矩陣的列代表各與自變量有關(guān)的因變量。正如這章前面所述，MATLAB函數(shù)trapz計(jì)算在某有限區(qū)間的梯形積分。這里我們尋找的積分是自變量為x的函數(shù)。即如果y=f(x)，我們尋找：式中的x1是向量x的第一個元素。用梯形規(guī)則，這個積分近似為：且S(x1)=0這樣，第k個數(shù)據(jù)點(diǎn)的積分是上述梯形面積的累加和。函數(shù)mmintgrl實(shí)現(xiàn)的這個算法如下：functionz=mmintgrl(x,y)%MMINTgrlComputeIntegralusingTrapezoidalRule.%MMINTGRL(X,Y)computestheintegralofthefunctiony二f(x)giventhe%datainXandY.XmustbeavectorbutYmaybeacolumnoriented%datamatrix.ThelengthofXmustequalthelengthofYifYisa%vector,oritmustequalthenumberofrowsinYifYisamatrix.%%Xneednotbeequallyspaced.Thetrapezoidalalgorithmisused.%%Seealsommderiv%Copyrigth(c)1996byPrentice-Hall,Inc.flag=0;%falgisTrueifyisarowx=x(:);nx=length(x);%makexacolumnflag=0;[ry,cy]二size(y);ifry==1&cy==nx,y=y.' ;ry=cy;cy=1;flag=1;endifnx~=ry,error('XandYnottherightsize'),enddx=x(2:nx)-x(1:nxT); %widthofeachtrapezoiddx=dx(:,ones(1,cy)); %duplicateforeachcolumninyyave=(y(2:ry, :)+y(1:ryT, :))/2; %averageofheightsz=[zeros(1,cy);cumsum(dx.*yave)]; %Usecumsumtofindareaifflag,z=z';end %ifywasarow,returnarow在介紹上述函數(shù)的使用之前，考慮微分。在這種情況下，人們感興趣的就是剛給定數(shù)據(jù)點(diǎn)的近似斜率。這里介紹一種下述的中心差分圖13.10加權(quán)中心差分方法從圖13.11可知，在第k個點(diǎn)的近似微分是：式中，并且Mk是連接yk-1到y(tǒng)k的直線的斜率。這樣，第k點(diǎn)的微分是相鄰兩點(diǎn)間斜率的加權(quán)平均，離該點(diǎn)越近的點(diǎn)權(quán)越重。在第一個和最后一個數(shù)據(jù)點(diǎn)上，不能簡單按照上述方法進(jìn)行處理，因?yàn)檫@兩個點(diǎn)都沒有伴隨的直線段。對于這些數(shù)據(jù)點(diǎn)，需要用另外的方法。這里所采取的方法是用二次多項(xiàng)式擬合前3個點(diǎn)(或最后3個點(diǎn))，并且計(jì)算這個多項(xiàng)式第一個(或最后一個)點(diǎn)的微分。函數(shù)mmderiv實(shí)現(xiàn)的這個算法如下：functionz=mmderiv(x,y)%MMDERIVComputeDerivativeUsingWeightedCentralDifferences.%MMDERIV(X,Y)computesthederivativeofthefunctiony二f(x)giventhe%datainXandY.Xmustbeavector,butYmaybeacolumnoriented%datamatrix.ThelengthofXmustequalthelengthofYifYisa%vector,oritmustequalthenumberofrowsinYifYisamatrix.

%Xneednotbeequallyspaced.Weightedcentraldifferenceareused.%Quadraticapproximationisusedattheendpoints.%%Seealsommintgrl%Copyrigth(c)1996byPrentice-Hall,Inc.flag=0; %flagisTrueifyisarowx=x(:);nx=length(x);%makexacolumn[ry,cy]二size(y);ifry==1&cy==nx,y=yry=cy;cy=1;flag=1;endifnx~=ry,error('XandYnottherightsize'),endifnx<3,error('XandYmusthavestleastthreeelements'),enddx=x(2:nx)-x(1dx=x(2:nx)-x(1:nxT);xdx=dx+(dx==0)*eps;dxx=x(3:nx)-x(1:nx-2);dxx=dxx+(dxx==0)*eps;%makeinfiniteslopesfinite%seconddifferenceinx%makeinfiniteslopesfinite

alpha二dx(l:nx-2)./dxxweight%centraldifferencealpha二dx(l:nx-2)./dxxweight%centraldifferencealpha二alpha(:,ones(l,cy));duplicateforeachcolumninydy=y(2:ry, :)-y(1:ry-1,:);%firstdifferenceinydx=dx(:,ones(1,cy));duplicatedxforeachcolumniny%nowapplyweightingtodyz二alpha.*dy(2:ry-1,:)./dx(2:nxT,:)+(1-alpha).*dy(1:ry-2, :)./dx(1:nx-2,:);z1=zeros(1,cy)>=z1;fori=1:cy%fitquadraticatendpointsofeachcolumnpl二polyfit(x(1:3),y(1:3,i),2); %quadraticatfirstpointz1(i)=2*p1(1)*x(1)+p1(2); %evalutepolyderivativepn=polyfit(x(nx-2:nx),y(ry-2:ry,i),2);%quadraticatlastpoint

%evaluatepoly%evaluatepolyderivativeendz=[z1;z;zn];ifflag,z=z';end %ifywasarow最后，給出一個例子：〉〉x=linspace(0,2*pi,30)〉〉y二sin(x); %createdata>>yi=mmintgrl(x,y);%findintegral>>yd=mmderiv(x,y);%findderivative〉〉plot(x,y,x,yi，‘—‘,x,yd,‘：‘)%注意這個積分定性地證明了等式：而微分定性地證明了等式：圖13.11y=sin(x)極其積分、微分曲線13.8小結(jié)表13.1總結(jié)了本章所討論的函數(shù)。,returnarowplotresults表13.1數(shù)值分析函數(shù)fplot(‘fname‘,[lbub]),returnarowplotresultsfmin(‘fname‘,[lbub]) 尋找上下限內(nèi)的標(biāo)量最小值fimis(‘fname‘,xo) 尋找xo附近的向量最小值fzero(‘fname‘,xo) 尋找xo附近的標(biāo)量函數(shù)的零點(diǎn)trapz(x,y)給定數(shù)據(jù)點(diǎn)x和y,計(jì)算y=f(x)下的梯形面積積分。diff(x) 數(shù)組元素間的差分[t,y]=ode23(‘fname‘,to,tf,yo) 用2階/3階龍格-庫塔算法解微分方程組[t,y]=ode45(‘fname‘,to,tf,yo) 用4階/5階龍格-庫塔算法解微分方程組本節(jié)介紹的是關(guān)于函數(shù)的函數(shù)(functionfunctions),這些函數(shù)是用來處理函數(shù)而非數(shù)值的；類別函數(shù)描述繪圖優(yōu)化求解fplot畫出函數(shù)fminbnd由一有范圍限制的變量找出函數(shù)的最小值fminsearch由幾個變量找出函數(shù)的最小值fzero找出函數(shù)的解(零值)數(shù)值積分quad低階數(shù)值估計(jì)積分quad8高階數(shù)值估計(jì)積分dblquad二重積分?jǐn)?shù)值微分見下一章MatLab中的函數(shù)表達(dá)MatLab中用M文件來表示函數(shù)，設(shè)有如下函數(shù)：他別表示為一稱為hump.m的文件中：functiony二humps(x)y=1./((x—0.3).八2+0.01)+1./((x—0.9).八2+0.04)—6;這個函數(shù)文件可用于數(shù)值分析的函數(shù)中.第二種方法就是創(chuàng)造一個行內(nèi)對象(inlineO)，方法如下：f=inline(‘l./((x—0.3).八2+0.01)+l./((x—0.9).八2+0.04)-6');用了上面的方法創(chuàng)造了函數(shù)文件，我們就可以找出函數(shù)在2的值：f(2.0)ans=—4.8552用創(chuàng)造行內(nèi)對象的方法還可以創(chuàng)造多參數(shù)的函數(shù)，如下：f=inline('y*sin(x)+x*cos(y)','x','y')ans=3.1416把函數(shù)畫出來fplot()可畫出在給定范圍內(nèi)的函數(shù)值，如下fplot('humps',[—55])gridon可通過限制y軸來放大圖形fplot('humps',[-55-1025])gridon你也可直接在fplot()中傳遞表達(dá)式，如：fplot('2*sin(x+3)',[—11])更可在一附圖中畫多個函數(shù)，如下fplot('[2*sin(x+3),humps(x)]',[-11])式中，[2*sin(x+3),humps(x)]組成了一個矩陣，每一列都是對應(yīng)于x的函數(shù)函數(shù)的最小值與解找出一變量的函數(shù)的極值x=fminbnd('humps',0.3,1)

0.6370你可通過向fminbnd()函數(shù)傳遞一個函數(shù)optimset()作為參數(shù)來把此過程顯示為列表形式：x=fminbnd('humps',0.3,l,optimset('Display','iter'))Func—countxf(x)Procedure10.56737612.9098initial20.73262413.7746golden30.46524825.1714golden40.64441611.2693parabolic50.641311.2583parabolic60.63761811.2529parabolic70.63698511.2528parabolic80.63701911.2528parabolic90.63705211.2528parabolicx二0.6370多變量函數(shù)極值先創(chuàng)造一個m文件，three_var.m:functionb=three_var(v)TOC\o"1-5"\h\zx= v(1);y= v(2)；z= v(3);b= x.八2+2.5*sin(y)-zJ*x八2*yJ;現(xiàn)在，以x=-0.6,y=-1.2,z=0.135為起始點(diǎn)找出函數(shù)的極值：v=[-0.6-1.20.135];a=fminsearch('three_var',v)a二0.0000-1.57080.1803設(shè)置尋找極值的參數(shù)x=fminbnd(fun,x1,x2,options)或x=fminsearch(fun,x0,options)其中，options是優(yōu)化工具箱中(OptimizationToolbox)中的函數(shù)所用的一個結(jié)構(gòu)，可如下設(shè)置options=optimset('Display','iter');options.DisplAy用來設(shè)置是否顯示中間過程,如為:〃iter〃則顯示，為〃off〃則不顯示，為〃finAl〃則只顯示最后結(jié)果；options.TolX設(shè)置結(jié)果的誤差范圍，默認(rèn)值是：l.e-4.options.MAxFunEvAl設(shè)置函數(shù)運(yùn)行次數(shù)的上限，默認(rèn)fminbnd()是500次,fmin