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

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

2、unctions ofFunctions!,一維數(shù)學(xué)函數(shù)的範(fàn)例,MATLAB 是以 M 檔案（副檔名為 m）來(lái)表示一個(gè)函數(shù)例如，內(nèi)建於MATLAB目錄的 humps.m 可用來(lái)計(jì)算下列函數(shù)：更多資訊：欲顯示此檔案的位置 ? which humps欲顯示此檔案的內(nèi)容 ? type humps,,提示,MATLAB 常被用到的測(cè)試函數(shù)humps：?jiǎn)屋斎牒瘮?shù)peaks：雙輸入函數(shù)「函式」和「函數(shù)」都代表「function

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

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

5、plot02.mx 的區(qū)間為[0, 1]y 的區(qū)間為[5, 25],fplot('humps', [0, 1, 5, 25]);grid on% 畫(huà)出格線,匿名函式,fplot 也接受匿名函式（當(dāng)場(chǎng)指定的函式）範(fàn)例：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])% 使用函式握把來(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]),,帶有參

7、數(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]),,Function handle is more flexible!,產(chǎn)生 X、Y 座標(biāo)點(diǎn),fplot 可進(jìn)行描點(diǎn)作圖，類(lèi)似 plot(x, y)，但x 座標(biāo)點(diǎn)的密度根據(jù)函數(shù)值的變化決定

8、我們顯示 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),若欲產(chǎn)生更密的 x 座標(biāo)點(diǎn)，可在 fplot 指令加入另一個(gè)輸入引數(shù)，已指定相對(duì)容忍度(Tolerance)範(fàn)例：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);,產(chǎn)生更密的 X 座標(biāo)點(diǎn) (2),在第一圖中，fplot 指令使用預(yù)設(shè)相對(duì)容忍度，其值為 0.002。在第二圖中，相對(duì)容忍度被設(shè)為 0.0001，可得到更準(zhǔn)確的圖形，但也要花更多計(jì)算及作圖時(shí)間。,,ezplot 指令,ezplot指

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

11、的參數(shù)式曲線,ezplot3 可畫(huà)出空間中的參數(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)例可以畫(huà)出範(fàn)例8-9：ezplot03.m,ezplot(@(x,y)x^3+2*x^2-3*x-y^2+15);,,函數(shù)的求根,fzero

12、 指令用於單變數(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)作若 x0 為一個(gè)起始點(diǎn) fzero 會(huì)自動(dòng)找出附近包含零點(diǎn)（即根，或函數(shù)變號(hào)點(diǎn)）的區(qū)間逐步縮小此區(qū)間以找出零點(diǎn)若 fzero 無(wú)法找出此區(qū)間，傳回 NaN若已知使函數(shù)值不同號(hào)的兩點(diǎn)由 x0 直接指定

13、尋根的區(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); % 帶入求值fprint

14、f('humps(%f) = %f\n', x , y);,humps(1.299550) = 0.000000,求根範(fàn)例 (2),若已知 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]

15、的根y = humps(x);% 帶入求值fprintf('humps(%f) = %f\n', x , y);,humps(-0.131618) = 0.000000,求根範(fàn)例 (3),若要畫(huà)出以上兩個(gè)零點(diǎn)，請(qǐng)見(jiàn)下列範(fàn)例範(fàn)例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');,顯示求解過(guò)程的中間結(jié)果 (1),MATLAB 可以顯示求解過(guò)程的中間結(jié)果使用 optimset

17、 指令來(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),顯示求解過(guò)程的中間結(jié)果 (2),,求零點(diǎn)過(guò)程中，

18、找下一點(diǎn)的兩個(gè)方法顯示在第四個(gè)欄位（Procedure 欄位）二分法(Bisection)內(nèi)插法(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,數(shù)值積分,MATLAB 可用於計(jì)算單變函數(shù)定積分quad：適應(yīng)式 Simpson 積分法（Adaptive Simpson Quadrature）quadl：適應(yīng)式 Lobatto 積分法（Adaptive Lobatto Quadrature）,定積分,計(jì)算 humps 在 [0, 1] 的定積分>> q = quad(@humps, 0, 1)

23、q = 29.8583quad 及 quad8 都應(yīng)用遞迴程序若遞迴次數(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.m,,t = 0:0.1:3*pi;plot3(sin(

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

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

26、, xMax, yMin, yMax)其中xMin：內(nèi)迴圈積分的下界值xMax：內(nèi)迴圈積分的上界值yMin：外迴圈積分的下界值yMax：外迴圈積分的上界值,雙重積分 (3),範(fàn)例：dblquad01.m一般的情況下dblquad 會(huì)呼叫 quad 計(jì)算定積分。若須呼叫更為精確的 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,函數(shù)的優(yōu)化,MATLAB 提供了數(shù)個(gè)基本指令來(lái)進(jìn)行數(shù)學(xué)函數(shù)的優(yōu)化，本節(jié)將介紹：?jiǎn)巫償?shù)函數(shù)的最小化: fminbnd

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

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

30、ps, 0.3, 1, opt),尋求最小值的中間過(guò)程 (2),,左表列出x值的變化及相對(duì)的函數(shù)值 f(x) 最後一欄位列出求極小值的方法，包含黃金分割搜尋(Golden Section Search)拋物線內(nèi)插法(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 提早傳回計(jì)算結(jié)果由 optimset 達(dá)成下例將 x 的誤差範(fàn)圍提高為

34、 0.1範(fàn)例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,多變數(shù)函數(shù)的極小值:fminsearch,fminsearch 指令求多變數(shù)函數(shù)的極小值必須指定一個(gè)起始點(diǎn)，讓 fminsearch 求出在起始點(diǎn)附近的局部最小值（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!,多變數(shù)函數(shù)的極小值範(fàn)例,若目標(biāo)函數(shù)為我們必須先產(chǎn)生一個(gè) MATLAB 的函數(shù) o

40、bjective.m 若起始點(diǎn)為範(fàn)例：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,最佳化選項(xiàng),MATLAB 最佳化的選項(xiàng)經(jīng)由另一個(gè)輸入引數(shù)（

41、Input Argument）來(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 ) o

42、ptions此結(jié)構(gòu)變數(shù)可代表各種最佳化的選項(xiàng)（或參數(shù)）,Extra parameters,Extra parameters,設(shè)定最佳化選項(xiàng) (1),如何設(shè)定最佳化選項(xiàng)用 optimset 指令options = optimset(prop1, value1, prop2, value2, …)prop1、prop2 ：欄位名稱(chēng) value1、value2 ：對(duì)應(yīng)的欄位值,設(shè)定最佳化選項(xiàng) (2),印出最佳化步驟的中間結(jié)果，並放鬆

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

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

45、 …,設(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è)

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

47、04FunValCheck: 'off',最佳化選項(xiàng)說(shuō)明 (1),Display 若為 0 （預(yù)設(shè)值），不顯示中間運(yùn)算結(jié)果若不為 0，則顯示運(yùn)算過(guò)程的中間結(jié)果MaxFunEvals 函數(shù)求值運(yùn)算（Function Evaluation）的最高次數(shù)對(duì) fminbnd 的預(yù)設(shè)值是 500對(duì) fminsearch 的預(yù)設(shè)值是 200 乘上 x0 的長(zhǎng)度MaxIter 最大疊代次數(shù)對(duì) fminbnd 的預(yù)設(shè)值

48、是 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)加快最佳化收斂的速度提高找到全

49、域最佳值(Global Optimum)的機(jī)會(huì),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),