《用 Python 學微積分》原文見參考資料 1。
16、優化
用一個給定邊長 4 的正方形來折一個沒有蓋的紙盒，設紙盒的底部邊長為 l，則紙盒的高為 (4-l)/2，那么紙盒的體積為：
$$V(l)=l^2\frac{4-l}{2}$$
怎樣才能使紙盒的容積最大？也就是在 l>0,4-l>0 的限制條件下，函數 V(l) 的最大值是多少？
優化問題關心的就是這樣的問題，在滿足限制條件的前提下，怎么才能使目標函數最大（或最小）？
先來看下 V(l) 的函數圖形：
import numpy as np import matplotlib.pyplot as plt l = np.linspace(0,4,100) V = lambda l: 0.5*l**2*(4-l) plt.plot(l,V(l)) plt.show()
看得出來，V(l) 在大約 2.5 處最大。
如果給定一個函數，知道它在點 x=a 處的函數導數為 0，或者該點的導數不存在，則稱該點為關鍵點。要想知道 f(a) 是局部最大值還是局部最小值，可以使用二次導數測試。
如果 f''(a)>0，則函數 f 在 a 處的函數值是局部最小值。
如果 f''(a)<0，則函數 f 在 a 處的函數值是局部最大值。
如果 f''(a)=0，則測試無法告訴我們結論。
上述二次導數測試可以從泰勒級數中推導出來。f(x) 在 x=a 處的泰勒級數為：
$$f(x)=f(a)+f'(a)(x-a)+\frac{1}{2}f''(a)(x-a)^2+\dots$$
因為 a 為關鍵點，f'(a)=0，所以：
$$f(x)=f(a)+\frac{1}{2}f''(a)(x-a)^2+O(x^3)$$
當 f''(a)≠0 時，f(x) 在 x=a 附近的表現近似於二次函數，二次項的系數 0.5f''(a) 決定了拋物線的開口朝向，因而決定了函數值在該點是怎樣的。
回到之前求最大盒子體積的問題，解法如下：
import sympy from sympy.abc import l V = 0.5*l**2*(4-l) # 看看一次導函數： print V.diff(l) # output is : -0.5*l**2 + 1.0*l*(-l + 4) # 一次導函數的定義域為(-oo,oo),因此關鍵點為V'(l)=0的根 cp = sympy.solve(V.diff(l),l) print cp # output is: [0.0, 2.66666666666667] # 找到關鍵點后，使用二次導數測試： for p in cp: print V.diff(l,2).subs(l,p) # output is: 4, -4 # 因此知道在l=2.666666處時，紙盒的體積最大
17、線性回歸
二維平面上有 n 個數據點，pi=(xi,yi)，現嘗試找到一條經過原點的直線，y=ax，使得所有數據點到該直線的殘差的平方和最小。
import numpy as np import matplotlib.pyplot as plt # 設定好隨機函數種子，確保模擬數據的可重現性 np.random.seed(123) # 隨機生成一些帶誤差的數據 x = np.linspace(0,10,10) res = np.random.randint(-5,5,10) y = 3*x + res # 求解回歸線的系數 a = sum(x*y)/sum(x**2) # 繪圖 plt.plot(x,y,'o') plt.plot(x,a*x,'red') for i in range(len(x)): plt.axvline(x[i],min((a*x[i]+5)/35.0,(y[i]+5)/35.0),\ max((a*x[i]+5)/35.0,(y[i]+5)/35.0),linestyle = '--',\ color = 'black') plt.show()
要找到這樣一條直線，實際上是一個優化問題：
$$\underset{a}{min}Err(a)=\sum_{i}{(y_i-ax_i)}^2$$
要找出函數 Err(a) 的最小值，首先計算一次導函數：
$$\frac{dErr}{da}=\sum_{i}2(y_i-ax_i)(-x_i)$$
$$\qquad = -2\sum_{i}x_iy_i + 2a\sum_{i}x_i^2$$
令該函數為 0，求解出關鍵點：
$$a = \frac{\sum_{i}x_iy_i}{\sum_{i}x_i^2}$$
使用二次導數測試：
$$\frac{d^2Err}{da^2}=2\sum_ix_i^2>0$$
因此
$$a = \frac{\sum_{i}x_iy_i}{\sum_{i}x_i^2}$$
是能夠使函數值最小。這也是上面 Python 代碼中，求解回歸線斜率所用的計算公式。
如果不限定直線一定經過原點，即公式為 y=ax+b，則同樣還是一個優化問題，只不過涉及的變量變成兩個。
$$\underset{a}{min}Err(a,b)=\sum_i{(y_i-ax_i-b)}^2$$
這個問題是多元微積分里要分析的問題，這里先給出 Python 中的求解方法。紅線為經過原點的直線，藍線為不限定經過原點的直線。
import numpy as np import matplotlib.pyplot as plt # 設定好隨機函數種子，確保模擬數據的可重現性 np.random.seed(123) # 隨機生成一些帶誤差的數據 x = np.linspace(0,10,10) res = np.random.randint(-5,5,10) y = 3*x + res # 求解回歸線的系數 a = sum(x*y)/sum(x**2) slope, intercept = np.polyfit(x,y,1) # 繪圖 plt.plot(x,y,'o') plt.plot(x,a*x,'red') plt.plot(x,slope*x+intercept, 'blue') for i in range(len(x)): plt.axvline(x[i],min((a*x[i]+5)/35.0,(y[i]+5)/35.0),\ max((a*x[i]+5)/35.0,(y[i]+5)/35.0),linestyle = '--',\ color = 'black') plt.show()
18、不定積分
根據加速度 a(t) 求速度 v(t)，可得：
$$\frac{dv}{dt}=a(t)$$
一旦找到了 v(t)，那么
$$\forall C\in\mathbb{R}, v(t)+C$$
也都是方程的解，因此常微分方程的解是 v(t)+C 這樣一系列函數。得出這一系列函數后，只需知道任一時刻汽車的速度，便可求出常數項 C。
Python 中積分的方法：
import sympy t = sympy.Symbol('t') v = t**3-3*t-6 a = v.diff() print a.integrate() # result is : t**3 - 3*t print sympy.integrate(sympy.E**t+3*t**2) # result is : t**3 + exp(t)
參考資料：
[1] https://ryancheunggit.gitbooks.io/calculus-with-python/content/