MATLAB 程式設計進階篇 一般數學函數的處理與分析 張智星 jang@cs.nthu.edu.tw http://www.cs.nthu.edu.tw/~jang 清大資工系 多媒體檢索實驗室

8-1 MATLAB 的函數表示法 函數(Functions) MATLAB 指令須以字串形式輸入函數或是函式握把（Function Handles）的方式輸入此函數 MATLAB 指令可對數學函數進行運算與分析 數值積分 最佳化（Optimization，求函數的極大極小值） 求解非線性方程式的根 求解微分方程式

MATLAB 的數學函數範例 MATLAB 數學函數是以 M 檔案（副檔名為 m）表示 顯示其位置，可輸入 Ex.以下的數學式子，寫成檔案名稱為 humps.m 的 MATLAB 數學函數 顯示其位置，可輸入 >> which humps C:\Program Files\MATLAB71\toolbox\matlab\demos\humps.m 上述的函數已內建於 MATLAB 目錄之下

提示 顯示 humps 函數的內容，可輸入 humps 是單輸入函數，peaks 是雙輸入函數，兩者都是 MATLAB 常被用到的測試函數 >> type humps humps 是單輸入函數，peaks 是雙輸入函數，兩者都是 MATLAB 常被用到的測試函數 「函式」和「函數」都代表「functions」，通常「函數」表「mathematic function」，「函式」表「subroutines or functions in a programming language」，兩者有時候會混用

8-2數學函數的作圖 使用函式握把（Function Handle）來表示：這是比較新的方式，適用於 MATLAB 6.x & 7.x，例如使用 @humps 來代表 humps.m 函式。 用 fplot 指令進行數學函數作圖 畫出 humps 函數在 [0,2] 間的曲線 範例8-1：fplot01.m subplot(2,1,1); fplot('humps', [0,2]); % 使用字串指定函式 subplot(2,1,2); fplot(@humps, [0 2]); % 使用函式握把來指定函式

改變 x、y 的區間 可同時改變 x 和 y 的區間 x 的區間為[0, 1] 範例8-2：fplot02.m y 的區間為[5, 25] grid on 指令用於畫出 x、y 平面的格線(Grid) fplot('humps', [0 1 5 25]); grid on % 畫出格線

字串的數學方程式作圖 fplot 也接受「當場表示」的數學函數 範例8-3：fplot021.m subplot(2,1,1); fplot('sin(2*x)+cos(x)', [-10, 10]) % 使用字串指定函式 subplot(2,1,2); fplot(@(x)sin(2*x)+cos(x), [-10, 10]) % 使用函式握把來指定函式

對多個函數作圖 fplot 也可同時對多個函數作圖 函數須放入一個向量 範例8-4：fplot022.m x 是行向量 (Column Vector) [sin(x), exp(-x)] 是二個行向量 每個行向量代表一個函 數（即一條曲線） fplot(@(x)[sin(x), exp(-x)], [0, 10])

產生 X 座標點 x 座標點的密度根據函數值的變化決定 fplot 是描點作圖，類似 plot(x, y) 顯示 fplot 所產生的 x 座標點 範例8-5：fplot03.m 函數變化平緩處，產生 稀疏的點 函數變化劇烈處，產生 緊密的點 [x, y] = fplot(@humps, [-1,2]); plot(x, y, '-o');

產生更密的 X 座標點 (I) 欲產生更密的 x 座標點 fplot指令加入相對容忍度(Tolerance)的輸入引數 範例8-6：fplot04.m subplot (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);

產生更密的 X 座標點 (II) 第一圖，fplot 指令使用預設相對容忍度，其值為 0.002 第二圖，相對容忍度被設為 0.0001 可得更準確的圖形，但也要花更多計算及作圖時間

ezplot 指令 ezplot指令和fplot指令類似，但使用更為簡便 若要畫出圖形 範例8-7：ezplot01.m 參數式函數的參數預設範圍為 -2*pi 到 2*pi 之間 ezplot(@(x)x^3-x^2+x);

參數式曲線 ezplot 也可畫出參數式的曲線 範例8-8：ezplot02.m 參數式函數的參數預設範圍仍是 -2*pi 到 2*pi 之間 ezplot(@(t)sin(3*t), @(t)cos(5*t));

隱函數作圖 ezplot 指令可用於隱函數作圖 範例8-9：ezplot03.m ezplot(@(x,y)x^3+2*x^2-3*x+5-y^2+10);

8-3 函數的求根 fzero 指令 用於單變數函數的求根 使用語法 x = fzero(fun, x0)

x0 不同對fzero的影響 fzero 指令根據 x0 不同而執行下列動作 若 x0 為一個起始點 逐步縮小此區間以找出零點 若 fzero 無法找出此區間，傳回 NaN 若已知使函數值不同號的兩點 由 x0 直接指定尋根的區間 fzero 更快速找到位於此區間內的根

求根範例 -1 找出humps在 x = 1.5 附近的根，並驗算 範例8-10：fzero01.m x = fzero(@humps, 1.5); % 求靠近 1.5 附近的根 y = humps(x); % 帶入求值 fprintf('humps(%f) = %f\n', x , y); humps(1.299550) = 0.000000

