MATLAB程式设计进阶篇一般数学函数的处理与分析

1、MATLAB 程式設計進階篇一般數學函數的處理與分析張智星/jang台大資工系 多媒體檢索實驗室MATLAB 程式設計進階篇：一般數學函數的處理與分析函數的函數nMATLAB 可對數學函數進行各種運算與分析，例如：n作圖n求根n優化：求函數的極大或極小值n數值積分n求解微分方程式n如何表示此種被分析的函數？n字串n函數握把 (function handles)n匿名函數 (anonymous function)Functions ofFunctions!MATLAB 程式設計進階篇：一般數學函數的處理與分析一維數學函數的範例nMATL

2、AB 是以 M 檔案（副檔名為 m）來表示一個函數n例如，內建於MATLAB目錄的 humps.m 可用來計算下列函數：n更多資訊：n欲顯示此檔案的位置 which humpsn欲顯示此檔案的內容 type humps604. 0)9 . 0(101. 0)3 . 0(1)(22xxxfMATLAB 程式設計進階篇：一般數學函數的處理與分析提示nMATLAB 常被用到的測試函數nhumps：單輸入函數npeaks：雙輸入函數n函式和函數都代表functions，兩者常會混用，若要正名，可區分如下：n函數：通常用來表示mathematic functionsn函式：通常用來表示subroutin

3、es or functions in a programming languageWe use 函數 to represent both.MATLAB 程式設計進階篇：一般數學函數的處理與分析數學函數的作圖 n表示函數的方式n函數握把：使用 humps 來代表 humps.mn字串：使用 humps 來代表 humps.mn用 fplot 指令進行數學函數作圖n畫出 humps 函數在 0,2 間的曲線n範例：fplot01.msubplot(2,1,1);fplot(humps, 0,2);% 使用字串指定函式subplot(2,1,2);fplot(humps, 0 2);% 使用函式握把

4、來指定函式00.81.82-5005010000.81.82-50050100Less flexible!MATLAB 程式設計進階篇：一般數學函數的處理與分析同時改變 x、y 的區間n我們可同時改變 x 和 y 的區間n範例：fplot02.mnx 的區間為0, 1ny 的區間為5, 25 fplot(humps, 0, 1, 5, 25);grid on% 畫出格線01510152025MATLAB 程式設計進階篇：一般數學函數的處理與分析匿名函式nfplot

5、也接受匿名函式（當場指定的函式）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)% 使用函式握把來指定函式-10-8-6-4-20246810-2-1012-10-8-6-4-20246810-2-1012MATLAB 程式設計進階篇：一般數學函數的處理與分析對多個函數作圖nfplot 也可同時對多個函數作圖，其中每個函數須以一個行向量來表示n範例：fplot022.mnx 是行向量（MATLAB 預設值）n

6、sin(x), exp(-x) 是二個行向量n每個行向量代表一個函數（即一條曲線）fplot(x)sin(x), exp(-x), 0, 10)012345678910-1-0.8-0.6-0.4-0.60.81MATLAB 程式設計進階篇：一般數學函數的處理與分析帶有參數的函數n匿名函式也可以帶有參數n範例：fplot023.mn此時 “(x)” 不可省略，以便指定自變數a=1; b=1.1; c=1.2;fplot(x)sin(a*x), sin(b*x), sin(c*x), 0, 10)012345678910-1-0.8-0.6-0.4-0.6

7、0.81Function handle is more flexible!MATLAB 程式設計進階篇：一般數學函數的處理與分析產生 X、Y 座標點nfplot 可進行描點作圖，類似 plot(x, y)，但x 座標點的密度根據函數值的變化決定n我們顯示 fplot 所產生的 x 座標點n範例：fplot03.mn函數變化平緩處，產生稀疏的點n函數變化劇烈處，產生緊密的點 x, y = fplot(humps, -1,2);plot(x, y, -o);-1-0.500.511.52-20020406080100MATLAB 程式設計進階篇：一般數學函數的處理與分析產生更密的 X 座標點 (1

8、)n若欲產生更密的 x 座標點，可在 fplot 指令加入另一個輸入引數，已指定相對容忍度(Tolerance)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);MATLAB 程式設計進階篇：一般數學函數的處理與分析產生更密的 X 座標點 (2)n在第一圖中，fplot 指令使用預設相對容忍度，其值為 0.002。n在第二圖中，相對容忍度被設為 0.0001，可得到更準確的圖形，但也要花更多計算及作圖時間。0.010.020

