matlab程式設計入門篇一般數學函數的處理與分析-mirlab

1、MATLAB 程式設計進階篇一般數學函數的處理與分析,張智星jang@mirlab.orghttp://mirlab.org/jang臺大資工系 多媒體檢索實驗室,函數的函數,MATLAB 可對數學函數進行各種運算與分析，例如：作圖求根優(yōu)化：求函數的極大或極小值數值積分求解微分方程式,如何表示此種被分析的函數？字串函數握把 (function handles)匿名函數 (anonymous function),F

2、unctions ofFunctions!,一維數學函數的範例,MATLAB 是以 M 檔案（副檔名為 m）來表示一個函數例如，內建於MATLAB目錄的 humps.m 可用來計算下列函數：更多資訊：欲顯示此檔案的位置 ? which humps欲顯示此檔案的內容 ? type humps,,提示,MATLAB 常被用到的測試函數humps：單輸入函數peaks：雙輸入函數「函式」和「函數」都代表「function

3、s」，兩者常會混用，若要正名，可區(qū)分如下：函數：通常用來表示「mathematic functions」函式：通常用來表示「subroutines or functions in a programming language」,We use 函數 to represent both.,數學函數的作圖,表示函數的方式函數握把：使用 @humps 來代表 humps.m字串：使用 'humps' 來代表 humps.

4、m用 fplot 指令進行數學函數作圖畫出 humps 函數在 [0,2] 間的曲線範例：fplot01.m,subplot(2,1,1);fplot('humps', [0,2]);% 使用字串指定函式subplot(2,1,2);fplot(@humps, [0 2]);% 使用函式握把來指定函式,,,Less flexible!,同時改變 x、y 的區(qū)間,我們可同時改變 x 和 y 的區(qū)間範例：f

5、plot02.mx 的區(qū)間為[0, 1]y 的區(qū)間為[5, 25],fplot('humps', [0, 1, 5, 25]);grid on% 畫出格線,匿名函式,fplot 也接受匿名函式（當場指定的函式）範例：fplot021.m,subplot(2,1,1);fplot('sin(2*x)+cos(x)', [-10, 10])% 使用字串指定函式subplot(2,1,2

6、);fplot(@(x)sin(2*x)+cos(x), [-10, 10])% 使用函式握把來指定函式,,對多個函數作圖,fplot 也可同時對多個函數作圖，其中每個函數須以一個行向量來表示範例：fplot022.mx 是行向量（MATLAB 預設值）[sin(x), exp(-x)] 是二個行向量每個行向量代表一個函數（即一條曲線）,fplot(@(x)[sin(x), exp(-x)], [0, 10]),,帶有參

7、數的函數,匿名函式也可以帶有參數範例：fplot023.m此時 “@(x)” 不可省略，以便指定自變數,a=1; b=1.1; c=1.2;fplot(@(x)[sin(a*x), sin(b*x), sin(c*x)], [0, 10]),,Function handle is more flexible!,產生 X、Y 座標點,fplot 可進行描點作圖，類似 plot(x, y)，但x 座標點的密度根據函數值的變化決定

8、我們顯示 fplot 所產生的 x 座標點範例：fplot03.m函數變化平緩處，產生稀疏的點函數變化劇烈處，產生緊密的點,[x, y] = fplot(@humps, [-1,2]);plot(x, y, '-o');,,產生更密的 X 座標點 (1),若欲產生更密的 x 座標點，可在 fplot 指令加入另一個輸入引數，已指定相對容忍度(Tolerance)範例：fplot04.m,subplo

