兩種近似圓周率的方法

兩種近似圓周率的方法

黃見利

本文謹獻給作者的恩師 Dar´ıo Castellanos 教授 (生於1937年12月4日; 逝於1995年11 月 23日)。 作者永遠懷念他的教誨。

介紹

要獲得近似圓周率的公式並不是困難的事情。 舉例來說, 我們已有下列級數, sin x = x − x³

6 + x⁵

120 − · · · , cos x = 1 −x²

2 + x⁴

24− · · · , 1

2sin 2x = sin x cos x = x − 2

3x³ + 2

15x⁵− · · · 如果只考慮前兩式, 然後忽略五次或五次以上的項, 我們便可得到

3 sin x − x cos x ≈ 2x, 因而

x≈ 3 sin x

2 + cos x, (1)

這是首先由 Nicolaus Cusanus 給出的公式; 以後荷蘭的數學家兼物理學家 Willebrord Snel- lius [1] 也獨立地發現了此公式。 Snellius亦即發現光學上著名的折射原理 (Snell’s law) 的 Snell。 兩個名字的差異起因於拉丁文和英文的拼法不同。

現在, 如果考慮所有三式, 然後忽略七次或七次以上的項, 我們便可得到 14 sin x − 6x cos x + sin x cos x ≈ 9x

或

x≈ sin x14 + cos x

9 + 6 cos x, (2)

77

這是首先由數學巨人 Isaac Newton 所發現。

由這兩條公式所得到的近似值會越來越好, 如果所取的角度越來越小。 若取角度為 π

12, 公式 (1) 給了圓周率值 π ≈ 3.141509994, 但是 Newton 的公式給了 π ≈ 3.141592169。 當然, 真 正的數值是 3.141592653589793238...

大師的足跡

在此得到 (1) 和 (2) 的方法是由 Newton 所導出 [2]。 現在就讓我們遵循著大師的路徑 來擴展他的過程。

首先考慮下列方程式

x≈ A¹sin x−A²xcos x+A3sin 2x−A⁴xcos 2x+· · ·+A^2s−1sin sx−A^2sxcos sx, (3) 在此 A_k 為有待決定的常數。 接著利用適當的Maclaurin 級數展開右手邊的項, 我們可以發現 下列系統

A₁− A₂+ 2A3− A₄+ · · · + sA_2s−1− A_2s= 1 A₁− 3A2+ 2³A₃− 3 · 2²A₄+ · · · + s³A_2s−1− 3s²A_2s= 0

· · · ·

A₁− (4s−1)A²+ 2^4s−1A₃− (4s−1)2^4s−2A₄+ · · · + s^4s−1A_2s−1− (4s−1)s^4s−2A_2s= 0 有 2s 條方程式及 2s 未知數可用來解 A_k。 解(3) 中的 x 我們發現

x≈ A₁sin x + A3sin 2x + A5sin 3x + · · · + A^2s−1sin sx 1 + A₂cos x + A₄cos 2x + A₆cos 3x + · · · + A2scos sx. 取 s = 1, 則 A1 = 3

2, A2 = 1

2, 我們得到了公式(1)。

為了列出其他 s 值所對應的公式, 我們利用加法公式 sin 2x = 2 sin x cos x, cos 2x = cos²x− sin²x, 等等, 以避免三角函數多倍角度上複雜的計算。 取 s = 2, 得到公式

x≈ 5

3sin x 16 + 5 cos x

17 + 16 cos x + 2 cos²x, (4) 取 s = 3, 得到公式

x≈ 7

5sin x 92 + 66 cos x + 7 cos²x

82 + 111 cos x + 36 cos²x+ 2 cos³x, (5) 取 s = 4, 得到公式

x≈ 1

35sin x 91648 + 103511 cos x + 28544 cos²x+ 1522 cos³x

1667 + 2944 cos x + 1560 cos²x+ 256 cos³x+ 8 cos⁴x. (6)

讓 x = π

6, 我們得到一些有趣的圓周率近似值。 公式(4), (5) 和 (6) 分別給出 5 32 + 5√

3 37 + 16√

3

= 3.141592229 . . . , (7)

21 5

389 + 132√ 3 436 + 225√

3

= 3.14159265346 . . . , (8)

3 70

45, 224 + 209, 305√ 3 5, 683 + 3, 136√

3

= 3.141592653589754 . . . ; (9)

再說一次, 真正的數值當然是 3.141592653589793238...

這裡有一個有趣的地方值得注意。(8) 和 (9) 改善了下列 Ramanujan 的公式 [5]

63 25

17 + 15√ 5 7 + 15√

5

= 3.14159265380 . . . (10)

雖然以相似大小的係數觀點來看, (10) 確實比 (7) 來得近似圓周率, 但是我們要強調的是 Ra- manujan 導出公式 (10) 所使用的方法是奠基於橢圓模函數的理論 (elliptic modular func- tions), 它所需要考慮到的和牽連到的內容和事情要比我們在此使用的簡單過程要複雜和困難 得多太多了。