9、.030.040.050.060.070.080.090.1-1-0.500.510.010.020.030.040.050.060.070.080.090.1-1-0.500.51MATLAB 程式設計進階篇：一般數學函數的處理與分析ezplot 指令nezplot指令和fplot指令類似，可進行描點作圖，但使用更為簡便，預設的作圖範圍為n範例8-7：ezplot01.mezplot(x)x3-x2+x);-6-4-20246-300-250-200-150-100-50050100150200 xx3-x2+x2,2MATLAB 程式設計進階篇：一般數學函數的處理與分析平面中的參數式曲線n

10、ezplot 也可畫出平面中的參數式曲線 n範例8-8：ezplot02.mn參數式函數的參數預設範圍仍是 ezplot(t)sin(3*t), (t)cos(5*t);-1-0.500.51-0.8-0.6-0.4-0.60.81xyx = sin(3 t), y = cos(5 t)2,2)5cos()3sin(tytx利薩如圖形 (Lissajous Figures)MATLAB 程式設計進階篇：一般數學函數的處理與分析空間中的參數式曲線nezplot3 可畫出空間中的參數式曲線 n範例8-8：ezplot021.mn參數式函數的參數預設範圍仍是 ezplot3(t)

11、sin(3*t), (t)cos(5*t), (t)t)2,2tztytx)5cos()3sin(3D利薩如圖形MATLAB 程式設計進階篇：一般數學函數的處理與分析隱函數作圖nezplot 指令可用於隱函數作圖n下列範例可以畫出n範例8-9：ezplot03.mezplot(x,y)x3+2*x2-3*x-y2+15);-6-4-20246-6-4-20246xyx3+2 x2-3 x+5-y2+10 = 001532223yxxxMATLAB 程式設計進階篇：一般數學函數的處理與分析函數的求根nfzero 指令n用於單變數函數的求根n語法x = fzero(fun, x0)nfun 是欲求

12、根的函數（以字串或函數握把來表示）nx0 是一個起始點或起始區間MATLAB 程式設計進階篇：一般數學函數的處理與分析X0 對 fzero 的影響nfzero 指令根據 x0 不同而執行下列動作n若 x0 為一個起始點n fzero 會自動找出附近包含零點（即根，或函數變號點）的區間n逐步縮小此區間以找出零點n若 fzero 無法找出此區間，傳回 NaNn若已知使函數值不同號的兩點n由 x0 直接指定尋根的區間n fzero 更快速找到位於此區間內的根 MATLAB 程式設計進階篇：一般數學函數的處理與分析求根範例 (1)n找出humps在 x = 1.5 附近的根，並驗算n範例8-10：fz

13、ero01.mnfzero 先找到在 1.5 附近變號的兩點（即 1.26 及 1.6697），然後再找出 humps 的零點 x = fzero(humps, 1.5);% 求靠近 1.5 附近的根y = humps(x); % 帶入求值fprintf(humps(%f) = %fn, x , y);humps(1.299550) = 0.000000MATLAB 程式設計進階篇：一般數學函數的處理與分析求根範例 (2)n若已知 humps 在 x = -1 及 1 間為異號n令 x0 = -1, 1 為起始區間來找出 humps 的零點n範例8-11：fzero02.mn此時 fzero

14、找到的是另一個零點x = fzero(humps, -1, 1);% 求落於區間 -1, 1 的根y = humps(x);% 帶入求值fprintf(humps(%f) = %fn, x , y);humps(-0.131618) = 0.000000MATLAB 程式設計進階篇：一般數學函數的處理與分析求根範例 (3)n若要畫出以上兩個零點，請見下列範例n範例8-12：fzero03.mfplot(humps, -1, 2); grid onz1 = fzero(humps, 1.5);z2 = fzero(humps, -1, 1);line(z1, humps(z1), marker,

15、 o, color, r);line(z2, humps(z2), marker, o, color, r);-1-0.500.511.52-20020406080100MATLAB 程式設計進階篇：一般數學函數的處理與分析顯示求解過程的中間結果 (1)nMATLAB 可以顯示求解過程的中間結果n使用 optimset 指令來設定顯示選項n再將 optimset 傳回結構變數送入 fzeron範例8-13：fzero04.mnoptimset 常用於設定最佳化的選項，下一節會有比較完整的介紹opt = optimset(disp, iter);% 顯示每個 iteration 的結果 a =