9、t (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 座標點 (2),在第一圖中，fplot 指令使用預設相對容忍度，其值為 0.002。在第二圖中，相對容忍度被設為 0.0001，可得到更準確的圖形，但也要花更多計算及作圖時間。,,ezplot 指令,ezplot指

10、令和fplot指令類似，可進行描點作圖，但使用更為簡便，預設的作圖範圍為範例8-7：ezplot01.m,ezplot(@(x)x^3-x^2+x);,,平面中的參數式曲線,ezplot 也可畫出平面中的參數式曲線 範例8-8：ezplot02.m參數式函數的參數預設範圍仍是,ezplot(@(t)sin(3*t), @(t)cos(5*t));,,?利薩如圖形 (Lissajous Figures),空間中

11、的參數式曲線,ezplot3 可畫出空間中的參數式曲線 範例8-8：ezplot021.m參數式函數的參數預設範圍仍是,ezplot3(@(t)sin(3*t), @(t)cos(5*t), @(t)t),,?3D利薩如圖形,隱函數作圖,ezplot 指令可用於隱函數作圖下列範例可以畫出範例8-9：ezplot03.m,ezplot(@(x,y)x^3+2*x^2-3*x-y^2+15);,,函數的求根,fzero

12、 指令用於單變數函數的求根語法x = fzero(fun, x0)fun 是欲求根的函數（以字串或函數握把來表示）x0 是一個起始點或起始區(qū)間,X0 對 fzero 的影響,fzero 指令根據 x0 不同而執(zhí)行下列動作若 x0 為一個起始點 fzero 會自動找出附近包含零點（即根，或函數變號點）的區(qū)間逐步縮小此區(qū)間以找出零點若 fzero 無法找出此區(qū)間，傳回 NaN若已知使函數值不同號的兩點由 x0 直接指定

13、尋根的區(qū)間 fzero 更快速找到位於此區(qū)間內的根,求根範例 (1),找出humps在 x = 1.5 附近的根，並驗算範例8-10：fzero01.mfzero 先找到在 1.5 附近變號的兩點（即 1.26 及 1.6697），然後再找出 humps 的零點,x = fzero(@humps, 1.5);% 求靠近 1.5 附近的根y = humps(x); % 帶入求值fprint

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

15、的根y = humps(x);% 帶入求值fprintf('humps(%f) = %f\n', x , y);,humps(-0.131618) = 0.000000,求根範例 (3),若要畫出以上兩個零點，請見下列範例範例8-12：fzero03.m,fplot(@humps, [-1, 2]); grid onz1 = fzero(@humps, 1.5);z2 = fzero(@humps, [-

16、1, 1]);line(z1, humps(z1), 'marker', 'o', 'color', 'r');line(z2, humps(z2), 'marker', 'o', 'color', 'r');,顯示求解過程的中間結果 (1),MATLAB 可以顯示求解過程的中間結果使用 optimset

17、 指令來設定顯示選項再將 optimset 傳回結構變數送入 fzero範例8-13：fzero04.moptimset 常用於設定最佳化的選項，下一節(jié)會有比較完整的介紹,opt = optimset('disp', 'iter');% 顯示每個 iteration 的結果 a = fzero(@humps, [-1, 1], opt),顯示求解過程的中間結果 (2),,求零點過程中，

18、找下一點的兩個方法顯示在第四個欄位（Procedure 欄位）二分法(Bisection)內插法(Interpolation)可由doc fzero找到所使用的演算法,Func-count x f(x) Procedure 1 -1 -5.13779 initial 2 1

19、 16 I initial 3 -0.513876 -4.02235 interpolation 4 0.243062 71.6382 bisection 5 -0.473635 -3.83767 interpolation 6 -0.115287 0

20、.414441 bisection 7 -0.150214 -0.423446 interpolation 8 -0.132562 -0.0226907 interpolation 9 -0.131666 -0.0011492 interpolation 10 -0.131618 1

21、.88371e-007 interpolation 11 -0.131618 -2.7935e-011 interpolation 12 -0.131618 8.88178e-016 interpolation 13 -0.131618 -9.76996e-015 interpolationZero found in th

22、e interval: [-1, 1].a = -0.1316,數值積分,MATLAB 可用於計算單變函數定積分quad：適應式 Simpson 積分法（Adaptive Simpson Quadrature）quadl：適應式 Lobatto 積分法（Adaptive Lobatto Quadrature）,定積分,計算 humps 在 [0, 1] 的定積分>> q = quad(@humps, 0, 1)

23、q = 29.8583quad 及 quad8 都應用遞迴程序若遞迴次數達 10 次，兩種方法均會傳回 Inf表示所計算之定積分可能不存在quad 及 quad8第四個引數用來指定積分的相對誤差容忍度,曲線的長度 (1),quad 及 quadl 計算曲線的長度一曲線是由下列參數化的方程式來表示 t 的範圍為 [0, 3*pi]範例：plotCurve.m,,t = 0:0.1:3*pi;plot3(sin(

24、2*t), cos(t), t);,曲線的長度 (2),此曲線的長度等於,,,曲線的長度 (3),先定義函數 curveLength.m>> type curveLength.mfunction out = curveLength(t)out = sqrt(4*cos(2*t).^2+sin(t).^2+1);曲線長度可計算如下>> len = quad(@curveLength, 0, 3*pi)l

25、en = 17.2220,雙重積分 (1),dblquad 指令用來計算雙重積分欲計算其中先建立被積分的函數 integrand.m>> type integarnd.mfunction out = integrand(x, y)out = y*sin(x) + x*cos(y);,,雙重積分 (2),計算雙重積分result = dblquad( 'integrand', xMin

26、, xMax, yMin, yMax)其中xMin：內迴圈積分的下界值xMax：內迴圈積分的上界值yMin：外迴圈積分的下界值yMax：外迴圈積分的上界值,雙重積分 (3),範例：dblquad01.m一般的情況下dblquad 會呼叫 quad 計算定積分。若須呼叫更為精確的 quadl，可執(zhí)行下列指令>>result = dblquad(@integrand, xMin, xMax, yMi

27、n, yMax, 'quadl')result = -9.8696,xMin = pi;xMax = 2*pi;yMin = 0;yMax = pi;result = dblquad(@integrand, xMin, xMax, yMin, yMax),result = -9.8698,函數的優(yōu)化,MATLAB 提供了數個基本指令來進行數學函數的優(yōu)化，本節(jié)將介紹：單變數函數的最小化: fminbnd

28、多變數函數的最小化: fminsearch設定最佳化的選項若讀者有興趣使用較複雜的方法，可以使用「最佳化工具箱」（Optimization Toolbox）,單變函數的最小化,fminbnd 指令尋求 humps 在 [0.3, 1] 中的最小值範例：fminbnd01.m最小值發(fā)生在 x = 0.637，且最小值為 11.2528,[x, minValue] = fminbnd(@humps, 0.3, 1),x

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

30、ps, 0.3, 1, opt),尋求最小值的中間過程 (2),,左表列出x值的變化及相對的函數值 f(x) 最後一欄位列出求極小值的方法，包含黃金分割搜尋(Golden Section Search)拋物線內插法(Parabolic Interpolation)x 值誤差小於 10^-4,Func-count x f(x) Procedure 1 0.567