求根範例 -2 若已知 humps 在 x = -1 及 1 間為異號 令 x0 = [-1, 1] 為起始區找出 humps 的零點 範例8-11：fzero02.m 此時 fzero 找到的是另一個零點 x = fzero(@humps, [-1, 1]); % 求落於區間 [-1, 1] 的根 y = humps(x); % 帶入求值 fprintf('humps(%f) = %f\n', x , y); humps(-0.131618) = 0.000000

求根範例 –3 (I) 若要畫出以上兩個零點 範例8-12：fzero03.m fplot(@humps, [-1, 2]); z1 = fzero(@humps, 1.5); z2 = fzero(@humps, [-1, 1]); line(z1, humps(z1), 'marker', 'o'); % 畫出第一個零點的位置 line(z2, humps(z2), 'marker', 'o'); % 畫出第二個零點的位置 grid on % 畫出格線

求根範例 –3 (II)

求解過程的中間結果 (I) 顯示求解過程的中間結果 使用 optimset 指令，再將 optimset 傳回的結構變數送入fzero 範例8-13：fzero04.m opt = optimset('disp', 'iter'); % 顯示每個 iteration 的結果 a = fzero(@humps, [-1, 1], opt)

求解過程的中間結果 (II) 求零點過程中，用到二分法(Bisection)或內插法(Interpolation) 若給定一個起始點 Func-count x f(x) Procedure 1 -1 -5.13779 initial 2 1 16 initial 3 -0.513876 -4.02235 interpolation 4 0.243062 71.6382 bisection 5 -0.473635 -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 -9.76996e-015 interpolation Zero found in the interval: [-1, 1]. a = -0.1316 求零點過程中，用到二分法(Bisection)或內插法(Interpolation) 方法在左表Procedure的第四個欄位中 若給定一個起始點 fzero 最初的步驟中會先找出一個適當的區間來搜尋零點

提示 optimset 常用於設定最佳化的選項

8-4 函數的極小值 MATLAB 指令進行數學函數最佳化，將介紹 單變數函數的最小化 多變數函數的最小化 設定最佳化的選項

單變函數的最小化 fminbnd 指令 尋求 humps 在 [0.3, 1] 中的最小值 範例8-14：fminbnd01.m 最小值發生在 x = 0.637，且最小值為 11.2528 [x, minValue] = fminbnd(@humps, 0.3, 1) x = 0.6370 minValue = 11.2528

尋求最小值的中間過程 (I) 尋求最小值的中間過程 使用 optimset 指令來設定顯示選項 再將 optimset 傳回結構變數送入 fminbnd 範例8-15：fminbnd02.m opt = optimset('disp', 'iter'); % 顯示每個步驟的結果 [x, minValue] = fminbnd(@humps, 0.3, 1, opt)

尋求最小值的中間過程 (II) 左表列出x值的變化及相對的函數值 f(x) 最後一欄位列出求極小值的方法 通常是 x 值誤差小於 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.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 parabolic Optimization terminated successfully: the current x satisfies the termination criteria using OPTIONS.TolX of 1.000000e-004 x = 0.6370 minValue = 11.2528 左表列出x值的變化及相對的函數值 f(x) 最後一欄位列出求極小值的方法 通常是 黃金分割搜尋 (Golden Section Search) 拋物線內插法 (Parabolic Interpolation) x 值誤差小於 10^-4

放鬆誤差管制 (I) 放鬆誤差管制 使 fminbnd 提早傳回計算結果 由 optimset 達成 範例8-16將 x 的誤差範圍提高為 0.1 範例8-16：fminbnd03.m opt = optimset( 'disp', 'iter', 'TolX', 0.1); [x, minValue] = fminbnd(@humps, 0.3, 1, opt)

放鬆誤差管制 (II) 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.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.6444 minValue = 11.2693

多變數函數的極小值 fminsearch 指令 求多變數函數的極小值 必須指定一個起始點，讓 fminsearch 求出在起始點附近的局部最小值（Local Minima） 目標函數： 先產生一個 MATLAB 的函數 objective.m >> type objective.m function y = objective(x) y = (x(1)-1).^2 +(x(2)-2).^2 + (x(3)-3).^2;

多變數函數的極小值 範例 若起始點為 ， ， 範例8-17：fminsearch01.m 若起始點為 ， ， 範例8-17：fminsearch01.m fminsearch 使用的方法是下坡式 Simplex 搜尋（Downhill Simplex Search） x0 = [0, 1, 2]; x = fminsearch(@objective, x0) x = 1.0000 2.0000 3.0000

最佳化選項 MATLAB 最佳化的選項 經由另一個輸入引數（Input Argument）來進入 fminbnd 或 fminsearch 使用語法 x = fminbnd(@function, x1, x2, options) 或 x = fminsearch(@function, x0, options ) options 結構變數 代表各種最佳化的選項（或參數）

