Lesson 4: 數值與符號為積分計算

Objective

學習 Python 在微分與積分的能力，包括數值（Numerical）與符號（Symbolic）方法。
學習使用數值與符號方法的時機。
從機率與統計學裡的積分找到應用問題。

Prerequisite:

基本的 Python 語法與繪圖技巧。
迴圈的應用。

範例 1：數值微分（Numerical differentiation）
For the function $f(x) = e^{2x}$ , compute $f'(x), f''(x), f'''(x)$ at $x=a$ .

注意事項：

使用 scipy 套件計算數值微分。
數值導函數的定義 $f'(x) = \frac{f(x+dx)-f(x)}{dx}$
數值導函數 $f'(x)$ 的準確度與 $dx$ 的大小有關，試著讓 $dx$ 逐漸變小來觀察準確度.

import numpy as np
import matplotlib.pyplot as plt
from scipy.misc import derivative

f = lambda x: np.exp(2*x)
a = 0 # evaluate at a

fp = derivative(f, a, dx = 0.001, n = 1) # nth derivative
fpp = derivative(f, a, dx = 0.001, n = 2)
fppp = derivative(f, a, dx = 0.001, n = 3, order = 5) # use 5 points to calculate

print(fp, fpp, fppp)

練習： Observe the value of $f'(a)$ by decreasing “dx” from $10^{-1}, 10^{-2}, \cdots, 10^{-16}$ , where $f(x) = -x^3+6x^2$ and a=1.

注意事項：

想想如何建立 $[10^{-1}, 10^{-2}, \cdots, 10^{-16}]$ 向量。再來對照numpy 的指令 logspace(-1, -16, num=16)。
熟練用 print 列印出迴圈過程的結果，如下左圖。必須顧及數字的長度與適合的格式（format），才能呈現準確的觀察值。
觀察哪一個 dx 值產生最準確的結果（9），並不是如理論定義的越小越好。
繪製如下右圖的折線圖，呈現 $f'(1)$ vs. “dx”，其中為了配合 dx 的大小，可以使用 matplotlib.pyplot.semilogx() 或調整 X 軸的格式為 ax.set_xscale(‘log’)，否則會擠成一團。

練習： Draw $f(x)$ and $f'(x)$ , where $f(x) = e^{-x} \sin(x)$

注意事項：

當呈現兩條線在同一個座標區域時，最好除了顏色外，也可以加上線條型態（line style）以示區別。
另外，也可以加入 legend，讓文字來說明線條型態。
這個問題不需要使用迴圈。

範例 2：符號微分（Symbolic differentiation）

Let $f(x) = e^{-x} \sin(x)$ .

以符號方式計算 $f'(x)$ 與 $f''(x)$ 。
計算 $f'(a), a$ 為自訂的數字。
依上述的方式，繪製 $f'(x), -3 \leq x \leq 3$ ，結果如下圖。

注意事項：

使用 sympy 套件計算符號微分 symbolic differentiation.
當使用 sympy 套件最符號運算時，必須使用專屬的數學函數，譬如，sympy.exp(), sympy.sin()，不再使用原來的 np.exp(), np.sin().
符號運算除了看到符號來表達的結果外，常常必須轉換為數值模式，提供後續的計算之用。
sym.lambdify 將符號函數轉為數值型的 lambda 函數，供後續使用，譬如繪圖。

import sympy as sym
import matplotlib.pyplot as plt

x = sym.Symbol('x')
f = sym.exp(-x) * sym.sin(x)
fp = f.diff(x) # the first derivative
fpp = f.diff(x, 2) # the second derivative 
print(fp) # print out to see the symbolic function

# substitution; need to convert to float if necessary
print(float(fp.subs(x, 3))) # evaluate f'(3)
# evaluate symbolic function at more than one points
feval = sym.lambdify(x, fp) # transfer from symbolic function to lambda function
a = np.linspace(-3, 3, 1000)

plt.plot(a, feval(a)) # use symbolic computation to draw a figure
plt.show()
print(feval(3)) # compare with the previous one

範例 3：符號繪圖（Symbolic plot）