31、376 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

32、 parabolic 6 0.637618 11.2529 parabolic 7 0.636985 11.2528 parabolic 8 0.637019 11.2528 parabolic 9 0.637052 11.2528 parabolicOptim

33、ization terminated successfully: the current x satisfies the termination criteria using OPTIONS.TolX of 1.000000e-004 x = 0.6370minValue = 11.2528,放鬆誤差管制 (1),放鬆誤差管制使 fminbnd 提早傳回計算結果由 optimset 達成下例將 x 的誤差範圍提高為

34、 0.1範例8-16：fminbnd03.m,opt = optimset( 'disp', 'iter', 'TolX', 0.1); [x, minValue] = fminbnd(@humps, 0.3, 1, opt),放鬆誤差管制 (2),,Func-count x f(x) Procedure 1 0.56737

35、6 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 pa

36、rabolic 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 =

37、 11.2693,多變數函數的極小值:fminsearch,fminsearch 指令求多變數函數的極小值必須指定一個起始點，讓 fminsearch 求出在起始點附近的局部最小值（Local Minima）Derivative free ? Less efficientMethod: Downhill Simplex Search (DSS), akaNelder-Mead methodAmoeba method,,

38、Downhill Simplex Search,DSS in WikipediaMany variationsBasic StepsUse a simplex to explore the objective function, with the operations:ReflectionExpansionContractionShrink,,,Downhill Simplex Search,About DSSStren

