香港中文大學數學系

2026年5月16日

數值方法透過小步長迭代，從初始條件計算出離散的近似解來求解微分方程，其準確度與效率取決於演算法和時間步長；常見的方法包括歐拉法 (Euler’s method) 和龍格－庫塔法 (Runge-Kutta method)。

MATLAB 提供了多種 ODE 求解器：ode45（通用型龍格－庫塔法）以及四種剛性求解器 (ode15s, ode23s, ode23t, ode2)。

剛性 (Stiff) 常微分方程涉及時間尺度變化差異極大的解（包含快/慢分量），除非使用極小的時間步長，否則會導致標準數值方法變得不穩定；剛性不僅取決於方程式本身，也取決於所使用的數值方法。

要使用 MATLAB 的 ODE 求解器，使用者需定義 ODE 函數、初始條件和時間跨度；求解器會回傳用於繪圖的時間與解向量。

MATLAB 的 dsolve 是一個符號（解析）求解器，與數值求解器不同，它無法處理複雜、非線性或沒有閉式解 (non-closed-form) 的常微分方程。

下表列出了 MATLAB ODE 求解器與其等效的 Python 方法（使用 \(\texttt{scipy.integrate.solve\_ivp}\) 進行數值求解，以及 \(\texttt{sympy.dsolve}\) 進行符號求解）之間的直接對應關係。

\[ \begin{array}{ccc} \hline \textbf{MATLAB 求解器} & \textbf{Python} ~\texttt{solve\_ivp} ~\textbf{方法} & \textbf{用途} \\ \hline \texttt{ode45} & \texttt{RK45} & \textrm{通用型（預設）} \\ \texttt{ode15s} & \texttt{Radau} & \textrm{剛性系統 (Stiff systems)} \\ \texttt{ode23s} & \texttt{BDF} & \textrm{剛性系統} \\ \texttt{ode23t} & \texttt{LSODA} & \textrm{剛性/非剛性自動切換} \\ \texttt{ode23tb} & \texttt{Radau} & \textrm{剛性系統} \\ \texttt{dsolve} & \texttt{sympy.dsolve} & \textrm{解析（符號）解} \\ \hline \end{array} \]

1 範例 1

利用 Python 繪製以下初值問題的解： \[ y' = 4x, \quad y(0) = 1, \quad x \in [0, 2.5] \]

# --------------------------------------------
# 尋找以下問題解析解的 Python 程式碼
# y' = 4x,  y(0) = 1,  x ∈ [0, 2.5]
# 使用符號數學 (SymPy) 來求得精確解
# --------------------------------------------
import sympy as sp
import numpy as np

# 定義符號變數
x = sp.Symbol('x')         # 自變數 x
y = sp.Function('y')(x)    # 應變數 y(x)

# 定義 ODE: y' = 4x
ode = sp.Eq(sp.diff(y, x), 4*x)

# 利用初始條件 y(0) = 1 分析求解 ODE
analytical_sol = sp.dsolve(ode, y, ics={y.subs(x, 0): 1})

# 印出精確的解析解
print("===== Analytical Solution (Exact) =====")

## ===== Analytical Solution (Exact) =====

print(analytical_sol)

## Eq(y(x), 2*x**2 + 1)

# 提取右側部分 (y(x) 的公式)
y_exact = analytical_sol.rhs
print("\nExplicit formula: y(x) =", y_exact)

## 
## Explicit formula: y(x) = 2*x**2 + 1

# --------------------------------------------
# 在任意 x 處評估解析解
# --------------------------------------------
# 將符號公式轉換為數值函數
y_func = sp.lambdify(x, y_exact, 'numpy')

# 在特定點進行測試 (對應您的問題)
x_test1 = 1.3125
x_test2 = 0.5625
y_test1 = y_func(x_test1)
y_test2 = y_func(x_test2)

print("\n===== Test Points =====")

## 
## ===== Test Points =====

print(f"At x = {x_test1}:  y = {y_test1}")

## At x = 1.3125:  y = 4.4453125

print(f"At x = {x_test2}:  y = {y_test2}")

## At x = 0.5625:  y = 1.6328125

# 在整個區間 [0, 2.5] 上進行評估
x_vals = np.linspace(0, 2.5, 10)
y_vals = y_func(x_vals)

print("\n===== Analytical Solution over x ∈ [0, 2.5] =====")

## 
## ===== Analytical Solution over x ∈ [0, 2.5] =====

for xi, yi in zip(x_vals, y_vals):
    print(f"x = {xi:6.3f}    y = {yi:8.4f}")

## x =  0.000    y =   1.0000
## x =  0.278    y =   1.1543
## x =  0.556    y =   1.6173
## x =  0.833    y =   2.3889
## x =  1.111    y =   3.4691
## x =  1.389    y =   4.8580
## x =  1.667    y =   6.5556
## x =  1.944    y =   8.5617
## x =  2.222    y =  10.8765
## x =  2.500    y =  13.5000

# ---------------------------
# Python 實作：ODE 求解與繪圖
# 問題: y' = 4x, y(0) = 1, x ∈ [0, 2.5]
# ---------------------------
import numpy as np
from scipy.integrate import solve_ivp
import matplotlib.pyplot as plt

# ---------------------------
# 1. 定義 ODE 函數
# ---------------------------
def ode_func(x, y):
    """ODE 的右側: dy/dx = 4x"""
    return 4 * x

# ---------------------------
# 2. 設定問題參數
# ---------------------------
x0 = 0                  # 初始 x 值
y0 = [1]                # 初始條件 y(x0) = 1
xspan = [0, 2.5]        # 求解區間

# ---------------------------
# 3. 使用固定的評估點求解 ODE (穩定，無索引錯誤)
# ---------------------------
x_eval = np.linspace(0, 2.5, 100)  # 定義 100 個穩定的點
sol = solve_ivp(
    fun=ode_func,
    t_span=xspan,
    y0=y0,
    method='RK45',
    t_eval=x_eval       # 使用固定點以保證有足夠的數據
)

# 提取求解結果
xsol = sol.t
ysol = sol.y[0]

# ---------------------------
# 4. 用於比較的解析解
# ---------------------------
def analytical_solution(x):
    """解析解: y(x) = 2x² + 1"""
    return 2 * x**2 + 1

y_analytical = analytical_solution(xsol)

# ---------------------------
# 5. 繪製結果
# ---------------------------
plt.figure(figsize=(8, 5))
plt.plot(xsol, ysol, 'b-', linewidth=1.5, label='Numerical Solution (RK45)')
plt.plot(xsol, y_analytical, 'r--', linewidth=1.5, label='Analytical Solution: $y=2x^2+1$')
plt.grid(True, linestyle='--', alpha=0.7)
plt.xlabel('x')
plt.ylabel('y(x)')
plt.title(r'Solution of $y^\prime = 4x, \ y(0)=1$')
plt.legend()
plt.tight_layout()
# plt.savefig('ode_solution_plot.png', dpi=300, bbox_inches='tight')
# plt.close()  # 避免在 R Markdown 中出現顯示停滯

# ---------------------------
# 6. 安全地印出結果 (無索引錯誤)
# ---------------------------
print("=== Numerical Solution Points (First 10) ===")

## === Numerical Solution Points (First 10) ===

print("x          y_numerical  y_analytical  error")

## x          y_numerical  y_analytical  error

print("-" * 45)

## ---------------------------------------------

# 使用 min 以避免索引超出範圍
n_print = min(10, len(xsol))
for i in range(n_print):
    x = xsol[i]
    y_num = ysol[i]
    y_an = y_analytical[i]
    err = abs(y_num - y_an)
    print(f"{x:<10.6f} {y_num:<12.6f} {y_an:<12.6f} {err:.2e}")

## 0.000000   1.000000     1.000000     0.00e+00
## 0.025253   1.001275     1.001275     2.22e-16
## 0.050505   1.005102     1.005102     2.22e-16
## 0.075758   1.011478     1.011478     0.00e+00
## 0.101010   1.020406     1.020406     2.22e-16
## 0.126263   1.031885     1.031885     0.00e+00
## 0.151515   1.045914     1.045914     0.00e+00
## 0.176768   1.062494     1.062494     0.00e+00
## 0.202020   1.081624     1.081624     0.00e+00
## 0.227273   1.103306     1.103306     0.00e+00

# 步長資訊
print("\n=== Step Size Information ===")

## 
## === Step Size Information ===

dx = np.diff(xsol)
if len(dx) >= 5:
    print(f"First 5 step sizes: {dx[:5]}")
else:
    print(f"Step sizes: {dx}")

## First 5 step sizes: [0.02525253 0.02525253 0.02525253 0.02525253 0.02525253]

print(f"Average step size: {np.mean(dx):.6f}")

## Average step size: 0.025253

# 測試點
print("\n=== Key Test Point Verification ===")

## 
## === Key Test Point Verification ===

for x_test in [1.3125, 0.5625]:
    y_test = analytical_solution(x_test)
    print(f"x = {x_test:.6f}, Analytical y = {y_test:.6f}")

## x = 1.312500, Analytical y = 4.445312
## x = 0.562500, Analytical y = 1.632812

2 範例 2

利用 Python 繪製以下初值問題的解： \[ y' + 3y = x e^{2x}, \quad y(2) = -4, \quad x \in [2, 2.7]. \]

精確解析解（來自 SymPy） \[ y(x) = \left( \frac{x}{5} - \frac{1}{25} \right) e^{2x} + \left( -4 e^{6} - \frac{9}{25} e^{4} \right) e^{-3x}. \]

簡化形式 \[ y(x) = \frac{5x - 1}{25} \, e^{2x} + C e^{-3x}, \] 其中 \[ C = -4 e^{6} - \frac{9}{25} e^{4}. \]

# ---------------------------
# Python 程式碼：求下列問題的解析解
# y' + 3y = x*exp(2x), y(2) = -4, x ∈ [2, 2.7]
# ---------------------------
import sympy as sp
import numpy as np

# 定義符號變數
x = sp.Symbol('x')
y = sp.Function('y')(x)

# 定義 ODE: y' + 3y = x*exp(2x)
ode = sp.Eq(sp.diff(y, x) + 3*y, x * sp.exp(2*x))

# 以初始條件 y(2) = -4 求解 ODE
solution = sp.dsolve(ode, y, ics={y.subs(x, 2): -4})

# 印出精確的解析解
print("=== Exact Analytical Solution ===")

## === Exact Analytical Solution ===

print(solution)

## Eq(y(x), x*exp(2*x)/5 - exp(2*x)/25 + (-9*exp(10)/25 - 4*exp(6))*exp(-3*x))

# 提取 y(x) 的顯式公式
y_exact = solution.rhs
print("\nExplicit formula:")

## 
## Explicit formula:

print(f"y(x) = {y_exact}")

## y(x) = x*exp(2*x)/5 - exp(2*x)/25 + (-9*exp(10)/25 - 4*exp(6))*exp(-3*x)

# 將符號解轉換為數值函數
y_func = sp.lambdify(x, y_exact, 'numpy')

# 在關鍵點上進行評估
x_points = np.array([2, 2.2, 2.5, 2.7])
y_values = y_func(x_points)

print("\n=== Solution at key points ===")

## 
## === Solution at key points ===

for xi, yi in zip(x_points, y_values):
    print(f"x = {xi:.2f} → y = {yi:.4f}")

## x = 2.00 → y = -4.0000
## x = 2.20 → y = 19.5980
## x = 2.50 → y = 62.9918
## x = 2.70 → y = 107.8065

# 選擇性：在整個區間上進行評估
x_interval = np.linspace(2, 2.7, 10)
y_interval = y_func(x_interval)

print("\n=== Solution over [2, 2.7] ===")

## 
## === Solution over [2, 2.7] ===

print("x\t\t y(x)")

## x         y(x)

print("---------------------------")

## ---------------------------

for xi, yi in zip(x_interval, y_interval):
    print(f"{xi:.2f}\t\t{yi:.4f}")

## 2.00     -4.0000
## 2.08     5.2233
## 2.16     14.3129
## 2.23     23.6600
## 2.31     33.6461
## 2.39     44.6591
## 2.47     57.1082
## 2.54     71.4389
## 2.62     88.1493
## 2.70     107.8065

# ---------------------------
# 範例 2 的 Python 實作
# 問題: y' + 3y = x*exp(2x), y(2) = -4, x ∈ [2, 2.7]
# ---------------------------
import numpy as np
from scipy.integrate import solve_ivp
import matplotlib.pyplot as plt
import sympy as sp

# ---------------------------
# 1. 定義 ODE 的右側
# ---------------------------
def ode_func(x, y):
    """
    ODE 的右側，經重組以進行數值求解：
    y' = -3y + x*exp(2x)
    """
    return -3 * y + x * np.exp(2 * x)

# ---------------------------
# 2. 設定問題參數
# ---------------------------
x0 = 2                  # 初始 x 值
y0 = [-4]               # 初始條件 y(2) = -4
xspan = [2, 2.7]        # 求解區間

# ---------------------------
# 3. 數值求解 ODE (RK45)
# ---------------------------
# 使用固定點以保證有穩定的輸出
x_eval = np.linspace(2, 2.7, 100)
sol = solve_ivp(
    fun=ode_func,
    t_span=xspan,
    y0=y0,
    method='RK45',
    t_eval=x_eval
)

# 安全地提取結果
xsol = sol.t
ysol = sol.y[0] if sol.y.ndim > 1 else sol.y

# ---------------------------
# 4. 解析解 (SymPy)
# ---------------------------
x_sym = sp.Symbol('x')
y_sym = sp.Function('y')(x_sym)
ode_sym = sp.Eq(sp.diff(y_sym, x_sym) + 3 * y_sym, x_sym * sp.exp(2 * x_sym))
solution_sym = sp.dsolve(ode_sym, y_sym, ics={y_sym.subs(x_sym, 2): -4})

# 將符號解轉換為數值函數
y_analytical_func = sp.lambdify(x_sym, solution_sym.rhs, 'numpy')
y_analytical = y_analytical_func(xsol)

# ---------------------------
# 5. 繪製解答圖
# ---------------------------
plt.figure(figsize=(8, 5))
plt.plot(xsol, ysol, 'b-', linewidth=1.5, label='Numerical Solution (RK45)')
plt.plot(xsol, y_analytical, 'r--', linewidth=1.5, label='Analytical Solution')
plt.grid(True, linestyle='--', alpha=0.7)
plt.xlabel('x')
plt.ylabel('y(x)')
plt.title(r'Solution of $y^\prime + 3y = x e^{2x}, \ y(2)=-4$')
plt.legend()
# plt.tight_layout()
# plt.close()  # 防止 R Markdown 顯示停滯

# ---------------------------
# 6. 安全驗證 (修正：無索引錯誤)
# ---------------------------
print("=== Exact Analytical Solution ===")

## === Exact Analytical Solution ===

print(solution_sym)

## Eq(y(x), x*exp(2*x)/5 - exp(2*x)/25 + (-9*exp(10)/25 - 4*exp(6))*exp(-3*x))

print("\n=== Verification at Key Points ===")

## 
## === Verification at Key Points ===

test_points = [2.0, 2.2, 2.5, 2.7]

# 為了安全起見使用預先計算的解答 (迴圈內不進行求解)
for x_test in test_points:
    # 取得解析值
    y_an = y_analytical_func(x_test)
    
    # 內插數值解 (安全，無索引風險)
    y_num = np.interp(x_test, xsol, ysol)
    
    print(f"x = {x_test:.2f}: Numerical y = {y_num:.4f}, Analytical y = {y_an:.4f}")

## x = 2.00: Numerical y = -4.0000, Analytical y = -4.0000
## x = 2.20: Numerical y = 19.5986, Analytical y = 19.5980
## x = 2.50: Numerical y = 62.9424, Analytical y = 62.9918
## x = 2.70: Numerical y = 107.8438, Analytical y = 107.8065

3 範例 3 一階常微分方程組

利用 Python 求解並繪製以下初值問題的解： \[ \begin{cases} \dfrac{dx}{dt} = y^2 - x \\[6pt] \dfrac{dy}{dt} = 0.5x^2 - y \end{cases} \quad \begin{bmatrix} x(0) \\ y(0) \end{bmatrix} = \begin{bmatrix} 1 \\ 0 \end{bmatrix}, \quad t \in [0, 3]. \]

# ---------------------------
# 範例 2.7 的 Python 實作
# ODE 系統:
# dx/dt = y^2 - x
# dy/dt = 0.5*x^2 - y
# 初始條件: x(0)=1, y(0)=0, t ∈ [0, 3]
# ---------------------------
import numpy as np
from scipy.integrate import solve_ivp
import matplotlib.pyplot as plt

# ---------------------------
# 1. 定義 ODE 系統
# ---------------------------
def ode_system(t, state):
    """
    定義一階 ODE 系統。
    輸入:
        t: 自變數 (時間)
        state: 包含當前狀態變數值的列表 [x, y]
    輸出:
        dstate_dt: 包含導數的列表 [dx/dt, dy/dt]
    """
    x, y = state
    dxdt = y**2 - x
    dydt = 0.5 * x**2 - y
    return [dxdt, dydt]

# ---------------------------
# 2. 設定問題參數
# ---------------------------
t0 = 0                  # 初始時間
u0 = [1, 0]             # 初始狀態 [x(0), y(0)]
tspan = [t0, 3]         # 時間區間 [0, 3]

# ---------------------------
# 3. 數值求解系統 (RK45，相當於 MATLAB 的 ode45)
# ---------------------------
# 定義用於穩定繪圖的固定評估點
t_eval = np.linspace(0, 3, 100)
sol = solve_ivp(
    fun=ode_system,
    t_span=tspan,
    y0=u0,
    method='RK45',      # 自適應 Runge-Kutta 方法 (對應 MATLAB ode45)
    t_eval=t_eval       # 使用固定點以保證輸出
)

# 提取結果
tsol = sol.t            # 時間點
xsol = sol.y[0]         # 解 x(t)
ysol = sol.y[1]         # 解 y(t)

# ---------------------------
# 4. 繪製解答圖
# ---------------------------
plt.figure(figsize=(8, 5))
plt.plot(tsol, xsol, 'b-', linewidth=1.5, label='x(t)')
plt.plot(tsol, ysol, 'r-', linewidth=1.5, label='y(t)')
plt.grid(True, linestyle='--', alpha=0.7)
plt.xlabel('t')
plt.ylabel('x(t), y(t)')
plt.title(r'Solutions of the System: $\frac{dx}{dt}=y^2-x$, $\frac{dy}{dt}=0.5x^2-y$')
plt.legend()
plt.tight_layout()

# 儲存圖片供 LaTeX/R Markdown 使用
# plt.savefig('ode_system_plot.png', dpi=300, bbox_inches='tight')
# plt.close()  # 防止在 R Markdown 中出現顯示停滯

# ---------------------------
# 5. 印出驗證資訊
# ---------------------------
print("=== System ODE Solution (RK45) ===")

## === System ODE Solution (RK45) ===

print("Initial conditions: x(0)=1, y(0)=0")

## Initial conditions: x(0)=1, y(0)=0

print("Time interval: [0, 3]")

## Time interval: [0, 3]

print("\n=== Values at Key Time Points ===")

## 
## === Values at Key Time Points ===

test_times = [0, 1, 2, 3]
for t in test_times:
    # 進行內插以取得測試時間點的數值
    x_val = np.interp(t, tsol, xsol)
    y_val = np.interp(t, tsol, ysol)
    print(f"t = {t:.1f}: x(t) = {x_val:.4f}, y(t) = {y_val:.4f}")

## t = 0.0: x(t) = 1.0000, y(t) = 0.0000
## t = 1.0: x(t) = 0.3756, y(t) = 0.1175
## t = 2.0: x(t) = 0.1428, y(t) = 0.0601
## t = 3.0: x(t) = 0.0535, y(t) = 0.0245

4 範例 4 具固定時間步長的一階常微分方程組

利用 Python 求解並繪製以下初值問題的解： \[ \begin{cases} \dfrac{dy_2}{dt} + 3y_2 = 6y_1 \\[6pt] \dfrac{dy_1}{dt} + 5y_1 = \sin(2t) \end{cases} \quad \begin{bmatrix} y_1(0) \\ y_2(0) \end{bmatrix} = \begin{bmatrix} 5 \\ -2 \end{bmatrix}, \quad t \in [0, 0.5]. \] 重組為標準形式： \[ \begin{cases} \dfrac{dy_1}{dt} = -5y_1 + \sin(2t) \\[6pt] \dfrac{dy_2}{dt} = 6y_1 - 3y_2 \end{cases} \]

# ---------------------------
# 範例 2.8 的 Python 實作
# ODE 系統:
# dy1/dt = -5*y1 + sin(2*t)
# dy2/dt = 6*y1 - 3*y2
# 初始條件: y1(0)=5, y2(0)=-2, t ∈ [0, 0.5]
# 固定時間步長: 0.001
# ---------------------------
import numpy as np
from scipy.integrate import solve_ivp
import matplotlib.pyplot as plt

# ---------------------------
# 1. 定義 ODE 系統
# ---------------------------
def ode_system(t, state):
    """
    定義標準形式的一階 ODE 系統。
    輸入:
        t: 自變數 (時間)
        state: 包含當前狀態變數值的列表 [y1, y2]
    輸出:
        dstate_dt: 包含導數的列表 [dy1/dt, dy2/dt]
    """
    y1, y2 = state
    dy1dt = -5 * y1 + np.sin(2 * t)
    dy2dt = 6 * y1 - 3 * y2
    return [dy1dt, dy2dt]

# ---------------------------
# 2. 設定問題參數
# ---------------------------
t0 = 0                  # 初始時間
y0 = [5, -2]            # 初始狀態 [y1(0), y2(0)]
tspan = [t0, 0.5]       # 時間區間 [0, 0.5]
dt = 0.001              # 固定時間步長 (如問題所指定)

# ---------------------------
# 3. 數值求解系統 (RK45，相當於 MATLAB 的 ode45)
# ---------------------------
# 建立固定步長 0.001 的時間點 (對應 MATLAB 的 tspan)
t_eval = np.arange(t0, 0.5 + dt, dt)
sol = solve_ivp(
    fun=ode_system,
    t_span=tspan,
    y0=y0,
    method='RK45',      # 自適應 Runge-Kutta 方法 (對應 MATLAB ode45)
    t_eval=t_eval       # 使用固定時間步長點以保證輸出
)

# 提取結果
tsol = sol.t            # 時間點 (步長為 0.001)
y1sol = sol.y[0]        # 解 y1(t)
y2sol = sol.y[1]        # 解 y2(t)

# ---------------------------
# 4. 繪製解答圖
# ---------------------------
plt.figure(figsize=(8, 5))
plt.plot(tsol, y1sol, 'b-', linewidth=1.5, label='y1(t)')
plt.plot(tsol, y2sol, 'r-', linewidth=1.5, label='y2(t)')
plt.grid(True, linestyle='--', alpha=0.7)
plt.xlabel('t')
plt.ylabel('y1(t), y2(t)')
plt.title(r'Solutions of the System: $\frac{dy_1}{dt}=-5y_1+\sin(2t)$, $\frac{dy_2}{dt}=6y_1-3y_2$')
plt.legend()
plt.tight_layout()

# 儲存圖片供 LaTeX/R Markdown 使用
# plt.savefig('ode_system_plot_2.png', dpi=300, bbox_inches='tight')
# plt.close()  # 防止在 R Markdown 中出現顯示停滯

# ---------------------------
# 5. 印出驗證資訊
# ---------------------------
print("=== System ODE Solution (RK45) ===")

## === System ODE Solution (RK45) ===

print("Initial conditions: y1(0)=5, y2(0)=-2")

## Initial conditions: y1(0)=5, y2(0)=-2

print("Time interval: [0, 0.5]")

## Time interval: [0, 0.5]

print(f"Fixed time step: {dt}")

## Fixed time step: 0.001

print("\n=== Values at Key Time Points ===")

## 
## === Values at Key Time Points ===

test_times = [0, 0.1, 0.2, 0.3, 0.4, 0.5]
for t in test_times:
    # 進行內插以取得測試時間點的數值 (安全，無索引風險)
    y1_val = np.interp(t, tsol, y1sol)
    y2_val = np.interp(t, tsol, y2sol)
    print(f"t = {t:.2f}: y1(t) = {y1_val:.4f}, y2(t) = {y2_val:.4f}")

## t = 0.00: y1(t) = 5.0000, y2(t) = -2.0000
## t = 0.10: y1(t) = 3.0408, y2(t) = 0.5352
## t = 0.20: y1(t) = 1.8685, y2(t) = 1.6270
## t = 0.30: y1(t) = 1.1716, y2(t) = 1.9683
## t = 0.40: y1(t) = 0.7616, y2(t) = 1.9447
## t = 0.50: y1(t) = 0.5240, y2(t) = 1.7648

5 範例 5 具固定時間步長的 3 變量常微分方程組

利用 Python 求解並繪製以下初值問題的解： \[ \begin{cases} \dfrac{dx}{dt} = y - z \\[6pt] \dfrac{dy}{dt} = 0.45z - e^{-t} \\[6pt] \dfrac{dz}{dt} = -0.25xy + t^2 \end{cases} \quad \begin{bmatrix} x(0) \\ y(0) \\ z(0) \end{bmatrix} = \begin{bmatrix} -2 \\ -5 \\ 7 \end{bmatrix}, \quad t \in [0, 4]. \]

# ---------------------------
# 範例 2.9 的 Python 實作
# 3 變量 ODE 系統:
# dx/dt = y - z
# dy/dt = 0.45*z - exp(-t)
# dz/dt = -0.25*x*y + t^2
# 初始條件: x(0)=-2, y(0)=-5, z(0)=7
# 時間區間: t ∈ [0, 4]
# 固定時間步長: 0.05
# ---------------------------
import numpy as np
from scipy.integrate import solve_ivp
import matplotlib.pyplot as plt

# ---------------------------
# 1. 定義 ODE 系統
# ---------------------------
def ode_system(t, state):
    """
    定義包含三個一階 ODE 的系統。
    輸入:
        t: 自變數 (時間)
        state: 包含當前狀態變數值的列表 [x, y, z]
    輸出:
        dstate_dt: 包含導數的列表 [dx/dt, dy/dt, dz/dt]
    """
    x, y, z = state
    dxdt = y - z
    dydt = 0.45 * z - np.exp(-t)
    dzdt = -0.25 * x * y + t**2
    return [dxdt, dydt, dzdt]

# ---------------------------
# 2. 設定問題參數
# ---------------------------
t0 = 0                  # 初始時間
u0 = [-2, -5, 7]       # 初始狀態 [x(0), y(0), z(0)]
tspan = [t0, 4]         # 時間區間 [0, 4]
dt = 0.05               # 固定時間步長 (如問題所指定)

# ---------------------------
# 3. 數值求解系統 (RK45，相當於 MATLAB 的 ode45)
# ---------------------------
# 建立固定步長 0.05 的時間點 (對應 MATLAB 的 tspan)
t_eval = np.arange(t0, 4 + dt, dt)
sol = solve_ivp(
    fun=ode_system,
    t_span=tspan,
    y0=u0,
    method='RK45',      # 自適應 Runge-Kutta 方法 (對應 MATLAB ode45)
    t_eval=t_eval       # 使用固定時間步長點以保證輸出
)

# 提取結果
tsol = sol.t            # 時間點 (步長為 0.05)
xsol = sol.y[0]         # 解 x(t)
ysol = sol.y[1]         # 解 y(t)
zsol = sol.y[2]         # 解 z(t)

# ---------------------------
# 4. 繪製解答圖
# ---------------------------
plt.figure(figsize=(8, 5))
plt.plot(tsol, xsol, 'b-', linewidth=1.5, label='x(t)')
plt.plot(tsol, ysol, 'r-', linewidth=1.5, label='y(t)')
plt.plot(tsol, zsol, 'y-', linewidth=1.5, label='z(t)')
plt.grid(True, linestyle='--', alpha=0.7)
plt.xlabel('t')
plt.ylabel('x(t), y(t), z(t)')
plt.title(r'Solutions of the 3-Variable ODE System')
plt.legend()
plt.tight_layout()

# 儲存圖片供 LaTeX/R Markdown 使用
# plt.savefig('ode_system_plot_3.png', dpi=300, bbox_inches='tight')
# plt.close()  # 防止在 R Markdown 中出現顯示停滯

# ---------------------------
# 5. 印出驗證資訊
# ---------------------------
print("=== 3-Variable System ODE Solution (RK45) ===")

## === 3-Variable System ODE Solution (RK45) ===

print("Initial conditions: x(0)=-2, y(0)=-5, z(0)=7")

## Initial conditions: x(0)=-2, y(0)=-5, z(0)=7

print("Time interval: [0, 4]")

## Time interval: [0, 4]

print(f"Fixed time step: {dt}")

## Fixed time step: 0.05

print("\n=== Values at Key Time Points ===")

## 
## === Values at Key Time Points ===

test_times = [0, 1, 2, 3, 4]
for t in test_times:
    # 進行內插以取得測試時間點的數值 (安全，無索引風險)
    x_val = np.interp(t, tsol, xsol)
    y_val = np.interp(t, tsol, ysol)
    z_val = np.interp(t, tsol, zsol)
    print(f"t = {t:.1f}: x(t) = {x_val:.4f}, y(t) = {y_val:.4f}, z(t) = {z_val:.4f}")

## t = 0.0: x(t) = -2.0000, y(t) = -5.0000, z(t) = 7.0000
## t = 1.0: x(t) = -10.3561, y(t) = -3.7487, z(t) = 0.3960
## t = 2.0: x(t) = -10.4528, y(t) = -5.9611, z(t) = -9.7007
## t = 3.0: x(t) = -4.2575, y(t) = -12.9089, z(t) = -19.7099
## t = 4.0: x(t) = -2.4123, y(t) = -21.5900, z(t) = -17.1283

6 範例 6 剛性 3 變量常微分方程組 (羅伯遜反應動力學 / Robertson Reaction Kinetics)

利用 Python 求解並繪製以下初值問題的解： \[ \begin{cases} \dfrac{dy_1}{dt} = -k_1 y_1 + k_3 y_2 y_3 \\[6pt] \dfrac{dy_2}{dt} = k_1 y_1 - k_2 y_2^2 - k_3 y_2 y_3 \\[6pt] \dfrac{dy_3}{dt} = k_2 y_2^2 \end{cases} \quad \begin{bmatrix} y_1(0) \\ y_2(0) \\ y_3(0) \end{bmatrix} = \begin{bmatrix} 1 \\ 0 \\ 0 \end{bmatrix}, \] 其中速率常數為 \(k_1 = 0.04\)、\(k_2 = 10^4\)、\(k_3 = 3 \times 10^7\)，區間為 \(t \in [0, 4 \times 10^6]\)。

# ---------------------------
# 範例 2.13 的 Python 實作
# 剛性 3 變量 ODE 系統 (羅伯遜反應動力學)
# dy1/dt = -k1*y1 + k3*y2*y3
# dy2/dt = k1*y1 - k2*y2^2 - k3*y2*y3
# dy3/dt = k2*y2^2
# 初始條件: y1(0)=1, y2(0)=0, y3(0)=0
# 速率常數: k1=0.04, k2=1e4, k3=3e7
# 時間區間: t ∈ [0, 4e6]
# ---------------------------
import numpy as np
from scipy.integrate import solve_ivp
import matplotlib.pyplot as plt

# ---------------------------
# 1. 定義參數與 ODE 系統
# ---------------------------
k1 = 0.04
k2 = 1e4
k3 = 3e7

def ode_system(t, state):
    """
    定義包含三個一階 ODE 的剛性系統 (羅伯遜動力學)。
    輸入:
        t: 自變數 (時間)
        state: 包含當前狀態變數值的列表 [y1, y2, y3]
    輸出:
        dstate_dt: 包含導數的列表 [dy1/dt, dy2/dt, dy3/dt]
    """
    y1, y2, y3 = state
    dy1dt = -k1 * y1 + k3 * y2 * y3
    dy2dt = k1 * y1 - k2 * y2**2 - k3 * y2 * y3
    dy3dt = k2 * y2**2
    return [dy1dt, dy2dt, dy3dt]

# ---------------------------
# 2. 設定問題參數
# ---------------------------
t0 = 0
u0 = [1, 0, 0]
tspan = [t0, 4e6]

# ---------------------------
# 3. 求解 ODE (修正: 無 t_eval 錯誤)
# ---------------------------
# 解法: 使用線性空間而非對數空間 (避免 0 → 1 的跳躍)
t_eval = np.linspace(0, 4e6, 200)

sol = solve_ivp(
    fun=ode_system,
    t_span=tspan,
    y0=u0,
    method='Radau',    # 剛性求解器 (類似 MATLAB 的 ode15s)
    t_eval=t_eval
)

# 安全地提取結果
tsol = sol.t
y1sol = sol.y[0]
y2sol = sol.y[1]
y3sol = sol.y[2]

# ---------------------------
# 4. 繪製解答圖
# ---------------------------
plt.figure(figsize=(8, 10))

plt.subplot(311)
plt.plot(tsol, y1sol, 'b-', linewidth=1.5)
plt.grid(True, alpha=0.7)
plt.ylabel('$y_1(t)$')
plt.title('Stiff Robertson ODE System Solution')

plt.subplot(312)
plt.plot(tsol, y2sol, 'r-', linewidth=1.5)
plt.grid(True, alpha=0.7)
plt.ylabel('$y_2(t)$')

plt.subplot(313)
plt.plot(tsol, y3sol, 'g-', linewidth=1.5)
plt.grid(True, alpha=0.7)
plt.xlabel('t')
plt.ylabel('$y_3(t)$')

# plt.tight_layout()
# plt.close()  # 對於 R Markdown 至關重要

# ---------------------------
# 5. 印出驗證資訊 (安全: 無索引/超出範圍錯誤)
# ---------------------------
print("=== Stiff 3-Variable ODE System Solution ===")

## === Stiff 3-Variable ODE System Solution ===

print("Initial conditions: y1(0)=1, y2(0)=0, y3(0)=0")

## Initial conditions: y1(0)=1, y2(0)=0, y3(0)=0

print("Time interval: [0, 4e6]")

## Time interval: [0, 4e6]

print("\n=== Key Time Points ===")

## 
## === Key Time Points ===

test_times = [0, 1e6, 2e6, 3e6, 4e6]

for t in test_times:
    y1_val = np.interp(t, tsol, y1sol)
    y2_val = np.interp(t, tsol, y2sol)
    y3_val = np.interp(t, tsol, y3sol)
    print(f"t = {t:.1e} | y1 = {y1_val:.6f} | y2 = {y2_val:.6e} | y3 = {y3_val:.6f}")

## t = 0.0e+00 | y1 = 1.000000 | y2 = 0.000000e+00 | y3 = 0.000000
## t = 1.0e+06 | y1 = 0.996242 | y2 = 3.535549e-07 | y3 = 0.003757
## t = 2.0e+06 | y1 = 0.995268 | y2 = 2.804462e-07 | y3 = 0.004731
## t = 3.0e+06 | y1 = 0.994586 | y2 = 2.449147e-07 | y3 = 0.005414
## t = 4.0e+06 | y1 = 0.994042 | y2 = 2.224892e-07 | y3 = 0.005958