16、fzero(humps, -1, 1, opt)MATLAB 程式設計進階篇：一般數學函數的處理與分析顯示求解過程的中間結果 (2) n求零點過程中，找下一點的兩個方法顯示在第四個欄位（Procedure 欄位）n二分法(Bisection)n內插法(Interpolation)n可由doc fzero找到所使用的演算法Func-count x f(x) Procedure 1 -1 -5.13779 initial 2 1 16 I initial 3 -0.513876 -4.02235 interpolation 4 0.243062 71.6382 bisection 5 -0.473

17、635 -3.83767 interpolation 6 -0.115287 0.414441 bisection 7 -0.150214 -0.423446 interpolation 8 -0.132562 -0.0226907 interpolation 9 -0.131666 -0.0011492 interpolation 10 -0.131618 1.88371e-007 interpolation 11 -0.131618 -2.7935e-011 interpolation 12 -0.131618 8.88178e-016 interpolation 13 -0.131618

18、 -9.76996e-015 interpolationZero found in the interval: -1, 1.a = -0.1316MATLAB 程式設計進階篇：一般數學函數的處理與分析數值積分nMATLAB 可用於計算單變函數定積分nquad：適應式 Simpson 積分法（Adaptive Simpson Quadrature）nquadl：適應式 Lobatto 積分法（Adaptive Lobatto Quadrature）MATLAB 程式設計進階篇：一般數學函數的處理與分析定積分n計算 humps 在 0, 1 的定積分 q = quad(humps, 0, 1)q

