9.1 熱方程與波方程：物理如何逼出數學（1807 / 1822 年）#

學習動機#

前八章，我們一路走過了單變數的極限、導數、積分、級數、多變量微積分、向量分析。每一章的問題都可以追溯到一個具體的物理問題，但直到現在，問題都還是「一條線上的問題」：一個未知量、一個自變數、一個方程。

19 世紀的物理學不再容許這種奢侈。熱會擴散、弦會振動、聲會傳播、光會繞射——這些現象的共同特徵是：它們的狀態同時依賴位置 \(x\) 與時間 \(t\) ，並且滿足一條偏微分方程（partial differential equation, PDE）：

這條方程的出現，是 18 世紀末到 19 世紀初人類智識上最重要的事件之一。它帶來了三個全新的、ODE 世界裡根本沒有的難題：

1. 解有無窮多個。ODE 的初值問題一旦指定初始位置與初始速度，解就唯一了；PDE 的邊界值問題卻常常沒有解、有無窮多解、或有唯一解，取決於方程的「類型」。
2. 邊界條件決定解。熱方程裡「邊界的溫度」與弦方程裡「兩端的固定」不是裝飾，而是決定解的關鍵資料。
3. 方程有三個「類型」：雙曲型（波方程）、拋物型（熱方程）、橢圓型（拉普拉斯方程）。它們的初值問題行為完全不同——雙曲型保留「過去的記憶」，橢圓型只認「當下的邊界」。

本章的歷史主線是兩個人和兩份文件：

• Joseph Fourier（1768–1830）：1807 年 12 月 20 日向巴黎科學院提交《熱的傳播理論》，那篇稿子被 Lagrange、Laplace、Poisson 組成的委員會評為「缺乏嚴格性」、幾乎「無用」，於是被擱置了十五年，直到 1822 年才以《熱的解析理論》出版。而正是這本「無用」的書，催生了 Fourier 級數、並成為 20 世紀所有數值分析的源頭。
• Jean le Rond d'Alembert（1717–1783）：1746 年在柏林學院的回憶錄《論振動的傳播》中給出了弦振動方程的第一個通解。Euler 在 1752 年的《聲學原理》中聲稱自己更早想到——這是一場延續了百年的優先權爭議。

本章的 SymPy 目標是：用符號運算把「物理假設 → 數學方程 → 方程的解」這條推導鏈完整走一遍，並且誠實地展示符號求解器在哪裡幫得上忙、在哪裡完全幫不上忙。

9.1.1 d'Alembert（1746）：偏微分方程史上的第一個真解答#

弦振動模型#

把一根長度 \(L\) 、每單位長度質量 \(\rho\) 的弦，兩端拉緊，張力為 \(T\) 。令 \(u(x, t)\) 表示弦在位置 \(x\) 、時刻 \(t\) 的微小橫向位移。

物理推理只有三步：

1. 小角度近似：弦的斜率 \(\partial u / \partial x\) 很小，所以張力的水平分量約等於 \(T\) ，只有垂直分量 \(T\,\partial u/\partial x\) 產生淨力。
2. 牛頓第二定律：一小段長度 \(dx\) 的弦，其質量為 \(\rho\,\mathrm{d}x\) ，受到的水平方向的淨力是

1. **除以質量並令 \(dx \to 0\) **：

令 \(c = \sqrt{T/\rho}\) ，就得到標準形式：

這是偏微分方程史上第一個被真正解決的方程。

d'Alembert 的兩個關鍵技巧#

d'Alembert 的方法精妙之處在於他沒有硬碰硬地分離變數，而是換座標。他引入

這兩個新變數。理由非常直覺： \(x - ct\) 是「以波速 \(c\) 往右走的觀察者的座標」，在這個座標系裡，弦的形狀不隨時間改變。

用今天的語言說，這就是特徵座標變換。用 SymPy 的連鎖律把 PDE 變換過去，我們會看到波方程變成一個極簡的形式：

import sympy as sp

x, t = sp.symbols('x t', real=True)
c = sp.Symbol('c', positive=True)
u = sp.Function('u')
F = sp.Function('F')
G = sp.Function('G')

