MATLAB-程式設(shè)計(jì)入門篇-一般數(shù)學(xué)函數(shù)的處理與分析

MATLAB程式設(shè)計(jì)進(jìn)階篇

一般數(shù)學(xué)函數(shù)的處理與分析張智星:///jang臺(tái)大資工系多媒體檢索實(shí)驗(yàn)室函數(shù)的函數(shù)MATLAB可對(duì)數(shù)學(xué)函數(shù)進(jìn)行各種運(yùn)算與分析，例如：作圖求根優(yōu)化：求函數(shù)的極大或極小值數(shù)值積分求解微分方程式如何表示此種被分析的函數(shù)？字串函數(shù)握把(functionhandles)匿名函數(shù)(anonymousfunction)FunctionsofFunctions!一維數(shù)學(xué)函數(shù)的範(fàn)例MATLAB是以M檔案〔副檔名為m〕來(lái)表示一個(gè)函數(shù)例如，內(nèi)建於MATLAB目錄的humps.m可用來(lái)計(jì)算以下函數(shù)：更多資訊：欲顯示此檔案的位置whichhumps欲顯示此檔案的內(nèi)容typehumps欲顯示此檔案位置，可輸入>>whichhumpsC:\ProgramFiles\MATLAB71\toolbox\matlab\demos\humps.m提示MATLAB常被用到的測(cè)試函數(shù)humps：?jiǎn)屋斎牒瘮?shù)peaks：雙輸入函數(shù)「函式」和「函數(shù)」都代表「functions」，兩者常會(huì)混用，假設(shè)要正名，可區(qū)分如下：函數(shù)：通常用來(lái)表示「mathematicfunctions」函式：通常用來(lái)表示「subroutinesorfunctionsinaprogramminglanguage」Weuse函數(shù)torepresentboth.數(shù)學(xué)函數(shù)的作圖

表示函數(shù)的方式函數(shù)握把：使用@humps來(lái)代表humps.m字串：使用'humps'來(lái)代表humps.m用fplot指令進(jìn)行數(shù)學(xué)函數(shù)作圖畫出humps函數(shù)在[0,2]間的曲線範(fàn)例：fplot01.msubplot(2,1,1);fplot('humps',[0,2]); %使用字串指定函式subplot(2,1,2);fplot(@humps,[02]); %使用函式握把來(lái)指定函式Lessflexible!同時(shí)改變x、y的區(qū)間我們可同時(shí)改變x和y的區(qū)間範(fàn)例：fplot02.mx的區(qū)間為[0,1]y的區(qū)間為[5,25]fplot('humps',[0,1,5,25]);gridon %畫出格線匿名函式fplot也接受匿名函式〔當(dāng)場(chǎng)指定的函式〕範(fàn)例：fplot021.msubplot(2,1,1);fplot('sin(2*x)+cos(x)',[-10,10]) %使用字串指定函式subplot(2,1,2);fplot(@(x)sin(2*x)+cos(x),[-10,10]) %使用函式握把來(lái)指定函式對(duì)多個(gè)函數(shù)作圖fplot也可同時(shí)對(duì)多個(gè)函數(shù)作圖，其中每個(gè)函數(shù)須以一個(gè)行向量來(lái)表示範(fàn)例：fplot022.mx是行向量〔MATLAB預(yù)設(shè)值〕[sin(x),exp(-x)]是二個(gè)行向量每個(gè)行向量代表一個(gè)函數(shù)〔即一條曲線〕fplot(@(x)[sin(x),exp(-x)],[0,10])帶有參數(shù)的函數(shù)匿名函式也可以帶有參數(shù)範(fàn)例：fplot023.m此時(shí)“@(x)”不可省略，以便指定自變數(shù)a=1;b=1.1;c=1.2;fplot(@(x)[sin(a*x),sin(b*x),sin(c*x)],[0,10])產(chǎn)生X、Y座標(biāo)點(diǎn)fplot可進(jìn)行描點(diǎn)作圖，類似plot(x,y)，但x座標(biāo)點(diǎn)的密度根據(jù)函數(shù)值的變化決定我們顯示fplot所產(chǎn)生的x座標(biāo)點(diǎn)範(fàn)例：fplot03.m函數(shù)變化平緩處，產(chǎn)生 稀疏的點(diǎn)函數(shù)變化劇烈處，產(chǎn)生 緊密的點(diǎn)[x,y]=fplot(@humps,[-1,2]);plot(x,y,'-o');產(chǎn)生更密的X座標(biāo)點(diǎn)(1)假設(shè)欲產(chǎn)生更密的x座標(biāo)點(diǎn)，可在fplot指令參加另一個(gè)輸入引數(shù)，已指定相對(duì)容忍度(Tolerance)範(fàn)例：fplot04.msubplot(2,1,1);fplot(@(x)sin(1./x),[0.01,0.1]);subplot(2,1,2);fplot(@(x)sin(1./x),[0.01,0.1],0.0001);產(chǎn)生更密的X座標(biāo)點(diǎn)(2)在第一圖中，fplot指令使用預(yù)設(shè)相對(duì)容忍度，其值為0.002。在第二圖中，相對(duì)容忍度被設(shè)為0.0001，可得到更準(zhǔn)確的圖形，但也要花更多計(jì)算及作圖時(shí)間。ezplot指令ezplot指令和fplot指令類似，可進(jìn)行描點(diǎn)作圖，但使用更為簡(jiǎn)便，預(yù)設(shè)的作圖範(fàn)圍為範(fàn)例8-7：ezplot01.mezplot(@(x)x^3-x^2+x);平面中的參數(shù)式曲線ezplot也可畫出平面中的參數(shù)式曲線範(fàn)例8-8：ezplot02.m參數(shù)式函數(shù)的參數(shù)預(yù)設(shè)範(fàn)圍仍是ezplot(@(t)sin(3*t),@(t)cos(5*t));