19、= 29.8583nquad 及 quad8 都應用遞迴程序n若遞迴次數達 10 次，兩種方法均會傳回 Infn表示所計算之定積分可能不存在nquad 及 quad8第四個引數用來指定積分的相對誤差容忍度MATLAB 程式設計進階篇：一般數學函數的處理與分析曲線的長度 (1)nquad 及 quadl 計算曲線的長度n一曲線是由下列參數化的方程式來表示n t 的範圍為 0, 3*pin範例：plotCurve.mttzttyttx)()cos()()2sin()(t = 0:0.1:3*pi;plot3(sin(2*t), cos(t), t);MATLAB 程式設計進階篇：一般數學函數的處理

20、與分析曲線的長度 (2)n此曲線的長度等於dtttdtdtdzdtdydtdx3022302221)(sin)2(cos4-1-0.500.51-1-0.500.510246810MATLAB 程式設計進階篇：一般數學函數的處理與分析曲線的長度 (3)n先定義函數 curveLength.m type curveLength.mfunction out = curveLength(t)out = sqrt(4*cos(2*t).2+sin(t).2+1);n曲線長度可計算如下 len = quad(curveLength, 0, 3*pi)len = 17.2220MATLAB 程式設計進階篇

21、：一般數學函數的處理與分析雙重積分 (1)ndblquad 指令n用來計算雙重積分n欲計算n其中n先建立被積分的函數 integrand.m type integarnd.mfunction out = integrand(x, y)out = y*sin(x) + x*cos(y); maxminmaxmin),(yyxxdxdyyxf)sin()sin(),(yxxyyxfMATLAB 程式設計進階篇：一般數學函數的處理與分析雙重積分 (2)n計算雙重積分result = dblquad( integrand, xMin, xMax, yMin, yMax)n其中nxMin：內迴圈積分的下

22、界值nxMax：內迴圈積分的上界值nyMin：外迴圈積分的下界值nyMax：外迴圈積分的上界值MATLAB 程式設計進階篇：一般數學函數的處理與分析雙重積分 (3)n範例：dblquad01.mn一般的情況下dblquad 會呼叫 quad 計算定積分。若須呼叫更為精確的 quadl，可執行下列指令result = dblquad(integrand, xMin, xMax, yMin, yMax, quadl)result = -9.8696xMin = pi;xMax = 2*pi;yMin = 0;yMax = pi;result = dblquad(integrand, xMin, x

23、Max, yMin, yMax)result = -9.8698MATLAB 程式設計進階篇：一般數學函數的處理與分析函數的優化nMATLAB 提供了數個基本指令來進行數學函數的優化，本節將介紹：n單變數函數的最小化: fminbndn多變數函數的最小化: fminsearchn設定最佳化的選項n若讀者有興趣使用較複雜的方法，可以使用最佳化工具箱（Optimization Toolbox）MATLAB 程式設計進階篇：一般數學函數的處理與分析單變函數的最小化n fminbnd 指令n尋求 humps 在 0.3, 1 中的最小值n範例：fminbnd01.mn最小值發生在 x = 0.637，

24、且最小值為 11.2528x, minValue = fminbnd(humps, 0.3, 1)x = 0.6370minValue = 11.2520.91102030405060708090100MATLAB 程式設計進階篇：一般數學函數的處理與分析尋求最小值的中間過程 (1)n尋求最小值的中間過程n使用 optimset 指令來設定顯示選項n再將 optimset 傳回結構變數送入 fminbndn範例8-15：fminbnd02.mopt = optimset(disp, iter);% 顯示每個步驟的結果x, minValue = fminbn

25、d(humps, 0.3, 1, opt)MATLAB 程式設計進階篇：一般數學函數的處理與分析尋求最小值的中間過程 (2) n左表列出x值的變化及相對的函數值 f(x) n最後一欄位列出求極小值的方法，包含n黃金分割搜尋(Golden Section Search)n拋物線內插法(Parabolic Interpolation)nx 值誤差小於 10-4 Func-count x f(x) Procedure 1 0.567376 12.9098 initial 2 0.732624 13.7746 golden 3 0.465248 25.1714 golden 4 0.644416 11

26、.2693 parabolic 5 0.6413 11.2583 parabolic 6 0.637618 11.2529 parabolic 7 0.636985 11.2528 parabolic 8 0.637019 11.2528 parabolic 9 0.637052 11.2528 parabolicOptimization terminated successfully: the current x satisfies the termination criteria using OPTIONS.TolX of 1.000000e-004 x = 0.6370minValue

27、= 11.2528MATLAB 程式設計進階篇：一般數學函數的處理與分析放鬆誤差管制 (1)n放鬆誤差管制n使 fminbnd 提早傳回計算結果n由 optimset 達成n下例將 x 的誤差範圍提高為 0.1n範例8-16：fminbnd03.mopt = optimset( disp, iter, TolX, 0.1); x, minValue = fminbnd(humps, 0.3, 1, opt)MATLAB 程式設計進階篇：一般數學函數的處理與分析放鬆誤差管制 (2) Func-count x f(x) Procedure 1 0.567376 12.9098 initial 2

28、0.732624 13.7746 golden 3 0.465248 25.1714 golden 4 0.644416 11.2693 parabolic 5 0.611083 11.4646 parabolic 6 0.677749 11.7353 parabolic Optimization terminated successfully: the current x satisfies the termination criteria using OPTIONS.TolX of 1.000000e-001 x = 0.6444minValue = 11.2693MATLAB 程式設計進

29、階篇：一般數學函數的處理與分析多變數函數的極小值:fminsearchnfminsearch 指令n求多變數函數的極小值n必須指定一個起始點，讓 fminsearch 求出在起始點附近的局部最小值（Local Minima）nDerivative free Less efficientnMethod: Downhill Simplex Search (DSS), akanNelder-Mead methodnAmoeba methodMATLAB 程式設計進階篇：一般數學函數的處理與分析Downhill Simplex SearchnDSS in WikipedianMany variatio

30、nsnBasic StepsnUse a simplex to explore the objective function, with the operations:nReflectionnExpansionnContractionnShrinkMATLAB 程式設計進階篇：一般數學函數的處理與分析Downhill Simplex SearchnAbout DSSnStrengthsnStraightforward conceptnEasy programmingnNo gradient or derivative needednWeaknessnSlownOnly good for con

31、tinuous objective functionnCould be trapped in local minimaQuiz!MATLAB 程式設計進階篇：一般數學函數的處理與分析多變數函數的極小值範例n若目標函數為我們必須先產生一個 MATLAB 的函數 objective.m n若起始點為n範例：fminsearch01.m 0 , 0 , 0,321xxxx0 = 0, 0, 0;x = fminsearch(objective, x0)332221321)3()2() 1(),(xxxxxxffunction y = objective(x)y = (x(1)-1)2 +(x(2)-

32、2)2 + (x(3)-3)2;x = 1.0000 2.0000 3.0000MATLAB 程式設計進階篇：一般數學函數的處理與分析最佳化選項nMATLAB 最佳化的選項n經由另一個輸入引數（Input Argument）來進入 fminbnd 或 fminsearchn使用語法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

33、 )n optionsn此結構變數可代表各種最佳化的選項（或參數）Extra parametersExtra parametersMATLAB 程式設計進階篇：一般數學函數的處理與分析設定最佳化選項 (1) n如何設定最佳化選項n用 optimset 指令options = optimset(prop1, value1, prop2, value2, )nprop1、prop2 ：欄位名稱 nvalue1、value2 ：對應的欄位值MATLAB 程式設計進階篇：一般數學函數的處理與分析設定最佳化選項 (2)n印出最佳化步驟的中間結果，並放鬆誤差範圍 options = optimset(Di

34、sp, iter, TolX, 0.1) Display: iter MaxFunEvals: MaxIter: TolFun: TolX: 0.1000 FunValCheck: OutputFcn: PlotFcns: ActiveConstrTol: MATLAB 程式設計進階篇：一般數學函數的處理與分析設定最佳化選項 (3)noptions 共有五十多個欄位n如果欄位值顯示是空矩陣，n使用此欄位的預設值來進行運算 options = optimset(fminbnd)n顯示非空矩陣的最佳化選項：nDisplay: notifynMaxFunEvals: 500nMaxIter: 500