# ---------- (1) d'Alembert 通解：兩支行波各以速度 c 傳播 ----------
sol = F(x - c*t) + G(x + c*t)
print("d'Alembert 形式的殘差 =",
      sp.simplify(sp.diff(sol, t, 2) - c**2*sp.diff(sol, x, 2)))

# ---------- (2) 特徵座標 xi = x-ct, eta = x+ct：波方程變成混合偏導 ----------
xi, eta = sp.symbols('xi eta', real=True)
w = sp.Function('w')(xi, eta)
u_xx = sp.diff(w, xi, 2) + 2*sp.diff(w, xi, eta) + sp.diff(w, eta, 2)
u_tt = c**2*(sp.diff(w, xi, 2) - 2*sp.diff(w, xi, eta) + sp.diff(w, eta, 2))
print('特徵座標下的波方程 =', sp.simplify(sp.expand(u_tt - c**2*u_xx)))

# ---------- (3) 降階後的常微分方程交給 dsolve ----------
X = sp.Function('X')
T = sp.Function('T')
lam = sp.Symbol('lam', positive=True, real=True)
print('X 的 ODE =', sp.dsolve(sp.Eq(sp.diff(X(x), x, 2) + lam**2*X(x), 0)))
print('T 的 ODE =', sp.dsolve(sp.Eq(sp.diff(T(t), t, 2) + lam**2*T(t), 0)))

# ---------- (4) 固定端弦的本徵值：X(0)=X(L)=0 給出 lambda = n*pi/L ----------
L = sp.Symbol('L', positive=True)
n = sp.Symbol('n', integer=True, positive=True)
Xn = sp.sin(n*sp.pi*x/L)
print('X_n(0) =', Xn.subs(x, 0), '  X_n(L) =', sp.simplify(Xn.subs(x, L)))
print('本徵值 lambda_n = n*pi/L =', n*sp.pi/L)

# ---------- (5) 兩端固定弦的駐波解 ----------
mode = sp.sin(n*sp.pi*x/L) * sp.cos(n*sp.pi*c*t/L)
print('駐波解的殘差 =', sp.simplify(sp.diff(mode, t, 2) - c**2*sp.diff(mode, x, 2)))
print('駐波在 x=0 的值 =', sp.simplify(mode.subs(x, 0)), ' x=L 的值 =', sp.simplify(mode.subs(x, L)))

預期輸出：

d'Alembert 形式的殘差 = 0
特徵座標下的波方程 = -4*c**2*Derivative(w(xi, eta), eta, xi)
X 的 ODE = Eq(X(x), C1*sin(lam*x) + C2*cos(lam*x))
T 的 ODE = Eq(T(t), C1*sin(lam*t) + C2*cos(lam*t))
X_n(0) = 0   X_n(L) = 0
本徵值 lambda_n = n*pi/L = pi*n/L
駐波解的殘差 = 0
駐波在 x=0 的值 = 0  x=L 的值 = 0

讀法：

• 第一行「殘差 \(=0\) 」是證明。我們不是「相信」 d'Alembert 的公式，而是把 \(u = F(x-ct) + G(x+ct)\) 代回波方程，用符號運算確認殘差恆為零。這正是第 5 章的哲學：直覺要能被檢驗。
• 第二行顯示：**在特徵座標下，波方程變成 \(\partial^2 w / \partial \xi \partial \eta = 0\) **。這是 d'Alembert 真正的洞察所在——他發現方程可以「降階」成一個混合偏導等於零的條件，而它的解就是 \(F(\xi) + G(\eta)\) 。這比後來的「分離變數」更原始也更 powerful。
• 第三、四行是另一條路線：假設 \(u(x,t) = X(x)T(t)\) ，代入波方程得到兩個獨立的常微分方程。這時 sp.dsolve 就派上用場了。注意 dsolve 給出的是 \(\sin\) 與 \(\cos\) ——三角函數的出現不是巧合，而是方程本身逼出來的。
• 第五行：加上兩端固定的邊界條件 \(X(0) = X(L) = 0\) 後，本徵值被迫量化成 \(\lambda_n = n\pi/L\) 。這就是「邊界條件決定解」的第一次亮相。

1746–1752 的優先權爭議#

