1768年 · 尤拉方法：微分追兇案#

卷宗編號：Case-1768-Euler｜偵探：Leonhard Euler｜委託人：牛頓留下的微分方程們

案發現場#

牛頓1687年立下 \(F = ma\) 之後，物理學家手裡全是微分方程，桌上卻沒有答案。除了少數幸運兒（拋體、單擺小角度、兩體問題），絕大多數方程根本寫不出初等函數解。

1760年代的聖彼得堡科学院，尤拉面對的報案堆積如山：彈道含空氣阻力怎麼算？月球攝動怎麼追？懸鏈線、船體搖擺、潮汐摩擦，哪個不是 \(y' = f(t,y)\) 的模樣？

案情的核心矛盾是：微分方程描述的是「每一瞬間的變化率」，人類卻想要「未來 whole 條軌跡」。從切線方向拼出整條路，這就是本案的委託。

當時沒有電腦，連像樣的函數表都不多。尤拉必須發明一種「瞎子摸象卻能走完全程」的方法：只用加減乘除，一步步往前挪。

偵查過程#

線索一：切線就是指南針#

尤拉的破案靈感樸素到令人心疼：既然方程告訴我們此刻斜率，那就沿著斜率走一小步。

考慮初值問題：

其中 \(f\) 為已知速率函數， \(y_0\) 為起始狀態。把時間切成步長 \(h\) 的小段，記 \(t_n = t_0 + nh\) ， \(y_n\) 為時刻 \(t_n\) 的近似值，則：

這就是尤拉方法：一階泰勒展開砍掉餘項，用切線代替曲線。幾何意義是沿著每一步的斜率 \(f\) 直行 \(h\) 距離。

線索二：誤差的兩張臉#

老練的偵探不會只給答案，還要給誤差口供。尤拉格式的誤差有兩種身份：

推導很直接。真解泰勒展開為 \(y(t+h) = y + h f + O(h^2)\) ，而格式只取前兩項，故單步丟掉 \(O(h^2)\) 。走 \(T/h\) 步，累積成 \(O(h)\) 。步長砍半，誤差約砍半——慢，但穩。

線索三：穩定性的陷阱#

案情在這裡出現反轉：步長太大，切線會把你帶進溝裡。考察試探方程 \(y' = \lambda y\) ，其中 \(\lambda\) 為常數（可為負實數或複數）。尤拉迭代給出：

要讓數值解不爆炸，必須有 \(|1 + h\lambda| \le 1\) 。若 \(\lambda = -10\) ，則步長須滿足 \(h \le 0.2\) 。超過此限，真解明明在衰減，數值解卻振盪發散——這就是顯式格式的「穩定性條件」。

剛性方程（如化學動力學快慢反應並存）正是拿這一點折磨後世，逼出了隱式尤拉與指數積分器。

理論定性：數值格式的原型#

尤拉方法確立了三條辦案規矩，後世所有 Runge-Kutta、Verlet、多步法都是其子孫：

1. 離散化：連續時間切成 \(h\) ，微分化為差分。
2. 一致性加穩定性等於收斂：格式要像真方程（局部誤差小），且誤差不被放大。
3. ** cost 精度權衡**： \(h\) 減半精度翻倍但計算量翻倍， \(O(h)\) 階是便宜但笨重的起點。

結案報告#

1768年尤拉在《積分學原理》中系統發表該法，微分追兇案形式上告破：任何 \(y' = f(t,y)\) 只要 \(f\) 夠光滑，都能用加減乘除逐步逼近。

此案遺產：

• 直接遺產：彈道、天體攝動首次能「算」而非「猜」；後世 Runge（1895）、Kutta（1901）把一階提升到四階。
• 警示遺產：穩定性概念埋下剛性方程、 CFL 條件、辛結構等整條支線。
• 哲學遺產：證明「不知道公式也能預言未來」，計算模擬從此與解析數學分庭抗禮。
• 未竟線索：一階精度太慢，實務上很快被高階法取代；但教學與理論分析至今仍以尤拉法為第一塊試金石。

尤拉沒見過電腦，卻為電腦寫好了第一行迴圈。

證據與工具#

• 核心格式： \(y_{n+1} = y_n + h f(t_n,y_n)\) ，由 \(y' = f(t,y)\) 導出。
• 誤差口供：局部 \(O(h^2)\) ，整體 \(O(h)\) ，一階收斂。
• 穩定性判據：對 \(y' = \lambda y\) 要求 \(|1 + h\lambda| \le 1\) 。
• 虛擬碼：輸入 \(y_0\) 、 \(h\) 、步數 \(N\) ；迴圈執行 \(y = y + h f(t,y)\) ， \(t = t + h\) 。
• 可重演實驗：解 \(y' = -y\) ， \(y_0 = 1\) ，分別取 \(h = 0.5\) 與 \(h = 0.1\) ，比較終點與真解 \(e^{-t}\) 的差距。
• 延伸閱讀：《Institutiones Calculi Integralis》卷一；後續案件直通傅立葉與 Navier-Stokes，因為偏微分方程正是常微分的無限維親戚。

補充：程式實作#

對應程式：1768-euler_method.py

本文核心理論是顯式尤拉格式：沿切線斜率走一小步，再以新位置斜率接續推進。 初值問題 \(y' = y\) 搭配 \(y(0) = 1\) 的真解為 \(y = e^t\) ，是檢驗精度的標準答案。 格式寫成 \(y_{n+1} = y_n + h y_n\) ，局部截斷誤差為 \(O(h^2)\) ，整體誤差為 \(O(h)\) 。 步長縮小十倍，終點誤差亦縮小約十倍，正是一階收斂的指紋特徵。 真兇 \(e\) 在終點 \(t = 1\) 露了馬腳，步長每砍十倍、誤差就乖乖縮水近十倍。

執行指令：

python3 _code/1768-euler_method.py

實測關鍵輸出（本機真實執行結果抄錄）：

h=0.1   steps=   10 y1=2.59374246 exact=2.71828183 err=1.245394e-01
h=0.01  steps=  100 y1=2.70481383 exact=2.71828183 err=1.346800e-02
h=0.001 steps= 1000 y1=2.71692393 exact=2.71828183 err=1.357896e-03
err ratio (0.1/0.01)   = 9.2471 (理論 $O(h)$ 比值約 $10$ )
err ratio (0.01/0.001) = 9.9183 (理論 $O(h)$ 比值約 $10$ )
VERIFY r1 = 9.2471（落在 $8$ 至 $12$ 區間，通過）
VERIFY r2 = 9.9183（落在 $8$ 至 $12$ 區間，通過）
PASS

讀者可改步長 \(h = 0.2\) 或換解 \(y' = -10 y\) ，觀察誤差比是否仍接近 \(10\) ，以及大步長何時觸發振盪發散。

上述終點值與誤差比均為程式逐步迭代所得，未經解析公式代入，呼應原文離散化與一致性等於收斂的主張。

完整程式如下：

# -*- coding: utf-8 -*-
"""1768 Euler 法（對應 wiki：計算模擬學 / Euler 顯式法解 ODE）
解 y'=y, y(0)=1 到 t=1（真解 y=e^t, y(1)=e），步長 h=0.1/0.01/0.001，
顯式 Euler：y_{k+1} = y_k + h·y_k。整體誤差為 O(h)，故 h 縮小 10 倍、
誤差亦縮小約 10 倍。驗證相鄰誤差比約 10 倍（8~12 內）。
只用 numpy，固定種子。
"""
import numpy as np

np.random.seed(0)

exact = np.e  # y(1) = e
hs = [0.1, 0.01, 0.001]
errs = []
for h in hs:
    n = int(round(1.0 / h))
    y = 1.0
    for _ in range(n):
        y = y + h * y
    err = abs(y - exact)
    errs.append(err)
    print(f"h={h:<6g} steps={n:>5d} y1={y:.8f} exact={exact:.8f} err={err:.6e}")

r1 = errs[0] / errs[1]
r2 = errs[1] / errs[2]
print(f"err ratio (0.1/0.01)   = {r1:.4f} (理論 O(h) 比值≈10)")
print(f"err ratio (0.01/0.001) = {r2:.4f} (理論 O(h) 比值≈10)")

# 驗證數字：理論值 vs 實測
ok1 = 8.0 < r1 < 12.0
ok2 = 8.0 < r2 < 12.0
print(f"VERIFY r1={r1:.4f} (8<r1<12? {ok1}), VERIFY r2={r2:.4f} (8<r2<12? {ok2})")
assert ok1, f"誤差比 r1 偏離 O(h) 預期: {r1}"
assert ok2, f"誤差比 r2 偏離 O(h) 預期: {r2}"
print("PASS")