設定最佳化選項 欲設定最佳化選項 用 optimset 指令 options = optimset(prop1, value1, prop2, value2, …) prop1、prop2 ：欄位名稱 value1、value2 ：對應的欄位值

設定最佳化選項 -1 印出最佳化步驟的中間結果，並放鬆誤差範圍 >> options = optimset('Disp', 'iter', 'TolX', 0.1) Display: 'iter' MaxFunEvals: [] MaxIter: [] TolFun: [] ………… TolRLPFun: [] TolXInteger: [] TypicalX: []

設定最佳化選項 -2 options 共有五十多個欄位 如果欄位值顯示是空矩陣， 使用此欄位的預設值來進行運算 顯示非空矩陣的最佳化選項： >> options = optimset('fminbnd') 顯示非空矩陣的最佳化選項： Display: 'notify MaxFunEvals: 500 MaxIter: 500 TolX: 1.0000e-004 FunValCheck: 'off'

設定最佳化選項 -3 若輸入 顯示非空矩陣的最佳化選項： Display: 'notify‘ >> options = optimset('fminsearch') 顯示非空矩陣的最佳化選項： Display: 'notify‘ MaxFunEvals: '200*numberofvariables' MaxIter: '200*numberofvariables‘ TolFun: 1.0000e-004 TolX: 1.0000e-004 FunValCheck: 'off'

最佳化選項說明 (I) Display MaxFunEvals MaxIter 若為 0 （預設值），不顯示中間運算結果 若不為 0，則顯示運算過程的中間結果 MaxFunEvals 函數求值運算（Function Evaluation）的最高次數 對 fminbnd 的預設值是 500 對 fminsearch 的預設值是 200 乘上 x0 的長度 MaxIter 最大疊代次數

最佳化選項說明 (II) TolFun TolX 決定終止搜尋的函數值容忍度 預設為 10^-4 此選項只被 fminsearch 用到，fminbnd 並不使用 TolX 終止搜尋的自變數值容忍度，預設為 10-4

提示 最佳化並非一蹴可及，通常一再重覆，最後才能收斂到最佳點 最佳化的結果和起始點的選定有很大的關聯，一個良好的起始點 加快最佳化收斂的速度 提高找到全域最佳值(Global Optimum)的機會

8-5 數值積分 MATLAB 可用於計算單變函數定積分 quad：適應式 Simpson 積分法（Adaptive Simpson Quadrature） quadl：適應式 Lobatto 積分法（Adaptive Lobatto Quadrature）

定積分 計算 humps 在 [0, 1] 的定積分 quad 及 quad8 都應用遞迴程序 >> q = quad(@humps, 0, 1) q = 29.8583 quad 及 quad8 都應用遞迴程序 若遞迴次數達 10 次，兩種方法均會傳回 Inf 表示所計算之定積分可能不存在 quad 及 quad8第四個引數用來指定積分的相對誤差容忍度

曲線的長度 (I) quad 及 quadl 計算曲線的長度 一曲線是由下列參數化的方程式來表示 t 的範圍為 [0, 3*pi] 範例8-18：plotCurve.m t = 0:0.1:3*pi; plot3(sin(2*t), cos(t), t);

曲線的長度 (II) 此曲線的長度等於

曲線的長度 (III) 先定義函數 curvleng.m 曲線長度可計算如下 >> type curvleng.m function out = curvleng(t) out = sqrt(4*cos(2*t).^2+sin(t).^2+1); 曲線長度可計算如下 >> len = quad(@curvleng, 0, 3*pi) len = 17.2220

雙重積分 (I) dblquad 指令 用來計算雙重積分 欲計算 其中 先建立被積分的函數 integrnd.m >> type integrnd.m function out = integrnd(x, y) out = y*sin(x) + x*cos(y);

雙重積分 (II) 計算雙重積分 其中 Xmin：內迴圈積分的下界值 Xmax：內迴圈積分的上界值 Ymin：外迴圈積分的下界值 result = dblquad( 'integrnd', Xmin, Xmax, Ymin, Ymax) 其中 Xmin：內迴圈積分的下界值 Xmax：內迴圈積分的上界值 Ymin：外迴圈積分的下界值 Ymax：外迴圈積分的上界值

雙重積分 (III) 範例8-19：dblquad01.m 一般的情況下dblquad 會呼叫 quad 計算定積分。若須呼叫更為精確的 quadl，可執行下列指令 >>result = dblquad(@integrnd,Xmin,Xmax,Ymin,Ymax,[],'quadl) result = -9.8696 Xmin = pi; Xmax = 2*pi; Ymin = 0; Ymax = pi; result = dblquad(@integrnd, Xmin, Xmax, Ymin, Ymax) result = -9.8698

8-6 本章指令彙整 類別 指令名稱 功能 作圖 ezplo fplot 簡便的函數作圖 一般的函數作圖 求根 fzero 單變數函數的求根 最佳化 fminbnd fminsearch 單變數函數的最小值多變數函數的最小值 數值積分(Quadrature) quad quad8 dblquad 低精度的數值積分 高精度的數值積分 雙重（二維）積分