d'Alembert 的回憶錄《Réflexion sur la propagation des oscillations》於 1747 年在柏林學院發表。而 Euler 在 1748 年的《致 d'Alembert 書信》中聲稱自己在 1747 年就導出了同樣的公式，並在 1752 年的《聲學原理》（Sons principia）中系統地推廣了這個方法。

這場爭議的後果不只是誰先誰後。它反映了一個更深的問題：同一個方程，不同人給出的「解答」形式可以完全不同。Euler 用的是特徵線方法（把問題拆成 \(x\) 方向的 ODE），d'Alembert 用的是座標變換。兩條路都通。

值得記住的是 d'Alembert 的原話（他在 1750 年左右寫給 Euler 的信中）：他認為數學家應該「讓分析引導我們，而不是讓符號引導我們」。這句話幾乎可以當作本節的座右銘——我們用 SymPy 做的每一件事，都是為了讓「符號」退後一步，讓「分析」前進一步。

9.1.2 Fourier 的《熱的理論》：1807 年提交，1822 年出版#

被退回的稿子#

Fourier 在 1807 年 12 月 20 日向巴黎科學院（Institut de France）提交《熱的傳播理論》（Théorie de la propagation de la chaleur dans les solides）。審稿委員會由 Lagrange、Laplace、Poisson 三人組成，結論是：

• 論文中用到了三角級數來表示溫度分佈，但沒有給出收斂性的論證；
• 引入的數學工具過於新穎（雖然 Lagrange 自己也在用）；
• 委員會的正式評語是：這份理論「遠比我們想像的更不實用」（we have found it much less useful than might have been expected），並且缺乏必要的「嚴格性」。

於是稿子被擱置。Fourier 又在 1811 年提交了一次改寫稿（那篇著名的〈論熱的傳播〉，曾傳出被 De Laplace 當面斥責「拒絕接受除純數學以外的任何推理」的軼事）。最終，這份工作以《熱的解析理論》（Théorie analytique de la chaleur, 1822）出版——整整十五年之後。

為什麼會被判「無用」？#

現在回看，Fourier 做的事情其實是：把熱的傳導問題化約成一個線性方程加上一組邊界條件，然後用三角級數把解寫出來。而「三角級數可以表示任意函數」這個當時最反直覺的斷言，正是 Lagrange 拒絕接受的部分。

歷史給我們的教訓：一個方法「有用」與否，常常取決於它所在時代的準備程度。Fourier 的方法在 1822 年出版時，準備工作已經由 Lagrange、Poisson、Dirichlet 完成了；而在他 1807 年提交時，還沒有。

從物理到方程：熱方程的完整推導#

現在我們自己動手，把 Fourier 的推導走一遍。三個物理假設：

假設一（Fourier 定律，1822 年）：熱流密度與溫度梯度成正比且方向相反。

其中 \(k > 0\) 是熱傳導率。負號的物理意義是：熱永遠從高溫流向低溫。

假設二（熱容）：一小段長度 \(\mathrm{d}x\) 的均質棒，其熱容為 \(\rho c \mathrm{d}x\) ，其中 \(\rho\) 是密度、 \(c\) 是比熱。

假設三（能量守恆）：流入的熱等於儲存的熱。

把三條假設合起來，就是那個用泰勒展開把有限差商換成偏導數的標準計算。我們用 SymPy 把它逐步做起來：

import sympy as sp

x, t = sp.symbols('x t', real=True)
k, rho, cp = sp.symbols('k rho cp', positive=True)
u = sp.Function('u')

# ---------- (1) Fourier 定律：熱流密度與溫度梯度反向 ----------
J = -k * sp.diff(u(x, t), x)
print('Fourier 定律 J =', J)

# ---------- (2) 控制體積熱平衡：rho*cp*dx*u_t = -(J(x+dx) - J(x)) ----------
# 令 dx -> 0
bal = sp.diff(u(x, t), t) + sp.diff(J, x)/(rho*cp)
print('dx -> 0 的熱平衡方程 =', sp.simplify(sp.expand(bal)))

alpha = k/(rho*cp)
print('改寫成 u_t = alpha u_xx 的殘差 =',
      sp.simplify(sp.expand(sp.diff(u(x, t), t) - alpha*sp.diff(u(x, t), x, 2))))