Draw $f(x) = e^{-x} \sin(x)$ in symbolic way.

注意事項：

Sympy 套件有自己的繪圖方式 sympy.plotting，有別於以數值方式繪圖的 matplotlib.pyplot 套件.
譬如，畫兩條函數線的方式便大不同，細節請參考下列程式碼。
請參考 Sympy 套件的使用手冊，了解繪圖指令的多種參數選擇與使用方式（有些與 matplotlib.pyplot 不一樣）。

import sympy as sym
from sympy.plotting import plot

x = sym.Symbol('x')
f = sym.exp(-x) * sym.sin(x)
fp = f.diff(x)
# draw two separate figures
# f1 = plot(f, xlim = [-4, 4], ylim = [-10, 5])
# f2 = plot(fp, xlim = [-4, 4], ylim = [-10, 5])
# draw two lines in a single figure
f1 = plot(f, xlim = [-4, 4], ylim = [-10, 6], \
    line_color = 'b', label = 'f(x)', show = False)
f2 = plot(fp, xlim = [-4, 4], ylim = [-10, 6], \
    line_color = 'g', label = "f'(x)", show = False)
f1.append(f2[0]) # append f2 to f1
f1.legend = True
f1.show()

練習： For $f(x) = \frac{\ln x}{x^3}$ , draw $f(x), f'(x)$ in a figure in symbolic way.

注意事項：

畫這張圖時，要小心選取適當的 x 值繪圖範圍與 Y 軸的成像範圍，方能同時看清楚兩個函數的極值。

範例 4：數值積分（Numerical Integration）

Compute the definite integral $\int_0^1 \sin x^2 \;dx$

注意事項：

數值積分當然是針對定積分（definite integral）。而定積分也可以使用符號運算的方式達成。
數值積分使用 Scipy.integrate 指令並選擇二次函數 quadratic rule 為定積分的計算方式。
Scipy.integrate 的定積分計算方式將產生兩個輸出值，請到官網查詢兩個輸出值代表的意思。
練習繪製定積分代表的面積（在函數下塗上顏色 fill_between）。

import numpy as np
import matplotlib.pyplot as plt
import scipy.integrate as integral

# define anonymous function
f = lambda x : np.sin(x**2)
lb, ub = 0, 1
result = integral.quad(f, lb, ub)
print("The integration is {:.4f}".format(result[0]))

# draw the function line -------------------
x = np.linspace(0, np.pi/2, 100)
y = f(x)
fig, ax = plt.subplots(1, 1) 
ax.plot(x, y)
# draw the area ----------------------------
x = np.linspace(lb, ub, 100)
y = f(x)
# fill color between y and 0
ax.fill_between(x, y, 0, alpha=0.3, color='b')
ax.set_xlabel('X')
ax.set_title("The area = {:.4f}".format(result[0]))
plt.grid(True, linestyle='--', which='major')
plt.show()

練習： 計算定積分 $\int_{-2}^6 \sqrt{x^2+4x+12}\;dx$ 並繪製所代表的不規則面積

注意事項：

順帶將上述定積分的數學式以 text 的方式秀在圖形空白處。

範例 4：數值積分（Numerical Integration）

Draw the function $F(x) = \int_{-\infty}^{x} \frac{1}{\sqrt{2\pi}} e^{-t^2/2}\; dt$ for $-5 \leq x \leq 5$ .

注意事項：

$F(x)$ 是標準常態分配的 CDF 函數。還沒畫圖前，已經先認出來，而且知道最後畫出來的樣子。我們常用這種熟悉的函數來幫助自己確認程式的對錯。
寫這支程式之初，也許第一個想法是利用迴圈的方式，來計算每個 x 值（[-5, 5] 之間）代入的定積分。那就依照第一個想法來寫，等完成任務了，再來問問有沒有更好的方法。
下列的程式碼採用「向量函數」的方式，避開傳統的迴圈的模式。也可以說是採用類似平行計算的概念，一次計算完所需要的 $F(x)$ 。這個方式值得記下來，常常用。
採用 Python 的副程式 def 負責計算定積分。
右下圖採用 step plot。
當縮小 step plot 的 step 時，就會變成一條平滑的曲線。

