工業数学完全攻略1: 常微分方程式（ODE）

工学生にとって避けられない科目が**工業数学**です。特に常微分方程式（ODE: Ordinary Differential Equation）は、電子回路解析、制御システム設計、機械振動解析など、ほぼすべての工学分野の数学的基盤となります。この記事では、ODEの基礎概念からラプラス変換まで、工学応用例を中心に体系的に説明します。

1. 微分方程式の基礎

微分方程式とは？

**微分方程式**は、未知関数とその導関数の関係を表す方程式です。例えば：

$$\frac{dy}{dx} = 2x, \quad y'' + 3y' + 2y = 0, \quad \frac{\partial u}{\partial t} = k\frac{\partial^2 u}{\partial x^2}$$

微分方程式は、自然現象や工学システムを数学的にモデル化する核心的なツールです。

分類：ODE vs PDE

**常微分方程式（ODE）**: 一つの独立変数に関する導関数のみを含みます。

$$\frac{d^2x}{dt^2} + 2\frac{dx}{dt} + x = 0$$

**偏微分方程式（PDE）**: 二つ以上の独立変数に関する偏微分を含みます。

$$\frac{\partial u}{\partial t} = c^2\frac{\partial^2 u}{\partial x^2}$$

この記事ではODEに焦点を当てます。

階数と次数

**階数**は方程式に現れる最高次の導関数の階数です。

- 1階ODE: $y' + y = x$ (最高導関数は $y'$)

- 2階ODE: $y'' + 2y' + y = 0$ (最高導関数は $y''$)

**次数**は最高階導関数の指数です。ほとんどの工学問題では次数は1です。

線形 vs 非線形

**線形ODE**: 未知関数とその導関数が1次の項のみで現れます。

$$a_n(x)y^{(n)} + a_{n-1}(x)y^{(n-1)} + \cdots + a_1(x)y' + a_0(x)y = f(x)$$

**非線形ODE**: $y^2$、$y \cdot y'$、$\sin(y)$ などの非線形項を含みます。

$$y'' + \sin(y) = 0 \quad \text{(振り子方程式)}$$

一般解、特殊解、特異解

- **一般解**: 任意定数を含む完全な解

- **特殊解**: 初期条件を適用して定数を決定した解

- **特異解**: 一般解から得られない特別な解

工学におけるODE

工学システムのほとんどはODEで記述されます。

**RC充電回路**:

$$RC\frac{dv}{dt} + v = V_s$$

**バネ-質量-ダンパーシステム**:

$$m\ddot{x} + c\dot{x} + kx = F(t)$$

**人口成長モデル**:

$$\frac{dP}{dt} = kP$$

**RL回路**:

$$L\frac{di}{dt} + Ri = E(t)$$

これらすべてがODEです。解法を学ぶことで、工学システムの動的挙動を完全に解析できます。

2. 1階常微分方程式

2.1 変数分離法

**形式**: $\frac{dy}{dx} = f(x)g(y)$

**解法手順**:

$$\frac{dy}{g(y)} = f(x)dx$$

両辺を積分します：

$$\int \frac{dy}{g(y)} = \int f(x)dx + C$$

**例題：RC充電回路**

回路方程式：

$$\frac{dv}{dt} = \frac{V_s - v}{RC}$$

変数分離：

$$\frac{dv}{V_s - v} = \frac{dt}{RC}$$

積分：

$$-\ln|V_s - v| = \frac{t}{RC} + C_1$$

初期条件 $v(0) = 0$ を適用：

$$V_s - v = V_s e^{-t/(RC)}$$

$$\boxed{v(t) = V_s\left(1 - e^{-t/(RC)}\right)}$$

この結果はRC回路の**過渡応答**を表します。時定数 $\tau = RC$ で決まる速度で電圧が $V_s$ に収束します。

**物理的意味**:

- $t = \tau$ のとき: $v \approx 0.632 V_s$ (63.2% 充電)

- $t = 5\tau$ のとき: $v \approx 0.993 V_s$ (ほぼ完全充電)

2.2 積分因子法

**1階線形ODEの標準形**:

$$\frac{dy}{dx} + P(x)y = Q(x)$$

**積分因子**: $\mu(x) = e^{\int P(x)dx}$

**導出過程**:

両辺に $\mu(x)$ を掛けます：