利薩如圖形(LissajousFigures)空間中的參數(shù)式曲線ezplot3可畫出空間中的參數(shù)式曲線範(fàn)例8-8：ezplot021.m參數(shù)式函數(shù)的參數(shù)預(yù)設(shè)範(fàn)圍仍是ezplot3(@(t)sin(3*t),@(t)cos(5*t),@(t)t)3D利薩如圖形隱函數(shù)作圖ezplot指令可用於隱函數(shù)作圖以下範(fàn)例可以畫出範(fàn)例8-9：ezplot03.mezplot(@(x,y)x^3+2*x^2-3*x-y^2+15);函數(shù)的求根fzero指令用於單變數(shù)函數(shù)的求根語(yǔ)法x=fzero(fun,x0)fun是欲求根的函數(shù)〔以字串或函數(shù)握把來(lái)表示〕x0是一個(gè)起始點(diǎn)或起始區(qū)間X0對(duì)fzero的影響fzero指令根據(jù)x0不同而執(zhí)行以下動(dòng)作假設(shè)x0為一個(gè)起始點(diǎn)fzero會(huì)自動(dòng)找出附近包含零點(diǎn)〔即根，或函數(shù)變號(hào)點(diǎn)〕的區(qū)間逐步縮小此區(qū)間以找出零點(diǎn)假設(shè)fzero無(wú)法找出此區(qū)間，傳回NaN假設(shè)使函數(shù)值不同號(hào)的兩點(diǎn)由x0直接指定尋根的區(qū)間fzero更快速找到位於此區(qū)間內(nèi)的根求根範(fàn)例(1)找出humps在x=1.5附近的根，並驗(yàn)算範(fàn)例8-10：fzero01.mfzero先找到在1.5附近變號(hào)的兩點(diǎn)〔即1.26及1.6697〕，然後再找出humps的零點(diǎn)x=fzero(@humps,1.5); %求靠近1.5附近的根y=humps(x);

