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

核心觀念

本題涵蓋四類常見微分方程：

• A：Cauchy–Euler 型齊次方程。
• B：Cauchy–Euler 型非齊次方程，利用算子因式分解。
• C：常係數線性微分方程，且右端與重根共振。
• D：含積分項的聯立微分方程，先消去積分項，再求解初值問題。

---

A. $\displaystyle \frac{d^2y}{dx^2}+\frac{1}{x}\frac{dy}{dx}+\frac{y}{x^2}=0$

解題方法

這是 Cauchy–Euler 方程，設

$$y=x^m,\qquad x\ne 0.$$

則

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

代入原方程：

$$m(m-1)x^{m-2}+\frac{1}{x}mx^{m-1} +\frac{1}{x^2}x^m=0.$$

整理得

$$\left[m(m-1)+m+1\right]x^{m-2}=0,$$

因此特徵方程為

$$m^2+1=0.$$

根為

$$m=\pm i.$$

對於複數根 $m=\alpha\pm i\beta$，實數通解為

$$y=x^\alpha\left[C_1\cos(\beta\ln|x|) +C_2\sin(\beta\ln|x|)\right].$$

本題 $\alpha=0,\beta=1$，所以

$$\boxed{y=C_1\cos(\ln|x|)+C_2\sin(\ln|x|)}.$$

解題技巧

看到

$$x^2y''+axy'+by=0$$

即可優先嘗試 $y=x^m$。本題乘上 $x^2$ 後即為

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

定義域需取在不包含 $x=0$ 的區間上。

【答案】

$$\boxed{y=C_1\cos(\ln|x|)+C_2\sin(\ln|x|)}$$

---

B. $\displaystyle x^2y''-2xy'+2y=x^2e^x$

核心觀念

左側 Cauchy–Euler 算子可因式分解為

$$x^2D^2-2xD+2=(xD-1)(xD-2),$$

其中 $D=\dfrac{d}{dx}$。

解題方法

先求齊次解。令 $y=x^m$，得到

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

即

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

因此

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

接著處理非齊次項。令

$$z=(xD-2)y=xy'-2y.$$

原方程化為

$$(xD-1)z=x^2e^x,$$

也就是

$$xz'-z=x^2e^x.$$

除以 $x^2$：

$$\frac{z'}{x}-\frac{z}{x^2}=e^x.$$

左側正好是

$$\left(\frac{z}{x}\right)'.$$

因此

$$\left(\frac{z}{x}\right)'=e^x,$$

積分得

$$\frac{z}{x}=e^x+C,$$

所以

$$z=xe^x+Cx.$$

回到 $y$：

$$xy'-2y=xe^x+Cx.$$

除以 $x$：

$$y'-\frac{2}{x}y=e^x+C.$$

這是一階線性方程，其積分因子為

$$\mu(x)=e^{\int -2/x\,dx}=x^{-2}.$$

乘上 $x^{-2}$：

$$\left(\frac{y}{x^2}\right)' =\frac{e^x}{x^2}+\frac{C}{x^2}.$$

利用

$$\int \frac{e^x}{x^2}\,dx =\operatorname{Ei}(x)-\frac{e^x}{x},$$

以及

$$\int \frac{C}{x^2}\,dx=-\frac{C}{x},$$

可得

$$\frac{y}{x^2} =\operatorname{Ei}(x)-\frac{e^x}{x} -\frac{C}{x}+C_2.$$

乘上 $x^2$：

$$y=x^2\operatorname{Ei}(x)-xe^x-Cx+C_2x^2.$$

將常數重新命名後，通解為

$$\boxed{ y=C_1x+C_2x^2+x^2\operatorname{Ei}(x)-xe^x }.$$

其中 $\operatorname{Ei}(x)$ 為指數積分函數，滿足

$$\frac{d}{dx}\operatorname{Ei}(x)=\frac{e^x}{x}.$$

解題技巧

算子因式分解是本題關鍵：

$$(xD-1)(xD-2)y =x^2y''-2xy'+2y.$$

分解後連續解兩個一階方程，比直接使用變數變換或變參數法更有效率。

【答案】

核心觀念

二階 Runge–Kutta 法的目標，是使數值解的一步增量與精確解的 Taylor 展開相同至 $h^2$ 項。因此：

• 單步截斷誤差為 $O(h^3)$；
• 累積的全域截斷誤差為 $O(h^2)$；
• 推導關鍵在於比較數值公式與
  $y(x+h)=y(x)+hy'(x)+\frac{h^2}{2}y''(x)+\frac{h^3}{6}y'''(x)+O(h^4).$

以下以最常見的二階顯式 Runge–Kutta 法——Heun 法為例：

$$k_1=f(x_n,y_n),$$ $$k_2=f(x_n+h,y_n+hk_1),$$ $$y_{n+1}=y_n+\frac{h}{2}(k_1+k_2).$$

其中 $h$ 為步長。

---

解題方法

令

$$x=x_n,\qquad y=y(x_n),\qquad f=f(x,y).$$

為了推導局部截斷誤差，將數值方法中的 $y_n$ 代換為精確值 $y(x_n)$。

1. 精確解的 Taylor 展開

由微分方程

$$y'=f(x,y),$$

可得

$$y(x+h) = y+hf+\frac{h^2}{2}y''+\frac{h^3}{6}y'''+O(h^4).$$

首先，

$$y''=\frac{d}{dx}f(x,y) =f_x+f_y y' =f_x+f_yf.$$

再對 $y''$ 微分：

$$y''' = \frac{d}{dx}(f_x+f_yf).$$

分別計算：

$$\frac{d}{dx}f_x=f_{xx}+f_{xy}f,$$

以及

$$\frac{d}{dx}(f_yf) = (f_{xy}+f_{yy}f)f+f_y(f_x+f_yf).$$

因此

$$y''' = f_{xx}+2ff_{xy}+f^2f_{yy}+f_xf_y+ff_y^2.$$

故精確解為

$$\begin{aligned} y(x+h) ={}&y+hf+\frac{h^2}{2}(f_x+ff_y)\\ &+\frac{h^3}{6} \left( f_{xx}+2ff_{xy}+f^2f_{yy}+f_xf_y+ff_y^2 \right) +O(h^4). \end{aligned}$$

---

2. 展開 Runge–Kutta 中的 $k_2$

由

$$k_2=f(x+h,y+hf),$$

對 $f(x+h,y+hf)$ 作二變數 Taylor 展開：

$$\begin{aligned} k_2 ={}&f+h(f_x+ff_y)\\ &+\frac{h^2}{2} \left( f_{xx}+2ff_{xy}+f^2f_{yy} \right) +O(h^3). \end{aligned}$$

因此數值方法的一步結果為

$$y_{n+1}^{\mathrm{RK2}} = y+\frac{h}{2}(f+k_2).$$

代入 $k_2$：

$$\begin{aligned} y_{n+1}^{\mathrm{RK2}} ={}&y+hf+\frac{h^2}{2}(f_x+ff_y)\\ &+\frac{h^3}{4} \left( f_{xx}+2ff_{xy}+f^2f_{yy} \right) +O(h^4). \end{aligned}$$

---

3. 比較精確解與數值解

定義局部截斷誤差為

$$\tau_{n+1} = y(x+h)-y_{n+1}^{\mathrm{RK2}}.$$

將兩個展開式相減，$h$ 與 $h^2$ 項完全抵消，得到

$$\begin{aligned} \tau_{n+1} ={}&h^3 \bigg[ \frac{1}{6} \left( f_{xx}+2ff_{xy}+f^2f_{yy}+f_xf_y+ff_y^2 \right)\\ &\qquad\qquad -\frac{1}{4} \left( f_{xx}+2ff_{xy}+f^2f_{yy} \right) \bigg] +O(h^4). \end{aligned}$$

整理後：

核心觀念

本題考查有限差分法（finite difference method）將二階常微分方程轉換為線性代數方程組。

微分方程為

$$y''+p(x)y'+q(x)y=r(x),$$

邊界條件為

$$y(0)=2,\qquad y'(1)=0.$$

網格點間距為

$$h=\frac{1-0}{4}=0.25.$$

各節點定義為

$$x_0=0,\quad x_1=0.25,\quad x_2=0.5,\quad x_3=0.75,\quad x_4=1,$$

其中

$$y_0=y(0)=2,$$

且未知向量為

$$Y= \begin{bmatrix} y_1\\y_2\\y_3\\y_4 \end{bmatrix}.$$

為簡化記號，令

$$p_i=p(x_i),\qquad q_i=q(x_i),\qquad r_i=r(x_i).$$

---

解題方法

1. 內部節點的中央差分

對於 $x_i$，二階導數與一階導數分別使用中央差分：

$$y''(x_i)\approx \frac{y_{i-1}-2y_i+y_{i+1}}{h^2},$$ $$y'(x_i)\approx \frac{y_{i+1}-y_{i-1}}{2h}.$$

代入原微分方程：

$$\frac{y_{i-1}-2y_i+y_{i+1}}{h^2} +p_i\frac{y_{i+1}-y_{i-1}}{2h} +q_i y_i=r_i.$$

整理後得到

$$\left(\frac{1}{h^2}-\frac{p_i}{2h}\right)y_{i-1} +\left(-\frac{2}{h^2}+q_i\right)y_i +\left(\frac{1}{h^2}+\frac{p_i}{2h}\right)y_{i+1} =r_i.$$

此式套用於內部節點 $x_1,x_2,x_3$。

---

2. 在 $x_1=0.25$ 的差分方程

$$\left(\frac{1}{h^2}-\frac{p_1}{2h}\right)y_0 +\left(-\frac{2}{h^2}+q_1\right)y_1 +\left(\frac{1}{h^2}+\frac{p_1}{2h}\right)y_2 =r_1.$$

由於 $y_0=2$，移至右側：

$$\left(-\frac{2}{h^2}+q_1\right)y_1 +\left(\frac{1}{h^2}+\frac{p_1}{2h}\right)y_2 = r_1-2\left(\frac{1}{h^2}-\frac{p_1}{2h}\right).$$

---

3. 在 $x_2=0.5$ 的差分方程

$$\left(\frac{1}{h^2}-\frac{p_2}{2h}\right)y_1 +\left(-\frac{2}{h^2}+q_2\right)y_2 +\left(\frac{1}{h^2}+\frac{p_2}{2h}\right)y_3 =r_2.$$

---

4. 在 $x_3=0.75$ 的差分方程

$$\left(\frac{1}{h^2}-\frac{p_3}{2h}\right)y_2 +\left(-\frac{2}{h^2}+q_3\right)y_3 +\left(\frac{1}{h^2}+\frac{p_3}{2h}\right)y_4 =r_3.$$

---

5. 在 $x_4=1$ 的 Neumann 邊界條件

在右端點使用後向差分：

核心觀念

本題為二維拉普拉斯方程

$$u_{xx}+u_{yy}=0.$$

主要使用分離變數法。設

$$u(x,y)=X(x)Y(y),$$

代入方程式得

$$X''(x)Y(y)+X(x)Y''(y)=0.$$

除以 $X(x)Y(y)$：

$$\frac{X''}{X}=-\frac{Y''}{Y}.$$

令分離常數為 $\lambda^2$，則

$$X''-\lambda^2X=0,\qquad Y''+\lambda^2Y=0.$$

依據邊界條件選擇適當的特徵函數，再利用傅立葉級數展開非齊次邊界值。

---

A. 半無限長矩形區域

條件為

$$u(0,y)=1,\qquad \lim_{x\to\infty}u(x,y)=0,\qquad 0<y<1,$$

以及

$$u_y(x,0)=0,$$ $$u_y(x,1)+u(x,1)=0.$$

解題方法

由於 $x\to\infty$ 時 $u\to 0$，$x$ 方向必須選擇衰減解。令

$$Y''+\lambda^2Y=0,$$

則

$$X''-\lambda^2X=0.$$

其一般解為

$$X(x)=C_1e^{\lambda x}+C_2e^{-\lambda x},$$ $$Y(y)=A\cos(\lambda y)+B\sin(\lambda y).$$

由 $u_y(x,0)=0$ 得

$$Y'(0)=0.$$

而

$$Y'(y)=-A\lambda\sin(\lambda y)+B\lambda\cos(\lambda y),$$

所以 $Y'(0)=B\lambda=0$，取 $\lambda>0$ 時有 $B=0$。因此

$$Y(y)=\cos(\lambda y).$$

再由 $u_y(x,1)+u(x,1)=0$：

$$Y'(1)+Y(1)=0.$$

代入 $Y(y)=\cos(\lambda y)$：

$$-\lambda\sin\lambda+\cos\lambda=0.$$

因此特徵值 $\lambda_n$ 滿足

$$\lambda_n\sin\lambda_n=\cos\lambda_n,$$

也可寫成

$$\tan\lambda_n=\frac{1}{\lambda_n}.$$

這些 $\lambda_n$ 為正根，對應的特徵函數為

$$Y_n(y)=\cos(\lambda_n y).$$

同時，由 $x\to\infty$ 時解必須衰減，故取

$$X_n(x)=e^{-\lambda_n x}.$$

因此解的形式為

$$u(x,y)=\sum_{n=1}^{\infty}A_ne^{-\lambda_n x}\cos(\lambda_n y).$$

求傅立葉係數

在 $x=0$ 時，必須滿足

$$1=\sum_{n=1}^{\infty}A_n\cos(\lambda_n y).$$

這是一組滿足 Neumann–Robin 邊界條件的特徵函數展開。由正交性，

$$A_n= \frac{\displaystyle\int_0^1\cos(\lambda_n y)\,dy} {\displaystyle\int_0^1\cos^2(\lambda_n y)\,dy}.$$

分子為

$$\int_0^1\cos(\lambda_n y)\,dy = \frac{\sin\lambda_n}{\lambda_n}.$$

分母為

$$\int_0^1\cos^2(\lambda_n y)\,dy = \frac12+\frac{\sin(2\lambda_n)}{4\lambda_n}.$$

故

$$A_n= \frac{\dfrac{\sin\lambda_n}{\lambda_n}} {\dfrac12+\dfrac{\sin(2\lambda_n)}{4\lambda_n}}.$$

因此 A 題解為

$$\boxed{ u(x,y)= \sum_{n=1}^{\infty} \frac{\dfrac{\sin\lambda_n}{\lambda_n}} {\dfrac12+\dfrac{\sin(2\lambda_n)}{4\lambda_n}} e^{-\lambda_n x}\cos(\lambda_n y) }$$

其中 $\lambda_n$ 是方程式

$$\tan\lambda_n=\frac1{\lambda_n}$$

的所有正根。

---

B. 半無限長矩形區域

條件為

$$u(0,y)=0,\qquad u(\pi,y)=e^{-y},\qquad y>0,$$

以及

$$u_y(x,0)=0,\qquad 0<x<\pi.$$

題目未明確寫出 $y\to\infty$ 的條件。以下採用此類半無限區域的自然條件：

$$u(x,y)\to 0\qquad (y\to\infty),$$

亦即要求解保持有界並衰減。

解題方法

先取一個能滿足左右邊界條件的函數：

$$u_0(x,y)=\frac{x}{\pi}e^{-y}.$$

確實有

$$u_0(0,y)=0,\qquad u_0(\pi,y)=e^{-y}.$$

令

$$u(x,y)=u_0(x,y)+v(x,y).$$

則 $v$ 滿足

$$v(0,y)=v(\pi,y)=0.$$

因為 $u_y(x,0)=0$，而

本題為半無限域（$y \ge 0$）上的二維拋物型偏微分方程（PDE），相當於一維非穩態擴散方程（其中下游距離 $x$ 扮演時間 $t$ 之角色）。採用**傅立葉餘弦轉換（Fourier Cosine Transform）**推導如下：

一、控制方程與邊界條件

將原偏微分方程改寫為：
$\frac{\partial C}{\partial x} = \frac{k}{u}\frac{\partial^2 C}{\partial y^2}$

邊界與初始條件為：

1. 對稱邊界條件：$\left.\frac{\partial C}{\partial y}\right|_{y=0} = 0$
2. 遠端邊界條件：$C(x, \infty) = 0$
3. $x=0$ 處之線污染源強度條件：$C(0, y) = C_0 \delta(y)$

---

二、對空間變數 $y$ 取傅立葉餘弦轉換

定義傅立葉餘弦轉換 $\mathcal{F}_c\{C(x, y)\} = \tilde{C}(x, \omega) = \int_0^\infty C(x, y) \cos(\omega y) \, dy$。

對 PDE 兩邊取轉換，利用二次微分性質及邊界條件 $\left.\frac{\partial C}{\partial y}\right|_{y=0} = 0$：
$\mathcal{F}_c\left\{\frac{\partial^2 C}{\partial y^2}\right\} = -\omega^2 \tilde{C}(x, \omega) - \left.\frac{\partial C}{\partial y}\right|_{y=0} = -\omega^2 \tilde{C}(x, \omega)$

代入 PDE 得到關於 $x$ 的一階常微分方程（ODE）：
$\frac{d\tilde{C}}{dx} + \frac{k\omega^2}{u} \tilde{C} = 0$

---

三、求解常微分方程

解此 ODE 可得：
$\tilde{C}(x, \omega) = \tilde{C}(0, \omega) \exp\left(-\frac{k x}{u} \omega^2\right)$

由初始條件 $C(0, y) = C_0 \delta(y)$ 求得轉換值：