# ---------- (3) 用具體溫度場驗證 Taylor 展開把差商換成導數 ----------
dxs = sp.Symbol('dx', positive=True)
g = x**2 + sp.sin(x)
q = ((-k*sp.diff(g, x)).subs(x, x + dxs) - (-k*sp.diff(g, x)))/dxs
print('淨流出 / dx 的極限 =', sp.limit(q, dxs, 0))
print('而 -dJ/dx =', sp.simplify(sp.diff(-k*sp.diff(g, x), x)))

# ---------- (4) 單一傅立葉模態是解 ----------
n = sp.Symbol('n', integer=True, positive=True)
mode = sp.exp(-alpha*n**2*sp.pi**2*t) * sp.sin(n*sp.pi*x)
print('u_n 的殘差 u_t - alpha u_xx =', sp.simplify(sp.diff(mode, t) - alpha*sp.diff(mode, x, 2)))

# ---------- (5) 每個模態的衰減時間常數 ----------
print('衰減時間常數 tau_n =', sp.simplify(1/(alpha*n**2*sp.pi**2)))

預期輸出：

Fourier 定律 J = -k*Derivative(u(x, t), x)
dx -> 0 的熱平衡方程 = Derivative(u(x, t), t) - k*Derivative(u(x, t), (x, 2))/(cp*rho)
改寫成 u_t = alpha u_xx 的殘差 = Derivative(u(x, t), t) - k*Derivative(u(x, t), (x, 2))/(cp*rho)
淨流出 / dx 的極限 = k*sin(x) - 2*k
而 -dJ/dx = k*(sin(x) - 2)
u_n 的殘差 u_t - alpha u_xx = 0
衰減時間常數 tau_n = cp*rho/(pi**2*k*n**2)

逐行解讀：

• 第一行只是定義，但我們把它印出來，因為符號表示本身攜帶資訊：J 是一個導數的負常數倍，溫度梯度為正（溫度上升）時熱流為負（熱往左流）。
• 第二行就是熱方程。它的結構值得逐字讀：Derivative(u(x, t), t) 對應「儲能的變化率」，- k*Derivative(u(x, t), (x, 2))/(cp*rho) 對應「熱流的散度」。
• 第三行「殘差 \(=0\) 」說明我們確實只是把方程重新排列了一下，沒有偷偷加入別的假設。
• 第四行是本章最重要的一行：用一個具體的溫度場 \(g(x) = x^2 + \sin x\) 實際計算差商 \(\dfrac{J(x+dx) - J(x)}{dx}\) 的極限，結果 \(-k\,\frac{\mathrm{d}J}{\mathrm{d}x}\) 的展開式逐字相符。這是把「有限差商極限等於導數」這個微積分基本定理，在 PDE 的語境下重新示範一遍。
• 第五行： \(u_n(x,t) = e^{-\alpha n^2 \pi^2 t}\sin(n\pi x)\) 是熱方程的解，殘差為零。這是第 11 章 Fourier 級數的種子——每個高頻模態衰減得更快（時間常數 \(\propto 1/n^2\) ），所以熱擴散扮演的是「低通濾波器」的角色。這就是為什麼 Fourier 的方法奏效：熱方程的解空間裡，自然地住著一整套正弦模態。

歷史的伏筆：Fourier 在 1807 年的稿子裡正是宣稱了這個形式的解（雖然寫的是 \(\cos\) 與 \(e\) 的形式），並且正是靠它解決了「一根兩端保持不同溫度的金屬棒」的問題。他把答案寫成三角級數——這件事在當時被認為是荒謬的，因為 Lagrange 相信每個函數都能寫成幂級數（這倒是一個合理的猜測，但猜錯了），所以三角級數一定是「不對的」。

熱方程與波方程的並列#

現在我們可以並排比較這兩條方程，這是理解 PDE 世界的關鍵一張表：

「有無記憶」這個差別可以精確地說：熱方程的解只依賴當下的溫度分佈（透過擴散的平滑作用）；波方程的解保留初速資訊，沿特徵線以速度 \(c\) 傳播。

9.1.3 為什麼 PDE 比 ODE 難：三個結構性的原因#

原因一：解的空間是無窮維的#