%帶入求值fprintf('humps(%f)=%f\n',x,y);humps(1.299550)=0.000000求根範(fàn)例(2)假設(shè)humps在x=-1及1間為異號(hào)令x0=[-1,1]為起始區(qū)間來(lái)找出humps的零點(diǎn)範(fàn)例8-11：fzero02.m此時(shí)fzero找到的是另一個(gè)零點(diǎn)x=fzero(@humps,[-1,1]); %求落於區(qū)間[-1,1]的根y=humps(x); %帶入求值fprintf('humps(%f)=%f\n',x,y);humps(-0.131618)=0.000000求根範(fàn)例(3)假設(shè)要畫出以上兩個(gè)零點(diǎn)，請(qǐng)見以下範(fàn)例範(fàn)例8-12：fzero03.mfplot(@humps,[-1,2]);gridonz1=fzero(@humps,1.5);z2=fzero(@humps,[-1,1]);line(z1,humps(z1),'marker','o','color','r');line(z2,humps(z2),'marker','o','color','r');顯示求解過程的中間結(jié)果(1)MATLAB可以顯示求解過程的中間結(jié)果使用optimset指令來(lái)設(shè)定顯示選項(xiàng)再將optimset傳回結(jié)構(gòu)變數(shù)送入fzero範(fàn)例8-13：fzero04.moptimset常用於設(shè)定最正確化的選項(xiàng)，下一節(jié)會(huì)有比較完整的介紹opt=optimset('disp','iter'); %顯示每個(gè)iteration的結(jié)果

a=fzero(@humps,[-1,1],opt)顯示求解過程的中間結(jié)果(2)

求零點(diǎn)過程中，找下一點(diǎn)的兩個(gè)方法顯示在第四個(gè)欄位〔Procedure欄位〕二分法(Bisection)內(nèi)插法(Interpolation)可由docfzero找到所使用的演算法Func-countxf(x)Procedure1-1-5.13779initial2116Iinitial3-0.513876-4.02235interpolation40.24306271.6382bisection5-0.473635-3.83767interpolation6-0.1152870.414441bisection7-0.150214-0.423446interpolation8-0.132562-0.0226907interpolation9-0.131666-0.0011492interpolation10-0.1316181.88371e-007interpolation11-0.131618-2.7935e-011interpolation12-0.1316188.88178e-016interpolation13-0.131618-9.76996e-015interpolationZerofoundintheinterval:[-1,1].a=-0.1316數(shù)值積分MATLAB可用於計(jì)算單變函數(shù)定積分quad：適應(yīng)式Simpson積分法〔AdaptiveSimpsonQuadrature〕quadl：適應(yīng)式Lobatto積分法〔AdaptiveLobattoQuadrature〕定積分計(jì)算humps在[0,1]的定積分>>q=quad(@humps,0,1)q=29.8583quad及quad8都應(yīng)用遞迴程序假設(shè)遞迴次數(shù)達(dá)10次，兩種方法均會(huì)傳回Inf表示所計(jì)算之定積分可能不存在quad及quad8第四個(gè)引數(shù)用來(lái)指定積分的相對(duì)誤差容忍度曲線的長(zhǎng)度(1)quad及quadl計(jì)算曲線的長(zhǎng)度一曲線是由以下參數(shù)化的方程式來(lái)表示t的範(fàn)圍為[0,3*pi]範(fàn)例：plotCurve.mt=0:0.1:3*pi;plot3(sin(2*t),cos(t),t);曲線的長(zhǎng)度(2)此曲線的長(zhǎng)度等於曲線的長(zhǎng)度(3)先定義函數(shù)

curveLength.m>>typecurveLength.mfunctionout=curveLength(t)out=sqrt(4*cos(2*t).^2+sin(t).^2+1);曲線長(zhǎng)度可計(jì)算如下>>len=quad(@curveLength,0,3*pi)len=17.2220雙重積分(1)dblquad指令用來(lái)計(jì)算雙重積分欲計(jì)算其中先建立被積分的函數(shù)integrand.m>>typeintegarnd.mfunctionout=integrand(x,y)out=y*sin(x)+x*cos(y);雙重積分(2)計(jì)算雙重積分result=dblquad('integrand',xMin,xMax,yMin,yMax)其中xMin：內(nèi)迴圈積分的下界值xMax：內(nèi)迴圈積分的上界值yMin：外迴圈積分的下界值yMax：外迴圈積分的上界值雙重積分(3)範(fàn)例：dblquad01.m一般的情況下dblquad會(huì)呼叫quad計(jì)算定積分。假設(shè)須呼叫更為精確的quadl，可執(zhí)行以下指令>>result=dblquad(@integrand,xMin,xMax,yMin,yMax,'quadl)result=-9.8696xMin=pi;xMax=2*pi;yMin=0;yMax=pi;result=dblquad(@integrand,xMin,xMax,yMin,yMax)result=

-9.8698函數(shù)的優(yōu)化MATLAB提供了數(shù)個(gè)根本指令來(lái)進(jìn)行數(shù)學(xué)函數(shù)的優(yōu)化，本節(jié)將介紹：?jiǎn)巫償?shù)函數(shù)的最小化多變數(shù)函數(shù)的最小化設(shè)定最正確化的選項(xiàng)假設(shè)讀者有興趣使用較複雜的方法，可以使用「最正確化工具箱」〔OptimizationToolbox〕單變函數(shù)的最小化