這裡還有一個問題尚待回答。 我們如何知道上述系統經常有解? 也就是說, 上述系統的係 數矩陣經常是非奇異的 (non-singular)? 這個問題的答案很有趣, 但是卻不太明顯。 它等同 於去找到 n 個函數使得它們的線性組合 k(x) 和某一被給的函數 f (x) 在下列的觀點下一致:

f^(j)(0) = k^(j)(0), j = 0, 1, . . . , n − 1。

為了證明這是經常可能的, 讓我們考慮一具有 n 個獨立解的 n 階齊次線性微分方程。 讓 k(x) 為此 n 個函數的線性組合, 且具有下列的初始條件: f^(j)(0) = k^(j)(0), j = 0, 1, . . . , n−1。

由於存在和唯一定理, 我們知道如此的 k 值是存在和唯一的。 更甚的是, 若 n = 2s 且 s 為此 微分方程獨立解的個數且為奇數, 則假使 f 為奇函數 (偶函數), 它便可以被解釋為 s 個奇函 數 (偶函數) 的線性組合。 特別地, 下列微分方程

(D²+ 1²)²(D²+ 2²)²· · · (D²+ s²)²y= 0

為 4s 階而且有 2s 個奇函數解如 sin x, x cos x, sin 2x, x cos 2x, . . ., sin sx, x cos sx.

因此, 選擇奇函數 y = x 我們可以寫成

x≈ A¹sin x − A²xcos x + A3sin 2x − A⁴xcos 2x + · · · + A^2s−1sin sx − A^2sxcos sx, 此即 (3)。

但是, 如果我們用矩陣來表達這個系統並且用行列式來證明這個矩陣經常有唯一解, 則可 能會更加有趣。 定義一 (2s) × (2s) 行列式

H_s(u1, . . . , u_s) =             

u₁ −1 u₂ −1 · · · u^s −1

u³₁ −3u²1 u³₂ −3u²2 · · · u³s −3u²s

u⁵₁ −5u⁴1 u⁵₂ −5u⁴2 · · · u⁵s −5u⁴s

... ... ... ... ... ... ... u^4s−1₁ (1 − 4s)u^4s−21 u^4s−1₂ (1 − 4s)u^4s−22 · · · u^4s−1s (1 − 4s)u^4s−2s

             我們可以證明

H_s(u1, . . . , u_s) = −2u³1(u²₂− u²1)⁴· · · (u²s− u²1)⁴H_s−1(u2, . . . , u_s), 而且

H_s(u1, . . . , u_s) = (−2)^s(u1· · · u^s)³ Y

1≤i<j≤s

(u²_j − u²i)⁴.

對所有的 j 讓 u_j = j, 我們看到了這個係數矩陣是非奇異的。 此種型式的行列式我們叫它 confluent Vandermonde 行列式。

超越近似圓周率

現在讓我們考慮下列方程式

x≈ A¹sin x − A²sin 2x + A3sin 3x − · · · + (−1)^s+1A_ssin sx (11)

在等號右手邊的項, 使用適當的正弦級數經過 Maclaurin 級數展開後, 我們可以發現下列系統 A₁− 2A2+ 3A3− · · · + (−1)^s−1sA_s= 1

A₁− 2³A₂+ 3³A₃− · · · + (−1)^s−1s³A_s= 0

... ... ... ... ... ...

A₁− 2^2s−1A₂+ 3^2s−1A₃− · · · + (−1)^s−1s^2s−1A_s= 0 有 2s 條方程式及 2s 未知數可用來解 A_k。

欲證明這個係數矩陣是非奇異的, 我們可以考慮微分方程 (D²+ 1²)(D²+ 2²) · · · (D²+ s²)y = 0。 它有 2s 階並且有 s 個奇函數解如 sin x, sin 2x, · · · , sin sx。

可是, 就如前面已敘述過的一樣, 我們可以使用行列式來證明這個矩陣經常有唯一解。 定 義

G(v1, . . . , v_s) =             

v₁ v₂ v₃ · · · v^s v₁³ v₂³ v₃³ · · · vs³

v₁⁵ v₂⁵ v₃⁵ · · · vs⁵

... ... ... ... ...

v₁^2s−1v₂^2s−1v₃^2s−1· · · vs^2s−1

             我們可以證明

G(v1, . . . , v_s) = v1· · · v^s Y

1≤i<j≤s

(v_j²− v²i).

由於這個行列式本質上也是屬於 Vandermonde 行列式, 而且對所有的 j 讓 vj = (−1)^j−1j, 我們看到了這個係數矩陣是非奇異的。

現在假設 P 是近似於圓周率的一個數值。 讓 x = P − π 表示近似值和真實值之間的誤 差, 同時注意到 sin (k(P − π)) = (−1)^ksin kP 。 因此, 在 (11) 式中, 我們用 P − π 取代 x, 得到

