工業数学完全攻略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)$

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

物理的意味のまとめ:

3.3 非斉次方程式

形式: $ay'' + by' + cy = f(x)$

完全解: $y = y_h + y_p$

• $y_h$: 斉次方程式の一般解（余関数）
• $y_p$: 特殊解

未定係数法

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

注意: $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 主要変換対

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の活用

import numpy as np
from scipy.integrate import solve_ivp
import matplotlib.pyplot as plt

# 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で解析解を求める

import sympy as sp

# シンボル定義
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 各種積分法の比較

import numpy as np
import matplotlib.pyplot as plt

# オイラー法 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で自動化

import sympy as sp

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ソルバー公式ドキュメント
• SymPy ラプラス変換ドキュメント

練習問題