fminbnd指令尋求humps在[0.3,1]中的最小值範(fàn)例：fminbnd01.m最小值發(fā)生在x=0.637，且最小值為11.2528[x,minValue]=fminbnd(@humps,0.3,1)x=0.6370minValue=11.2528尋求最小值的中間過程(1)尋求最小值的中間過程使用optimset指令來(lái)設(shè)定顯示選項(xiàng)再將optimset傳回結(jié)構(gòu)變數(shù)送入fminbnd範(fàn)例8-15：fminbnd02.mopt=optimset('disp','iter'); %顯示每個(gè)步驟的結(jié)果[x,minValue]=fminbnd(@humps,0.3,1,opt)尋求最小值的中間過程(2)

左表列出x值的變化及相對(duì)的函數(shù)值f(x)最後一欄位列出求極小值的方法，包含黃金分割搜尋

(GoldenSectionSearch)拋物線內(nèi)插法

(ParabolicInterpolation)x值誤差小於10^-4

Func-countxf(x)Procedure10.56737612.9098initial20.73262413.7746golden30.46524825.1714golden40.64441611.2693parabolic50.641311.2583parabolic60.63761811.2529parabolic70.63698511.2528parabolic80.63701911.2528parabolic90.63705211.2528parabolicOptimizationterminatedsuccessfully:thecurrentxsatisfiestheterminationcriteriausingOPTIONS.TolXof1.000000e-004x=0.6370minValue=11.2528放鬆誤差管制(1)放鬆誤差管制使fminbnd提早傳回計(jì)算結(jié)果由optimset達(dá)成下例將x的誤差範(fàn)圍提高為0.1範(fàn)例8-16：fminbnd03.mopt=optimset('disp','iter','TolX',0.1);[x,minValue]=fminbnd(@humps,0.3,1,opt)放鬆誤差管制(2)

Func-countxf(x)Procedure10.56737612.9098initial20.73262413.7746golden30.46524825.1714golden40.64441611.2693parabolic50.61108311.4646parabolic60.67774911.7353parabolic

Optimizationterminatedsuccessfully:thecurrentxsatisfiestheterminationcriteriausingOPTIONS.TolXof1.000000e-001