P − π ≈

s

X

k=1

(−1)^k−1A_ksin(k(P − π)). (11.1) 緊接著, 使用 sin (k(P − π)) = (−1)^ksin kP 在 (11.1) 中, 我們得到

P − π ≈ −

s

X

k=1

A_ksin kP , (11.3) 經過移項後, 我們得到了

π≈ P + A¹sin P + A2sin 2P + A3sin 3P + · · · + A^ssin sP . (11.5) 解上面的系統後並令 s = 1 我們發現

π ≈ P + sin P, (12)

這是由 D. Shanks 所給出的 [3]。 讀者知道他是誰嗎? 他在 1961年的七月和作者的恩師 J. W. Wrench, Jr., 利用位於紐約 IBM 資料處理中心的一台超級電腦 IBM 7090, 以 C.

St¨ormer 發現的公式 π

4 = 6 tan⁻¹1 8

+ 2 tan⁻¹ 1 57

+ tan⁻¹ 1 239

 .

作為計算主程式, 而用 Gauss 發現的公式 π

4 = 12 tan⁻¹ 1 18

+ 8 tan⁻¹ 1 57

− 5 tan⁻¹ 1 239

 .

作為核對副程式, 使用 8 小時 43 分鐘的時間計算了 100,265位圓周率。 因此, 間接來說, 他也是 作者的師祖輩呢, 哈哈。

假設 P 為圓周率的 n 位小數近似值, 則 (12) 式將可得到 3n 位小數近似值。 這是因為 我們注意到 (12) 式也可以寫成 π + x − (x − x³

3! −x⁵

5! − · · · ) = π + x³ 3! − x⁵

5! + · · · , 則誤 差 x 可以減少至 x³

3!。

解上面的系統後並令 s 為更大的正整數, 則我們就可得到更佳的近似值:

π ≈ P + 4

3sin P + 1

6sin 2P (13)

得到 5n 位小數,

π ≈ P + 3

2sin P + 3

10sin 2P + 1

30sin 3P (14)

得到 7n 位小數,

π ≈ P + 8

5sin P + 2

5sin 2P + 8

105sin 3P + 1

140sin 4P (15) 得到 9n 位小數, 等等。 例如, 讓 P = 3.1, (12) 給了 π = 3.1415806, (13) 給了 π = 3.1415926494, (14) 給了 π = 3.141592653588, (15) 給了 π = 3.1415926535897922。

若將上面所敘述過的概念和下列的 Fourier 展開作一比較, 我們會發現另一個有趣之處:

x= 2

∞

X

n=1

(−1)ⁿ⁻¹sin nx

n , −π < x < π.

用 P − π 取代 x, 而且如同前面所敘述過的, 利用 sin(k(P − π)) = (−1)^ksin kP , 得到 π= P + 2

∞

X

n=1

sin nP

n , 0 < P < 2π. (16) 在研究過 Fourier 級數後, 我們知道 [4], 當一個系統 S = {φ0, φ₁, φ₂, . . .} 在 [a, b] 上 為單範正交 (orthonormal) 時, 則在此系統的元素所有可能線性組合之中, 若以最小均方誤差 的觀點來看, 會有一個在 [a, b] 上為黎曼可積分的 (Riemann-integrable) 函數 f 其 Fourier 級數的部份和 sn 將會產生對於 f 的最佳可能近似值。 從 (13) 式到 (16) 式都是由正弦函數 所線性組合而成的, 而且對於奇函數 f (x) = π, 0 < x < π, 這些周期為 2π 的正弦函數會產 生越來越好的近似值。

感謝

作者在此希望對其恩師, 臺灣大學數學系程舜仁教授, 在 Vandermonde 行列式方面的教 導表達內心由衷的感謝。

參考文獻

1. A. G. Kaestner, Geschichte der Mathematik, G¨ottingen, 1 (1796), 415.

2. Isaac Newton, Treatise on the Method of Flexions and Infinite Series, London, 1737.

3. D. Shanks, Improving an Approximation for pi, Amer. Math. Monthly 99 (1992), no.

3, 263.

4. T. M. Apostol, Mathematical Analysis, 2nd ed., Addison-Wesley, Reading, MA, (1974), 464-465.

5. S. Ramanujan, Modular equations and approximations to pi, Quart. J. Math. Oxford Ser. 45(1914), 350-372.

—本文作者現就讀於國立臺灣大學數學研究所碩士班三年級_—

∼∼ 跨越二世紀 50年剪報談數學 ∼∼

石厚高先生 (建中退休教師) 將自民國 41 年起所陸續收集之數學剪報已全面電 子化, 放置在網站上供數學愛好者點閱及觀看。 剪報網頁請點選 <建國中學-駝 客記事-教學研究會-數學科-教師團隊-教師陣容-石厚高>或在 yahoo 搜尋裡直 接點取網址<http://math1.ck.tp.edu.tw/石厚高/index.htm>