Mathematica

Mathematica 戴忠淵 [email protected] http://chungyuandye.blogspot.com

July 30, 2009

Contents

1

2

微積分的基本操作 . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .

1

1.1

極限 . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .

1

1.2

微分 . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .

2

1.3

全微分 . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .

5

1.4

積分 . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .

6

1.5

級數 . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .

18

非線性最佳化 . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 23 2.1

2.2

方程式求解 . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .

23

2.1.1 Solve . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .

23

2.1.2 FindRoot . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .

26

非線性最佳化 . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .

34

2.2.1 無限制條件下最佳化求解 . . . . . . . . . . . . . . . . . . . . . . . . . . . .

36

第1章 微積分的基本操作

1.1 極限 Mathematica計算極限的指令為Limit[f, x]。該指令的用法有以下幾種 1. Limit[f, x->a]：計算函數f在x → a的極限值 2. Limit[f, x->a, Direction -> 1]：計算函數f在x → a的右極限值 3. Limit[f, x->a, Direction ->-1]：計算函數f在x → a的左極限值

範例 1-1. sin(x) x→0 x

計算 lim

In[1]: Limit[Sin[x]/x, x -> 0] Out[1]: 1

計算 lim + x→−2

|x + 2| x+2

In[2]: Limit[Abs[x + 2]/(x + 2), x -> -2, Direction -> 1] Out[2]: −1

計算 lim − x→−2

|x + 2| x+2

In[3]: Limit[Abs[x + 2]/(x + 2), x -> -2, Direction -> -1] Out[3]: 1

1.2 微分

計算 lim

x→−2

2

|x + 2| x+2

In[4]: Limit[Abs[x + 2]/(x + 2), x -> -2] Out[4]: 1

在上式中很明顯發現Mathematica在計算極限值時，並不會考慮到極限的存在性。所以在函數極限發 散或不存在時，使用Limit指令必須更為小心。

1.2 微分 在Mathematica中計算函數的微分或積分是非常方便的，其指令為D[f,x]，表示函數f對x做或偏微 分。該指令的用法有以下幾種： 1. D[f, x]：函數f對x微分，

∂f ∂x

2. D[f, x1 , x2 , . . . , xn ]： 函數f對x1 , x2 , . . . 偏微分， 3. D[f,{x,n}]：函數f對x做n次微分，

∂n f ∂xn