$$\mu(x)\frac{dy}{dx} + \mu(x)P(x)y = \mu(x)Q(x)$$

左辺は積の微分規則により：

$$\frac{d}{dx}[\mu(x)y] = \mu(x)Q(x)$$

積分すると：

$$\mu(x)y = \int \mu(x)Q(x)dx + C$$

$$\boxed{y = \frac{1}{\mu(x)}\left[\int \mu(x)Q(x)dx + C\right]}$$

**例題：RL回路**

$$L\frac{di}{dt} + Ri = E$$

標準形に変換：

$$\frac{di}{dt} + \frac{R}{L}i = \frac{E}{L}$$

積分因子: $\mu(t) = e^{(R/L)t}$

$$\frac{d}{dt}\left[e^{(R/L)t} \cdot i\right] = \frac{E}{L}e^{(R/L)t}$$

積分：

$$e^{(R/L)t} \cdot i = \frac{E}{R}e^{(R/L)t} + C$$

初期条件 $i(0) = 0$ を適用：

$$\boxed{i(t) = \frac{E}{R}\left(1 - e^{-(R/L)t}\right)}$$

時定数 $\tau = L/R$ で電流が定常状態 $E/R$ に収束します。

2.3 完全微分方程式

**形式**: $M(x,y)dx + N(x,y)dy = 0$

**完全条件**:

$$\frac{\partial M}{\partial y} = \frac{\partial N}{\partial x}$$

この条件を満たせば、$F(x,y) = C$ となるポテンシャル関数 $F$ が存在します。

**解法**:

1. $\frac{\partial F}{\partial x} = M$ から $F$ を $x$ について積分

2. $\frac{\partial F}{\partial y} = N$ の条件で積分定数を決定

3. 一般解: $F(x,y) = C$

**例題**:

$(2xy + y^2)dx + (x^2 + 2xy)dy = 0$

完全性の確認：

$$\frac{\partial M}{\partial y} = 2x + 2y = \frac{\partial N}{\partial x}$$

完全方程式なので：

$$F = \int M\,dx = \int (2xy + y^2)dx = x^2y + xy^2 + g(y)$$

$\frac{\partial F}{\partial y} = x^2 + 2xy + g'(y) = N = x^2 + 2xy$ より $g'(y) = 0$、$g(y) = C$

一般解: $x^2y + xy^2 = C$

2.4 ベルヌーイ方程式

**形式**:

$$\frac{dy}{dx} + P(x)y = Q(x)y^n \quad (n \neq 0, 1)$$

**置換**: $v = y^{1-n}$ とすると線形ODEに変換されます。

$v' = (1-n)y^{-n}y'$ なので、元の方程式を $y^n$ で割ると：

$$y^{-n}\frac{dy}{dx} + P(x)y^{1-n} = Q(x)$$

$$\frac{1}{1-n}\frac{dv}{dx} + P(x)v = Q(x)$$

$$\frac{dv}{dx} + (1-n)P(x)v = (1-n)Q(x)$$

これは $v$ に関する1階線形ODEで、積分因子法で解けます。

**応用**: ロジスティック成長モデル $\frac{dP}{dt} = rP\left(1 - \frac{P}{K}\right)$ は $n=2$ のベルヌーイ方程式です。

3. 2階線形常微分方程式

3.1 斉次方程式

**標準形**:

$$ay'' + by' + cy = 0$$

**特性方程式**の導出: $y = e^{rx}$ を代入します。

$$a r^2 e^{rx} + b r e^{rx} + c e^{rx} = 0$$

$$\boxed{ar^2 + br + c = 0}$$

判別式 $\Delta = b^2 - 4ac$ の符号により3つの場合が生じます。

**場合1：異なる2つの実根** ($\Delta > 0$)

$r_1 \neq r_2$ のとき：

$$y = C_1 e^{r_1 x} + C_2 e^{r_2 x}$$

**場合2：重実根** ($\Delta = 0$)

$r_1 = r_2 = r = -b/(2a)$ のとき：

$$y = (C_1 + C_2 x)e^{rx}$$

$C_2 x e^{rx}$ を追加する理由は、**ロンスキー行列式**が0にならないように独立な第2の解を構成する必要があるためです。

**場合3：共役複素根** ($\Delta < 0$)