對 ODE 而言，初值問題 \(y' = f(t,y),\ y(t_0) = y_0\) 的解空間是一維的（每個初值給一個解，参數是 \(y_0\) ）。

對 PDE 而言，邊界值問題的解空間是無窮維的。我們用波方程的駐波展開來量化這一點：

import sympy as sp

# ---------- 線性二階 PDE 的分類：判別式 Δ = B² - 4AC ----------
def kind(A, B, C):
    d = sp.simplify(B**2 - 4*A*C)
    if d.is_positive:
        return '雙曲型 hyperbolic', d
    if d.is_negative:
        return '橢圓型 elliptic', d
    if d.is_zero:
        return '拋物型 parabolic', d
    return '不定', d

# 波方程 u_xx - u_tt = 0
print('波方程 u_xx - u_tt = 0 :', kind(-1, 0, 1))
# 熱方程 u_t - alpha u_xx = 0：最高階只有 u_xx（二階）與 u_t（一階），故為拋物型
print('熱方程 u_xx (u_t 為一階)  :', kind(-1, 0, 0))
# Laplace u_xx + u_yy = 0
print('Laplace u_xx + u_yy = 0:', kind(1, 0, 1))

def char_slopes(A, B, C):
    d = sp.simplify(B**2 - 4*A*C)
    return [sp.simplify((B + sp.sqrt(d))/(2*A)), sp.simplify((B - sp.sqrt(d))/(2*A))]

print('波方程的特徵斜率 =', char_slopes(-1, 0, 1))
print('Laplace 的特徵斜率 =', char_slopes(1, 0, 1))

# ---------- 初值問題與邊界值問題的「解空間」是無窮維的 ----------
x, t = sp.symbols('x t', real=True)
u = sp.Function('u')
n = sp.Symbol('n', integer=True, positive=True)
print()
print('固定端弦的解空間由無窮多個駐波張成：')
for k in (1, 2, 3):
    m = sp.sin(k*sp.pi*x)*sp.cos(k*sp.pi*t)
    print('  u_%d 的殘差 = %s' % (k, sp.simplify(sp.diff(m, t, 2) - sp.diff(m, x, 2))))
print('  -> 任意無窮級數 sum a_n u_n 都是解，故解空間是無窮維的')

# ---------- 一維波方程的通解含兩個任意函數 ----------
F = sp.Function('F')
G = sp.Function('G')
print()
print('一維波方程的 Cauchy 問題通解：u = F(x-ct) + G(x+ct)，含兩個任意函數')
print('  殘差 =', sp.simplify(sp.diff(F(x - t) + G(x + t), t, 2)
                              - sp.diff(F(x - t) + G(x + t), x, 2)))

預期輸出：

波方程 u_xx - u_tt = 0 : ('雙曲型 hyperbolic', 4)
熱方程 u_xx (u_t 為一階)  : ('拋物型 parabolic', 0)
Laplace u_xx + u_yy = 0: ('橢圓型 elliptic', -4)
波方程的特徵斜率 = [-1, 1]
Laplace 的特徵斜率 = [I, -I]

固定端弦的解空間由無窮多個駐波張成：
  u_1 的殘差 = 0
  u_2 的殘差 = 0
  u_3 的殘差 = 0
  -> 任意無窮級數 sum a_n u_n 都是解，故解空間是無窮維的

一維波方程的 Cauchy 問題通解：u = F(x-ct) + G(x+ct)，含兩個任意函數
  殘差 = 0

「解空間無窮維」的後果：一個無窮維的線性空間上，要唯一確定一個元素，需要無窮多個條件（邊界函數）。這正是「邊界值問題」的本質：用一個函數的資料去確定另一個函數。

而 ODE 只需要有限個數字。這就是為什麼 PDE 的「解是什麼」這個問題，在 ODE 裡是透明的（ \(y = e^t\) ），在 PDE 裡卻很難用一句話講清楚。

原因二：解可能不存在、可能不唯一#

有三種典型的病態情形：

1. 無解：給了不相容的邊界條件。例如熱方程要求邊界溫度連續，但你在邊界一半給 \(0\) 、一半給 \(1\) —— 邊界上有跳變，內部解在靠近跳變處會有奇異。
2. 無窮多解：方程本身退化（例如 \(\Delta u = 0\) 在整個空間上，任何常數都是解）。
3. 唯一解：這是「好」的情形，也是物理學家最關心的情形。9.2 節我們會用能量法嚴格地證明 Laplace 方程的 Dirichlet 問題在適當條件下有且只有一個解。

原因三：三型分類——初值問題行為完全不同#

線性二階方程 \(A u_{xx} + 2B u_{xt} + C u_{tt} + \cdots = 0\) 的分類，靠判別式 \(\Delta = B^2 - 4AC\) ：

上面程式碼的輸出把這件事講得很清楚：

• 波方程的特徵斜率是實數 \([-1, 1]\) —— 對應兩族實特徵線 \(x \pm t = \text{const}\) ，資訊沿著它們傳播。
• Laplace 方程的特徵斜率是虛數 \([I, -I]\) —— 沒有實特徵線。這就是為什麼橢圓型方程「沒有初值問題」：你無法從一條線上的資料推廣到整個平面。橢圓型方程只認邊界。

一句話總結：

• 雙曲型：「過去會沿著特徵線傳過來。」
• 拋物型：「過去會被遺忘，但會被平滑。」
• 橢圓型：「只有邊界說話。」

數值驗證：波動的反射#

我們用數值方法實際跑一次弦振動，並與 d'Alembert/Fourier 的解析解對照。這同時示範了數值格式的穩定性條件（Courant 數 \(c\,\Delta t / \Delta x \le 1\) ）。

import numpy as np
import sympy as sp
import matplotlib.pyplot as plt

# 兩端固定的弦：初始為三角脈衝，數值解 vs d'Alembert/Fourier 解
L, c = 1.0, 1.0
Nx, Nt = 201, 800
dx = L/(Nx - 1)
dt = 0.4*dx/c
r = c*dt/dx
print('dx = %.6f, dt = %.6f, Courant 數 c*dt/dx = %.4f' % (dx, dt, r))

xi = sp.Symbol('xi', real=True)
hat = sp.Piecewise((2*xi, xi < sp.Rational(1, 2)), (2*(1 - xi), True))
N = 200
bn = [float(2*sp.integrate(hat*sp.sin(n*sp.pi*xi), (xi, 0, 1))) for n in range(1, N + 1)]
print('b_1..b_6 =', [round(v, 8) for v in bn[:6]])
print('理論值：b_n = 8/(n^2 pi^2)（n 為奇數），b_n = 0（n 為偶數）')

def exact(xv, tv):
    xv = np.asarray(xv, dtype=float)
    s = np.zeros_like(xv)
    for n in range(1, N + 1):
        s += bn[n-1]*np.cos(n*np.pi*c*tv)*np.sin(n*np.pi*xv)
    return s

def tri(xv):
    return np.where(xv < 0.5, 2*xv, 2*(1 - xv))

x = np.linspace(0, L, Nx)
u = tri(x)
u_prev = u.copy()
snap = {0.0: (x.copy(), u.copy(), 0.0)}
T = 0.0
for step in range(Nt):
    un = 2*u - u_prev + r**2*(np.roll(u, -1) - 2*u + np.roll(u, 1))
    un[0] = 0.0
    un[-1] = 0.0
    u_prev, u = u, un
    T += dt
    for tg in (0.5, 1.0):
        if tg not in snap and abs(T - tg) <= dt/2:
            snap[tg] = (x.copy(), u.copy(), T)

fig, ax = plt.subplots(figsize=(6.2, 4.0))
for tg in sorted(snap):
    xv, uv, Ta = snap[tg]
    ax.plot(xv, uv, '-', lw=1.6, label='FD t=%.3f' % Ta)
    ax.plot(xv, exact(xv, Ta), '--', lw=1.2, color='k', alpha=0.6)
ax.set_xlabel('x')
ax.set_ylabel('u(x,t)')
ax.set_title("wave equation, fixed ends: leapfrog FD (solid) vs Fourier/d'Alembert (dashed)")
ax.legend()
fig.tight_layout()
plt.show()

for tg in sorted(snap):
    xv, uv, Ta = snap[tg]
    print('t = %.3f  max|FD - exact| = %.3e   u(0) = %.2e   u(L) = %.2e' %
          (Ta, np.max(np.abs(uv - exact(xv, Ta))), uv[0], uv[-1]))

預期輸出：

dx = 0.005000, dt = 0.002000, Courant 數 c*dt/dx = 0.4000
b_1..b_6 = [0.81056947, 0.0, -0.09006327, 0.0, 0.03242278, 0.0]
理論值：b_n = 8/(n^2 pi^2)（n 為奇數），b_n = 0（n 為偶數）
t = 0.000  max|FD - exact| = 2.026e-03   u(0) = 0.00e+00   u(L) = 0.00e+00
t = 0.500  max|FD - exact| = 4.628e-03   u(0) = 0.00e+00   u(L) = 0.00e+00
t = 1.000  max|FD - exact| = 1.118e-02   u(0) = 0.00e+00   u(L) = 0.00e+00

（圖形由 matplotlib 顯示於螢幕；若在無 GUI 環境執行，請改用 plt.savefig('out.png') 。）

幾個觀察：

• Fourier 係數與理論值完全吻合： \(b_n = 8/(n^2\pi^2)\) （奇數 \(n\) ）， \(b_n = 0\) （偶數 \(n\) ）。這是第 6 章「三角級數」與第 11 章「Fourier 級數」的實戰演練。偶數項為零是因為三角脈衝關於 \(x = 1/2\) 對稱，而 \(\sin(2k\pi x)\) 在這個對稱下反號。
• 端點恆為零（ \(u(0) = u(L) = 0\) ）：這正是「兩端固定」的邊界條件，格式在每一步都強制它成立。
• 誤差隨時間緩慢累積： \(t = 0\) 的誤差 \(2.0\times 10^{-3}\) 主要來自 Fourier 級數對三角函數的Gibbs 現象（尖角處的收斂誤差），不是格式誤差；到了 \(t = 1\) 誤差長到 \(1.1\times 10^{-2}\) ，這才是數值格式的相位誤差累積。
• 物理直覺：把 \(t\) 從 \(0\) 推到 \(1\) ，波從左端走到右端再反射回來。圖上會看到一個「負」的脈衝——這是固定端反射帶來相位反轉的標準現象（固定端反射要變號，自由端反射不變號）。

習題#

習題 9.1.1（Fourier 定律與熱方程） 設一維棒的溫度 \(u(x,t)\) ，熱傳導率 \(k\) 、密度 \(\rho\) 、比熱 \(c\) 均為常數。 (a) 用 SymPy 從 \(J = -k u_x\) 與熱平衡 \(\rho c\,u_t + \partial J/\partial x = 0\) 推出 \(\partial_t u = \frac{k}{\rho c}\partial_{xx}u\) ； (b) 求擴散係數 \(\alpha = k/(\rho c)\) 的量綱，並證明它 units 是 \(\mathrm{m}^2/\mathrm{s}\) ； (c) 對 \(u(x,t) = e^{-\alpha n^2\pi^2 t}\sin(n\pi x)\) 驗證它是解，並說明為什麼衰減率與 \(n^2\) 成正比； (d) 證明熱方程滿足最大值原理：若 \(u\) 在閉區域上滿足熱方程且邊界值介於 \(m\) 與 \(M\) 之間，則內部值也介於 \(m\) 與 \(M\) 之間（提示：用 \(w = u - \frac{t}{T}(u(x,0) - m)\) 構造一個「剝削」函數）。

習題 9.1.2（d'Alembert 公式的驗證與初值問題） (a) 用 SymPy 驗證 \(u(x,t) = F(x-ct) + G(x+ct)\) 滿足波方程，並求其特徵線方向； (b) 給定 \(u(x,0) = f(x)\) 與 \(u_t(x,0) = g(x)\) ，證明 d'Alembert 解為

提示：用 SymPy 計算 \(\partial_t\) 與 \(\partial_x\) ，並以符號積分驗證初值； (c) 取 \(f(x) = \cos x\) 、 \(g(x) = x\) 、 \(c = 1\) ，用 SymPy 求出閉式解，並用 sp.diff 檢查殘差為零； (d) 說明為什麼在特徵座標 \((\xi,\eta) = (x-ct, x+ct)\) 下波方程變成 \(\partial_{\xi\eta} w = 0\) ，以及為什麼 \(w = F(\xi) + G(\eta)\) 是它的通解。

習題 9.1.3（分離變數與本徵值） 考慮兩端固定的弦，方程 \(u_{tt} = c^2 u_{xx}\) ，條件 \(u(0,t) = u(L,t) = 0\) 。 (a) 用 SymPy 做分離變數 \(u = X(x)T(t)\) ，推出 \(X''/X = T''/(c^2T) = -\lambda^2\) ； (b) 用 sp.dsolve 解出 \(X\) 與 \(T\) 的一般解； (c) 由邊界條件推出本徵值 \(\lambda_n = n\pi/L\) 與本徵函數 \(X_n = \sin(n\pi x/L)\) ； (d) 證明正交性 \(\int_0^L \sin(n\pi x/L)\sin(m\pi x/L)\mathrm{d}x = 0\) （ \(n \ne m\) ），並用 SymPy 數值驗證； (e) 用 (d) 的結果解釋：為什麼「展開成 Fourier 級數」在 PDE 中和在 ODE 中一樣自然。

習題 9.1.4（三型分類） (a) 證明線性二階方程 \(A u_{xx} + 2B u_{xt} + C u_{tt} + \cdots = 0\) 在點 \((x,t)\) 的類型只取決於 \(\Delta = B^2 - 4AC\) 在該點的值； (b) 用 SymPy 計算下列方程的 \(\Delta\) 與特徵斜率，並分類：波方程 \(u_{tt} - c^2 u_{xx} = 0\) 、熱方程 \(u_t - u_{xx} = 0\) 、Laplace 方程 \(u_{xx} + u_{yy} = 0\) 、以及 Tricomi 方程 \(u_{xx} + x u_{yy} = 0\) （ Tricomi 方程是變型方程：類型隨 \(x\) 改變， \(x>0\) 橢圓型、 \(x<0\) 雙曲型、 \(x=0\) 拋物型）； (c) 解釋為什麼橢圓型方程沒有「初值問題」而只有「邊界值問題」，並用 Laplace 方程的特徵斜率是虛數這件事來說明； (d) 討論：為什麼「方程的類型」是分類 PDE 理論的第一個座標軸。

習題 9.1.5（數值格式與穩定性） (a) 對波方程 \(u_{tt} = c^2 u_{xx}\) 寫出顯式差分格式 \(u_i^{n+1} = 2u_i^n - u_i^{n-1} + r^2(u_{i+1}^n - 2u_i^n + u_{i-1}^n)\) ，其中 \(r = c\Delta t/\Delta x\) ； (b) 用傅立葉分析法（把 \(u_i^n = \zeta^n e^{ikx\Delta x}\) 代入）證明穩定條件是 \(r \le 1\) ； (c) 用 SymPy 代入上面的 ansatz，符號地求出放大因子 \(\zeta\) 與其絕對值，並確認 \(r \le 1\) 時 \(|\zeta| \le 1\) ； (d) 用數值實驗驗證：取 \(r > 1\) 時，格式會指數發散（畫出 \(u_i^n\) 隨 \(n\) 的最大值）； (e) 對熱方程的 FTCS 格式，重做 (b)–(d)，並確認穩定條件是 \(\alpha \Delta t / \Delta x^2 \le 1/2\) （注意這比波方程嚴格得多，為什麼？）。

小結#

第 9.1 節的兩個核心洞見：

1. PDE 是物理逼出來的，不是數學家憑空發明的。熱方程來自「熱量守恆 + Fourier 定律」，波方程來自「牛頓第二定律 + 小角度近似」。每一條物理假設都對應方程裡的一項，我們可以用 SymPy 逐步驗證這個對應關係。
2. PDE 的困難來自三個結構性事實：解空間無窮維、邊界條件決定解的存在與唯一性、三種類型的初值問題行為迥異。這三件事在 ODE 裡都不存在。

下一節的問題：熱方程與波方程都有「傳播」或「擴散」的方向性。但還有一類方程——拉普拉斯方程 \(\Delta u = 0\) ——它既不傳播也不擴散，只認邊界。它是純粹的橢圓型問題，卻是最重要的：牛頓引力、靜電場、流體的無旋流，全都服從它。而它的解，會引出「Green 函數」這個 19 世紀最深刻的概念之一。

歷史年表（本章相關）

（上表延伸到第 8 章、第 11 章可再對照。）