《用 Python 學微積分》原文見參考資料 1。
13、大 O 記法
比較兩個函數時，我們會想知道，隨着輸入值 x 的增長或減小，兩個函數的輸出值增長或減小的速度究竟誰快誰慢。通過繪制函數圖，我們可以獲得一些客觀的感受。
比較 x!、ex、x3 和 log(x) 的變化趨勢。
import numpy as np import sympy import matplotlib.pyplot as plt x = range(1,7) factorial = [np.math.factorial(i) for i in x] exponential = [np.e**i for i in x] polynomial = [i**3 for i in x] logarithmic = [np.log(i) for i in x] plt.plot(x,factorial,'black',\ x,exponential, 'blue',\ x,polynomial, 'green',\ x,logarithmic, 'red') plt.show()
根據上圖，當 x—>∞ 時，x!>ex>x3>ln(x)。要想證明的話，可以取極限，用洛必達法則，例如：
$$\lim_{x\rightarrow \infty}\frac{e^x}{x^3}=\infty$$
表明當 x—>∞ 時，雖然分子分母都在趨向無窮大，但分子遠遠凌駕於分母之上。類似地，也可以這樣看：
$$\lim_{x\rightarrow \infty}\frac{ln(x)}{x^3}=0$$
表明分母遠遠凌駕於分子之上。
SymPy 是 Python 的數學符號計算庫，用它可以進行數學公式的符號推導。下面代碼用 SymPy 來推導上面兩式。
import sympy from sympy.abc import x # sympy中無限infty用oo表示 print ((sympy.E**x)/(x**3)).limit(x,sympy.oo) # result is: oo print (sympy.ln(x)/x**3).limit(x,sympy.oo) # result is 0
為了描述這種隨着 x—>∞ 或 x—>0 時函數的表現，我們定義如下大 O 記法：
若我們稱函數 f(x) 在 x—>0 時是 O(g(x))，則需要找到一個常數 C，對於所有足夠小的 x，均有 |f(x)|<C|g(x)|。
若我們稱函數 f(x) 在 x—>∞ 時是 O(g(x))，則需要找到一個常數 C，對於所有足夠大的 x，均有 |f(x)|<C|g(x)|。
之所以叫大 O 記法，是因為函數的增長速率很多時候被稱為函數的階（Order）。
例如，當 x—>∞ 時，x(1+x2)1/2 是 O(x2)，下面來個直觀感受。
下圖是兩個函數的變化趨勢，紅線是 x(1+x2)1/2 ，藍線是 2x2 。
import sympy from sympy.abc import x import numpy as np import matplotlib.pyplot as plt xvals = np.linspace(0,100,1000) f = x*sympy.sqrt(1+x**2) g = 2*x**2 y1 = [f.evalf(subs={x:xval}) for xval in xvals] y2 = [g.evalf(subs={x:xval}) for xval in xvals] plt.plot(xvals[:10],y1[:10],'r',xvals[:10],y2[:10],'b') #plt.plot(xvals,y1,'r',xvals,y2,'b') plt.show()
Sympy 可以幫助我們分析函數的階，如下面求 x(1+x2)1/2 的階。
import sympy from sympy.abc import x f = x*sympy.sqrt(1+x**2) print sympy.O(f, (x, sympy.oo)) # result is : O(x**2, (x, oo))
計算機中使用大 O 記法，通常是分析當輸入數據 —>∞ 時，程序在時間或空間上的表現。然而，從上面的介紹，我們知道這個位置可以是 0，甚至可以是任何有意義的位置。
import sympy from sympy.abc import x f = x*sympy.sqrt(1+x**2) print sympy.O(f, (x, 0)) # result is : O(x)
在前面泰勒級數一節，利用 Sympy 取函數泰勒級數的前幾項時，代碼是這樣：
import sympy from sympy.abc import x exp = sympy.E**x sum15 = exp.series(x,0,15).removeO() print sum15
其中 removeO() 的作用是讓 sympy 忽略掉級數展開后的大 O 表示項。不然結果如下：
import sympy from sympy.abc import x exp = sympy.E**x print exp.series(x, 0, 3) # result is: 1 + x + x**2/2 + O(x**3)
這表示從泰勒級數的第 4 項起，剩余所有項在 x—>0 時是 O(x3)。這表明，當 x—>0 時，用 1+x+0.5x2 來近似 ex ，我們得到的誤差上限將是 Cx3，其中 C 是一個常數。也就是說，大 O 記法能用來描述我們使用多項式近似時的誤差。
另外大 O 記法也可以直接參與計算中去，例如要計算：
$$cos(x^2)\sqrt{(x)}$$
在 x—>0 時階 O(x5) 以內的多項式近似，可以這樣：
$$cos(x^2)\sqrt{(x)}=(1-\frac{1}{2}x^4+O(x^6))x^{\frac{1}{2}}$$
$$\qquad = x^{\frac{1}{2}}- \frac{1}{2}x^{\frac{9}{2}} + O(x^{\frac{13}{2}})$$
import sympy from sympy.abc import x print (sympy.cos(x**2)*sympy.sqrt(x)).series(x,0,5) # result is: sqrt(x) - x**(9/2)/2 + O(x**5)
14、導數
對函數某一點求導，得到的是函數在該點處切線的斜率。選中函數圖像中某一點，不斷放大，最后會發現函數圖像一條直線，這條直線就是切線。下面獲得一些直觀的感受。
import numpy as np from sympy.abc import x import matplotlib.pyplot as plt # 函數 f = x**3-2*x-6 # 在x=6處正切於函數的切線 line = 106*x-438 d1 = np.linspace(2,10,1000) d2 = np.linspace(4,8,1000) d3 = np.linspace(5,7,1000) d4 = np.linspace(5.8,6.2,100) domains = [d1,d2,d3,d4] # 畫圖的函數 def makeplot(f,l,d): plt.plot(d,[f.evalf(subs={x:xval}) for xval in d],'b',\ d,[l.evalf(subs={x:xval}) for xval in d],'r') for i in range(len(domains)): # 繪制包含多個子圖的圖表 plt.subplot(2, 2, i+1) makeplot(f,line,domains[i]) plt.show()
導數定義 1：
$$f'(a)=\frac{df}{dx}\bigg|_{x=a}=\lim_{x\rightarrow a}\frac{f(x)-f(a)}{x-a}$$
若該極限不存在，則函數在 x=a 處的導數不存在。
導數定義 2：
$$f'(a)=\frac{df}{dx}\bigg|_{x=a}=\lim_{h\rightarrow 0}\frac{f(a+h)-f(a)}{h}$$
若該極限不存在，則函數在 x=a 處的導數不存在。
導數定義 3：
函數 f(x) 在 x=a 處的導數 f'(a) 是滿足如下條件的常數 C，對於在 a 附近輸入值的微小變化 h，有 f(a+h)=f(a)+Ch+O(h2) 始終成立。也就是說導數 C 輸入值變化中一階項的系數。上式稍加變化，兩邊同時除以 h，並同時取極限可得：
$$\lim_{h\rightarrow 0}\frac{f(a+h)-f(a)}{h}=\lim_{h\rightarrow 0}C+O(h)=C$$
便與上面定義 2 相一致了。
例如求 cos(x) 在 x=a 處的導數：
$$cos(a+h)=cos(a)cos(h)-sin(a)sin(h)$$
$$\qquad = cos(a)(1+O(h^2))-sin(a)(h+O(h^3))$$
$$\qquad = cos(a)-sin(a)h +O(h^2)$$
所以：
$$\frac{d}{dx}{cos(x)}\bigg|_{x=a}=-sin(a)$$
我們可以自己定義求導的函數：
import numpy as np from sympy.abc import x f = lambda x: x**3-2*x-6 # 我們設定參數h的默認值，如果調用函數時沒有指明參數h的值，便會使用默認值 def derivative(f,h=0.00001): return lambda x: float(f(x+h)-f(x))/h fprime = derivative(f) print fprime(6) # result is：106.000179994
Sympy 也提供求導的方法：
from sympy.abc import x f = x**3-2*x-6 print f.diff() # result is :3*x**2 - 2 print f.diff().evalf(subs={x:6}) # result is : 106.0000000000
依據導數的定義 3，有 f(a+h)=f(a)+f'(a)h+O(h2)，如果將高階項丟掉，就獲得了 f(a+h) 的線性近似式子：f(a+h)≈f(a)+f'(a)h。
例如，用線性近似的方法估算 2551/2：
$$\sqrt{256-1}\approx \sqrt{256}+\frac{1}{2\sqrt{256}}(-1)$$
$$\qquad = 16 - \frac{1}{32}$$
$$\qquad = 15\frac{31}{32}$$
15、牛頓迭代法
如何在不使用 x1/2 前提下求 C 的正二次根。
上述問題等價於求 f(x)=x2+C=0 的解，根據上面介紹的線性近似：f(x+h)≈f(x)+f'(x)h 。如果 f(x+h)=0，那么：
$$h\approx -\frac{f(x)}{f'(x)}$$
$$x+h\approx x - \frac{f(x)}{f'(x)}$$
如果我們對 f(x)=0 的解有一個初始估計 x0，便可以用上面的近似不斷獲取更加准確的估計值，方法為：
$$x_{n+1} = x_{n} - \frac{f(x_n)}{f'_{x_n}}$$
將 f(x)=x2+C 代入上式，即得 xn 的更新規則：
$$x_{n+1} = \frac{1}{2}(x_{n}+\frac{C}{x_{n}})$$
from sympy.abc import x def mysqrt(c, x = 1, maxiter = 10, prt_step = False): for i in range(maxiter): x = 0.5*(x+ c/x) if prt_step == True: # 在輸出時，{0}和{1}將被i+1和x所替代 print "After {0} iteration, the root value is updated to {1}".format(i+1,x) return x print mysqrt(2,maxiter =4,prt_step = True) # After 1 iteration, the root value is updated to 1.5 # After 2 iteration, the root value is updated to 1.41666666667 # After 3 iteration, the root value is updated to 1.41421568627 # After 4 iteration, the root value is updated to 1.41421356237 # 1.41421356237
通過繪圖進一步了解這個方法，例如，我們要猜 f(x)=x2-2x-4=0 的解，從 x0=2 開始，找到 f(x) 在 x=x0 處的切線 y=2x-8，找到其與 y=0 的交點 (4,0)，將該交點作為新的猜測解 x1=4，如此循環。
import numpy as np import matplotlib.pyplot as plt f = lambda x: x**2-2*x-4 l1 = lambda x: 2*x-8 l2 = lambda x: 6*x-20 x = np.linspace(0,5,100) plt.plot(x,f(x),'black') plt.plot(x[30:80],l1(x[30:80]),'blue', linestyle = '--') plt.plot(x[66:],l2(x[66:]),'blue', linestyle = '--') l = plt.axhline(y=0,xmin=0,xmax=1,color = 'black') l = plt.axvline(x=2,ymin=2.0/18,ymax=6.0/18, linestyle = '--') l = plt.axvline(x=4,ymin=6.0/18,ymax=10.0/18, linestyle = '--') plt.text(1.9,0.5,r"$x_0$", fontsize = 18) plt.text(3.9,-1.5,r"$x_1$", fontsize = 18) plt.text(3.1,1.3,r"$x_2$", fontsize = 18) plt.plot(2,0,marker = 'o', color = 'r' ) plt.plot(2,-4,marker = 'o', color = 'r' ) plt.plot(4,0,marker = 'o', color = 'r' ) plt.plot(4,4,marker = 'o', color = 'r' ) plt.plot(10.0/3,0,marker = 'o', color = 'r' ) plt.show()
如下定義牛頓迭代法：
def NewTon(f, s = 1, maxiter = 100, prt_step = False): for i in range(maxiter): # 相較於f.evalf(subs={x:s}),subs()是更好的將值帶入並計算的方法。 s = s - f.subs(x,s)/f.diff().subs(x,s) if prt_step == True: print "After {0} iteration, the solution is updated to {1}".format(i+1,s) return s from sympy.abc import x f = x**2-2*x-4 print NewTon(f, s = 2, maxiter = 4, prt_step = True) # After 1 iteration, the solution is updated to 4 # After 2 iteration, the solution is updated to 10/3 # After 3 iteration, the solution is updated to 68/21 # After 4 iteration, the solution is updated to 3194/987 # 3194/987
Sympy 可以幫助我們求解方程：
import sympy from sympy.abc import x f = x**2-2*x-4 print sympy.solve(f,x) # result is：[1 + sqrt(5), -sqrt(5) + 1]
參考資料：
[1] https://ryancheunggit.gitbooks.io/calculus-with-python/content/