$r = \alpha \pm i\beta$ (ただし $\alpha = -b/(2a)$、$\beta = \sqrt{4ac - b^2}/(2a)$) のとき：

$$y = e^{\alpha x}(C_1\cos\beta x + C_2\sin\beta x)$$

オイラーの公式 $e^{i\beta x} = \cos\beta x + i\sin\beta x$ を使って実数解に変換した形です。

3.2 RLC回路の完全解析

**回路方程式**（電荷 $q(t)$ 基準）：

$$L\frac{d^2q}{dt^2} + R\frac{dq}{dt} + \frac{q}{C} = 0$$

標準形に整理：

$$q'' + \frac{R}{L}q' + \frac{1}{LC}q = 0$$

**特性方程式**：

$$r^2 + \frac{R}{L}r + \frac{1}{LC} = 0$$

$$r = -\frac{R}{2L} \pm \sqrt{\left(\frac{R}{2L}\right)^2 - \frac{1}{LC}}$$

減衰係数 $\alpha = R/(2L)$、共振周波数 $\omega_0 = 1/\sqrt{LC}$ と定義すると：

$$r = -\alpha \pm \sqrt{\alpha^2 - \omega_0^2}$$

**過減衰（Overdamped）**: $\alpha > \omega_0$ ($R > 2\sqrt{L/C}$)

$$q(t) = C_1 e^{r_1 t} + C_2 e^{r_2 t} \quad (r_1, r_2 < 0)$$

2つの指数項の和で、振動なく単調減少します。

**臨界減衰（Critically Damped）**: $\alpha = \omega_0$ ($R = 2\sqrt{L/C}$)

$$q(t) = (C_1 + C_2 t)e^{-\alpha t}$$

振動なく最も速く平衡状態に戻ります。制御システム設計で非常に重要です。

**不足減衰（Underdamped）**: $\alpha < \omega_0$ ($R < 2\sqrt{L/C}$)

減衰振動周波数 $\omega_d = \sqrt{\omega_0^2 - \alpha^2}$ と定義すると：

$$q(t) = e^{-\alpha t}(C_1\cos\omega_d t + C_2\sin\omega_d t)$$

振動しながら減衰します。通信回路、振動システムでよく見られます。

**物理的意味のまとめ**:

| 条件 | 名称 | 応答特性 |

| ----------------- | -------- | ---------------- |

| $R > 2\sqrt{L/C}$ | 過減衰 | 振動なく遅く減衰 |

| $R = 2\sqrt{L/C}$ | 臨界減衰 | 振動なく最速減衰 |

| $R < 2\sqrt{L/C}$ | 不足減衰 | 振動しながら減衰 |

| $R = 0$ | 無減衰 | 持続振動 |

3.3 非斉次方程式

**形式**:

$$ay'' + by' + cy = f(x)$$

**完全解**: $y = y_h + y_p$

- $y_h$: 斉次方程式の一般解（余関数）

- $y_p$: 特殊解

未定係数法

$f(x)$ の形式に応じて $y_p$ の形式を推測します。

| $f(x)$ の形式 | $y_p$ の推測形式 |

| -------------------------------------- | --------------------------------------- |

| $ke^{ax}$ | $Ae^{ax}$ |

| $kx^n$ | $A_n x^n + \cdots + A_1 x + A_0$ |

| $k\cos\omega x$ または $k\sin\omega x$ | $A\cos\omega x + B\sin\omega x$ |

| $ke^{ax}\cos\omega x$ | $e^{ax}(A\cos\omega x + B\sin\omega x)$ |

**注意**: $y_p$ の形式が $y_h$ の項と重なる場合は $x$ を掛けて修正します。

**例題**: $y'' + 4y = 3\cos 2x$

斉次解: $y_h = C_1\cos 2x + C_2\sin 2x$

$y_p$ の推測形が $y_h$ と重なるため修正：

$$y_p = x(A\cos 2x + B\sin 2x)$$

代入後に計算すると: $A = 0$、$B = 3/4$

$$y_p = \frac{3}{4}x\sin 2x$$

これは**共振（resonance）**現象で、振幅が時間に比例して増加します。

**強制振動応答**:

$$m\ddot{x} + c\dot{x} + kx = F_0\cos\omega t$$

定常状態の特殊解：

$$x_p(t) = \frac{F_0/k}{\sqrt{(1-r^2)^2 + (2\zeta r)^2}}\cos(\omega t - \phi)$$

ここで振動数比 $r = \omega/\omega_n$、減衰比 $\zeta = c/(2\sqrt{mk})$。

$r = 1$（共振条件）で振幅が最大になります。

パラメータ変換法

未定係数法を適用できない一般的な $f(x)$ に使います。

$y_h = C_1 y_1 + C_2 y_2$ が与えられるとき：

$$y_p = -y_1 \int \frac{y_2 f}{W}dx + y_2 \int \frac{y_1 f}{W}dx$$

ロンスキー行列式: $W = y_1 y_2' - y_1' y_2$

この公式はすべての連続関数 $f(x)$ に適用できます。

4. 連立微分方程式

4.1 行列表現

$n$ 個の1階ODEからなる連立方程式：

$$\frac{d\mathbf{X}}{dt} = A\mathbf{X}$$

ここで $\mathbf{X} = [x_1, x_2, \ldots, x_n]^T$、$A$ は $n \times n$ 係数行列。

4.2 固有値・固有ベクトルによる解法

解の形式: $\mathbf{X} = \mathbf{v}e^{\lambda t}$

代入すると: $\lambda \mathbf{v} = A\mathbf{v}$

すなわち $\lambda$ は $A$ の**固有値**、$\mathbf{v}$ は対応する**固有ベクトル**です。

特性方程式: $\det(A - \lambda I) = 0$

**例題**: 2つのタンクが接続されたシステム

$$\frac{d}{dt}\begin{pmatrix} x_1 \\ x_2 \end{pmatrix} = \begin{pmatrix} -2 & 1 \\ 1 & -2 \end{pmatrix}\begin{pmatrix} x_1 \\ x_2 \end{pmatrix}$$

特性方程式: $(\lambda+2)^2 - 1 = 0$

$$\lambda_1 = -1, \quad \lambda_2 = -3$$

固有ベクトル: $\mathbf{v}_1 = [1, 1]^T$、$\mathbf{v}_2 = [1, -1]^T$

一般解：

$$\mathbf{X} = C_1 \begin{pmatrix} 1 \\ 1 \end{pmatrix} e^{-t} + C_2 \begin{pmatrix} 1 \\ -1 \end{pmatrix} e^{-3t}$$

4.3 回路網解析への応用

2ループ回路（キルヒホッフの電圧則）：

$$L_1\frac{di_1}{dt} + R_1 i_1 - M\frac{di_2}{dt} = E_1(t)$$

$$L_2\frac{di_2}{dt} + R_2 i_2 - M\frac{di_1}{dt} = E_2(t)$$

ここで $M$ は相互インダクタンス。行列形式に整理して固有値法またはラプラス変換で解けます。

5. ラプラス変換

5.1 定義と収束

**定義**:

$$\mathcal{L}\{f(t)\} = F(s) = \int_0^\infty f(t)e^{-st}dt$$

$s = \sigma + j\omega$ は複素変数で、積分が収束するには $\text{Re}(s) > \sigma_c$ を満たす必要があります。

**なぜラプラス変換を使うのか？**

- 微分演算が代数演算に変換されます

- 初期条件を自動的に処理します

- 複雑なODEを簡単な代数方程式として解けます

- 伝達関数（Transfer Function）概念の基盤です

5.2 主要変換対

| $f(t)$ | $F(s) = \mathcal{L}\{f(t)\}$ |

| ------------------------ | ---------------------------------- |

| $1$ (単位ステップ) | $\dfrac{1}{s}$ |

| $t$ | $\dfrac{1}{s^2}$ |

| $t^n$ | $\dfrac{n!}{s^{n+1}}$ |

| $e^{at}$ | $\dfrac{1}{s-a}$ |

| $\sin\omega t$ | $\dfrac{\omega}{s^2+\omega^2}$ |

| $\cos\omega t$ | $\dfrac{s}{s^2+\omega^2}$ |

| $e^{at}\sin\omega t$ | $\dfrac{\omega}{(s-a)^2+\omega^2}$ |

| $e^{at}\cos\omega t$ | $\dfrac{s-a}{(s-a)^2+\omega^2}$ |

| $\delta(t)$ (インパルス) | $1$ |

| $u(t-a)$ (遅延ステップ) | $\dfrac{e^{-as}}{s}$ |

5.3 主要性質

**線形性**:

$$\mathcal{L}\{af(t) + bg(t)\} = aF(s) + bG(s)$$

**時間移動（第1移動定理）**:

$$\mathcal{L}\{e^{at}f(t)\} = F(s-a)$$

**周波数移動（第2移動定理）**:

$$\mathcal{L}\{f(t-a)u(t-a)\} = e^{-as}F(s)$$

**微分変換**（最も重要な性質）:

$$\mathcal{L}\{f'(t)\} = sF(s) - f(0)$$

$$\mathcal{L}\{f''(t)\} = s^2F(s) - sf(0) - f'(0)$$

$$\mathcal{L}\{f^{(n)}(t)\} = s^nF(s) - s^{n-1}f(0) - s^{n-2}f'(0) - \cdots - f^{(n-1)}(0)$$

初期条件がs領域で自動的に含まれます。

**積分変換**:

$$\mathcal{L}\left\{\int_0^t f(\tau)d\tau\right\} = \frac{F(s)}{s}$$

**畳み込み定理**:

$$\mathcal{L}\{(f*g)(t)\} = F(s)G(s)$$

ここで $(f*g)(t) = \int_0^t f(\tau)g(t-\tau)d\tau$

**最終値定理**:

$$\lim_{t\to\infty} f(t) = \lim_{s\to 0} sF(s)$$

**初期値定理**:

$$\lim_{t\to 0^+} f(t) = \lim_{s\to\infty} sF(s)$$

5.4 逆ラプラス変換

**部分分数分解**:

$F(s) = P(s)/Q(s)$ で分母を因数分解し、変換対の和に分解します。

**場合1：異なる実数極**

$$F(s) = \frac{N(s)}{(s-p_1)(s-p_2)\cdots(s-p_n)} = \frac{A_1}{s-p_1} + \frac{A_2}{s-p_2} + \cdots$$

$$A_k = \lim_{s\to p_k}(s-p_k)F(s)$$

**場合2：重複極** ($s = p$ が $m$ 重根)

$$F(s) = \cdots + \frac{A_m}{(s-p)^m} + \frac{A_{m-1}}{(s-p)^{m-1}} + \cdots + \frac{A_1}{s-p}$$

$$A_k = \frac{1}{(m-k)!}\left[\frac{d^{m-k}}{ds^{m-k}}(s-p)^m F(s)\right]_{s=p}$$

**場合3：共役複素極**

$$\frac{As + B}{s^2 + 2\alpha s + (\alpha^2 + \beta^2)} \Rightarrow e^{-\alpha t}(C_1\cos\beta t + C_2\sin\beta t)$$

5.5 ラプラス変換でODEを解く

**例題：RLC直列回路の過渡応答**

$$L\frac{d^2i}{dt^2} + R\frac{di}{dt} + \frac{i}{C} = E_0\delta(t)$$

初期条件: $i(0) = 0$、$i'(0) = 0$（キャパシタに電荷なし）

**解法**:

ラプラス変換を適用：

$$I(s) = \frac{E_0/L}{s^2 + (R/L)s + 1/(LC)}$$

$\alpha = R/(2L)$、$\omega_0 = 1/\sqrt{LC}$、$\omega_d = \sqrt{\omega_0^2 - \alpha^2}$ (不足減衰の場合) と定義すると：

$$I(s) = \frac{E_0/L}{(s+\alpha)^2 + \omega_d^2}$$

逆ラプラス変換：

$$i(t) = \frac{E_0}{L\omega_d}e^{-\alpha t}\sin\omega_d t$$

これがRLC回路の**インパルス応答**であり、**グリーン関数**でもあります。

6. PythonでODEを数値的に解く

6.1 scipy.integrate.solve_ivpの活用

from scipy.integrate import solve_ivp

RLC回路シミュレーション（2階ODE -> 1階ODE連立）

L * q'' + R * q' + q/C = Vs（ステップ入力）

状態変数: y[0] = q（電荷）、y[1] = i = dq/dt（電流）

def rlc_circuit(t, y, R, L, C, Vs):

q, dq_dt = y

q'' = (Vs - R*i - q/C) / L

d2q_dt2 = (Vs - R * dq_dt - q / C) / L

return [dq_dt, d2q_dt2]

回路パラメータ

R = 100 # 抵抗 (Ohm)

L = 0.1 # インダクタンス (H)

C = 1e-6 # キャパシタンス (F)

Vs = 10.0 # 電源電圧 (V)

共振周波数と減衰比の計算

omega_0 = 1.0 / np.sqrt(L * C)

alpha = R / (2 * L)

zeta = alpha / omega_0

print(f"共振周波数: {omega_0/(2*np.pi):.1f} Hz")

print(f"減衰比 zeta: {zeta:.3f}")

if zeta > 1:

print("-> 過減衰")

elif zeta == 1:

print("-> 臨界減衰")

else:

print("-> 不足減衰")

数値解法

y0 = [0.0, 0.0] # 初期条件: q(0) = 0, i(0) = 0

t_span = (0, 0.01)

t_eval = np.linspace(0, 0.01, 2000)

sol = solve_ivp(

rlc_circuit, t_span, y0,

t_eval=t_eval,

args=(R, L, C, Vs),

method='RK45',

rtol=1e-8,

atol=1e-10

)

電圧計算: v_C = q/C

v_C = sol.y[0] / C

i = sol.y[1]

可視化

fig, axes = plt.subplots(2, 1, figsize=(10, 8))

axes[0].plot(sol.t * 1e3, v_C, 'b-', linewidth=2, label='キャパシタ電圧')

axes[0].axhline(y=Vs, color='r', linestyle='--', label=f'定常状態 {Vs}V')

axes[0].set_xlabel('時間 (ms)')

axes[0].set_ylabel('電圧 (V)')

axes[0].set_title('RLC直列回路 - ステップ応答')

axes[0].legend()

axes[0].grid(True, alpha=0.3)

axes[1].plot(sol.t * 1e3, i * 1e3, 'g-', linewidth=2, label='電流')

axes[1].set_xlabel('時間 (ms)')

axes[1].set_ylabel('電流 (mA)')

axes[1].set_title('RLC回路電流')

axes[1].legend()

axes[1].grid(True, alpha=0.3)

plt.tight_layout()

plt.savefig('rlc_response.png', dpi=150)

plt.show()

6.2 SymPyで解析解を求める

シンボル定義

t, s = sp.symbols('t s')

R_sym, L_sym, C_sym = sp.symbols('R L C', positive=True)

Vs_sym = sp.Symbol('Vs')

RC回路の解析解（sympyで微分方程式を解く）

v = sp.Function('v')

tau = sp.Symbol('tau', positive=True)

dv/dt = (Vs - v) / (R*C)

ode = sp.Eq(v(t).diff(t), (Vs_sym - v(t)) / (R_sym * C_sym))

初期条件 v(0) = 0

solution = sp.dsolve(ode, v(t), ics={v(0): 0})

print("RC回路の解:")

print(solution)

sp.pprint(solution)

ラプラス変換の計算

f = sp.exp(-t) * sp.sin(2*t)

F = sp.laplace_transform(f, t, s, noconds=True)

print(f"\nL{{e^(-t) sin(2t)}} = {F}")

逆ラプラス変換

F_inv = sp.inverse_laplace_transform(

1 / (s**2 + 2*s + 5), s, t

)

print(f"L^(-1){{1/(s^2+2s+5)}} = {F_inv}")

6.3 各種積分法の比較

オイラー法 vs RK4比較

例題: y' = -y, y(0) = 1, 厳密解 y = e^(-t)

def f(t, y):

return -y

def euler_method(f, t0, y0, h, n):

t = np.zeros(n+1)

y = np.zeros(n+1)

t[0], y[0] = t0, y0

for i in range(n):

y[i+1] = y[i] + h * f(t[i], y[i])

t[i+1] = t[i] + h

return t, y

def rk4_method(f, t0, y0, h, n):

t = np.zeros(n+1)

y = np.zeros(n+1)

t[0], y[0] = t0, y0

for i in range(n):

k1 = h * f(t[i], y[i])

k2 = h * f(t[i] + h/2, y[i] + k1/2)

k3 = h * f(t[i] + h/2, y[i] + k2/2)

k4 = h * f(t[i] + h, y[i] + k3)

y[i+1] = y[i] + (k1 + 2*k2 + 2*k3 + k4) / 6

t[i+1] = t[i] + h

return t, y

h = 0.5

n = 10

t0, y0 = 0.0, 1.0

t_euler, y_euler = euler_method(f, t0, y0, h, n)

t_rk4, y_rk4 = rk4_method(f, t0, y0, h, n)

t_exact = np.linspace(0, n*h, 200)

y_exact = np.exp(-t_exact)

plt.figure(figsize=(10, 5))

plt.plot(t_exact, y_exact, 'k-', linewidth=2, label='厳密解')

plt.plot(t_euler, y_euler, 'r--o', label=f'オイラー法 (h={h})')

plt.plot(t_rk4, y_rk4, 'b--s', label=f'RK4 (h={h})')

plt.xlabel('t')

plt.ylabel("y(t)")

plt.title("オイラー法 vs RK4比較: y' = -y")

plt.legend()

plt.grid(True, alpha=0.3)

plt.tight_layout()

plt.show()

誤差比較

print(f"t=5でのオイラー法誤差: {abs(y_euler[-1] - np.exp(-5)):.6f}")

print(f"t=5でのRK4誤差: {abs(y_rk4[-1] - np.exp(-5)):.6f}")

6.4 ラプラス変換をPythonで自動化

s, t = sp.symbols('s t', real=True)

a, omega = sp.symbols('a omega', positive=True)

ラプラス変換表の自動生成

functions = {

'1': sp.Integer(1),

't': t,

't^2': t**2,

't^n': t**3,

'e^(at)': sp.exp(a*t),

'sin(wt)': sp.sin(omega*t),

'cos(wt)': sp.cos(omega*t),

't*e^(at)': t * sp.exp(a*t),

'e^(at)*sin(wt)': sp.exp(a*t) * sp.sin(omega*t),

}

print("ラプラス変換表:")

print("-" * 50)

for name, func in functions.items():

F = sp.laplace_transform(func, t, s, noconds=True)

print(f"L{{{name}}} = {sp.simplify(F)}")

7. 工学応用総合例題

7.1 電源付きRLC回路の過渡解析

**問題**: $R = 50\,\Omega$、$L = 0.2\,\text{H}$、$C = 100\,\mu\text{F}$ の直列RLC回路に $t=0$ で $V_s = 100\,\text{V}$ DC電源が接続される。初期条件は $v_C(0) = 0$、$i(0) = 0$。

**ラプラス変換で解く**:

キルヒホッフの電圧則：

$$L\frac{di}{dt} + Ri + \frac{1}{C}\int i\,dt = V_s u(t)$$

$i = C dv_C/dt$ で置換して電荷 $q$ で整理：

$$L\ddot{q} + R\dot{q} + \frac{q}{C} = V_s$$

ラプラス変換：

$$\left(Ls^2 + Rs + \frac{1}{C}\right)Q(s) = \frac{V_s}{s}$$

$$Q(s) = \frac{V_s/L}{s\left(s^2 + \frac{R}{L}s + \frac{1}{LC}\right)}$$

数値代入: $\alpha = R/(2L) = 125\,\text{s}^{-1}$、$\omega_0 = 1/\sqrt{LC} \approx 223.6\,\text{s}^{-1}$

$\omega_0 > \alpha$ なので不足減衰です。

$$\omega_d = \sqrt{\omega_0^2 - \alpha^2} \approx 185.0\,\text{s}^{-1}$$

部分分数分解後、逆ラプラス変換：

$$v_C(t) = V_s\left[1 - e^{-\alpha t}\left(\cos\omega_d t + \frac{\alpha}{\omega_d}\sin\omega_d t\right)\right]$$

$$= 100\left[1 - e^{-125t}\left(\cos 185t + 0.676\sin 185t\right)\right]\,\text{V}$$

7.2 バネ-質量-ダンパーの共振解析

**問題**: $m = 1\,\text{kg}$、$k = 100\,\mathrm{N/m}$、$c = 4\,\mathrm{N \cdot s/m}$、$F(t) = 10\cos\omega t\,\mathrm{N}$

固有角振動数: $\omega_n = \sqrt{k/m} = 10\,\text{rad/s}$

減衰比: $\zeta = c/(2\sqrt{mk}) = 0.2$ (不足減衰)

定常状態応答の振幅：

$$X = \frac{F_0/k}{\sqrt{(1-r^2)^2 + (2\zeta r)^2}}$$

$r = \omega/\omega_n = 1$（共振）のとき：

$$X_{\mathrm{res}} = \frac{F_0/k}{2\zeta} = \frac{10/100}{2 \times 0.2} = 0.25\,\text{m}$$

これは共振なし時の静的たわみ $F_0/k = 0.1\,\text{m}$ の**2.5倍**です。これが共振の危険性です。

まとめと次のステップ

この記事で扱った内容：

1. **微分方程式の基礎**: 分類、階数、線形性、解の種類

2. **1階ODEの解法**: 変数分離法、積分因子法、完全ODE、ベルヌーイ方程式

3. **2階線形ODE**: 特性方程式、3つの場合、RLC回路の完全解析

4. **非斉次方程式**: 未定係数法、パラメータ変換法、強制振動

5. **連立ODE**: 行列表現、固有値・固有ベクトルによる解法

6. **ラプラス変換**: 定義、変換対、性質、逆変換、ODEの解法

7. **Python実装**: scipy、sympyによる数値・記号解法

次回は**フーリエ級数・変換と偏微分方程式（PDE）**を扱います。信号処理と電磁気学の核心的な数学ツールです。

参考文献

- Kreyszig, E. "Advanced Engineering Mathematics", 10th Edition, Wiley

- Boyce, W. & DiPrima, R. "Elementary Differential Equations", 11th Edition, Wiley

- Simmons, G. "Differential Equations with Applications and Historical Notes", 3rd Edition

- [SciPy ODEソルバー公式ドキュメント](https://docs.scipy.org/doc/scipy/reference/integrate.html)

- [SymPy ラプラス変換ドキュメント](https://docs.sympy.org/latest/modules/integrals/integrals.html)

練習問題

**答え**: 積分因子法は1階線形ODE dy/dx + P(x)y = Q(x) の標準形を解くために使います。積分因子は mu(x) = exp(P(x)のインテグラル) です。

**解説**: 両辺に mu(x) を掛けると左辺が mu(x)\*y の導関数になり、直接積分できます。これはすべての1階線形ODEに適用できる系統的な手法です。

**答え**: 3つの場合は減衰係数 alpha = R/(2L) と共振周波数 omega_0 = 1/sqrt(LC) の大小関係で決まります：

- 過減衰: alpha が omega_0 より大きい (R が 2\*sqrt(L/C) より大きい) — 振動なく緩やかに減衰

- 臨界減衰: alpha が omega_0 に等しい (R が 2\*sqrt(L/C) に等しい) — 振動なく最速で減衰

- 不足減衰: alpha が omega_0 より小さい (R が 2\*sqrt(L/C) より小さい) — 振動しながら減衰

**解説**: 特性方程式 r^2 + (R/L)r + 1/(LC) = 0 の根 r = -alpha ± sqrt(alpha^2 - omega_0^2) において、判別式の符号が根が実数か複素数かを決定します。

**答え**: 畳み込み定理は、2つの関数の畳み込みのラプラス変換が、それぞれの関数のラプラス変換の積に等しいことを述べます: L\{(f*g)(t)\} = F(s)*G(s)

**解説**: この定理は、線形時不変システムの出力が入力信号とシステムのインパルス応答の畳み込みであることを証明する工学の基本定理です。s領域では単純な乗算になるため、システム解析が大幅に簡単になります。

**答え**: 共振は振動数比 r = omega/omega_n が1になったとき、すなわち駆動周波数がシステムの固有振動数と等しいときに発生します。共振条件では定常状態の振幅が X = (F_0/k) / (2\*zeta) となり、静的たわみの何倍にもなります。

**解説**: 未定係数法で y'' + 4y = 3\*cos(2x) を解くとき、加振周波数が固有振動数と一致するため特殊解に x の係数が現れます (永年項)。実際のシステムでは減衰により振幅は無限大にはなりませんが、共振点で最大値をとります。構造・機械工学では共振回避が重要な設計課題です。

**答え**: dX/dt = AX の形の連立ODEに対して X = v*exp(lambda*t) の形の解を仮定します。代入すると lambda*v = A*v となり、lambda は行列 A の固有値、v は固有ベクトルです。一般解はすべての固有値・固有ベクトル対の線形結合です: X = 合計(C*k * v*k * exp(lambda_k \* t))

**解説**: 固有値の実部が負であればその固有モードは減衰します。すべての固有値の実部が負なら系は安定です。固有ベクトルはシステムの各固有モードの方向を定義します。