汪群超 Chun-Chao Wang

Dept. of Statistics, National Taipei University, Taiwan

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

Objective

  1. 學習 Python 在微分與積分的能力,包括數值(Numerical)與符號(Symbolic)方法。
  2. 學習使用數值與符號方法的時機。
  3. 從機率與統計學裡的積分找到應用問題。

Prerequisite:

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

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

注意事項:

  1. 使用 scipy 套件計算數值微分。

  2. 數值導函數的定義 f'(x) = \frac{f(x+dx)-f(x)}{dx}

  3. 數值導函數 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.

注意事項:

  1. 想想如何建立 [10^{-1}, 10^{-2}, \cdots, 10^{-16}] 向量。

  2. 熟練用 print 列印出迴圈過程的結果,如下左圖。必須顧及數字的長度與適合的格式(format),才能呈現準確的觀察值。

  3. 觀察哪一個 dx 值產生最準確的結果(9),並不是如理論定義的越小越好。

  4. 繪製如下右圖的折線圖,呈現 f'(1) vs. “dx”,其中為了配合 dx 的大小,必須使用 matplotlib.pyplot.semilogx(),否則會擠成一團。


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

注意事項:

  1. 當呈現兩條線在同一個座標區域時,最好除了顏色外,也可以加上線條型態(line style)以示區別。

  2. 另外,也可以加入 legend,讓文字來說明線條型態。

  3. 這個問題不需要使用迴圈。


範例 2:符號微分(Symbolic differentiation)

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

  1. 以符號方式計算 f'(x)f''(x)

  2. 計算 f'(a), a為自訂的數字。

  3. 依上述的方式,繪製 f'(x), -3 \leq x \leq 3



注意事項:

  1. 使用 sympy 套件計算符號微分 symbolic differentiation.

  2. 當使用 sympy 套件最符號運算時,必須使用專屬的數學函數,譬如,sympy.exp(), sympy.sin(),不再使用原來的 np.exp(), np.sin().

  3. 符號運算除了看到符號來表達的結果外,常常必須轉換為數值模式,提供後續的計算之用。

    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.

注意事項:

  1. Sympy 套件有自己的繪圖方式 sympy.plotting,有別於以數值方式繪圖的 matplotlib.pyplot 套件.

  2. 譬如,畫兩條函數線的方式便大不同,細節請參考下列程式碼。

  3. 請參考 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.

注意事項:

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


範例 4:數值積分(Numerical Integration)

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

注意事項:

  1. 數值積分當然是針對定積分(definite integral)。而定積分也可以使用符號運算的方式達成。

  2. 數值積分使用 Scipy.integrate 指令並選擇二次函數 quadratic rule 為定積分的計算方式。

  3. Scipy.integrate 的定積分計算方式將產生兩個輸出值,請到官網查詢兩個輸出值代表的意思。

  4. 練習繪製定積分代表的面積(在函數下塗上顏色 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 並繪製所代表的不規則面積

注意事項:

  1. 順帶將上述定積分的數學式以 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.

注意事項:

  1. F(x) 是標準常態分配的 CDF 函數。 還沒畫圖前,已經先認出來,而且知道最後畫出來的樣子。我們常用這種熟悉的函數來幫助自己確認程式的對錯。

  2. 寫這支程式之初,也許第一個想法是利用迴圈的方式,來計算每個 x 值([-5, 5] 之間)代入的定積分。那就依照第一個想法來寫,等完成任務了,再來問問有沒有更好的方法。

  3. 下列的程式碼採用「向量函數」的方式,避開傳統的迴圈的模式。也可以說是採用類似平行計算的概念,一次計算完所需要的 F(x)。 這個方式值得記下來,常常用。

  4. 採用 Python 的副程式 def 負責計算定積分。

  5. 右下圖採用 step plot。

  6. 當縮小 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

注意事項:

  1. 在本範例,令黎曼和的 N 區間為等寬距離即 \Delta_k = (8-1)/N.

  2. 學習如何為長方形塗上顏色(patch a rectangle)。

  3. 當附檔名為 py 時,可以加入時間暫留的動畫技巧。

  4. 慢慢提高 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()    

練習:

  1. 前一個範例的黎曼和使用 right endpoints 的函數值,做為長方形的高。試試 left endpoints 與 midpoints 並比較三種方式在相同 N 個區間的條件下,面積的差異(哪一個比較準確)。

  2. 將黎曼和的長方形換成矩形,比較面積的差異。

  3. Scipy.integrate 套件也有採用梯形法的選項,試試這個做法的定積分值。並與自己做的梯形法比較看看。


範例 5:符號積分(Symbolic Integration)

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

  1. Compute \int f(x) \; dx
  2. Evaluate \int_{-2}^6 f(x) \; dx.

    注意事項:

    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:

  1. \int_0^5 100 e^x\; dx
  2. \int_0^7 \frac{x^2}{2} e^{-x}\;dx (This is Gamma(3,1))
  3. \int_{-1.96}^{1.96} \frac{1}{\sqrt{2\pi}} e^{-x^2/2}\;dx
  4. 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}


注意事項:

  1. 盡量採用多種方式計算同一題,看看不同方式是否產生差異。做為未來選擇積分計算方式的參考。

  2. 前三題也可以加入圖形繪製,繪出定積分所代表的面積,將積分值秀在 title,並將整個積分式打在圖中的空白處。


練習題:

繪製下列兩個函數(或隱函數)所圍成的不規則面積(塗色),並計算面積。

  1. f(x) = 2x - x^2, \;\; g(x) = x^2.

  2. f(x) = x^3-3x^2+1, \;\; g(x) = -x^2+2x.

  3. x+2y=4,\;\;  x-y^2=1.

%d bloggers like this: