香港中文大學數學系

2026年5月16日

一階常微分方程（First-order ordinary differential equation）只包含未知函數對其自變量的一階導數。這些方程式廣泛應用於數學、物理、工程及其他許多科學領域。本文探討兩種線性一階微分方程經典應用的數學形式與解法：

人口增長與衰減 (Population growth and decay)
混合問題 (Mixing problems)

1 人口增長與衰減

1.1 指數變化

在許多現實生活的情況中，某個量 \(y(t)\) 隨時間 \(t\) 變化的速度取決於其當前的值。例子包括放射性物質質量的衰減和人口規模的變化。這些過程被稱為指數變化。

我們希望求出能解出以下初值問題 (IVP) 的函數 \(y(t)\)： \[\begin{align} \frac{dy}{dt} &= ky & \quad \quad (1) \\ y(0) &= y_0 & \quad \quad (2) \\ \end{align}\] 其中 \(k\) 為常數。

若 \(y>0\) 且該量在增長，則 \(k > 0\)。
若 \(y>0\) 且該量在縮減，則 \(k < 0\)。
若 \(y_0 = 0\)，則對所有時間 \(t\) 皆有 \(y(t) = 0\)（方程式 (1) 的一個簡單解）。

當 \(y \neq 0\) 時，我們重組方程式 (1) 進行分離變量： \[ \frac{1}{y} \frac{dy}{dt} = k. \]

兩邊對 \(t\) 積分： \[ \int \frac{1}{y} \frac{dy}{dt} \, dt = \int k \, dt. \]

在左側使用 \(dy = \frac{dy}{dt} dt\)： \[ \int \frac{1}{y} \, dy = \int k \, dt. \]

計算積分得出： \[ \ln|y| = kt + C, \] 其中 \(C\) 為積分常數。我們使用指數來求解 \(y\)： \[\begin{align*} |y| &= e^{kt + C} \\ y &= \pm e^{kt} e^{C}. \end{align*}\]

設 \(A = \pm e^{C}\)。通解為： \[ y(t) = A e^{kt}. \] 若我們設定 \(A=0\)，此公式亦包含了簡單解 \(y(t)=0\)。

現在我們利用初始條件 (2) 來求 \(A\)。在 \(t=0\) 時： \[ y_0 = A e^{0} = A. \]

代入 \(A = y_0\) 得到最終解： \[ y(t) = y_0 e^{kt}. \quad \quad (3) \]

當 \(k > 0\) 時，方程式 (3) 表示指數增長；當 \(k < 0\) 時，表示指數衰減。常數 \(k\) 稱為速率常數 (rate constant)。

放射性衰變的一個關鍵特徵是半衰期 (half-life) \(T\)：物質衰變至初始量一半所需的時間。對於衰減 (\(k<0\))，在方程式 (3) 中設定 \(y(T) = \frac{1}{2}y_0\)： \[ \frac{1}{2}y_0 = y_0 e^{-kT}. \]

兩邊除以 \(y_0\) (\(y_0 \neq 0\))： \[ \frac{1}{2} = e^{-kT}. \]

兩邊取自然對數： \[\begin{align*} \ln\left(\frac{1}{2}\right) &= -kT \\ -\ln 2 &= -kT. \end{align*}\]

解出 \(T\)： \[ T = \frac{\ln 2}{k}. \quad \quad (4) \]

1.1.1 前向歐拉公式 (Forward Euler Formula)

前向歐拉法（前向格式）是一種簡單的數值技巧，用來求解常微分方程 (ODE)： \[ \frac{dy}{dt} = ky,\quad y(0) = y_0 \]

前向歐拉公式 \[ y_{n+1} = y_n + dt \cdot k \cdot y_n \] 其中

\(y_n\) = 當前值
\(dt\) = 時間步長
\(y_{n+1}\) = 下一個值

此程式碼：

實作前向歐拉格式
比較數值結果與精確解析解 (\(y(t) = y_0 e^{kt}\))
繪製兩種解的圖形以便視覺化
適用於增長 (\(k>0\)) 及衰減 (\(k<0\))

import numpy as np
import matplotlib.pyplot as plt

# --------------------------
# 使用者定義參數
# --------------------------
y0 = 100        # 初始值 (例如：人口、質量)
k = 0.1         # 速率常數：k>0 = 增長，k<0 = 衰減
t_end = 20      # 總模擬時間
dt = 0.5        # 時間步長 (越小越準確)

# --------------------------
# 步驟 1：建立時間陣列
# --------------------------
t = np.arange(0, t_end + dt, dt)
n_steps = len(t)

# --------------------------
# 步驟 2：前向歐拉法 (前向格式)
# --------------------------
y_euler = np.zeros(n_steps)
y_euler[0] = y0  # 設定初始條件

# 利用前向格式進行迭代
for i in range(n_steps - 1):
    y_euler[i+1] = y_euler[i] + dt * k * y_euler[i]

# --------------------------
# 步驟 3：精確解析解
# --------------------------
y_exact = y0 * np.exp(k * t)

# --------------------------
# 步驟 4：計算半衰期 (僅適用於衰減)
# --------------------------
if k < 0:
    half_life = np.log(2) / abs(k)
    print(f"=== Radioactive Decay ===")
    print(f"Half-life T = {half_life:.2f} time units")
else:
    print(f"=== Exponential Growth ===")

## === Exponential Growth ===

# --------------------------
# 步驟 5：繪製結果圖形
# --------------------------
plt.figure(figsize=(10, 5))
plt.plot(t, y_euler, 'o-', label='Forward Euler (Numerical)', linewidth=2, markersize=4)
plt.plot(t, y_exact, '-', label='Exact Solution (Analytical)', linewidth=2)
plt.xlabel('Time t')
plt.ylabel('Quantity y(t)')
plt.title(f'Exponential Change (k={k}, y0={y0})')
plt.legend()
plt.grid(True)
plt.show()

# 驗證用，印出前五個數值
print("\nFirst 5 time points:")

## 
## First 5 time points:

print("t\tEuler\t\tExact")

## t    Euler       Exact

for i in range(5):
    print(f"{t[i]:.1f}\t{y_euler[i]:.2f}\t\t{y_exact[i]:.2f}")

## 0.0  100.00      100.00
## 0.5  105.00      105.13
## 1.0  110.25      110.52
## 1.5  115.76      116.18
## 2.0  121.55      122.14

1.1.2 後向歐拉公式 (Backward Euler Formula)

後向歐拉法（隱式/後向格式）是一種求解常微分方程的隱式數值方法： \[ \frac{dy}{dt} = ky,\quad y(0) = y_0 \]

後向歐拉公式 \[ y_{n+1} = y_n + dt\cdot k\cdot y_{n+1} \]

我們將其重組以得到顯式的更新規則（解出 \(y_{n+1}\)）： \[ y_{n+1} = \frac{y_n}{1 - k\cdot dt} \]