import numpy as np
import matplotlib.pyplot as plt
import scipy.integrate as integral

n = 100
f = lambda x : np.exp(-x**2/2) / np.sqrt(2*np.pi)
xlim = [-5, 5]
x = np.linspace(xlim[0], xlim[1], n)

def cdf(x):
    return integral.quad(f, -np.inf, x)[0]

vec_F = np.vectorize(cdf) # define a vectorized function
F = vec_F(x) # evalaute cdf for all x
plt.plot(x, F, drawstyle='steps-pre') # stairs plot
plt.title('A stair plot of CDF function by integration')
plt.show()

範例 4：數值積分（Numerical Integration）

Express the following definite integration as a Riemann sum.

$\int_1^8 x^2 + 3x + 5 \; dx \approx \sum_{k=1}^{N} f(x_k)\Delta_k$

注意事項：

在本範例，令黎曼和的 N 區間為等寬距離即 $\Delta_k = (8-1)/N$ .
學習如何為長方形塗上顏色（patch a rectangle）。
當附檔名為 py 時，可以加入時間暫留的動畫技巧。
慢慢提高 N，即增加長方形數量，看看圖形變成甚麼？

import numpy as np
import matplotlib.pyplot as plt
import scipy.integrate as integral

N = 10 # how many intervals
n = 100
lb, ub = 1, 8
f = lambda x : x**2 + 3*x +5
xlim = [0, 10]
x = np.linspace(xlim[0], xlim[1], n)
plt.plot(x, f(x))

rec = np.linspace(lb, ub, N+1)
for i in range(N) :
    plt.fill_between([rec[i], rec[i+1]], 
        [f(rec[i+1]), f(rec[i+1])], 0, alpha=0.3, color='b')
    # plt.pause(0.5) # animation not work in Jupyter notbook 
plt.title('Integration by Riemann Sum')    
plt.show()

練習：

前一個範例的黎曼和使用 right endpoints 的函數值，做為長方形的高。試試 left endpoints 與 midpoints 並比較三種方式在相同 N 個區間的條件下，面積的差異（哪一個比較準確）。
將黎曼和的長方形換成矩形，比較面積的差異。
Scipy.integrate 套件也有採用梯形法的選項，試試這個做法的定積分值。並與自己做的梯形法比較看看。

範例 5：符號積分（Symbolic Integration）

For the function $f(x) = \int \sqrt{x^2 + 4x + 12}$

Compute $\int f(x) \; dx$
Evaluate .

注意事項：
1. 下列程式碼中，符號積分的用法採 f.integrate(x)，另一個方法是 sym.integrate(f)。也就是將 integrate 當作 function 或是物件的 method。
2. 用符號積分計算不定積分時，可以用 print 顯示出符號的結果。但是當積分式比較複雜時，便呈現出來，譬如本題的函數。試試比較簡單的函數，如 $f(x) = \sin x$ or $f(x) = \sqrt{1+x}$ 並用 print 來觀察積分的結果。

import sympy as sym

x = sym.Symbol('x')
f = sym.sqrt(x**2 + 4*x + 12)
# f = sym.sin(x) # check this simple function
# F = sym.integrate(f, x)
F = f.integrate(x) # indefinite integral
print(F)
# Definite integral with a symbolic result
I = sym.integrate(f, (x, -2, 6)) 

print("The definite integral is {:.4f}".format(float(I)))

練習題：
Evaluate the following definite integrals:

$\int_0^5 100 e^x\; dx$
$\int_0^7 \frac{x^2}{2} e^{-x}\;dx$ (This is Gamma(3,1))
$\int_{-1.96}^{1.96} \frac{1}{\sqrt{2\pi}} e^{-x^2/2}\;dx$
$E[\theta] =\frac{\int_0^1 (2+\theta)^{125}(1-\theta)^{38}\theta^{35}\;d\theta}{\int_0^1 (2+\theta)^{125}(1-\theta)^{38}\theta^{34}\;d\theta}$