x=0.6444minValue=11.2693多變數(shù)函數(shù)的極小值fminsearch指令求多變數(shù)函數(shù)的極小值必須指定一個(gè)起始點(diǎn)，讓fminsearch求出在起始點(diǎn)附近的局部最小值〔LocalMinima〕假設(shè)目標(biāo)函數(shù)為 我們必須先產(chǎn)生一個(gè)MATLAB的函數(shù)objective.m>>typeobjective.mfunctiony=objective(x)y=(x(1)-1)^2+(x(2)-2)^2+(x(3)-3)^2;多變數(shù)函數(shù)的極小值範(fàn)例假設(shè)起始點(diǎn)為範(fàn)例：fminsearch01.mfminsearch使用的方法是下坡式Simplex搜尋〔DownhillSimplexSearch,簡(jiǎn)稱DSS〕x0=[0,0,0];x=fminsearch(@objective,x0)x=1.00002.00003.0000最正確化選項(xiàng)MATLAB最正確化的選項(xiàng)經(jīng)由另一個(gè)輸入引數(shù)〔InputArgument〕來(lái)進(jìn)入fminbnd或fminsearch使用語(yǔ)法x=fminbnd(@objFun,x1,x2,options)x=fminbnd(@(x)objFun(x,a),x1,x2,options) 或x=fminsearch(@objFun,x0,options)x=fminsearch(@(x)objFun(x,a),x0,options)options此結(jié)構(gòu)變數(shù)可代表各種最正確化的選項(xiàng)〔或參數(shù)〕設(shè)定最正確化選項(xiàng)(1)如何設(shè)定最正確化選項(xiàng)用optimset指令options=optimset(prop1,value1,prop2,value2,…)prop1、prop2：欄位名稱value1、value2：對(duì)應(yīng)的欄位值設(shè)定最正確化選項(xiàng)(2)印出最正確化步驟的中間結(jié)果，並放鬆誤差範(fàn)圍>>options=optimset('Disp','iter','TolX',0.1)Display:'iter'MaxFunEvals:[]MaxIter:[]TolFun:[]TolX:0.1000FunValCheck:[]OutputFcn:[]PlotFcns:[]ActiveConstrTol:[]…設(shè)定最正確化選項(xiàng)(3)options共有五十多個(gè)欄位如果欄位值顯示是空矩陣，使用此欄位的預(yù)設(shè)值來(lái)進(jìn)行運(yùn)算>>options=optimset('fminbnd')顯示非空矩陣的最正確化選項(xiàng)：Display:'notifyMaxFunEvals:500MaxIter:500TolX:1.0000e-004FunValCheck:'off'設(shè)定最正確化選項(xiàng)(4)假設(shè)輸入>>options=optimset('fminsearch')顯示非空矩陣的最正確化選項(xiàng)：Display:'notify‘MaxFunEvals:'200*numberofvariables'MaxIter:'200*numberofvariables‘TolFun:1.0000e-004TolX:1.0000e-004FunValCheck:'off'最正確化選項(xiàng)說(shuō)明(1)Display假設(shè)為0〔預(yù)設(shè)值〕，不顯示中間運(yùn)算結(jié)果假設(shè)不為0，則顯示運(yùn)算過程的中間結(jié)果MaxFunEvals函數(shù)求值運(yùn)算〔FunctionEvaluation〕的最高次數(shù)對(duì)fminbnd的預(yù)設(shè)值是500對(duì)fminsearch的預(yù)設(shè)值是200乘上x0的長(zhǎng)度MaxIter最大疊代次數(shù)對(duì)fminbnd的預(yù)設(shè)值是500對(duì)fminsearch的預(yù)設(shè)值是200乘上x0的長(zhǎng)度最正確化選項(xiàng)說(shuō)明(2)TolFun

決定終止搜尋的函數(shù)值容忍度預(yù)設(shè)為10^-4此選項(xiàng)只被fminsearch用到，fminbnd並不使用TolX終止搜尋的自變數(shù)值容忍度，預(yù)設(shè)為10-4

提示最正確化並非一蹴可及，通常一再重覆，最後才能收斂到最正確點(diǎn)最正確化的結(jié)果和起始點(diǎn)的選定有很大的關(guān)聯(lián)，一個(gè)良好的起始點(diǎn)加快最正確化收斂的速度提高找到全域最正確值(GlobalOptimum)的機(jī)會(huì)StepsforDown-hillSimplexsearch(DSS)DefineanobjectivefunctionmyFun(x)Setinitialguess“x0”Startthesearchviafminsearchx=fminsearch(@myFun,x0);VariantsObjectivefunctionmyFun(x,prm)Startthesearchx=fminsearch(@myFun,x0,[],prm);x=fminsearch(@(x)myFun(x,prm),x0);Optimizationvia

DownhillSimplexSearch

DSS:Example1(1/3)Properties!DefinitionofTheFermatpoint!ProblemdefinitionFindapointPinaplanesuchthatPA+PB+PCisminimized.SolutionAnalyticsolutionPistheFermatpointsatisfyingAFB=BFC=CFA=120oNumericalsolutionAllkindsofoptimizationmethodsDSS:Example1(2/3)

Objectivefunctiondist2ABC.mMainprogram:goMinDist2ABC.mfunctionsumDistance=dist2ABC(x)%dist2ABC:ThedistancesumtopointsA,B,CA=[0,0];B=[3,0];C=[0,4];sumDistance=norm(x-A)+norm(x-B)+norm(x-C);p0=[5,5]; %Initialguessp=fminsearch(@dist2ABC,p0);…BadideatohardwirethedatapointsDSS:Example1(3/3)

DSS:Example2(1/3)

ProblemdefinitionFindapointPinaplanesuchthatthetotaldistancetoasetofpoint