6.1 迭代求圓周率與超越函數#

概念說明#

圓周率的計算史就是一部迭代法演化史，四個里程碑各有教學價值：

阿基米德夾逼：從正六邊形出發，反覆倍增多邊形邊數，用周長（或面積）從內外兩側夾住 \(\pi\) 。每次倍增使誤差約縮為 1/4（線性收斂， \(L = 1/4\) ），到 96 邊形只能得兩位小數——歷史悠久但效率低。
萊布尼茲級數： \(\dfrac{\pi}{4} = 1 - \frac13 + \frac15 - \frac17 + \cdots\) 。形式優美，但誤差 \(\sim 1/N\) ，要 100 萬項才有 6 位——教科書經典、實務災難，恰好示範「收斂速度才是王道」。
Machin 型公式： \(\dfrac{\pi}{4} = 4\arctan\frac15 - \arctan\frac{1}{239}\) ，配合 \(\arctan\) 的泰勒級數（每多一項誤差縮 \(1/25\) 以上），線性收斂但速率高，19 世紀靠它算到數百位。
AGM 迭代（Gauss-Legendre）：算術-幾何平均迭代 \(a_{k+1} = \frac{a_k + b_k}{2}\) 、 \(b_{k+1} = \sqrt{a_k b_k}\) ，配上 \(c\) 與 \(S\) 的更新， \(n\) 步後 \(\pi \approx (a+b)^2 / (4S)\) 。平方收斂：每步正確位數翻倍，25 次迭代即達十億位——這是記錄級 \(\pi\) 計算的基礎演算法。

超越函數的計算同理。 \(\sin x\) 、 \(e^x\) 用泰勒級數逐項累加（迭代加總，直到項小於機器精度即停）； \(\arctan\) 用級數或範圍縮減。核心觀念一致：把無窮過程截斷成有限迭代，並用餘項估計決定何時停。交錯級數（如萊布尼茲）的好處是截斷誤差被下一項控制；正項級數則需比較餘項與幾何級數。

生活類比#

量一棵樹的高度：阿基米德法像爬上旁邊更高的樹去比；級數法像用無數張小紙條逐步貼近樹頂；AGM 法像拿一根能「自動對摺精修」的神奇尺，量兩三次就鎖定答案。方法不同，逼近本質相同。

演算法步驟#

以 AGM 求 \(\pi\) 為例：

初始化 \(a_0 = 1\) ， \(b_0 = 1/\sqrt{2}\) ， \(S = \frac12\) 。
迭代： \(a_{k+1} = \frac{a_k + b_k}{2}\) ， \(b_{k+1} = \sqrt{a_k b_k}\) ， \(c_k = a_k - b_k\) ， \(S \leftarrow S - 2^k c_k^2\) 。
當 \(|c_k|\) 足夠小（位數已翻倍足夠多次）停止。
\(\pi \approx \dfrac{(a_{k+1} + b_{k+1})^2}{4S}\) 。

Python 範例#

用 Python 內建的任意精度整數（mpmath 若無則用浮點）示範三法比較：

import math

def archimedes(n_sides0=6, rounds=15):
    n, s = n_sides0, 1.0          # 內接正 n 邊形邊長（半徑 1）
    for _ in range(rounds):
        s = math.sqrt(2 - 2 * math.sqrt(1 - (s/2)**2))  # 倍增邊長
        n *= 2
        yield n, n * s / 2        # 面積逼近 pi

def leibniz(n):
    return 4 * sum((-1)**k / (2*k + 1) for k in range(n))

def agm_pi(rounds=5):
    a, b, S = 1.0, 1/math.sqrt(2), 0.5
    for k in range(1, rounds + 1):
        a, b = (a + b) / 2, math.sqrt(a * b)
        S -= (2**k) * (a - b)**2
    return (a + b)**2 / (4 * S)

print(f"萊布尼茲 10^6 項誤差: {abs(leibniz(10**6) - math.pi):.2e}")
for n, area in archimedes():
    pass
print(f"阿基米德 {n} 邊形誤差: {abs(area - math.pi):.2e}")
print(f"AGM 5 步誤差      : {abs(agm_pi() - math.pi):.2e}")

# 迭代計算超越函數：sin x 之泰勒級數
def my_sin(x, eps=1e-16):
    term, s, k = x, 0.0, 0
    while abs(term) > eps:
        s += term
        k += 1
        term *= -x*x / ((2*k) * (2*k + 1))
    return s, k

val, n = my_sin(1.0)
print(f"sin(1) ≈ {val:.15f}（{n} 項），誤差 {abs(val - math.sin(1)):.2e}")

複雜度與收斂性分析#

萊布尼茲誤差 \(\sim \frac{1}{2N}\) （線性、極慢）；阿基米德每輪位數乘 \(\log_{10} 4 \approx 0.6\) 位；AGM 平方收斂， \(k\) 步後約 \(2^k\) 位正確，故達 \(D\) 位需 \(O(\log D)\) 次迭代，配合快速乘法總成本近 \(O(D \log^2 D)\) 。級數法計算超越函數時，對小 \(x\) 用泰勒級數快，對大 \(x\) 需先做範圍縮減（如 \(\sin(x) = \sin(x \bmod 2\pi)\) ），否則消去誤差吃掉精度。教訓：選迭代法先看收斂階，再談實作。

練習#

實作 Machin 公式求 \(\pi\) ，統計達 12 位正確所需的級數項數，與萊布尼茲比較。
寫出迭代版 \(e^x\) （泰勒級數），與 math.exp 比較，並測 \(x = 50\) 時為何應先算 \(e^{50} = (e^{25})^2\) 。
用 mpmath（pip install mpmath）以 50 位精度重跑 AGM，確認每步正確位數翻倍。