108 年 國立成功大學環境工程學系碩士班乙組《工程數學》

1. 齊次解（Homogeneous Solution）

對應之齊次微分方程式為：
$y''' + 2y'' + 4y' + 8y = 0$

寫出特徵方程式（Characteristic Equation）：
$m^3 + 2m^2 + 4m + 8 = 0$

因式分解：
$m^2(m + 2) + 4(m + 2) = 0 \implies (m + 2)(m^2 + 4) = 0$

求得特徵根：
$m_1 = -2, \quad m_{2,3} = \pm 2i$

故齊次解為：
$y_h(x) = c_1 e^{-2x} + c_2 \cos(2x) + c_3 \sin(2x)$
（其中 $c_1, c_2, c_3$ 為任意常數）

---

2. 特解（Particular Solution）

非齊次項為 $r(x) = 6xe^{-2x}$，採用逆運算子法（Inverse Operator Method）：
$y_p(x) = \frac{1}{D^3 + 2D^2 + 4D + 8} (6xe^{-2x}) = \frac{1}{(D+2)(D^2+4)} (6xe^{-2x})$

應用微分位移定理（Exponential Shift Theorem）：
$y_p(x) = e^{-2x} \frac{1}{[(D-2)+2][(D-2)^2 + 4]} (6x) = e^{-2x} \frac{1}{D(D^2 - 4D + 8)} (6x)$

核心觀念

本題是二階線性 Euler–Cauchy 微分方程：

$$x^2y''-2xy'+2y=x^4e^x.$$

齊次方程可利用試誤解 $y=x^m$，因為各項的次方一致。非齊次部分則可利用齊次解中的因子 $x$ 作代換，降低計算難度。

需注意 $x=0$ 為方程的奇異點，因此一般在 $x>0$ 或 $x<0$ 的區間上討論解。

---

解題方法一：先求齊次解

先考慮齊次方程

$$x^2y''-2xy'+2y=0.$$

令

$$y=x^m,$$

則

$$y'=mx^{m-1},\qquad y''=m(m-1)x^{m-2}.$$

代入得

$$x^2m(m-1)x^{m-2}-2x(mx^{m-1})+2x^m=0,$$

整理為

$$\bigl[m(m-1)-2m+2\bigr]x^m=0.$$

因此特徵方程為

$$m(m-1)-2m+2=0,$$

即

$$m^2-3m+2=0 =(m-1)(m-2)=0.$$

所以

$$m=1,\qquad m=2.$$

齊次解為

$$y_h=C_1x+C_2x^2.$$

---

求非齊次特解

觀察齊次解含有因子 $x$，令

$$y=xv(x).$$

則

$$y'=v+xv',$$ $$y''=2v'+xv''.$$

代回原方程左側：

\begin{align*}
x^2y''-2xy'+2y
&=x^2(2v'+xv'')-2x(v+xv')+2xv\
&=2x^2v'+x^3v''-2xv-2x^2v'+2xv\
&=x^3v''.
\end{align*}

因此原方程化為

$$x^3v''=x^4e^x,$$

所以

$$v''=xe^x.$$

第一次積分：

$$v'=\int xe^x\,dx.$$

利用分部積分，

$$\int xe^x\,dx =xe^x-\int e^x\,dx =e^x(x-1).$$

故

$$v'=e^x(x-1).$$

再次積分：

利用冪級數解法（Power Series Method），設通解為 $y(t) = \sum_{n=0}^{\infty} a_n t^n$。

微分後代入原微分方程 $y'' - 4ty' + 4y = 0$：
$\sum_{n=2}^{\infty} n(n-1)a_n t^{n-2} - 4t \sum_{n=1}^{\infty} n a_n t^{n-1} + 4 \sum_{n=0}^{\infty} a_n t^n = 0$

平移指標並合併次項：
$\sum_{n=0}^{\infty} \left[ (n+2)(n+1)a_{n+2} - 4(n-1)a_n \right] t^n = 0$

令各次項係數為零，得到遞迴關係式（Recurrence Relation）：

觀念重點

本題屬於可透過變數變換化簡的微分方程式。觀察到方程式中含有 $y \frac{dy}{dx}$ 與 $y^2$ 項，令 $u = y^2$ 即可將原式轉化為一階線性微分方程式，再以積分因子法（Integrating Factor Method）求解。

---

詳細推導

步驟一：變數代換
令 $u = y^2$，對 $x$ 微分得：
$\frac{du}{dx} = 2y \frac{dy}{dx} \implies y \frac{dy}{dx} = \frac{1}{2} \frac{du}{dx}$

將 $u$ 與 $y \frac{dy}{dx}$ 代回原微分方程式 $y \frac{dy}{dx} = x^2 + y^2 - x$：
$\frac{1}{2} \frac{du}{dx} = x^2 + u - x$

移項整理得一階線性微分方程式標準式：
$\frac{du}{dx} - 2u = 2x^2 - 2x$

---

本題缺少 $y=1$ 與 $x=0$ 處之邊界條件，在此假設 $u(x,1)=0$ 及 $u(0,y)=f(y)$ 進行求解。

解題步驟

本題為半無限帶狀區域 $0 < x < \infty, 0 < y < 1$ 之拉普拉斯方程式（Laplace's Equation）邊界值問題，採用分離變數法（Separation of Variables）求解。

一、 分離變數

設 $u(x,y) = X(x)Y(y)$，代入偏微分方程式：
$X''(x)Y(y) + X(x)Y''(y) = 0 \implies \frac{X''(x)}{X(x)} = -\frac{Y''(y)}{Y(y)} = \lambda$

二、 邊界條件分析與求解

1. 對 $x$ 方向：
   由條件 $\lim_{x \to \infty} u(x,y) = 0$，解在無窮遠處衰減，故分離常數須為正數 $\lambda = k^2 > 0$（$k > 0$）。
   $X''(x) - k^2 X(x) = 0 \implies X(x) = c_1 e^{-kx} + c_2 e^{kx}$
   由 $\lim_{x \to \infty} X(x) = 0$ 得 $c_2 = 0$，故 $X(x) = c_1 e^{-kx}$。

2. 對 $y$ 方向：
   $Y''(y) + k^2 Y(y) = 0 \implies Y(y) = A \cos(ky) + B \sin(ky)$
   微商 $Y'(y) = -Ak \sin(ky) + Bk \cos(ky)$。

對變數 $y$（$0 < y < \infty$）取傅立葉餘弦轉換（Fourier Cosine Transform）：
$U(x, \omega) = \mathcal{F}_c \{u(x, y)\} = \int_0^\infty u(x, y) \cos(\omega y) \, dy$

利用邊界條件 $\left.\frac{\partial u}{\partial y}\right|_{y=0} = 0$，將原 PDE 轉換為 $x$ 的二階常微分方程式：
$\mathcal{F}_c \left\{\frac{\partial^2 u}{\partial y^2}\right\} = -\omega^2 U(x, \omega)$
$\frac{d^2 U(x, \omega)}{dx^2} - \omega^2 U(x, \omega) = 0$

將 $x$ 方向的邊界條件進行轉換：

1. $u(0, y) = 0 \implies U(0, \omega) = 0$
2. $u(\pi, y) = e^{-y} \implies U(\pi, \omega) = \int_0^\infty e^{-y} \cos(\omega y) \, dy = \frac{1}{1 + \omega^2}$

求解 ODE 通解：
$U(x, \omega) = A(\omega) \cosh(\omega x) + B(\omega) \sinh(\omega x)$

本題為非齊次邊界條件與非齊次源項（Source term）之偏微分方程式，採用穩態解（Steady-state solution）與暫態解（Transient solution）拆解法求解。

令 $u(x,t) = v(x) + w(x,t)$，其中 $v(x)$ 為不隨時間改變的穩態解，$w(x,t)$ 為滿足齊次邊界條件的暫態解。

---

步驟一：求穩態解 $v(x)$

穩態時 $\frac{\partial v}{\partial t} = 0$，微分方程式轉化為常微分方程式：
$v''(x) + \sin x = 0$

邊界條件為：
$v(0) = 400, \quad v(\pi) = 200$

對 $v''(x) = -\sin x$ 連續積分兩次：
$v'(x) = \cos x + C_1$
$v(x) = \sin x + C_1 x + C_2$

帶入邊界條件求解常數 $C_1, C_2$：
$v(0) = \sin(0) + C_1(0) + C_2 = 400 \implies C_2 = 400$
$v(\pi) = \sin(\pi) + C_1 \pi + 400 = 200 \implies C_1 = -\frac{200}{\pi}$

故穩態解為：
$v(x) = 400 - \frac{200}{\pi} x + \sin x$

---

步驟二：求暫態解 $w(x,t)$

將 $u(x,t) = v(x) + w(x,t)$ 代回原偏微分方程式，得到 $w(x,t)$ 的齊次熱傳導方程式：
$\frac{\partial w}{\partial t} = \frac{\partial^2 w}{\partial x^2}$

核心觀念

本題使用前向時間、中心空間差分法（FTCS）近似熱傳導方程式

$$\frac{\partial T}{\partial t} = \frac{\partial^2 T}{\partial x^2}.$$

差分格式為

$$\frac{T_i^{n+1}-T_i^n}{\Delta t} = \frac{T_{i+1}^n-2T_i^n+T_{i-1}^n}{\Delta x^2}.$$

需要分析：

1. 將精確解代入差分格式，求出截斷誤差。
2. 對數值誤差進行傅立葉模態分析，求出穩定條件。

令

$$r=\frac{\Delta t}{\Delta x^2}.$$

---

（A）截斷誤差

1. 時間方向泰勒展開

在 $(x_i,t_n)$ 對精確解 $T(x,t)$ 作泰勒展開：

$$T(x_i,t_n+\Delta t) = T + \Delta t\,T_t + \frac{\Delta t^2}{2}T_{tt} + \frac{\Delta t^3}{6}T_{ttt} + O(\Delta t^4).$$

因此

$$\frac{T_i^{n+1}-T_i^n}{\Delta t} = T_t + \frac{\Delta t}{2}T_{tt} + \frac{\Delta t^2}{6}T_{ttt} + O(\Delta t^3).$$

---

2. 空間方向泰勒展開

分別對 $T(x_i+\Delta x,t_n)$ 與 $T(x_i-\Delta x,t_n)$ 展開：

$$T_{i+1}^n = T + \Delta x\,T_x + \frac{\Delta x^2}{2}T_{xx} + \frac{\Delta x^3}{6}T_{xxx} + \frac{\Delta x^4}{24}T_{xxxx} + O(\Delta x^5),$$ $$T_{i-1}^n = T - \Delta x\,T_x + \frac{\Delta x^2}{2}T_{xx} - \frac{\Delta x^3}{6}T_{xxx} + \frac{\Delta x^4}{24}T_{xxxx} + O(\Delta x^5).$$

相加後，奇數階導數抵消：

$$T_{i+1}^n-2T_i^n+T_{i-1}^n = \Delta x^2 T_{xx} + \frac{\Delta x^4}{12}T_{xxxx} + O(\Delta x^6).$$

所以

$$\frac{T_{i+1}^n-2T_i^n+T_{i-1}^n}{\Delta x^2} = T_{xx} + \frac{\Delta x^2}{12}T_{xxxx} + O(\Delta x^4).$$

---

3. 求差分格式的殘差

定義局部截斷誤差為

$$\tau_i^n = \frac{T_i^{n+1}-T_i^n}{\Delta t} - \frac{T_{i+1}^n-2T_i^n+T_{i-1}^n}{\Delta x^2}.$$

代入上述展開式：

$$\tau_i^n = \left( T_t+\frac{\Delta t}{2}T_{tt}+O(\Delta t^2) \right) - \left( T_{xx}+\frac{\Delta x^2}{12}T_{xxxx}+O(\Delta x^4) \right).$$

由原微分方程式 $T_t=T_{xx}$，得到

$$\boxed{ \tau_i^n = \frac{\Delta t}{2}T_{tt} - \frac{\Delta x^2}{12}T_{xxxx} + O(\Delta t^2+\Delta x^4) }$$

因此此差分格式的截斷誤差階數為

$$\boxed{\tau_i^n=O(\Delta t+\Delta x^2)}.$$

由於方程式滿足 $T_t=T_{xx}$，再對時間微分可得

$$T_{tt} = \frac{\partial}{\partial t}(T_{xx}) = \frac{\partial^2}{\partial x^2}(T_t) = T_{xxxx}.$$

因此也可寫成

$$\boxed{ \tau_i^n = \left( \frac{\Delta t}{2} - \frac{\Delta x^2}{12} \right)T_{xxxx} + O(\Delta t^2+\Delta x^4) }.$$

其中時間為一階精度，空間為二階精度。

---

（B）穩定條件

1. 誤差方程式

令數值解與精確解之間的誤差為

$$e_i^n=T_i^n-\widetilde{T}_i^n.$$

因為差分格式為線性格式，誤差滿足相同的差分關係：

$$\frac{e_i^{n+1}-e_i^n}{\Delta t} = \frac{e_{i+1}^n-2e_i^n+e_{i-1}^n}{\Delta x^2}.$$

一、 臭氧濃度與假一級速率常數計算

假設大氣環境為標準條件（$25^\circ\text{C}$, $1\text{ atm}$），空氣分子數密度為：
$N_{\text{air}} = \frac{P N_A}{R T} \approx 2.46 \times 10^{19} \text{ molecules/cm}^3$

臭氧分子數濃度 $[\text{O}_3]$ 為：
$[\text{O}_3] = 100 \text{ ppbv} = 100 \times 10^{-9} \times 2.46 \times 10^{19} = 2.46 \times 10^{12} \text{ molecules/cm}^3$

因背景臭氧濃度保持不變，反應視為假一級（pseudo-first order）化學反應，假一級速率常數 $k'$ 為：
$k' = k [\text{O}_3] = (2 \times 10^{-14} \text{ cm}^3/\text{molecule}\cdot\text{sec}) \times (2.46 \times 10^{12} \text{ molecules/cm}^3) \approx 0.0492 \text{ sec}^{-1}$

---

二、 一級反應動力學與濃度比值推導

設總排放量為 $[\text{NO}_x]$，初始條件（$t = 0$）為：