39、gthsStraightforward conceptEasy programmingNo gradient or derivative neededWeaknessSlowOnly good for continuous objective functionCould be trapped in local minima,,Quiz!,多變數函數的極小值範例,若目標函數為我們必須先產生一個 MATLAB 的函數 o

40、bjective.m 若起始點為範例：fminsearch01.m,,,,x0 = [0, 0, 0];x = fminsearch(@objective, x0),function y = objective(x)y = (x(1)-1)^2 +(x(2)-2)^2 + (x(3)-3)^2;,x = 1.0000 2.0000 3.0000,最佳化選項,MATLAB 最佳化的選項經由另一個輸入引數（

41、Input Argument）來進入 fminbnd 或 fminsearch使用語法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 ) o

42、ptions此結構變數可代表各種最佳化的選項（或參數）,Extra parameters,Extra parameters,設定最佳化選項 (1),如何設定最佳化選項用 optimset 指令options = optimset(prop1, value1, prop2, value2, …)prop1、prop2 ：欄位名稱 value1、value2 ：對應的欄位值,設定最佳化選項 (2),印出最佳化步驟的中間結果，並放鬆

43、誤差範圍>> options = optimset('Disp', 'iter', 'TolX', 0.1) Display: 'iter' MaxFunEvals: [] MaxIter: []

44、 TolFun: [] TolX: 0.1000 FunValCheck: [] OutputFcn: [] PlotFcns: [] ActiveConstrTol: []

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

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

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

48、是 500對 fminsearch 的預設值是 200 乘上 x0 的長度,最佳化選項說明 (2),TolFun 決定終止搜尋的函數值容忍度預設為 10^-4此選項只被 fminsearch 用到，fminbnd 並不使用TolX 終止搜尋的自變數值容忍度，預設為 10-4,提示,最佳化並非一蹴可及，通常一再重覆，最後才能收斂到最佳點最佳化的結果和起始點的選定有很大的關聯，一個良好的起始點加快最佳化收斂的速度提高找到全

49、域最佳值(Global Optimum)的機會,Steps for Down-hill Simplex search (DSS)Define an objective function myFun(x)Set initial guess “x0”Start the search via fminsearchx=fminsearch(@myFun, x0);VariantsObjective function myFun(x,

50、 prm)Start the searchx=fminsearch(@(x) myFun(x, prm), x0);x=fminsearch(@myFun, x0, [], prm);,More Examples of DSS,Old usage,Extra parametersTo be passed to myFun(),DSS: Example 1 (1/3),Fermat pointA point P in a pla

51、ne such that PA+PB+PC is minimized.SolutionAnalytic solutionP satisfies AFB=BFC=CFA=120oNumerical solutionAll kinds of optimization methods,Quiz!,Properties!,DSS: Example 1 (2/3),Objective functiondist2ABC.mMai

52、n program:goMinDist2ABC.m,function 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(

53、@dist2ABC, p0);…,Bad idea tohardwire the data points,DSS: Example 1 (3/3),,DSS: Example 2 (1/3),Problem definitionFind a point P such that the total distance to a set of points is minimized.SolutionAnalytic solution

54、Not sure if it exists…Numerical solutionAll kinds of optimization methods,DSS: Example 2 (2/3),Objective functiondist2points.mMain program:goMinDist2points01.m,function sumDistance=dist2points(x, points)% dist2

55、points: The sum of distance sum to pointssumDistance=0;for i=1:size(points,2)sumDistance=sumDistance+norm(x-points(:,i));end,points=[-1 -1; 3 0; 0 4; 1 1; 2 5; 4 2]’;p0=[5 5]’;% Initial guessp=fminsearch(@(x) dis

56、t2points(x, points), p0);…,It’s Better to acceptthe data points externally,Data matrixof 2xn,DSS: Example 2 (3/3),,DSS: Example 3 (1/2),Objective functionminDist2points02.mMain program:goMinDist2points02.m,funct

57、ion p=minDist2points02(points, x0, showPlot)% minDist2points: Compute the point that has the min total …if nargin<2||isempty(x0), x0=mean(points, 2); end % Initial guessp=fminsearch(@(x)dist2points(x, points),