隨機採樣黎曼和 : 定積分數值算法
武國寧 · 孫 娜
引言
在一些工程數學模型求解中, 往往遇到下述問題: (1) 實際問題中有些函數沒有具體的表 達式, 僅僅知道函數在一些離散點的函數值; (2) 有些函數的原函數不能夠用基本初等函數表 示, 例如橢圓積分等; (3) 一些函數的原函數需要採用函數級數或者函數無限乘積來表示。 還有 的函數的原函數需要借助於特殊函數來表示, 例如含參變量積分中的伽馬、 貝塔函數等等。 以 上這些問題使得數值積分具有了重要的意義。 常見的數值積分方法有: (1) 利用函數在一些節 點的加權和來表示定積分。 例如數值積分中的梯形公式、 辛普森公式和高斯節點積分方法等等。
這些積分方法被廣泛應用於一元函數的數值積分。 若求解函數為多元函數, 計算量會按指數增 長。 (2) 另外一種數值積分方法為蒙地卡羅方法, 該方法採用隨機投點, 通過計算落入積分區域 與包含積分區域的一個規則區域內點的比例來計算定積分。 蒙地卡羅方法提高了多元函數數值 積分的計算速度, 應用較為廣泛。 定積分中對積分區間 (或區域) 分隔後, 黎曼和與每個小區間 上隨機取點相關。 若把每個小區間上的選點過程看作一個在該小區間上服從均勻分布隨機變量 的取值, 則黎曼和為隨機變量的函數。 大數定律表明: 對於獨立重複試驗, 頻率穩定於概率。 受 此啟發, 我們可以通過重複試驗, 研究黎曼和的平均值, 以此研究定積分。 數值試驗表明: (1) 當分割給定, 黎曼和平均值隨著試驗次數的增加收斂於函數的定積分, 且隨機黎曼和直方圖呈 常態分布; (2) 當試驗次數給定, 隨機黎曼和平均值隨著分割份數的增加收斂於積分值, 方差隨 著分割份數的增加趨於零。
隨機黎曼和 : 定積分數值算法
在微積分教材中, 定積分定義如下:
定義1: 設 f(x) 為定義在 [a, b] 上的有界函數, 在 [a, b] 中任取分點 {xi}n0, 做成一種劃分:
P : a = x0 < x1 < x2 <· · · < xn=b
並任意取點 ξi ∈ [xi−1, xi]。 記小區間的長度為 Δxi =xi− xi−1,令 λ(P ) = max1≤i≤nΔxi,
89
存在, 且極限值與分割 P 和 ξi 的取法無關, 則稱 f(x) 在 [a, b] 上黎曼可積分。 其極限值 I 為 f (x) 在 [a, b] 上的定積分, 記為: b
a
f (x)dx.
在上述定義中, 任意取點 ξi ∈ [xi−1, xi], 可以認為 [xi−1, xi] 上服從均勻分布隨機變量 Xi 的 取值。 當函數 f(x) 和劃分 P 給定時, 函數的黎曼和
σ(f ;ξ, P ) =
n i=1
f (ξi)Δxi,
可以看作隨機變量 Xi=ξi的表達式 (求一次黎曼和, 對應一次隨機試驗), 這裡 ξ =(ξ1, ξ2, . . ., ξn)。 為分析隨機黎曼和隨著試驗次數和分割份數增加的統計規律, 我們從兩個方面展開數值試 驗:
(1) 當分割給定, 通過增加重複試驗的次數研究黎曼和的穩定性 (頻率穩定於概率);
(2) 當分割加細, 隨機變量增多時, 黎曼和是否近似於常態分布 (中央極限定理)。
數值試驗設計如下: 討論積分
1 0
x2dx, 這裡我們把 [0, 1] 區間 n 等份:
P : 0 n < 1
n < 2
n <· · · < i− 1 n < i
n <· · · < n n = 1 假設 Xi 為第 i 個區間 i− 1
n , i n
上服從均勻分布的隨機變量: Xi ∼ Ui− 1 n , i
n
, 黎 曼和
σ(f ;X = ξ1, P ) =
n i=1
f (ξi1)Δxi 為一次隨機試驗, 這裡 ξ1 = (ξ11, ξ21, . . . , ξn1)。 我們定義平均值如下:
pm = 1 m
m j=1
σ(f ;ξj, P ) = 1 m
m j=1
n i=1
f (ξij)Δxi (1) 這裡 ξj = (ξ1j, ξ2j, . . . , ξnj) 表示隨機向量 X 第 j 次取值, pm 表示 [0, 1] 區間 n 等份後 m 次隨機試驗黎曼和的平均值。
圖 1 為 n = 10 (10等份), 102 次試驗, 黎曼和平均值 pm 的變化趨勢; 圖 2 為 n = 10 (10等份), 104 次試驗, 黎曼和平均值 pm 的變化趨勢。 從圖中可以看出: 隨著試驗次數的增加, 黎曼和平均值趨於 1/3, 這驗證了頻率穩定於概率這一規律。
0 20 40 60 80 100 0.33
0.335
Numbers
Values
0 2k 4k 6k 8k 10k
0.33 0.335 0.34 0.345
Numbers
Values
圖1: 黎曼和平均值隨著試驗次數的變化趨勢, 試驗次數 102 次。
圖2: 黎曼和平均值隨著試驗次數的變化趨勢, 試驗次數 104 次。
圖3 (a) n = 10 (10等份), 104 次隨機試驗, 黎曼和直方分佈圖; (b) n = 10 (10等 份), 104 次隨機試驗。 上曲線為每次隨機試驗黎曼和的變化趨勢, 下曲線為每次試驗黎曼和與 真實值之間的誤差; (c) n = 10 (10等份), 105 次隨機試驗, 黎曼和直方分佈圖; (d) n = 10 (10等份), 105 次隨機試驗。 上曲線為每次隨機試驗黎曼和的變化趨勢, 下曲線為每次試驗黎曼 和與真實值之間的誤差; (e) n = 100 (100等份), 104 次隨機試驗, 黎曼和直方分佈圖; (f) n = 100 (100等份), 104 次隨機試驗。 上曲線為每次隨機試驗黎曼和的變化趨勢, 下曲線為每 次試驗黎曼和與真實值之間的誤差。 從圖中可以看出隨機黎曼和直方圖近似於常態分布, 其中 心為積分的真實值。 當分割的份數增加為 100 份時, 其精準度得到了較大的提高。
下面分析隨機黎曼和平均值 pm 隨分割份數增加的變化趨勢。 圖 4 為等分份數分別為 10, 20, 30,. . ., 500份, 對於每個分割份數試驗次數為 104 次。 從圖中可以看出隨着分割份數的增 加, 黎曼和平均值收斂於真實的積分值。 圖 5 分析了分割份數為 10, 20, 30, . . ., 500 份, 104 次試驗隨機黎曼和的標準差。 從圖中可以看出, 隨著分割份數的增加, 黎曼和的方差趨於零。
0.3 0.31 0.32 0.33 0.34 0.35 0.36 0
0.5 1 1.5 2
2.5 0 2k 4k 6k 8k
0 2k 4k 6k 8k
−0.02 0 0.02
Value
Numbers
Numbers
Count Values
0.3 0.32 0.34 0.36
0 0.5 1 1.5 2
0 20k 40k 60k 80k 100k
0.3 0.32 0.34 0.36
0 20k 40k 60k 80k 100k
−0.04
−0.02 0 0.02 0.04
Value
Numbers
Numbers
Count ValuesValues
0.3325 0.333 0.3335 0.334
0 0.5 1 1.5 2 2.5
0 2k 4k 6k 8k
0.333 0.334
0 2k 4k 6k 8k
−0.001 0 0.001
Value
Numbers
Numbers
Count ValuesValues
(a) (b)
(c) (d)
(e) (f)
圖3: (a) 隨機黎曼和分布直方圖, 等分份數 n = 10, 試驗次數 104 次; (b) 上曲線為隨機黎曼 和, 下曲線為黎曼和與真實積分值的誤差, 等分份數 n = 10, 試驗次數 104 次; (c) 隨機黎曼 和分布直方圖, 等分份數 n = 10, 試驗次數 105 次; (d) 上曲線為隨機黎曼和, 下曲線為黎曼 和與真實積分值的誤差, 等分份數 n = 10, 試驗次數 105 次; (e) 隨機黎曼和分布直方圖, 等 分份數 n = 100, 試驗次數 104 次; (f) 上曲線為隨機黎曼和, 下曲線為黎曼和與真實積分值的 誤差, 等分份數 n = 100, 試驗次數 104 次。
0 100 200 300 400 500 0.33
0.34 0.35 0.36 0.37
Numbers
Values
0 100 200 300 400 500
0 0.005 0.01
Numbers
Values
圖4: 隨機黎曼和平均值的變化趨勢, 分割份數 n = 10, 20, 30, . . . , 500, 試驗次數 n = 104。
圖5: 隨機黎曼和標準差的變化趨勢, 分割份數 n = 10, 20, 30, . . . , 500, 試驗次數 n = 104。
圖 6 為採用蒙地卡羅方法計算積分
1 0
x2dx 得到的圖形。 這裡投點個數為 104 個, 藍 色點為落在函數和座標軸所圍成區域的內部, 黃色點表示落在函數和座標軸所圍成區域的外部。
該方法得到的積分近似值為 0.3300。
0 0.2 0.4 0.6 0.8 1
0 0.2 0.4 0.6 0.8
1 Below
Above func
圖6: 蒙地卡羅隨機投點試驗, 投點個數為 104 個, 積分值近似為 0.3300。
為了對比隨機黎曼和平均值與蒙地卡羅方法的收斂精準度, 我們做了以下對比試驗: 對於
定積分 0 x dx, 2 ln xdx, 0 x dx, 0 e−x dx 採用以上兩種方法結果對比圖。 圖中藍色 曲線代表蒙地卡羅方法, 紅色曲線代表隨機黎曼和方法。 試驗結果表明: 在相同點數的情況下, 隨機黎曼和方法具有較快的收斂速度。
0 2k 4k 6k 8k 10k
0.31 0.32 0.33 0.34 0.35
M. Carlo Riemann
Numbers
Values
0 2k 4k 6k 8k 10k
0.003 0.004 0.005 0.006 0.007 0.008
0.009 M. Carlo
Riemann
Numbers
Values
0 2k 4k 6k 8k 10k
0.78 0.8 0.82 0.84 0.86 0.88 0.9 0.92
0.94 M. Carlo
Riemann
Numbers
Values
0 2k 4k 6k 8k 10k
0.73 0.74 0.75 0.76 0.77 0.78
0.79 M. Carlo
Riemann
Numbers
Values
(a) (b)
(c) (d)
圖7: 數值積分, 從上到下, 從左到右函數依次為1
0 x2dx, 10 2
1
ln xdx, 1 0
sin x x dx, 1
0 e−x2dx。 積分方法: 隨機黎曼和 (紅色); 蒙地卡羅方法 (藍色) 。
結語
本文將定積分中黎曼和在每個小區間的取值看作均勻分布的隨機變量, 從數值試驗的角度 探討了黎曼和的分佈。 試驗表明: 當分割給定, 隨機黎曼和的平均值隨著試驗次數增加收斂於 函數的定積分, 其方差隨著分割份數的增加而減小。 對比蒙地卡羅方法, 說明了隨機黎曼和方法 具有較好的收斂速度。 不足之處在於: 沒有從理論上證明隨機黎曼和隨著分割份數的增加趨於 常態分布。
參考文獻
1. 王能超。 數值分析簡明教程[M], 2版。 武漢: 華中科技大學出版社, 43-59, 2002。
2. 克萊因。 古今數學思想[M]。 上海: 上海科學技術出版社, 1979。
3. Stewart, J. M., Python for Scientists [M], Cambridge University Press, 2017.
4. Robert, C. P., Casella, G., Monte Carlo Statistical Methods [M], Springer, 2009.
5. 陳紀修, 於崇華, 近路。 數學分析, 上冊[M], 2版。 北京: 高等教育出版社, 272-280, 2004。
—本文作者武國寧、 孫娜任教中國石油大學數學系—