∂n f ∂x1 ∂x2 · · · ∂xn {

} ∂f ∂f ∂f 4. D[f, {{x1 , x2 , . . . ,xn }]： 函數f分別對x1 , x2 , . . . 偏微分， , ,..., ，即傳回函 ∂x1 ∂x2 ∂xn 數f的梯度。必須注意的是此處所傳回的為一陣列 對於單變數函數的微分，Mathematica有個方便的技巧，就是直接以f’[x]表示即可；同理，f”[x]則表 示函數f(x)的二階導數，以此類推。 範例 1-2.

定義函數f(x) = xex

2

In[1]: f[x_]:=x*Exp[x^2] Out[1]:

傳回一階、二階和三階導數 In[2]: {f’[x],f’’[x],f’’’[x]} { 2 } 2 2 2 2 2 2 Out[2]: ex + 2ex x2 , 6ex x + 4ex x3 , 6ex + 24ex x2 + 8ex x4

CHAPTER 1. 微積分的基本操作

3

傳回f(x)在x = 0的5階泰勒展開式 In[3]: Plus@@Table[(D[f[x],{x,n}]/.x->0)/n!*x^n,{n,0,5}] 1 1 1 1 (5) f (0)x5 + f(4) (0)x4 + f(3) (0)x3 + x2 f ′′ (0) + xf ′ (0) + f(0) Out[3]: 120 24 6 2 2

若將函數f(x) = xex 帶入上式，可得f(x)在x = 0的5階泰勒展開式 In[4]: Plus@@Table[(D[f[x],{x,n}]/.x->0)/n!*x^n,{n,0,5}] x5 Out[4]: + x3 + x 2

以上為Mathematica對單變數函數微分的方式，多變數函數微分在Mathematica中也是非常簡單，僅 需要注意微分的順序即可。

範例 1-3. f(x, y) = xyexy

定義函數f(x, y) = xyexy In[1]: f[x_,y_]:=x*y*Exp[x*y] Out[1]:

函數f(x, y)對x做偏導數，即求

∂f ∂x

In[2]: D[f[x,y],x] Out[2]: exy xy2 + exy y

函數f(x, y)對x做二階偏微分，

∂2 f ∂x2

In[3]: D[f[x,y],x,x] Out[3]: exy xy3 + 2exy y2

以下指令同樣也是對x做二階偏微分 In[4]: D[f[x,y],{x,2}] Out[4]: exy xy3 + 2exy y2

1.2 微分

4

函數f(x, y)分別對x, y偏微分，傳回梯度 In[5]: D[f[x,y],{{x,y}}] { } Out[5]: exy xy2 + exy y, exy yx2 + exy x

上式結果分別再對x, y做偏微分即傳回海賽矩陣（Hessian Matrix），並以矩陣方式表示 In[6]: D[%, {{x, y}}] // MatrixForm   xy 3 xy 2 xy 2 2 xy xy e x y + 3e xy + e  e xy + 2e y  Out[6]:   exy x2 y2 + 3exy xy + exy exy yx3 + 2exy x2

運用函數D[f,{x1 , n}]和D[f, {{x1 ,x2 ,. . . ,xn }}]的概念，我們可以把程式寫得更簡潔 In[7]: D[f[x,y],{{x,y},2}]//MatrixForm   xy 3 xy 2 xy 2 2 xy xy Out[7]: e x y + 3e xy + e  e xy + 2e y    xy 2 2 xy xy xy 3 xy 2 e x y + 3e xy + e e yx + 2e x

產生f(x, y)的多變數麥克勞尼展開式前3階各項元素 In[8]: Table[Binomial[n,i]*(D[f[x,y],{x,i},{y,n-i}]/.{x->0,y->0})* x^i*y^(n-i)/n!,{n,0,3},{i,0,n}]//TableForm Out[8]: f(0, 0) yf(0,1) (0, 0) 1 2 (0,2) y f (0, 0) 2 1 3 (0,3) y f (0, 0) 6

xf(1,0) (0, 0) xyf(1,1) (0, 0) 1 2 (1,2) xy f (0, 0) 2

1 2 (2,0) x f (0, 0) 2 1 2 (2,1) x yf (0, 0) 2

1 3 (3,0) x f (0, 0) 6

將上式所有元素相加即為f(x, y)在(x, y) = (0, 0)的3階泰勒展開式 In[9]: Plus@@(%//Flatten) 1 1 Out[9]: f[0, 0] + yf(0,1) [0, 0] + y2 f(0,2) [0, 0] + y3 f(0,3) [0, 0] + xf(1,0) [0, 0] 2 6 1 2 (1,2) 1 2 (2,0) 1 1 (1,1) +xyf [0, 0] + xy f [0, 0] + x f [0, 0] + x2 yf(2,1) [0, 0] + x3 f(3,0) [0, 0] 2 2 2 6

CHAPTER 1. 微積分的基本操作

5

1.3 全微分 在上一小節中，我們介紹如何用Mathematica來求解y′ ，然而並非所有的y都可以表示成x的封閉解 dy （close-form）。當y無法表示成x的封閉解時，要求解 則必須倚賴隱函數微分。Mathematica中要進 dx 行隱含數微分則必須借助Dt這個全微分指令。全微分指令的用法有以下幾種： 1. Dt[f, x]：函數f對x做全微分 2. Dt[f]： 函數f對做全微分 3. Dt[f,{x,n}]：函數f對x做n階全微分 4. Dt[f, {{x1 , x2 , . . . ,xn }]： 函數f分別對x1 , x2 , . . . 做全微分

範例 1-4. 函數f(x, y) = 4x2 + 2xy − xy3

定義函數f(x, y) In[1]: f[x_, y_] := 4x^2 + 2x*y - x*y^3 Out[1]:

函數f(x, y)對x做全微分 In[2]: Dt[f[x,y],x] Out[2]: 8x + 2y − y3 + 2xDt[y, x] − 3xy2 Dt[y, x]

上式中Dt[y, x]即為

dy dy ，接下來我們用指令Solve來求解 dx dx

In[3]: Solve[%==0,Dt[y,x]] {{ }} 8x + 2y − y3 Out[3]: Dt[y, x] → x (−2 + 3y2 )

函數f(x, y)做全微分 In[4]: Dt[f[x,y]] Out[4]: 8xDt[x] + 2yDt[x] − y3 Dt[x] + 2xDt[y] − 3xy2 Dt[y]

1.4 積分

6

上式中Dt[x]和Dt[y]分別表示dx和dy，移項後也可傳回與Dt[y, x]相同之結果 In[5]: Solve[% == 0, Dt[y]] {{ }} 8xDt[x] + 2yDt[x] − y3 Dt[x] Out[5]: Dt[y] → x (−2 + 3y2 )

1.4 積分 在Mathematica中計算函數的積分指令為Integrate，該指令的用法有以下幾種： ∫ 1. Integrate[f, x]： 傳回函數f的不定積分，即 f(x)dx 2. Integrate[f, {x, a, b}]： 傳回函數f在區間(a, b)間的定積分值 3. Integrate[f, {x1 , a, b}, {x2 , c, d}, . . . ] 傳回函數f在區間(a, b) × (c, d) × · · · 的定積分值。要 注意的是x1 為最外面的積分變數，x2 為倒數第二個的積分變數，即 ∫ ∫b ∫d · · · f(x1 , x2 , . . . , xn )dxn · · · dx2 dx1 a

∫ 範例 1-5. 計算

計算

c

1 dx xp

1 的不定積分 xp

In[1]: Integrate[1/x^p,x] x1−p Out[1]: 1−p

由於當p > 1時積分會發散，所以Mathematica會給於提示 In[2]: Integrate[1/x^p, {x,0,1}] [ ] [ ] 1 Out[2]: If Re(p) < 1, , Integrate x−p , {x, 0, 1}, Assumptions → Re(p) > 1 1−p

如果要避免發生錯誤，可以在指令後面加上一個假設參數，Assumptions->{0{0
CHAPTER 1. 微積分的基本操作

7

另外一個方式為設定參數GenerateConditions為False，要求Mathematica不產生條件判斷式 In[4]: Integrate[1/x^p, {x,0,1}, GenerateConditions->True] 1 Out[4]: 1−p

若希望Mathematica在未來的計算過程中都不要自動產生判斷式，我們可以藉由改變Integrate函數的 屬性來達成。

要求Mathematica在未來所有的計算過程都不產生條件判斷式 In[6]: SetOptions[Integrate,GenerateConditions->False]

∫∫ 範例 1-6.

(x2 − y)dydx，其中Ω = {(x, y)| − x2 6 y 6 x2 , −1 6 x 6 1}

Ω

繪出Ω的範圍 In[1]: RegionPlot[y<=x^2&&y>=-x^2,{x,-1,1},{y,-1,1}] Out[1]: 1.0

0.5

0.0

-0.5

-1.0 -1.0

-0.5

0.0

0.5

1.0

∫ 1 ∫ x2 我們可以把積分寫為 x 2 − yd yd x −1 −x2

In[2]: Integrate[x^2-y,{x,-1,1},{y,-x^2,x^2}] 4 Out[2]: 5

1.4 積分

8

我們也可以配合布爾函數來計算積分值 In[3]: Integrate[(x^2-y)*Boole[-x^2<=y&&y<=x^2&&-1<=x&&x<=1], {y,-1,1},{x,-1,1}] 4 Out[3]: 5

此處的Boole為一指標函數，若(x, y)在Ω內則傳回1，否則傳為0，即   1, (x, y) ∈ Ω Boole[-x^2<=y&&y<=x^2&&-1<=x&&x<=1] =  0, (x, y) ∈ /Ω 因此，當積分範圍複雜時，搭配布爾函數則可不考慮積分順序。

範例 1-7. 求由f(x) = x3 − 3x及g(x) = 3x所夾之區域面積。

定義函數f(x)及g(x) In[1]: f[x_]:=x^3-x;g[x_]:=3x; Out[1]:

定義函數f(x)及g(x)並繪出f(x)及g(x)夾之區域 In[2]: f[x_]:=x^3-x;g[x_]:=3x; In[3]: Show[Graphics[Plot[{x^3-x,3x},{x,-2.5,2.5}]], Graphics[Plot[{x^3-x,3x},{x,-2,2},Filling->{1->{2}}]], PlotRange->All] Out[3]: 10

5

-2

1

-1 -5

-10

2

CHAPTER 1. 微積分的基本操作

9

使用布爾函數計算面積 In[4]: Integrate[Boole[-2<x<=2]*Abs[f[x]-g[x]],{x,-2,2}] Out[4]: 8

若不使用布爾函數，本範例積分範圍必須拆成兩部分計算 In[5]: Integrate[f[x]-g[x],{x,-2,0}]+Integrate[g[x]-f[x],{x,0,2}] Out[5]: 8

∫∫ 範例 1-8.

24xydydx，其中Ω = {(x, y)|0 6 x 6 1, 0 6 y 6 1, 1 6 x + y 6 1.5}

Ω

首先畫出函數f(x, y) = 24xy在Ω的3D立體圖 In[1]: Plot3D[24*x*y,{x,0,1},{y,0,1},Filling->Bottom, RegionFunction->Function[{x,y},x>=0&&x<=1&&y>=0&&y<=1&&x+y>=1&&x+y<=1.5]] Out[1]:

10

1.0

5 0 0.0

0.5

0.5

1.0 0.0

配合布爾函數來計算積分值 In[2]: Integrate[24*x*y*Boole[x>=0&&x<=1&&y>=0&&y<=1&&x+y>=1&&x+y<=1.5], {x,0,1},{y,0,1}] Out[2]: 2.9375

1.4 積分

10

若不使用布爾函數，本範例積分範圍必須拆成兩部分計算 In[3]: Integrate[24*x*y,{x,0,0.5},{y,1-x,1}]+ Integrate[24*x*y,{x,0.5,1},{y,1-x,1.5-x}] Out[3]: 2.9375

範例 1-9. 假設f(x, y) = xyex

2 +y2

定義函數f(x, y) In[1]: f[x_,y_]:=x*y*Exp[x^2+y^2] Out[1]:

∫ 2 2 f(x, y)對x做積分，即求 xyex +y dx In[2]: Integrate[f[x,y],x] 1 2 2 Out[2]: ex +y y 2

∫b 2 2 f(x, y)對x做積分，其中x ∈ (a, b)，即求 xyex +y dx a

In[3]: Integrate[f[x,y],{x,a,b}] ) 1 2( 2 2 Out[3]: ey −ea + eb y 2

∫b ∫d 2 2 f(x, y)對x做積分，其中x ∈ (a, b)及y ∈ (c, d)，即求 xyex +y dydx a

c

In[4]: Integrate[f[x,y],{x,a,b},{y,c,d}] )( 2 ) 1 ( a2 2 2 Out[4]: e − eb ec − ed 4 x2 1 當然並不是所有的函數都是可積，例如標準常態分配的機率密度函數，f(z) = √ e− 2 , −∞ < 2π z < ∞就無法求得其不定積分。此時，Mathematica也就無能為力。所以在求解此類問題時就必須借助數

值積分。數值積分是解決定積分的另一種有效的方法，它可以求出一個近似解。特別是對於用Integrate指 令無法求出的定積分，數值積分更是可以發揮巨大作用。

CHAPTER 1. 微積分的基本操作

11

數值積分與傳統定積分基本概念是一樣的，但是傳統定積分使用對應的積分公式，對於一般人來說不 容易推導高複雜度的非線性函數定積分公式，藉由電腦或計算機來快速計算積分值。數值積分的概念源自 於黎曼積分，將積分範圍在一個極小的區域內分割成[x1 , x2 , x3 , . . . , xn ]，然後使用電腦或電子計算機程 式進行累加積分值。 z2 1 範例 1-10. 假設f(z) = √ e− 2 ，求解f(z)在區間−∞ < z < 1.96之積分值 2π

定義函數 In[1]: f[z_]:=1/Sqrt[2Pi]*Exp[-z^2/2]; Out[1]:

在−5與1.96間隨機選取49個點，即將區間劃分成50個子區間 In[2]: pts=RandomReal[{-5,1.96},49]//Sort; Out[2]:

將50的個矩形繪出，可以發現此圖形面積近似原來函數面積 In[3]: Show[Graphics[Table[{RGBColor[i/50,0,0],AspectRatio->0.65, Rectangle[{pts[[i]],0},{pts[[i+1]],f[pts[[i+1]]]}]},{i,1,48}]], Plot[f[x],{x,-5,5},PlotStyle->Thickness[0.015]]] Out[3]:

在數值分析上，常用的數值積分的方法有梯形法則和辛普森法則，這兩種方法都屬於牛頓-寇次公式。 牛頓-寇次公式的原理是 以n + 1點進行插值，求得對應函數f(x)的Lagrange多項式來近似原來的函 數，再進行定積分。所以，當n = 1時，即為梯形法則；n = 2時，為辛普森法則；n = 3時，為辛普 森3/8法則；n = 4時，為保爾法則。在Mathematica中，以n個點建構n − 1次方的Lagrange多項式指令

1.4 積分

為InterpolatingPolynomial，該指令用法為： InterpolatingPolynomial[{{x1 ,f(x1 )},{x2 ,f(x2 )},. . . ,{xn ,f(xn )}},x] 範例 1-11. 承上例，計算f(z)在區間−∞ < z < 1.96之積分值

建立牛頓-寇次公式 In[1]: NC[n_]:=NC[n]=Block[{temp=n,x,h,a,b,f,i}, x[0]=a;x[temp]=b;x[i]=a+i*h;h=(b-a)/temp; Integrate[InterpolatingPolynomial[Table[{x[i],f[x[i]]}, {i,0,temp}],t],{t,a,b}]] Out[1]:

輸出牛頓-寇次公式的前五近似公式 In[2]: TableForm[NC/@Range[5]] Out[2]:

1 − (a − b)(f [a] + f [b]) 2 ( [ ]) 1 a+b − (a − b) f [a] + f [b] + 4f 6 ( ( [2 ] [ ])) 1 1 1 − (a − b) f [a] + f [b] + 3 f (2a + b) + f (a + 2b) 8 3[ ( (3 ( ) ] [ ])) 1 a+b 1 1 − (a − b) 7f [a] + 7f [b] + 4 3f + 8f (3a + b) + 8f (a + 3b) 90 2 4 4

將−5與1.96均分成1000等分 In[3]: pts=Table[-5+(1.96+5)/1000*i,{i,0,1000}]; Out[3]:

以梯形法則計算面積 In[4]: Sum[(a=pts[[i]];b=pts[[i+1]];NC[1]),{i,1,1000,1}] Out[4]: 0.975001

以辛普森法則計算面積 In[5]: Sum[(a=pts[[i]];b=pts[[i+2]];NC[2]),{i,1,999,2}] Out[5]: 0.975002

12

CHAPTER 1. 微積分的基本操作

13

分別以梯形法則、辛普森法則及保爾法則計算面積 In[6]: Sum[(a=pts[[i]];b=pts[[i+#]];NC[#]),{i,1,1001-#,#}]&/@{1,2,4} Out[6]: {0.975001, 0.975002, 0.975002}

上述三種計算法則即為牛頓-寇次公式的1,2,4次方近似式，由上式可發現並非次方越高，近似效果越 好。牛頓-寇次公式對於次數較高的多項式而有很大誤差，此即為龍格現象。

分別將−5與1.96均分成10, 20, . . . ,1000等分在以梯形法計算面積。由輸出圖形發以發現分割越 多，越趨近0.975 In[7]: ListPlot[(pts=Table[-5+(1.96+5)/#*i,{i,0,#}]; Sum[(a=pts[[i]];b=pts[[i+1]];-(1/2)(a-b)(f[a]+f[b])),{i,1,#,1}])& /@Table[10*i,{i,1,100}]] Out[7]: 0.975000 0.974995 0.974990 0.974985 0.974980 0.974975 20

40

60

80

100

上述的計算方式，雖然可以達到近似的效果，但效率並不高。為了提高收斂速度，減少計算次數，數 學家也開始尋求其他方法。我們在此介紹一個較常用的計算方式，Romberg積分法 。Romberg積分法是 在以牛頓-寇次公式為基礎所構造出的一種加速計算積分的方法，它在不增加計算次數的前提下提高了準 確度。有關Romberg演算法的計算公式如下：  1   R (0, 0) = (b − a) (f (a) + f (b))   2   ) 2n−1 (  b−a ∑ b−a 1 f a + (2k − 1) n R(n, 0) = R (n − 1, 0) + n  2 2 2   k=1   1   R (n, m) = R (n, m − 1) + [R (n, m − 1) − R (n − 1, m − 1)] 4m − 1 根據Romberg演算法，我們可以將上述公式以迴圈的方式撰寫程式如下：

1.4 積分

14

∫ 1.96 x2 1 √ e− 2 d x 利用Romberg演算法求解 2π −5 In[8]: Clear[R,n,m];f[x_]:=1/Sqrt[2Pi]*Exp[-x^2/2];a=-5;b=1.96;n=1;error=100; R[0,0]=((b-a)/2*(f[a]+f[b]))//N; While[error>10^-10, For[m=1,m<=n, R[n,0]=(R[n-1,0]/2+(b-a)/2^n*Sum[f[a+(2k-1)*(b-a)/2^n], {k,1,2^(n-1)}])//N; R[n,m]=(R[n,m-1]+(b-a)/(4^m-1)*(R[n,m-1]-R[n-1,m-1]))//N; error=Abs[R[n,m]-R[n,m-1]]; m++]; n++]; Out[8]: 0.975001

將計算過程列出 In[8]: Table[R[n,m],{n,0,5},{m,0,5}]//TableForm Out[8]: 0.20338 0.539003 0.950481 0.967869 0.973201 0.974551

0.650877 1.08764 1.11676 0.973665 0.966066 0.963674 0.974978 0.975065 0.975208 0.975001 0.975002 0.975001

0.975253 0.975 0.975

當然Mathematica中作數值積分不需要自己撰寫程式，Mathematica有自己的數值積分指令。 在Mathematica中數值積分的指令為NIntegrate，該指令組合了N以及Integrate的用途，該指令的用法有 以下幾種： 1. NIntegrate[f, {x, a, b}]： 傳回函數f的數值積分 2. NIntegrate[f, {x1 , a, b}, {x2 , c, d}, . . . ] 傳回函數f在區間(a, b) × (c, d) × · · · 的定積分值。要 注意的是x1 為最外面的積分變數，x2 為倒數第二個的積分變數

z2 1 範例 1-12. 假設f(z) = √ e− 2 ，求解f(z)在區間−∞ < z < 1.96之積分值，即求 2π

∫ 1.96 −∞

z2 1 √ e − 2 dz 2π

CHAPTER 1. 微積分的基本操作

15

繪出常態分配在−∞ < z < 1.96的圖形 In[1]: p1=Plot[1/(2*Pi)^0.5*Exp[-z^2/2],{z,-5,5},RegionFunction->(#<=1.96&), Filling->Bottom,PlotStyle->Thickness[0.01]]; p2=Plot[1/(2*Pi)^0.5*Exp[-z^2/2],{z,-5,5}, PlotStyle->Thickness[0.01],Filling->y]; Show[p1,p2] Out[1]:

0.4

0.3

0.2

0.1

-4

-2

2

4

以Integrate計算不定積分值 In[1]: Integrate[1/(2*Pi)^0.5*Exp[-z^2/2],z] [ ] z Out[1]: 0.5Erf √ 2

上式中Erf(x)為誤差函數， 2 Erf(x) = √ π

∫x

2

e−t dt 0

由於函數f(z)不可積，所以Mathematica並沒有真正傳回不定積分。若在上式中將x以數值取代，則可傳 回其近似值。

以Integrate計算不定積分值 ，在統計學上此數表示為Z0.05 In[2]: Integrate[1/(2*Pi)^0.5*Exp[-z^2/2],{z,-Infinity,1.96}] Out[2]: 0.975002

1.4 積分

16

我們也可以用NIntegrate來求解相同的問題。要注意的是Nintegrate指令中的積分範圍參數必須為一 數值，否則會無法運算並出現錯誤訊息。

以NIntegrate計算，由於x並沒有指定特定數值，所以無法計算 In[3]: NIntegrate[1/(2*Pi)^0.5*Exp[-z^2/2],{z, -Infinity, x}] NIntegrate::nintp : Encountered the non-number x at z = z. More. . . [ ] z2 e− 2 Out[3]: Nintegrate , {z, −∞, x} (2π)0.5

將x以1.96取代則可求得積分值為0.975002 In[4]: NIntegrate[1/(2*Pi)^0.5*Exp[-z^2/2],{z, -Infinity, 1.96}] Out[4]: 0.975002

若在模式建立的過程中需要將NIntegrate中將積分參數以符號表示，則必須以函數方式表示。否 則，Mathematica會在後續的計算中以Integrate來取代NIntegrate。

定義函數Φ(x)，其中?Number表示指定x為數值型態的變數 In[5]: Φ[x_?NumberQ]:=NIntegrate[1/(2*Pi)^0.5*Exp[-z^2/2], {z,-Infinity,x}]; Out[6]:

計算Φ(z), z = −3, −2, . . . , 3, 的常態機率值 In[6]: Table["P(Z<"<>ToString[i]<>")="<>ToString[Φ[i]], {i,-3,3,1}]//TableForm Out[6]: P(Z < −3) = 0.0013499 P(Z < −2) = 0.0227501 P(Z < −1) = 0.158655 P(Z < 0) = 0.5 P(Z < 1) = 0.841345 P(Z < 2) = 0.97725 P(Z < 3) = 0.99865

CHAPTER 1. 微積分的基本操作

17

利用Φ(x)產生1000個標準常態分配的隨機變數 In[7]: Table[x/.FindRoot[Φ[x]==Random[],{x,1}],{1000}]; Out[8]:

引入直方圖套件，並將上式得到的隨機變數以直方圖繪出 In[9]: Needs["Histograms`"];Histogram[%] Out[9]: 80

60

40

20

-3

-2

0

-1

1

2

此處所使用的FindRoot為方程式求根的指令，使用方式留待後面章節說明。

範例 1-13. 假設f(x, y) = e−

x2 +y2 2

，求解f(x, y)在區間0 < x < y < 1之雙重積分值，即求 ∫1 ∫y e− 0

0

x2 +y2 2

d xd y

以Integrate計算雙重積分值 In[1]: Integrate[Exp[-(x^2+y^2)/2],{y,0,1},{x,0,y}] [ ]2 1 1 √ Out[1]: πErf 4 2

將上式以數值方式輸出 In[2]: %//N Out[2]: 0.366047

1.5 級數

18

以NIntegrate計算雙重積分值 In[3]: NIntegrate[Exp[-(x^2+y^2)/2],{y,0,1},{x,0,y}] Out[3]: 0.366047

搭配布爾函數，就可以不必理會積分順序 In[4]: NIntegrate[Exp[-(x^2+y^2)/2]*Boole[0<=x&&x<=y&&y<=1],{y,0,1},{x,0,1}] Out[4]: 0.366047

1.5 級數 在Mathematica中計算函數的泰勒級數指令為Series[f,{x, a, n}]，表示函數f在x = a的n階的泰勒 展開式。同時，Series也可處理多變數函數的泰勒展開式。當a = 0時，泰勒級數也稱為麥克勞林級數。 該指令的用法有： 1. Series[f,{x, a, n}]：函數f在x = a的n階泰勒展開式，所傳回的運算式中O (xn )稱為第n項餘 式。 2. Series[f,{x1 , a1 , n1 }, {x2 , a2 , n2 }]： 函 數f先 對x1 在x1 = a1 做n1 階 展 開 後 再 對x2 在 x2 = a2 展開n2 階的泰勒展開式。 與級數較相關的指令有SeriesCoefficient與CoefficientList，其用法簡單介紹如下： 1. SeriesCoefficient[f[x],n]：傳回多項式函數f(x)中xn 項的係數。 2. CoefficientList[f[x],x]：傳回多項式函數f(x)中的係數。

CHAPTER 1. 微積分的基本操作

19

範例 1-14.

計算f(x)在x = a的5階泰勒展開式 In[1]: Series[f[x], {x, a, 5}] 1 1 1 1 (5) Out[1]: f(a) + f ′ (a)(x − a) + f ′′ (a)(x − a)2 + f(3) (a)(x − a)3 + f(4) (a)(x − a)4 + f (a)(x − a)5 + 2 6 24 120 ) ( O (x − a)6

計算f(x) = Sin(x)在x = 0的5階泰勒展開式 In[2]: Series[Sin[x],{x, 0, 10}] ( ) x3 x5 x7 x9 Out[2]: x − + − + + O x11 6 120 5040 362880

將上式中的第5項餘式省略 In[3]: Normal[%] x3 x5 x7 x9 Out[3]: x − + − + 6 120 5040 362880

將Sin(x)及上式同時繪出 In[4]: Plot[{Sin[x],%},{x,0,2Pi}] Out[4]: 3

2

1

1

2

3

4

5

6

-1

分別將Sin(x)在x = 0的5階泰勒展開式中的各項係數取出 In[5]: Table[SeriesCoefficient[%3, k], {k, 1, 10}] } { 1 1 1 1 , 0, − , 0, ,0 Out[5]: 1, 0, − , 0, 6 120 5040 362880

1.5 級數

20

覺得以Table方式取出太麻煩，用CoefficientList更方便 In[6]: CoefficientList[%3,x] { } 1 1 1 1 Out[6]: 1, 0, − , 0, , 0, − , 0, ,0 6 120 5040 362880

在泰勒展開式中我們是以多項式函數來逼近原式。因此，隨著n的增加，泰勒展開式就可以更逼近原 式。 在以下範例我們將以不同的n來近似Sin(x)做說明。 範例 1-15.

分別產生Sin(x)的2, 4, 6, 8, 10, . . . , 20階泰勒展開式 In[1]: y=Table[Normal[Series[Sin[x],{x,0,i}]],{i,2,20,2}]; Out[1]:

將上式產生的各展開式與Sin(x)同時繪出 In[2]: Plot[Evaluate[{y,Sin[x]}], {x, 0, 2Pi},PlotRange->{-1.5,1.5}, PlotStyle->Table[Hue[k], {k, 0, 1, 0.05}], Ticks->{Flatten[Table[{i},{i,0,2Pi,Pi/4}]],Automatic}] Out[2]:

1.5 1.0 0.5

Π

Π

3Π

4

2

4

Π

5Π

3Π

7Π

4

2

4

2Π

-0.5

-1.0

-1.5

由上圖可以發現，當n = 20時，泰勒展開式與Sin(x)在x ∈ (0, 2π)之間幾乎無差別。以上我們所探 討的是單變數下Mathematica中計算級數的方法；接下來我們介紹函數Series在多變數的情況下的應用。

CHAPTER 1. 微積分的基本操作

21

範例 1-16.

先將f(x, y) = Sin(x + y)以x = 0求3階泰勒展開式，接著在以y = 0展開3階泰勒展開式 In[1]: Series[Sin[x+y],{x,0,3},{y,0,3}] ( ) ( ) ( ) ( 4) ( 4) ( 4) y3 y2 y y3 2 Out[1]: y − +O y +x 1− +O y +x − + +O y 2 2 12 (6 ) 2 ( ) ( ) 1 y + O y4 + O x4 +x3 − + 6 12

我們必須注意上式並非f(x, y)以(x, y) = (0, 0)展開的泰勒展開式。而是先將f(x, y)以x = 0求3階展 開式，之後再對x次方項的係數y求3階的泰勒展開式，說明如下：

先將f(x, y) = Sin(x + y)以x = 0求3階泰勒展開式 In[2]: Series[Sin[x+y],{x,0,3}}] ( ) 1 1 Out[2]: sin(y) + x cos(y) − x2 sin(y) − x3 cos(y) + O x4 2 6

將上式的係數項取出 In[3]: mycoff=CoefficientList[Series[Sin[x+y],{x,0,3}],x] { } sin(y) cos(y) Out[3]: sin(y), cos(y), − ,− 2 6

分別將各係數項對y取3階泰勒展開式 In[4]: Series[mycoff[[#]],{y,0,3}]&/@{1,2,3,4} { } ( ) ( ) y y3 ( ) 1 y2 ( ) y3 y2 Out[4]: y − + O y4 , 1 − + O y4 , − + + O y4 , − + + O y4 6 2 2 12 6 12

我們在此可以發現，將%4的結果代回%2中可得到與%1相同的結果。所以，Mathematica中函 數Series在多變數中所求得的展開式並非泰勒展開式。若要求得Sin(x + y)之3階泰勒展開式語法可參考 範例3，語法如下：

建立一個以(x, y) = (a, b)展開n階泰勒展開式的公式 In[5]: mytaylor[a_,b_,k_]:=Plus@@(Table[Binomial[n,i]* (D[Sin[x+y],{x,i},{y,n-i}]/.{x->a,y->b})*(x-a)^i*(y-b)^(n-i)/n!, {n,0,k},{i,0,n}]//Flatten); Out[5]:

1.5 級數

傳回(x, y) = (0, 0)展開3階泰勒展開式 In[6]: mytaylor[0,0,3] x3 x2 y xy2 y3 Out[6]: − − − +x− +y 6 2 2 6

22

第2章 非線性最佳化

2.1 方程式求解 在Mathematica中求解方程式根的方式有很多。當方程式具有封閉解時，我們可以用Solve這個指令 來求解。要是函數的解無函數形式，只能藉助數值方法求其數值解，這類問題所得是精確的數值解。在此 分別將Mathematica求解方程式和聯立方程組的指令介紹如下。

2.1.1 Solve 在Mathematica中函數求解方程式或聯立方程組最基本的方法就是使用指令Solve，Solve的使用方法 有以下幾種：

1. Solve[lhs == rhs, x]：以x為變數求解方程式lhs == rhs 2. Solve[{eqn1, eqn2, . . . }, {x1 ,x2 . . . }]：以x1 , x2 , . . . 為變數求解聯立方程組

範例 2-1. 求解ax2 + bx + c = 0

以Solve求解 In[1]: Solve[a*x^2+b*x+c==0,x] {{ } { }} √ √ −b − b2 − 4ac −b + b2 − 4ac Out[1]: x→ , x→ 2a 2a

2.1 方程式求解

24

求解6x2 + 13x + 5 = 0 In[2]: Solve[#1*x^2+#2*x+#3==0,x]&@@{6,13,5} {{ } { }} 5 1 Out[2]: x→− , x→− 3 2

範例 2-2. 求解x使得x5 − 4x + 2 = 0

繪出函數圖，由圖形可知方程式有三組實數解 In[1]: Plot[x^5-4*x+2,{x,-2,2}] Out[1]: 10

5

-2

1

-1

2

-5

以Solve求解，由於並無解析解，因此Mathematica在此以Root表示之 In[2]: Solve[x^5-4*x+2==0,x] Out[2]:

{{

[ ]} { [ ]} x → Root 2 − 4#1 + #15 &, 1 , x → Root 2 − 4#1 + #15 &, 2 , { [ ]} { [ ]} x → Root 2 − 4#1 + #15 &, 3 , x → Root 2 − 4#1 + #15 &, 4 , { [ ]}} x → Root 2 − 4#1 + #15 &, 5

若要求得數值解，可以搭配N要求輸出數值解 In[3]: sol=x/.Solve[x^5-4*x+2==0,x]//N Out[3]: {−1.51851, 0.508499, 1.2436, −0.116792 − 1.43845i, −0.116792 + 1.43845i}

在Mathematica中Root指令為傳回多項式f(x) = 0的第k個根，其用法為Root[f,k]。所以要傳回%1的 第二個根有以下兩種方式。

CHAPTER 2. 非線性最佳化

25

以Solve方式傳回第2個根 In[4]: Solve[x^5-4*x+2==0,x][[2]]//N Out[4]: {x → 0.508499}

以Root方式傳回x5 − 4x + 2 = 0所有根 In[5]: (Root[x^5-4*x+2,x,#]//N)&/@{1,2,3,4,5} Out[5]: {−1.51851, 0.508499, 1.2436, −0.116792 − 1.43845i, −0.116792 + 1.43845i}

以Cases傳回實數根 In[6]: Cases[sol,_Real] Out[6]: {−1.51851, 0.508499, 1.2436}

以Select傳回實數根 In[7]: Select[sol,#∈Reals&] Out[7]: {−1.51851, 0.508499, 1.2436}

繪出函數與x軸交點圖形 In[8]: Plot[x^5-4x+2,{x,-2,2},Epilog->{Red,PointSize[Large], Point[{#,0}]&/@Cases[sol,_Real]}] Out[8]: 10

5

-2

1

-1

-5

2

2.1 方程式求解

26

  x 2 + y2 = 1 範例 2-3. 解聯立方程組  |x| + |y| = 1

以Solve求解，可得四組解依序為 In[1]: sol={x,y}/.Solve[{Abs[x]+Abs[y]==1,x^2+y^2==1},{x,y}] Out[1]: {{−1, 0}, {0, −1}, {0, 1}, {1, 0}}

以Select傳回x > 0及y > 0的解 In[2]: Select[sol,#[[1]]>=0&&#[[2]]>=0&] Out[2]: {{0, 1}, {1, 0}}

繪出聯立方程組交點 In[3]: ContourPlot[{Abs[x]+Abs[y]==1,x^2+y^2==1},{x,-1,1},{y,-1,1}, Epilog->{Red,PointSize[Large],Point[sol]}] Out[3]: 1.0

0.5

0.0

-0.5

-1.0 -1.0

-0.5

0.0

0.5

1.0

2.1.2 FindRoot 在Mathematica中函數求數值解的指令為FindRoot，該指令的用法有以下幾種： 1. FindRoot[lhs == rhs, {x, a}]：以x = a為起始值尋找左右兩方程式交點的數值解。

CHAPTER 2. 非線性最佳化

27

2. FindRoot[f, {x, a}]：以x = a為起始值尋找方程式f的數值解。必須注意的是，該指令預設為求 解f = 0，若要求解其他數值可以f == c定義。 3. FindRoot[{f1 , f2 , . . . }, {{x, a1 }, {y, a2 }, . . . }]：以(x, y, . . . )為起始值尋找聯立方程 組(f1 , f2 , . . . )的數值解。 FindRoot指令中預設以Newton Method做為預設的搜尋方法，假設聯立方程組 fi (X) = 0, i = 1, 2, . . . , m 在給定Xk 時，由泰勒展開式可知 ( ) ( )( ) fi (X) ≈ fi Xk + ∇fi Xk X − Xk , i = 1, 2, . . . , m 由於fi (X) = 0, i = 1, 2, . . . , m，故可知 ( ) ( )( ) fi Xk + ∇fi Xk X − Xk = 0 經移項整理後，可求得

[ ( )]−1 ( k ) X = Xk − ∇fi Xk fi X

所以可求得單變量與多變數Newton Metom的迭代計算過程間遞迴關係式如下： xk+1 = xk − 及

f(xk ) f ′ (xk )

[ ( )]−1 ( k ) Xk+1 = Xk − ∇fi Xk fi X

[ ( )] 其中 ∇fi Xk 為Jacobian。為避免迭代過程發散，建議在求解之前先繪製函數圖形以利設定起始值； 起始值可以單點設定，也可以一個區間，(x, min, max)，作為其設定值。 範例 2-4. 求ex − x2 = 0

定義函數 In[1]: f[x_]:=Exp[x]-x^2; Out[1]:

2.1 方程式求解

28

繪製圖形以方便設定起始值 In[2]: Plot[f[x], {x, -2, 2}]; Out[2]: 3 2 1

-2

-1

1

2

-1 -2 -3 -4

以While自行編寫程式 In[2]: init=2;x=init; err=10;sol={}; While[err>10^(-6),x=N[x-f[x]/f’[x]];err=Abs[f[x]];AppendTo[sol,x]]; sol Out[2]: {2.21298,1.31291,-0.513828,-0.71937,-0.703567,-0.703467}

上式所使用的While為一般程式所使用的控制流程，在使用方式與一般程式無太大差別，但上式程式 的寫法無法發揮Mathematica真正的功能。上式中我們可以發現x=N[x-f[x]/f’[x]]為一個遞迴函數， 因此，我們在這邊可以NestList來將程式寫得更簡潔。

以NestList傳回逐次迭代值 In[2]: NestList[N[#-f[#]/f’[#]]&,2,6] Out[2]: {3,2.21298,1.31291,-0.513828,-0.71937,-0.703567,-0.703467}

若想直接求得方程式ex − x2 = 0的根，那我們亦可以NestWhile來求得 In[3]: init=2;f[x_]:=Exp[x]-x^2;g[x_]:=N[x-f[x]/f’[x]]; NestWile[g,init,Unequal,All] Out[3]: -0.703467

CHAPTER 2. 非線性最佳化

29

根據Newton Method前六次迭代的傳回數值繪製收斂過程 In[3]: myplot=Plot[{f’[#]*(x-#)+f[#],f[x]},{x,-3,5},PlotRange->{-5,15}, PlotStyle->{Thickness[0.003],Thickness[0.01]}, Epilog->{PointSize[0.02],Black,Point[{#,N[f[#]]}], Red,Point[{#,0}]}]&/@NestList[N[#-f[#]/f’[#]]&,3,5];

將myplot所得到的六個圖形以GraphicsGrid繪製成2 × 3的矩陣圖 In[4]: GraphicsGrid[{myplot[[{1, 2, 3}]], myplot[[{4, 5, 6}]]}] Out[4]: 15

15

15

10

10

10

5

5

5

2

-2

4

2

-2

4

-2

-5

-5

-5

15

15

15

10

10

10

5

5

5

2

-2

4

-5

2

-2 -5

4

-2

2

4

2

4

-5

利用FindRoot進行求解 In[4]: FindRoot[f[x],{x,3}] Out[4]: {x → −0.703467}

在使用函數FindRoot時，我們也可藉由指定參數EvaluationMonitor來要求將收斂過程輸出至指定字 串。

定義一個sol的一維陣列，並利用AppendTo這個指令將迭代過程加入sol中。 In[5]: sol={};FindRoot[f[x],{x,2},EvaluationMonitor:>AppendTo[sol, x]];sol Out[5]: {3.,2.21298,1.31291,-0.51383,-0.71937,-0.70357,-0.70347,-0.70347}

2.1 方程式求解

30

Mathematica也有提供一個繪製FindRoot迭代步驟圖的指令FindRootPlot。FindRootPlot所傳回資 料陣列包含三個元素，依序為收斂值，求解資訊 ，收斂步驟圖。

引入最佳化套件繪製求解步驟圖 In[6]: <
傳回收斂步驟圖 In[7]: FindRootPlot[f[x], {x, 3}, PlotStyle->{Thickness[.01],RGBColor[0.4,0,0.4]}][[3]] Out[7]: 10

8

6

4

2

-0.5

0.5

1.0

1.5

2.0

2.5

3.0

由圖形可知，收斂過程為x0 = 3 ⇒ x1 = 2.21298 ⇒ · · · ⇒ x7 = −0.703467

傳回求解資訊 In[9]: FindRootPlot[f[x], {x, 2}, PlotStyle->{Thickness[.01],RGBColor[0.4,0,0.4]}][[{1,2}]] Out[9]: {{x → −0.703467}, {Steps → 8, Residual → 8, Jacobian → 7}}

由輸出資料可知經由8次迭代計算後，x = −0.703467

CHAPTER 2. 非線性最佳化

31

範例 2-5. 求解聯立方程組   x2 + y2 = 1  y − xex = 0

首先定義函數f(x, y)和g(x, y) In[1]: f[x_,y_]:=x^2+y^2-1;g[x_,y_]:=y-x*Exp[x]; Out[1]: 由圖形可了解交點有兩處 In[2]: ContourPlot[{f[x,y],g[x,y]},{x,-1.5,1.5},{y,-1.5,1.5}] Out[2]: 1.5

1.0

0.5

0.0

-0.5

-1.0

-1.5 -1.5

-1.0

-0.5

0.0

0.5

1.0

1.5

首先計算Jacobian In[3]: J=D[{f[x,y],g[x,y]},{{x,y}}] Out[3]: {{2x, 2y}, {−ex − ex x, 1}}

[ ( )]−1 ( k ) 定義 ∇fi Xk fi X ，其中"."在Mathematica中為矩陣相乘符號 In[4]: delta[x_,y_]=Inverse[J].{f[x,y],g[x,y]} { ( ( )) } −1 + x2 + 2ex xy − y2 2xy + ex −1 + y2 + x −1 + (−1 + x)x + y2 , Out[4]: 2 (x + ex (1 + x)y) 2 (x + ex (1 + x)y)

2.1 方程式求解

32

分別以(1, 0.5)和(−1, 0.5)為起始值利用NestWhile傳回數值解 In[5]: NestWhile[N[#-Apply[delta,#]]&,#,Unequal,All]&/@{{1,0.5},{-1,0.5}} Out[5]: {{0.513489,0.858096},{-0.930244,-0.366942}}

上式所求得的解之所以有差異，主要原因為起始值所導致。

繪出最佳解 In[6]: ContourPlot[{f[x,y],g[x,y]},{x,-1.5,1.5},{y,-1.5,1.5}, Epilog->{Red,PointSize[Large], Point[{x,y}/.FindRoot[{f[x,y],g[x,y]},{x,#},{y,#}]&/@{1,-1}]}] Out[6]: 1.5

1.0

0.5

0.0

-0.5

-1.0

-1.5 -1.5

-1.0

-0.5

0.0

0.5

1.0

1.5

定義一個sol的二維陣列，並利用AppendTo將以(1, 0.5)為起始值的迭代過程加入sol中。 In[7]: sol={};FindRoot[{f[x,y],g[x,y]},{x,1},{y,0.5}, EvaluationMonitor:>AppendTo[sol,{Length[sol],x,y}]];sol//TableForm Out[7]: 0 1 2 3 4 5 6 7

0.1 0.5 0.897833 1.080430 0.291842 0.639567 0.586557 0.900495 0.518046 0.859091 0.513507 0.858098 0.513489 0.858096 0.513489 0.858096

CHAPTER 2. 非線性最佳化

33

分別以(1, 0.5)和(−1, 0.5)為起始值利用FindRoot求解聯立方程組 In[8]: FindRoot[{f[x,y],g[x,y]},{x,#[[1]]},{y,#[[2]]}]&/@{{1,0.5},{-1,0.5}} Out[8]: {{x → 0.513489, y → 0.858096}, {x → −0.930244, y → −0.366942}}

引入最佳化套件繪製求解步驟圖 In[9]: <
繪製求解的步驟圖 In[10]: FindRootPlot[{f[x,y],g[x,y]},{{x,1},{y,0.5}}][[3]] Out[10]: 0.9

0.8

0.7

0.6

0.5 0.6

0.7

0.8

0.9

1

傳回求解資訊，由輸出資料可知經由7次迭代計算後，(x, y) = (0.513489, 0.858096) In[11]: FindRootPlot[{f[x,y],g[x,y]},{{x,1},{y,0.5}}][[{1,2}]] Out[11]: {{x->0.513489,y->0.858096},{Steps->7,Residual->7,Jacobian>6}}

2.2 非線性最佳化

34

2.2 非線性最佳化 在微積分中，我們知道目標函數f (X)若存在極值，則極值必發生在∇f (X) = 0或不存在的點。 當∇f (X) = 0有解析解，我們可以用Solve來求得臨界點；反之，若∇f (X) = 0只有數值解，則可 用FindRoot來解聯立方程組以求得臨界點。 範例 2-6. min f(x1 , x2 ) = 100(x2 − x21 )2 + (1 − x1 )2 x1 ,x2

定義函數f(x1 , x2 ) In[1]: f[x1_,x2_]:=100*(x2-x1^2)^2+(1-x1)^2

繪出函數3D立體圖 In[2]: Plot3D[f[x1,x2],{x1,-2,2},{x2,-2,2} Out[2]:

2000 1500

2

1000 1

500 0 -2

0 -1 0

-1 1 2 -2

計算梯度 In[3]: g=D[f[x1,x2],{{x1,x2}}] { ( ) ( )} Out[3]: −2(1 − x1) − 400x1 −x12 + x2 , 200 −x12 + x2

計算臨界點 In[4]: {x1,x2}/.FindRoot[D[f[x1,x2],{{x1,x2}}],{x1,1},{x2,1}] Out[4]: {1., 1.}

CHAPTER 2. 非線性最佳化

35

計算海賽矩陣 In[5]: h=D[f[x1,x2],{{x1,x2}},{{x1,x2}}] {{ ( ) } } Out[5]: 2 + 800x12 − 400 −x12 + x2 , −400x1 , {−400x1, 200}

計算海賽矩陣各階行列值 In[6]: Det[h[[Table[i,{i,1,#}],Table[i,{i,1,#}]]]]&/@{1,2}/. FindRoot[D[f[x1,x2],{{x1,x2}}],{x1,1},{x2,1}] Out[6]: {802., 400.}

上式中應用到Map這個映射指令，在此也可利用Table達到相同的目的 In[7]: Table[Det[h[[Table[i,{i,1,k}],Table[i,{i,1,k}]]]],{k,2}]/. FindRoot[D[f[x1,x2],{{x1,x2}}],{x1,1},{x2,1}] Out[7]: {802., 400.}

由上式傳回之海賽矩陣行列值可知f(x1 , x2 )在{x1, x2} = {1., 1.} 有極小值，極小值可計算如下： In[8]: f[x1,x2]/.FindRoot[D[f[x1,x2],{{x1,x2}}],{x1,1},{x2,1}] Out[8]: 0.

範例 2-7.

min 2x1 x2 x3 − 4x1 x3 − 2x2 x3 + x21 + x22 + x23 − 2x1 − 4x2 + 4x3

x1 ,x2 ,x3

定義函數f(x1 , x2 , x3 ) In[1]: f[x1_,x2_,x3_]:=2x1*x2*x3-4x1*x3-2x2*x3+x1^2+x2^2+x3^2-2x1-4x2+4x3;

由於臨界點有五個，若採用FindRoot無法將此五組臨界點一次求出，在此可採用Solve In[1]: {x1,x2,x3}/.Solve[D[f[x1,x2,x3],{{x1,x2,x3}}]==0,{x1,x2,x3}] Out[1]: {{0, 1, −1}, {0, 3, 1}, {1, 2, 0}, {2, 1, 1}, {2, 3, −1}}

2.2 非線性最佳化

36

計算海賽矩陣 In[2]: h=D[f[x1,x2,x3],{{x1,x2,x3}},{{x1,x2,x3}}]/. Solve[D[f[x1,x2,x3],{{x1,x2,x3}}]==0,{x1,x2,x3}];

計算海賽矩陣各階行列值。由海賽矩陣可知，上述解只有第三組為極小值，其餘皆為鞍點。 In[2]: Table[Det[h[[m,Table[i,{i,1,#}],Table[i,{i,1,#}]]]]&/@{1,2,3},{m,1,5}] //TableForm Out[2]: 2 2 2 2 2

0 0 4 0 0

−32 −32 8 −32 −32

計算極小值 In[3]: f[x1, x2, x3]/.Solve[D[f[x1,x2,x3],{{x1,x2,x3}}]==0,{x1,x2,x3}][[3]] Out[3]: −5

2.2.1 無限制條件下最佳化求解 然而當目標函數是包含多變數時，求解聯立方程組∇f (X) = 0就更為複雜。所以我們需要一些有效率 的演算法來幫助求解。假設函數f (X)，由泰勒展開式可知 1 f (X + △X) ≈ f (X) + (∇f (X))⊤ △X + △X⊤ H △ X 2 欲求f (X)極值，必要條件為 ∂f (X + △X) ≈ ∇f (X) +H △ X = 0 ∂△X 由於△X = Xk+1 − Xk ，Xk 為第k階段的迭代值 ( ) Xk+1 − Xk = △X = −H−1 ∇f Xk 所以第k + 1階段的迭代值為

( ) Xk+1 = Xk − 1 × H−1 ∇f Xk

CHAPTER 2. 非線性最佳化

37

以上的搜尋方式稱為牛頓法（Newton’s Method）。在這個推導中步長1，但實際演算法進行時容易造 成發散。為避免發生上述情況，我們以α取代1並求得使得f (X)有極大或極小的α值作為下一階段的迭代 值，即

( ) Xk+1 = Xk −αH−1 ∇f Xk

這種作法稱為修正牛頓法（Modified Newton’s Method）。 Mathematica提 供 了 一 些 求 解 最 佳 化 問 題 的 指 令 ， 例 如FindMinimum、FindMaximum 、Minimize、Maximize、Minimize和 NMaximize來求解有限制情況和無限制情況下的極大極小化問 題。以上指令使用方式介紹如下： 1. FindMinimum[f(x),{x,x0 }]：以x = x0 為起始值求解函數f(x)的區域極小值 2. FindMinimum[f(x1 , x2 , . . . ),{x1 ,x10 },{x1 ,x10 },

. . . ]： 以(x10 , x20 , . . . )為 起 始 值 求 解

f(x1 , x2 , . . . )的區域極小值 3. FindMinimum[{f(x), 限制式},x]：求解f(x)在限制下的區域極小值，其中限制式可以是等式、不等 式以及邏輯條件 4. FindMinimum[{f(x1 , x2 , . . . ), 限制式},{x1 ,x1 , . . . }]：求解f(x1 , x2 , . . . )在限制下的區域極小值， 其中限制式可以是等式、不等式以及邏輯條件 5. Minimize[f(x1 , x2 , . . . ),{x1 ,x1 , . . . }]：求解f(x1 , x2 , . . . )的全域極小值 6. Minimize[{f(x1 , x2 , . . . ), 限制式},{x1 ,x1 , . . . }]：求解f(x1 , x2 , . . . )在限制下的全域極小值，其中 限制式可以是等式、不等式以及邏輯條件 7. NMinimize[f(x1 , x2 , . . . ),{x1 ,x1 , . . . }]：求解f(x1 , x2 , . . . )的全域極小值 8. NMinimize[{f(x1 , x2 , . . . ), 限制式},{x1 ,x1 , . . . }]：求解f(x1 , x2 , . . . )在限制下的全域極小值，其 中限制式可以是等式、不等式以及邏輯條件 9. NMinimize[f(x1 , x2 , . . . ),{x1 ,x1 , . . . }]：求解f(x1 , x2 , . . . )的全域極小值 10. NMinimize[{f(x1 , x2 , . . . ), 限制式},{x1 ,x1 , . . . }]：求解f(x1 , x2 , . . . )在限制下的全域極小值，其 中限制式可以是等式、不等式以及邏輯條件

2.2 非線性最佳化

38

Mathematica在進行FindMinimum時，可以指定以下幾種搜尋方式進行求解：Newton, QuasiNewton, LevenbergMarquardt, ConjugateGradient 及 PrincipalAxis。一般Mathematica以準牛頓法（QuasiNewton Method）進行迭代求解，所以收斂的速度較牛頓法為慢。然而牛頓法需要計算海賽矩陣的 反矩陣，如果海賽矩陣是奇異矩陣（singular matrix）或狀況不佳的矩陣（ill-conditioned matrix）， 求反矩陣時會產生數值問題，應用上應該加以注意。若不指定迭代方式，當目標函數含有平方和 時，Mathematica則改以LevenbergMarquardt法進行迭代程序。當變數的起始值給定兩個時，則改以使 用主軸法。在這邊必須提醒的是，FindMinimum所求得的是區域極值，並不保證為全域極值。 若 以Minimize進 行 求 解 ， 則 必 須 是 有 解 析 解 的 前 提 下 才 可 使 用 ， 否 則 必 須 改 以NMinimize求 解。若以NMinimize進行求解，則可以指定以下幾種搜尋方式：NelderMead, DifferentialEvolution, SimulatedAnnealing, RandomSearch。在這邊必須提醒的是，Minimize或NMinimize所求得的通常是全 域極值，但並不保證一定是全域極值。 範例 2-8. min f(x) = 3x4 − 28x3 + 84x2 − 96x + 42 x

定義函數f(x) In[1]: f[x_]:=3x^4-28x^3+84x^2-96x+42

計算臨界點並繪出函數圖，由圖形可知函數f(x)的最小值為f(4) In[2]: Plot[f[x],{x,0,5},Epilog->{Red,PointSize[Large], Map[Point,{x,f[x]}/.Solve[f’[x]==0,x]]}] Out[2]: 40 30 20 10

1

2

3

4

5

-10 -20

使用FindMinimum以不同的起始值{0.5,1.5,3}求解 In[4]: FindMinimum[f[x],{x,#}][[1]]&/@{0.5,1.5,3} Out[4]: {5., 5., −22.}

CHAPTER 2. 非線性最佳化

39

由上式輸出可發現當起始值為0.5或1時，都只求得區域極值。

以NMinimize求解也只求得區域極小值 In[5]: NMinimize[f[x],x] Out[5]: {5., {x → 1.}}

以Minimize求解即可求得全域極值 In[6]: Minimize[f[x],x] Out[6]: {−22, {x → 4}}

範例 2-9. min f(x1 , x2 ) = 100(x2 − x21 )2 + (1 − x1 )2 x1 ,x2

定義函數f(x1 , x2 ) In[1]: f[x1_,x2_]:=100*(x2-x1^2)^2+(1-x1)^2 Out[1]:

計算梯度及海賽矩陣 In[2]: g=D[f[x1,x2],{{x1,x2}}];h=D[f[x1,x2],{{x1,x2},2}]; Out[2]:

以(0, 0)為起始值，利用修正牛頓法進行求解 In[3]: delta[x1_,x2_]=Inverse[h].g; sol=NestList[(#-a*delta@@#)/.FindRoot[D[f@@(#-a*delta@@#),a], {a,0},AccuracyGoal->6]&,{0,0},8];sol//TableForm Out[3]:

0 0.161262 0.367189 0.630157 0.807286 0.976581 0.997991 1.00003 1.

0 0 0.106009 0.380126 0.639082 0.951533 0.996209 1.00006 1.

2.2 非線性最佳化

40

利用上式所得資料繪出修正牛頓法求解的步驟圖 In[4]: ContourPlot[f[x1,x2],{x1,0,1},{x2,0,1},Epilog->{Green,Line[sol], Red,PointSize[Large],Point[sol]}] Out[4]: 1.0

0.8

0.6

0.4

0.2

0.0 0.0

0.2

0.4

0.6

0.8

1.0

以(0, 0)為起始值利用FindMinimum求解 In[5]: FindMinimum[f[x1,x2],{x1,0},{x2,0}] Out[5]: {0., {x1 → 1., x2 → 1.}}

FindMinimum使用各種方式進行求解，並傳回迭代次數 In[6]: (Clear[x1,x2];sol={};FindMinimum[f[x1,x2],{x1,0},{x2,0},Method->#, EvaluationMonitor:>AppendTo[sol,{Length[sol],x1,x2}]];Length[sol])& /@{Gradient,ConjugateGradient,InteriorPoint,QuasiNewton,Newton, LevenbergMarquardt,PrincipalAxis} Out[6]: {149, 57, 20, 30, 17, 15, 335}

引入最佳化套件繪製求解步驟圖 In[7]: <
CHAPTER 2. 非線性最佳化

41

由於目標函數包含平方和，故FindMinimum以LevenbergMarquardt法進行求解。由輸出資料可知 經由12次迭代計算後，(x∗ , y∗ ) = (1, 1) In[8]: FindMinimumPlot[f[x1,x2],{{x1,0},{x2,0}}][[{1,2}]] Out[8]: {{0., {x1 → 1., x2 → 1.}}, {Steps → 12, Function → 15}}

繪製求解的步驟圖 In[9]: FindMinimumPlot[f[x1,x2],{{x1,0},{x2,0}}][[3]] Out[9]: 1.0

0.8

0.6

0.4

0.2

0.0 0.0

0.2

0.4

0.6

0.8

1.0

由於在本範例中具有解析解，所以Minimize可順利求解 In[10]: Minimize[f[x1,x2],{x1,x2}] Out[10]: {0, {x1 → 1, x2 → 1}}

NMinimize使用各種方式進行求解，並傳回迭代次數 In[11]: (Clear[x1,x2];sol={};NMinimize[f[x1,x2],{x1,x2},Method->#, EvaluationMonitor:>AppendTo[sol,{Length[sol],x1,x2}]];Length[sol])& /@{NelderMead,DifferentialEvolution,SimulatedAnnealing,RandomSearch} Out[11]: {129, 2032, 264, 86}

2.2 非線性最佳化

42

( ) dp ∫ ∞ dK y 範例 2-10. T CU(y, R) = +h R+ −µ + (x − R)f(x)dx y 2 y R

定義成本函數 In[1]: TCU[y_,R_]:=d*K/y+h*(R+y/2-mu) +d*p/y*Integrate[(-R+x)*f[x],{x,R,Infinity}]; Out[1]:

定義常態分配 In[2]: dist=NormalDistribution[mu,sigma]; Out[2]:

定義常態分配機率密度函數 In[3]: f[x_]:=PDF[dist,x]; Out[3]:

定義各參數值 In[4]: {d,K,h,p,mu,sigma}={1000,100,2,10,50,2}; Out[4]:

以NMinimize求解，可求得期望成本最小值為640.304 In[5]: NMinimize[TCU[y,R], {y,R}] Out[5]: {640.304, {y → 317.098, R → 53.0534}}

最後，我們為參數d, K, h及p進行敏感性分析。敏感性分析的方式為針對每個參數進行-50%, -25%, 25%及50%的變動以觀察各參數對期望成本的影響。

建立敏感性分析所需參數資料 In[6]: w=Range[0.5,1.5,0.25];para={1000,100,2,10}; sen=Flatten[Table[para/.para[[i]]->para[[i]]*w[[j]],{i,3},{j,5}],1]; Out[6]:

CHAPTER 2. 非線性最佳化

43

傳回敏感性分析資料，並將資料資定給sendata In[7]: sendata=Block[{d=#[[1]],K=#[[2]],h=#[[3]],p=#[[4]],ans}, ans=NMinimize[{TCU[y,R],y>0&&R>0},{y,R}]; {ans[[1]],y/.ans[[2]],R/.ans[[2]]}]&/@sen; Out[7]:

建立各參數敏感性分析表格 In[8]: myexample[x_?IntegerQ]:=TableForm[Partition[data,5][[x]], TableHeadings->{If[#==3,ToString[0],ToString[(#-1)*25-50]<>"%"]&/@ Range[5],{"TCU","y","R"}}] Out[8]:

傳回d的敏感性分析 In[9]: myexample[1] Out[9]: -50% -25% 0 25% 50%

TCU y R 454.439 224.536 52.6838 555.317 274.755 52.9038 640.304 317.098 53.0534 715.147 354.407 53.1661 782.792 388.139 53.2562