後向歐拉法是隱式的（兩邊都使用到了 \(y_{n+1}\)），但對於線性常微分方程 \(y' = ky\)，我們可以直接解出它： \[\begin{align*} y_{n+1} &= y_n + k\,dt\,y_{n+1} \\ y_{n+1} - k\,dt\,y_{n+1} &= y_n \\ y_{n+1}(1 - k\,dt) &= y_n \\ y_{n+1} &= \frac{y_n}{1 - k\,dt} \end{align*}\]

此程式碼：

實作後向歐拉格式
與精確解析解進行比較
繪製結果圖形
適用於指數增長 (\(k>0\)) 及衰減 (\(k<0\))

import numpy as np
import matplotlib.pyplot as plt

# --------------------------
# 模擬參數
# --------------------------
y0 = 100        # 初始值
k = 0.1         # 速率常數 (k>0=增長，k<0=衰減)
t_end = 20      # 總模擬時間
dt = 0.5        # 時間步長

# --------------------------
# 時間陣列
# --------------------------
t = np.arange(0, t_end + dt, dt)
n_steps = len(t)

# --------------------------
# 後向歐拉法 (後向格式)
# --------------------------
y_backward = np.zeros(n_steps)
y_backward[0] = y0

# 利用後向歐拉公式進行迭代
for i in range(n_steps - 1):
    # 線性 ODE dy/dt = ky 重組後的顯式公式
    y_backward[i+1] = y_backward[i] / (1 - k * dt)

# --------------------------
# 精確解析解
# --------------------------
y_exact = y0 * np.exp(k * t)

# --------------------------
# 半衰期計算 (適用於衰減)
# --------------------------
if k < 0:
    half_life = np.log(2) / abs(k)
    print("=== Exponential Decay ===")
    print(f"Half-life T = {half_life:.2f} time units")
else:
    print("=== Exponential Growth ===")

## === Exponential Growth ===

# --------------------------
# 繪製結果
# --------------------------
plt.figure(figsize=(10, 5))
plt.plot(t, y_backward, 's-', label='Backward Euler (Numerical)', linewidth=2, markersize=5)
plt.plot(t, y_exact, '-', label='Exact Solution', linewidth=2)
plt.xlabel('Time t')
plt.ylabel('Quantity y(t)')
plt.title(f'Backward Euler - Exponential Change (k={k}, y0={y0})')
plt.legend()
plt.grid(True)
plt.show()

# --------------------------
# 印出部分樣本數值
# --------------------------
print("\nFirst 5 time points comparison:")

## 
## First 5 time points comparison:

print(f"t\tBackward Euler\tExact")

## t    Backward Euler  Exact

for i in range(5):
    print(f"{t[i]:.1f}\t{y_backward[i]:.2f}\t\t{y_exact[i]:.2f}")

## 0.0  100.00      100.00
## 0.5  105.26      105.13
## 1.0  110.80      110.52
## 1.5  116.64      116.18
## 2.0  122.77      122.14

1.1.3 中心差分格式 (Central Difference Scheme)

對於常微分方程 \[ \frac{dy}{dt} = ky, \] 中心差分格式（二階精度）使用以下公式： \[ \frac{y_{n+1} - y_{n-1}}{2dt} = k y_n. \]

重組後得到更新規則： \[ y_{n+1} = y_{n-1} + 2\cdot dt\cdot k\cdot y_n. \]

這是一種多步法（需要兩個初始值：\(y_0\) 和 \(y_1\)）。我們先使用前向歐拉法計算 \(y_1\) 以進行初始化。

import numpy as np
import matplotlib.pyplot as plt

# --------------------------
# 模擬參數
# --------------------------
y0 = 100        # 初始值 y(0)
k = 0.1         # 速率常數 (k>0=增長，k<0=衰減)
t_end = 20      # 總模擬時間
dt = 0.5        # 時間步長

# --------------------------
# 建立時間陣列
# --------------------------
t = np.arange(0, t_end + dt, dt)
n_steps = len(t)

# --------------------------
# 中心差分格式 (蛙跳法)
# --------------------------
y_central = np.zeros(n_steps)
y_central[0] = y0  # 初始條件

# 步驟 1：使用前向歐拉法計算 y1 (多步法的初始化)
y_central[1] = y_central[0] + dt * k * y_central[0]

# 步驟 2：利用中心差分法更新剩餘的步數
for n in range(1, n_steps - 1):
    y_central[n+1] = y_central[n-1] + 2 * dt * k * y_central[n]

# --------------------------
# 精確解析解
# --------------------------
y_exact = y0 * np.exp(k * t)

# --------------------------
# 半衰期 (僅適用於衰減)
# --------------------------
if k < 0:
    half_life = np.log(2) / abs(k)
    print("=== Exponential Decay ===")
    print(f"Half-life T = {half_life:.2f} time units")
else:
    print("=== Exponential Growth ===")

## === Exponential Growth ===

# --------------------------
# 繪製結果
# --------------------------
plt.figure(figsize=(10, 5))
plt.plot(t, y_central, 'D-', label='Central Difference Scheme', linewidth=2, markersize=5)
plt.plot(t, y_exact, '-', label='Exact Solution', linewidth=2)
plt.xlabel('Time t')
plt.ylabel('Quantity y(t)')
plt.title(f'Central Difference - Exponential Change (k={k}, y0={y0})')
plt.legend()
plt.grid(True)
plt.show()

# --------------------------
# 印出比較表
# --------------------------
print("\nFirst 5 time points:")

## 
## First 5 time points:

print(f"t\tCentral\t\tExact")

## t    Central     Exact

for i in range(5):
    print(f"{t[i]:.1f}\t{y_central[i]:.2f}\t\t{y_exact[i]:.2f}")

## 0.0  100.00      100.00
## 0.5  105.00      105.13
## 1.0  110.50      110.52
## 1.5  116.05      116.18
## 2.0  122.11      122.14

1.2 放射性衰變

考慮一種半衰期為 1720 年的放射性物質。在時間 \(t=0\) 時，該物質的質量為 6 克。我們的目標是求出：

860 年後剩餘的質量。
剩餘質量達到 2.5 克時所需的時間 \(t_1\)。

由初始條件可知，\(t_0 = 0\) 且 \(y_0 = 6\) 克。將這些值代入通解 (3)，得到隨時間變化的質量函數為： \[ y(t) = 6 e^{-kt}. \quad \quad (5) \] 請注意指數中的負號，這明確地表示了衰減 (\(k>0\))。

已知半衰期 \(T = 1720\) 年。使用半衰期公式 (4)，解出速率常數 \(k\)： \[\begin{align*} k &= \frac{\ln 2}{T} \\ &= \frac{\ln 2}{1720}. \end{align*}\]

將此 \(k\) 值代入方程式 (5) 中，得出該物質的特解： \[ y(t) = 6 e^{-\frac{t \ln 2}{1720}}. \quad \quad (6) \]

1.2.1 860年後的質量

為求出 860 年後剩餘的質量，我們計算 \(t=860\) 時的 \(y(t)\)： \[\begin{align*} y(860) &= 6 e^{-\frac{860 \ln 2}{1720}} \\ &= 6 e^{-\frac{\ln 2}{2}} \\ &= 6 \left( e^{\ln 2} \right)^{-1/2} \\ &= 6 (2)^{-1/2} \\ &= \frac{6}{\sqrt{2}} \\ &= 3\sqrt{2} \\ &\approx 4.24 \text{ 克}. \end{align*}\] 因此，860年後剩餘的質量約為 4.24 克。

1.2.2 達到 2.5 克所需的時間

為求出 \(y(t_1) = 2.5\) 克時的時間 \(t_1\)，我們將這些值代入方程式 (6)： \[ 2.5 = 6 e^{-\frac{t_1 \ln 2}{1720}}. \] 請注意 \(2.5 = \frac{5}{2}\)，因此方程式變為： \[ \frac{5}{2} = 6 e^{-\frac{t_1 \ln 2}{1720}}. \] 兩邊除以 \(6\)： \[ \frac{5}{12} = e^{-\frac{t_1 \ln 2}{1720}}. \] 兩邊取自然對數： \[ \ln\left(\frac{5}{12}\right) = -\frac{t_1 \ln 2}{1720}. \] 利用對數性質 \(\ln(a/b) = -\ln(b/a)\)，我們改寫左邊： \[ -\ln\left(\frac{12}{5}\right) = -\frac{t_1 \ln 2}{1720}. \] 消去負號並解出 \(t_1\)： \[\begin{align*} t_1 &= 1720 \cdot \frac{\ln\left(\frac{12}{5}\right)}{\ln 2} \\ &\approx 2172.42 \text{ 年}. \end{align*}\] 因此，大約在 2172.42 年後，該物質的質量將減少至 2.5 克。

1.3 人口增長

珠海市在 2010 年的人口為 1,108,000 人，2011 年為 1,138,000 人。假設人口以恆定速率呈指數增長，我們的目標是預測 2026 年的人口。

我們定義 2010 年為 \(t=0\)。初始人口為 \(y_0 = 1,108,000\)。使用指數增長的通解 (3)，在時間 \(t\) 的人口為： \[ y(t) = 1108000 e^{kt}. \quad \quad (7) \]

為了確定增長常數 \(k\)，我們使用 2011 年的人口數據，對應於 \(t=1\)： \[ y(1) = 1138000 = 1108000 e^{k \cdot 1}. \] 兩邊除以 1,108,000： \[ e^k = \frac{1138000}{1108000} = \frac{1138}{1108}. \] 兩邊取自然對數： \[ k = \ln\left(\frac{1138}{1108}\right) \approx 0.026716. \]

將此 \(k\) 值代回方程式 (7)，人口模型變為： \[ y(t) = 1108000 e^{0.026716 t}. \]

為了預測 2026 年的人口，請注意 2026 年對應於 \(t = 2026 - 2010 = 15\)。計算 \(t=15\) 時的模型： \[\begin{align*} y(15) &= 1108000 e^{0.026716 \cdot 15} \\ &\approx 1108000 e^{0.40074} \\ &\approx 1108000 \cdot 1.4924 \\ &\approx 1654165.39. \end{align*}\] 四捨五入至整數，預計 2026 年珠海市的人口約為 1,654,165 人。

2 混合問題

2.1 問題陳述

考慮一個混合均勻的水箱，初始裝有 \(V_0\) 升的溶液，並含有 \(A_0\) 克的已溶解化學物質。含有相同化學物質、濃度為 \(c_1\) 克/升的溶液以恆定速率 \(r_1\) 升/分鐘流入水箱中。均勻混合的溶液則以恆定速率 \(r_2\) 升/分鐘流出水箱。由於持續攪拌，假設水箱內化學物質的濃度在任何時刻都是均勻的。令 \(A(t)\) 為時間 \(t\) 時水箱內化學物質的量（克），\(V(t)\) 為溶液體積（升）。則流出水箱的化學物質濃度 \(c_2(t)\) 為： \[ c_2(t) = \frac{A(t)}{V(t)}. \quad \quad (7) \]

Mixing Problem.

2.2 數學建模

我們透過分析 \(V(t)\) 和 \(A(t)\) 在微小時間區間 \(\Delta t\) 內的變化，來推導控制它們的微分方程。

體積的變化量 \(\Delta V\) 等於流入體積減去流出體積： \[ \Delta V = r_1 \Delta t - r_2 \Delta t = (r_1 - r_2)\Delta t. \label{eq:vol-change} \quad \quad (8) \]

化學物質量的變化 \(\Delta A\) 等於流入量減去流出量。流入量為 \(c_1 r_1 \Delta t\)。流出量約為 \(c_2 r_2 \Delta t = \frac{A(t)}{V(t)} r_2 \Delta t\)。因此： \[ \Delta A \approx c_1 r_1 \Delta t - \frac{A(t)}{V(t)} r_2 \Delta t = \left(c_1 r_1 - \frac{A(t)}{V(t)} r_2\right)\Delta t. \label{eq:chem-change} \quad \quad (9) \]

將 (8) 和 (9) 式同除以 \(\Delta t\)，得到平均變化率： \[ \frac{\Delta V}{\Delta t} = r_1 - r_2, \quad \frac{\Delta A}{\Delta t} \approx c_1 r_1 - \frac{A(t)}{V(t)} r_2. \]

為求出瞬時變化率，我們取 \(\Delta t \to 0\) 的極限，將近似值轉換為等式： \[\begin{align} \frac{dV}{dt} &= r_1 - r_2 \label{eq:ode-vol} & \quad \quad (10) \\ \frac{dA}{dt} &= c_1 r_1 - \frac{A(t)}{V(t)} r_2. \label{eq:ode-chem} & \quad \quad (11) \end{align}\]

我們先求解 \(V(t)\)。由於 \(r_1\) 和 \(r_2\) 是常數，方程式 (10) 是一個簡單的常微分方程 (ODE)，可直接積分。其通解為： \[ V(t) = (r_1 - r_2)t + C, \] 其中 \(C\) 是積分常數。利用初始條件 \(V(0) = V_0\)，求得 \(C=V_0\)。因此，在任何時間 \(t\) 的體積為： \[ V(t) = (r_1 - r_2)t + V_0. \]

將此 \(V(t)\) 的表達式代入方程式 (11)，得到關於 \(A(t)\) 的線性一階常微分方程： \[ \frac{dA}{dt} = c_1 r_1 - \frac{r_2}{(r_1 - r_2)t + V_0} A(t). \] 將其重新整理為線性常微分方程的標準形式，我們得到： \[ \frac{dA}{dt} + \frac{r_2}{(r_1 - r_2)t + V_0} A(t) = c_1 r_1. \label{eq:linear-ode-chem} \quad \quad (12) \] 此方程式可利用積分因子求解，並受限於初始條件 \(A(0) = A_0\)。

2.3 線性混合問題 ODE 的通解

我們考慮由質量守恆定律推導得出，控制均勻混合水箱中化學物質量 \(A(t)\) 的一階線性常微分方程 (ODE)： \[ \frac{dA}{dt} + \frac{r_2}{(r_1 - r_2)t + V_0} A(t) = c_1 r_1. \label{main_ode} \quad \quad (13) \] 此 ODE 符合線性一階常微分方程的標準形式： \[ \frac{dA}{dt} + P(t)A(t) = Q(t), \] 我們在此明確定義其係數函數： \[ P(t) = \frac{r_2}{(r_1 - r_2)t + V_0}, \quad Q(t) = c_1 r_1. \] 此問題受限於初始條件 (initial condition)： \[ A(0) = A_0. \label{initial_condition} \quad \quad (14) \] 所有參數均為常數實數：

\(r_1\) = 體積流入率 (單位時間內的體積)
\(r_2\) = 體積流出率 (單位時間內的體積)
\(V_0\) = \(t=0\) 時水箱內溶液的初始體積
\(c_1\) = 化學物質的流入濃度 (單位體積內的質量)
\(A_0\) = \(t=0\) 時水箱內化學物質的初始質量

2.3.1 步驟 1：計算積分因子 (Integrating Factor) \(\mu(t)\)

線性 ODE 的積分因子定義為： \[ \mu(t) = \exp\left( \int P(t) \, dt \right). \] 將 \(P(t)\) 代入公式中： \[ \mu(t) = \exp\left( \int \frac{r_2}{(r_1 - r_2)t + V_0} dt \right). \] 我們使用 u-代換法 (u-substitution) 來計算積分：令 \[ u = (r_1 - r_2)t + V_0 \implies du = (r_1 - r_2)dt \implies dt = \frac{du}{r_1 - r_2}. \] 代入積分中： \[ \int \frac{r_2}{u} \cdot \frac{du}{r_1 - r_2} = \frac{r_2}{r_1 - r_2} \int \frac{1}{u} du. \] \(1/u\) 的反導數 (不定積分) 是自然對數： \[ \frac{r_2}{r_1 - r_2} \ln|u| + C = \frac{r_2}{r_1 - r_2} \ln\left[(r_1 - r_2)t + V_0\right]. \] 我們去掉絕對值，因為體積永遠是正數。

取指數以獲得積分因子： \[ \mu(t) = \exp\left( \frac{r_2}{r_1 - r_2} \ln\left[(r_1 - r_2)t + V_0\right] \right). \] 使用對數恆等式 \(a\ln b = \ln b^a\) 及 \(\exp(\ln x)=x\)： \[ \mu(t) = \left[(r_1 - r_2)t + V_0\right]^{\frac{r_2}{r_1 - r_2}}. \quad \quad (15) \label{integrating_factor} \]

2.3.2 步驟 2：將 ODE 乘以積分因子

將方程式 (13) 兩邊乘以 \(\mu(t)\)： \[ \mu(t)\frac{dA}{dt} + \mu(t)P(t)A(t) = \mu(t)Q(t). \] 基於積分因子的設計，左邊即為 \(\mu(t)A(t)\) 的乘積法則導數 (product rule derivative)： \[ \frac{d}{dt}\left[ \mu(t) A(t) \right] = c_1 r_1 \cdot \left[(r_1 - r_2)t + V_0\right]^{\frac{r_2}{r_1 - r_2}}. \label{product_form} \quad \quad (16) \]

2.3.3 步驟 3：兩邊積分

將方程式 (16) 對時間 \(t\) 積分： \[ \mu(t) A(t) = c_1 r_1 \int \left[(r_1 - r_2)t + V_0\right]^{\frac{r_2}{r_1 - r_2}} dt + C, \] 其中 \(C\) 為積分常數。

使用相同的代換法 \(u=(r_1 - r_2)t + V_0\)、\(du=(r_1 - r_2)dt\) 計算積分： \[ \int u^{\frac{r_2}{r_1 - r_2}} \cdot \frac{du}{r_1 - r_2} = \frac{1}{r_1 - r_2} \cdot \frac{u^{\frac{r_2}{r_1 - r_2}+1}}{\frac{r_2}{r_1 - r_2}+1} + C. \] 化簡指數： \[ \frac{r_2}{r_1 - r_2} + 1 = \frac{r_2 + (r_1 - r_2)}{r_1 - r_2} = \frac{r_1}{r_1 - r_2}. \] 因此積分化簡為： \[ \frac{1}{r_1 - r_2} \cdot \frac{u^{\frac{r_1}{r_1 - r_2}}}{\frac{r_1}{r_1 - r_2}} = \frac{1}{r_1} u^{\frac{r_1}{r_1 - r_2}}. \] 將 \(u\) 代回並置入方程式中： \[ \mu(t) A(t) = c_1 r_1 \cdot \frac{1}{r_1}\left[(r_1 - r_2)t + V_0\right]^{\frac{r_1}{r_1 - r_2}} + C. \] 消去 \(r_1\)： \[ \mu(t) A(t) = c_1 \left[(r_1 - r_2)t + V_0\right]^{\frac{r_1}{r_1 - r_2}} + C. \quad \quad (17) \label{integrated_eq} \]

2.3.4 步驟 4：解出 \(A(t)\) (通解)

將方程式 (15) 的 \(\mu(t)\) 代入 (17)： \[ \left[(r_1 - r_2)t + V_0\right]^{\frac{r_2}{r_1 - r_2}} A(t) = c_1 \left[(r_1 - r_2)t + V_0\right]^{\frac{r_1}{r_1 - r_2}} + C. \] 兩邊除以 \(\left[(r_1 - r_2)t + V_0\right]^{\frac{r_2}{r_1 - r_2}}\)： \[ A(t) = c_1 \left[(r_1 - r_2)t + V_0\right] + C \cdot \left[(r_1 - r_2)t + V_0\right]^{-\frac{r_2}{r_1 - r_2}}. \label{general_solution} \quad \quad (18) \] 這是此 ODE 的通解。

2.3.5 步驟 5：代入初始條件 \(A(0)=A_0\)

在 \(t=0\) 時計算方程式 (18)： \[ A(0) = c_1 V_0 + C \cdot V_0^{-\frac{r_2}{r_1 - r_2}} = A_0. \] 解出常數 \(C\)： \[ C \cdot V_0^{-\frac{r_2}{r_1 - r_2}} = A_0 - c_1 V_0, \] \[ C = (A_0 - c_1 V_0) \cdot V_0^{\frac{r_2}{r_1 - r_2}}. \label{constant_c} \quad \quad (19) \]

2.3.6 步驟 6：最終特解

將方程式 (19) 的 \(C\) 代入通解 (18) 中： \[ \begin{aligned} A(t) &= c_1\left[(r_1 - r_2)t + V_0\right] \\ &\quad + (A_0 - c_1 V_0) V_0^{\frac{r_2}{r_1 - r_2}} \left[(r_1 - r_2)t + V_0\right]^{-\frac{r_2}{r_1 - r_2}}. \end{aligned} \] 這就是該初值問題的唯一顯式解，在體積保持正數的所有時間 \(t\) 內皆有效。

2.3.7 鹽水混合問題

考慮一個水桶，初始裝有 100 加侖的鹽水，其中溶解了 50 磅的鹽。含有 2 磅/加侖鹽份的鹽水溶液以 5 加侖/分鐘的速率流入桶中。混合液透過攪拌保持均勻，並以 4 加侖/分鐘的速率流出桶外。我們的目標是求出：

鹽份進入水桶的速率。
時間 \(t\) 時桶內的鹽水體積。
鹽份流出水桶的速率。
描述此混合過程之初值問題的解。
25 分鐘後桶內的鹽份濃度。

我們整理出給定參數：

初始體積：\(V_0 = 100\) 加侖
初始鹽量：\(A_0 = 50\) 磅
流入濃度：\(c_1 = 2\) 磅/加侖
流入速率：\(r_1 = 5\) 加侖/分鐘
流出速率：\(r_2 = 4\) 加侖/分鐘

鹽份進入速率

鹽份進入水桶的速率為流入濃度與流入速率的乘積： \[ \text{Rate in} = c_1 r_1 = (2 \text{ 磅/加侖}) \cdot (5 \text{ 加侖/分鐘}) = 10 \text{ 磅/分鐘}. \]

時間 \(t\) 時的溶液體積

桶中溶液體積以恆定速率變化。使用我們先前推導的公式： \[\begin{align*} V(t) &= V_0 + (r_1 - r_2)t \\ &= 100 + (5 - 4)t \\ &= 100 + t \text{ 加侖}. \end{align*}\]

鹽份流出速率 時間 \(t\) 時桶內的鹽份濃度為 \(\frac{A(t)}{V(t)} = \frac{A(t)}{100 + t}\) 磅/加侖。鹽份流出水桶的速率為該濃度與流出速率的乘積： \[ \text{Rate out} = \left( \frac{A(t)}{100 + t} \text{ 磅/加侖} \right) \cdot (4 \text{ 加侖/分鐘}) = \frac{4A(t)}{100 + t} \text{ 磅/分鐘}. \]

微分方程模型

鹽量的變化率 \(\frac{dA}{dt}\) 為進入速率減去流出速率： \[ \frac{dA}{dt} = 10 - \frac{4A}{100 + t}. \] 重組為標準線性常微分方程形式： \[ \frac{dA}{dt} + \frac{4}{100 + t}A = 10. \] 這裡，\(P(t) = \frac{4}{100 + t}\) 且 \(Q(t) = 10\)。

求解線性常微分方程

我們使用積分因子 \(\mu(t)\) 來解此方程： \[\begin{align*} \mu(t) &= e^{\int P(t) dt} \\ &= e^{\int \frac{4}{100 + t} dt} \\ &= e^{4 \ln|100 + t|} \\ &= e^{\ln(100 + t)^4} \\ &= (100 + t)^4. \end{align*}\]

將標準形式的常微分方程兩邊乘以 \(\mu(t)\)： \[ (100 + t)^4 \frac{dA}{dt} + (100 + t)^4 \frac{4}{100 + t}A = 10(100 + t)^4. \] 左邊即為乘積 \(\mu(t)A(t)\) 的導數： \[ \frac{d}{dt} \left[ (100 + t)^4 A(t) \right] = 10(100 + t)^4. \] 兩邊對 \(t\) 積分： \[ (100 + t)^4 A(t) = \int 10(100 + t)^4 dt. \] 計算積分： \[ (100 + t)^4 A(t) = 10 \cdot \frac{(100 + t)^5}{5} + C = 2(100 + t)^5 + C. \] 解出 \(A(t)\) 得到通解： \[ A(t) = 2(100 + t) + \frac{C}{(100 + t)^4}. \]

我們套用初始條件 \(A(0) = 50\) 磅來找出常數 \(C\)： \[\begin{align*} 50 &= 2(100 + 0) + \frac{C}{(100 + 0)^4} \\ 50 &= 200 + \frac{C}{100^4} \\ C &= (50 - 200) \cdot 100^4 = -150 \cdot 100^4. \end{align*}\] 將 \(C\) 代回通解中得到特解： \[ A(t) = 2(100 + t) - \frac{150 \cdot 100^4}{(100 + t)^4}. \]

25 分鐘後的濃度

首先，我們求出 \(t=25\) 分鐘時桶內的鹽量： \[\begin{align*} A(25) &= 2(100 + 25) - \frac{150 \cdot 100^4}{(100 + 25)^4} \\ &= 250 - \frac{150 \cdot 100^4}{125^4} \\ &\approx 250 - 61.44 \\ &\approx 188.56 \text{ 磅}. \end{align*}\] \(t=25\) 分鐘時溶液的體積為： \[ V(25) = 100 + 25 = 125 \text{ 加侖}. \] 此時的鹽份濃度為： \[ \text{Concentration} = \frac{A(25)}{V(25)} \approx \frac{188.56}{125} \approx 1.51 \text{ 磅/加侖}. \]

中心差分 (Central Difference)

此程式碼實作了混合水槽常微分方程的中心差分（蛙跳）法： \[ \frac{dA}{dt} = 10 - \frac{4A}{100+t} \] 且初始條件為 \(A(0)=50\)。

中心差分法具有二階精度，且比歐拉法精確得多。我們使用前向歐拉法來初始化第一步（多步法的必要步驟）。

使用的中心差分公式 \[ A_{i+1} = A_{i-1} + 2\Delta t\cdot \frac{dA}{dt}\left(t_i,A_i\right) \]

多步法：需要 \(A_0\) 和 \(A_1\)
利用前向歐拉法計算 \(A_1\)
二階精度 \(\to\) 非常接近精確解

範例參數

\(A_0 = 50\) 磅
\(V_0 = 100\) 加侖
\(r_1 = 5,\ r_2 = 4\) 加侖/分鐘
\(c_1 = 2\) 磅/加侖
運算至 \(t=25\) 分鐘

import numpy as np
import matplotlib.pyplot as plt

# --------------------------
# 參數設定 (與您的問題相符)
# --------------------------
A0 = 50          # 初始鹽量 (lb)
V0 = 100         # 初始體積 (gal)
r1 = 5           # 流入速率 (gal/min)
r2 = 4           # 流出速率 (gal/min)
c1 = 2           # 流入濃度 (lb/gal)
t_end = 25       # 模擬至 25 分鐘
dt = 0.1         # 時間步長 (越小越準確)

# --------------------------
# 時間陣列
# --------------------------
t = np.arange(0, t_end + dt, dt)
n = len(t)

# --------------------------
# 體積函數 V(t) (解析解)
# --------------------------
def volume(t):
    return V0 + (r1 - r2) * t

# --------------------------
# ODE 函數：dA/dt = 10 - (4*A)/(100 + t)
# --------------------------
def dA_dt(t_val, A_val):
    V = volume(t_val)
    return c1*r1 - (r2 * A_val) / V

# --------------------------
# 中心差分 (蛙跳) 格式
# --------------------------
A_central = np.zeros(n)
A_central[0] = A0  # 初始條件

# 步驟 1：利用前向歐拉法計算 A1 (初始化多步法)
A_central[1] = A_central[0] + dt * dA_dt(t[0], A_central[0])

# 步驟 2：對剩餘的所有步驟進行中心差分
for i in range(1, n-1):
    # 蛙跳公式：A_{i+1} = A_{i-1} + 2*dt*dA/dt(t_i, A_i)
    A_central[i+1] = A_central[i-1] + 2 * dt * dA_dt(t[i], A_central[i])

# --------------------------
# 精確解析解 (供比較用)
# --------------------------
def exact_solution(t_val):
    V = volume(t_val)
    C = -150 * (100**4)
    return 2 * V + C / (V**4)

A_exact = exact_solution(t)

# --------------------------
# t=25 min 時的結果
# --------------------------
t_target = 25
idx = np.where(np.isclose(t, t_target))[0][0]

V25 = volume(t_target)
A25_central = A_central[idx]
A25_exact = exact_solution(t_target)
conc_central = A25_central / V25
conc_exact = A25_exact / V25

# --------------------------
# 列印輸出
# --------------------------
print("="*60)

## ============================================================

print("MIXING PROBLEM - CENTRAL DIFFERENCE vs EXACT SOLUTION")

## MIXING PROBLEM - CENTRAL DIFFERENCE vs EXACT SOLUTION

print("="*60)

## ============================================================

print(f"Time = {t_target} min")

## Time = 25 min

print(f"Volume V(t) = {V25:.0f} gal")

## Volume V(t) = 125 gal

print(f"\nSalt Amount A(t):")

## 
## Salt Amount A(t):

print(f"Central Difference: {A25_central:.2f} lb")

## Central Difference: 188.56 lb

print(f"Exact Solution:     {A25_exact:.2f} lb")

## Exact Solution:     188.56 lb

print(f"\nConcentration:")

## 
## Concentration:

print(f"Central Difference: {conc_central:.3f} lb/gal")

## Central Difference: 1.508 lb/gal

print(f"Exact Solution:     {conc_exact:.3f} lb/gal")

## Exact Solution:     1.508 lb/gal

print("="*60)

## ============================================================

# --------------------------
# 繪圖
# --------------------------
plt.figure(figsize=(10,5))
plt.plot(t, A_central, 'o-', label='Central Difference', markersize=2, linewidth=1.5)
plt.plot(t, A_exact, '-', label='Exact Solution', linewidth=2)
plt.xlabel('Time t (minutes)')
plt.ylabel('Amount of Salt A(t) (lb)')
plt.title('Mixing Tank Problem: Central Difference Scheme')
plt.legend()
plt.grid(True)
plt.show()

2.3.8 範例：具最大值的肥料混合問題

一個水桶初始裝有 100 加侖的純水。含有 1 磅/加侖肥料的溶液以 1 加侖/分鐘的速率流入桶內。混合液透過攪拌保持均勻，並以 3 加侖/分鐘的速率泵出桶外。我們的目標是求出桶內肥料的最大量，以及達到該最大量所需的時間。

令 \(A(t)\) 為時間 \(t\) 時桶內肥料的量（磅）。我們已知 \(A(0) = 0\)。

時間 \(t\) 時的溶液體積

桶中溶液體積以恆定速率變化： \[\begin{align*} V(t) &= V_0 + (r_1 - r_2)t \\ &= 100 + (1 - 3)t \\ &= 100 - 2t \text{ 加侖}. \end{align*}\] 當 \(V(t)=0\) 時，水桶將會清空，這發生在 \(t=50\) 分鐘。因此，此模型適用於 \(0 \le t < 50\)。

微分方程模型

肥料進入水桶的速率為： \[ \text{Rate in} = c_1 r_1 = (1 \text{ 磅/加侖}) \cdot (1 \text{ 加侖/分鐘}) = 1 \text{ 磅/分鐘}. \] 肥料流出水桶的速率為： \[ \text{Rate out} = \left( \frac{A(t)}{V(t)} \text{ 磅/加侖} \right) \cdot (3 \text{ 加侖/分鐘}) = \frac{3A(t)}{100 - 2t} \text{ 磅/分鐘}. \] 肥料量的變化率 \(\frac{dA}{dt}\) 為進入速率減去流出速率： \[ \frac{dA}{dt} = 1 - \frac{3A}{100 - 2t}. \] 重組為標準線性常微分方程形式： \[ \frac{dA}{dt} + \frac{3}{100 - 2t}A = 1. \] 這裡，\(P(t) = \frac{3}{100 - 2t}\) 且 \(Q(t) = 1\)。

求解線性常微分方程

我們使用積分因子 \(\mu(t)\) 來解此方程： \[\begin{align*} \mu(t) &= e^{\int P(t) dt} \\ &= e^{\int \frac{3}{100 - 2t} dt} \\ &= e^{-\frac{3}{2} \ln|100 - 2t|} \\ &= e^{\ln(100 - 2t)^{-3/2}} \\ &= (100 - 2t)^{-3/2}. \end{align*}\]

將標準形式的常微分方程兩邊乘以 \(\mu(t)\)： \[ (100 - 2t)^{-3/2} \frac{dA}{dt} + (100 - 2t)^{-3/2} \frac{3}{100 - 2t}A = (100 - 2t)^{-3/2}. \] 左邊即為乘積 \(\mu(t)A(t)\) 的導數： \[ \frac{d}{dt} \left[ (100 - 2t)^{-3/2} A(t) \right] = (100 - 2t)^{-3/2}. \] 兩邊對 \(t\) 積分： \[ (100 - 2t)^{-3/2} A(t) = \int (100 - 2t)^{-3/2} dt. \] 使用代換法 \(u=100-2t\)、\(du=-2dt\) 計算積分： \[\begin{align*} \int (100 - 2t)^{-3/2} dt &= -\frac{1}{2} \int u^{-3/2} du \\ &= -\frac{1}{2} \left( \frac{u^{-1/2}}{-1/2} \right) + C \\ &= u^{-1/2} + C \\ &= (100 - 2t)^{-1/2} + C. \end{align*}\] 因此，我們得到： \[ (100 - 2t)^{-3/2} A(t) = (100 - 2t)^{-1/2} + C. \] 解出 \(A(t)\) 得到通解： \[ A(t) = (100 - 2t) + C(100 - 2t)^{3/2}. \]

我們套用初始條件 \(A(0) = 0\) 來求出常數 \(C\)： \[\begin{align*} 0 &= (100 - 0) + C(100 - 0)^{3/2} \\ 0 &= 100 + C(100^{3/2}) \\ 0 &= 100 + C(1000) \\ C &= -\frac{100}{1000} = -\frac{1}{10}. \end{align*}\] 將 \(C\) 代回通解中得到特解： \[ A(t) = (100 - 2t) - \frac{1}{10}(100 - 2t)^{3/2}. \]

求取最大值

為求出肥料的最大量，我們將導數 \(\frac{dA}{dt}\) 設為零。首先，計算我們解答的導數： \[\begin{align*} \frac{dA}{dt} &= \frac{d}{dt} \left[ (100 - 2t) - \frac{1}{10}(100 - 2t)^{3/2} \right] \\ &= -2 - \frac{1}{10} \cdot \frac{3}{2}(100 - 2t)^{1/2}(-2) \\ &= -2 + \frac{3}{10}\sqrt{100 - 2t}. \end{align*}\]

設定 \(\frac{dA}{dt} = 0\)： \[\begin{align*} -2 + \frac{3}{10}\sqrt{100 - 2t} &= 0 \\ \frac{3}{10}\sqrt{100 - 2t} &= 2 \\ \sqrt{100 - 2t} &= \frac{20}{3} \\ 100 - 2t &= \left(\frac{20}{3}\right)^2 = \frac{400}{9} \\ 2t &= 100 - \frac{400}{9} = \frac{900 - 400}{9} = \frac{500}{9} \\ t &= \frac{250}{9} \approx 27.78 \text{ 分鐘}. \end{align*}\]

為了找出肥料的最大量，我們計算 \(t = \frac{250}{9}\) 時的 \(A(t)\)： \[\begin{align*} A\left(\frac{250}{9}\right) &= \left(100 - 2\left(\frac{250}{9}\right)\right) - \frac{1}{10}\left(100 - 2\left(\frac{250}{9}\right)\right)^{3/2} \\ &= \left(100 - \frac{500}{9}\right) - \frac{1}{10}\left(100 - \frac{500}{9}\right)^{3/2} \\ &= \left(\frac{400}{9}\right) - \frac{1}{10}\left(\frac{400}{9}\right)^{3/2} \\ &= \frac{400}{9} - \frac{1}{10} \left(\frac{20}{3}\right)^3 \\ &= \frac{400}{9} - \frac{1}{10} \left(\frac{8000}{27}\right) \\ &\approx 44.44 - 29.63 \\ &\approx 14.81 \text{ 磅}. \end{align*}\] 因此，桶內的最大肥料量約為 14.8 磅，且大約發生在 27.78 分鐘時。