35、nTolX: 1.0000e-004nFunValCheck: offMATLAB 程式設計進階篇：一般數學函數的處理與分析設定最佳化選項 (4)n若輸入 options = optimset(fminsearch) n顯示非空矩陣的最佳化選項：nDisplay: notifynMaxFunEvals: 200*numberofvariables nMaxIter: 200*numberofvariablesnTolFun: 1.0000e-004nTolX: 1.0000e-004nFunValCheck: offMATLAB 程式設計進階篇：一般數學函數的處理與分析最佳化選項說明 (1)

36、nDisplay n若為 0 （預設值），不顯示中間運算結果n若不為 0，則顯示運算過程的中間結果nMaxFunEvals n函數求值運算（Function Evaluation）的最高次數n對 fminbnd 的預設值是 500n對 fminsearch 的預設值是 200 乘上 x0 的長度nMaxIter n最大疊代次數n對 fminbnd 的預設值是 500n對 fminsearch 的預設值是 200 乘上 x0 的長度MATLAB 程式設計進階篇：一般數學函數的處理與分析最佳化選項說明 (2)nTolFun n決定終止搜尋的函數值容忍度n預設為 10-4n此選項只被 fminsea

37、rch 用到，fminbnd 並不使用nTolX n終止搜尋的自變數值容忍度，預設為 10-4 MATLAB 程式設計進階篇：一般數學函數的處理與分析提示n最佳化並非一蹴可及，通常一再重覆，最後才能收斂到最佳點n最佳化的結果和起始點的選定有很大的關聯，一個良好的起始點n加快最佳化收斂的速度n提高找到全域最佳值(Global Optimum)的機會MATLAB 程式設計進階篇：一般數學函數的處理與分析nSteps for Down-hill Simplex search (DSS)nDefine an objective function myFun(x)nSet initial guess “

38、x0”nStart the search via fminsearchnx=fminsearch(myFun, x0);nVariantsnObjective function myFun(x, prm)nStart the searchnx=fminsearch(x) myFun(x, prm), x0);nx=fminsearch(myFun, x0, , prm);More Examples of DSSOld usageExtra parametersTo be passed to myFun()MATLAB 程式設計進階篇：一般數學函數的處理與分析DSS: Example 1 (1/

39、3)nFermat pointnA point P in a plane such that PA+PB+PC is minimized.nSolutionnAnalytic solutionnP satisfies AFB=BFC=CFA=120onNumerical solutionnAll kinds of optimization methodsQuiz!Properties!MATLAB 程式設計進階篇：一般數學函數的處理與分析DSS: Example 1 (2/3) nObjective functionndist2ABC.mnMain program:ngoMinDist2ABC

40、.mfunction sumDistance=dist2ABC(x)% dist2ABC: The distance sum to points A, B, CA=0 0;B=3 0;C=0 4;sumDistance=norm(x-A)+norm(x-B)+norm(x-C);p0=5 5;% Initial guessp=fminsearch(dist2ABC, p0);Bad idea tohardwire the data pointsMATLAB 程式設計進階篇：一般數學函數的處理與分析DSS: Example 1 (3/3) MATLAB 程式設計進階篇：一般數學函數的處理與分析DSS: Example 2 (1/3) nProblem definitionnFind a point P such that the total