1階線形常微分方程式は、物理学・工学で最も頻繁に登場する微分方程式の1つです。RC回路の過渡応答、放射性崩壊、ニュートンの冷却法則、薬物動態モデルなど、「ある量の変化率が、その量自体と外部からの入力に比例する」現象はすべてこの形式で記述されます。
本記事では、1階線形常微分方程式の一般的な解法である 積分因子法（integrating factor method） を、動機づけから導出まで省略せずに解説します。さらに、具体的な物理応用とPython実装を通じて理解を深めます。
本記事の内容
- 1階線形ODEの標準形と分類
- 積分因子の動機づけと導出
- 一般解の公式の証明
- 定数係数・変数係数の具体例
- RC回路の過渡応答とニュートンの冷却法則への応用
- Pythonによる解析解と数値解の比較
前提知識
この記事を読む前に、以下の記事を読んでおくと理解が深まります。
前提記事では、微分方程式の基本概念（階数、線形性、一般解と特殊解、初期値問題）と変数分離形の解法を扱いました。本記事ではそれらの知識を前提として、より一般的な1階線形ODEの解法に進みます。
1階線形ODEの標準形
1階線形常微分方程式の 標準形 は次のように書かれます。
$$ \frac{dy}{dx} + P(x)\,y = Q(x) $$
ここで、$P(x)$ と $Q(x)$ は $x$ の既知の関数（連続関数を仮定）であり、未知関数 $y(x)$ とその1階導関数 $dy/dx$ が1次の項としてのみ現れています。このことが「線形」と呼ばれる理由です。
いくつかの例を見てみましょう。
| 方程式 | $P(x)$ | $Q(x)$ | 備考 |
|---|---|---|---|
| $\dfrac{dy}{dx} + 2y = 6$ | $2$ | $6$ | 定数係数 |
| $\dfrac{dy}{dx} + \dfrac{1}{x}y = x^2$ | $\dfrac{1}{x}$ | $x^2$ | 変数係数 |
| $\dfrac{dy}{dx} – 3y = e^x$ | $-3$ | $e^x$ | 定数係数、非斉次 |
$Q(x) = 0$ の場合を 斉次（homogeneous） 方程式、$Q(x) \neq 0$ の場合を 非斉次（nonhomogeneous） 方程式と呼びます。斉次の場合は変数分離法で解けますが、非斉次の場合は一般に変数分離法だけでは解けません。ここで登場するのが 積分因子法 です。
積分因子の動機づけ — なぜ積分因子が必要か
標準形の方程式
$$ \frac{dy}{dx} + P(x)\,y = Q(x) $$
の左辺をよく見てみましょう。もしこの左辺が何かの関数の微分（全微分）の形に書けたら、両辺を積分するだけで解が求まるはずです。
ここで思い出したいのが 積の微分法則 です。
$$ \frac{d}{dx}\bigl[\mu(x)\,y\bigr] = \mu(x)\,\frac{dy}{dx} + \mu'(x)\,y $$
もし適切な関数 $\mu(x)$ を見つけて、元の方程式の両辺に $\mu(x)$ を掛けたとき、左辺がちょうど $\dfrac{d}{dx}\bigl[\mu(x)\,y\bigr]$ の形になれば、解が簡単に求まります。このような $\mu(x)$ を 積分因子（integrating factor） と呼びます。
具体的に、元の方程式の両辺に $\mu(x)$ を掛けると
$$ \mu(x)\,\frac{dy}{dx} + \mu(x)\,P(x)\,y = \mu(x)\,Q(x) $$
となります。左辺が $\dfrac{d}{dx}\bigl[\mu(x)\,y\bigr]$ に等しくなるためには、積の微分法則と比較して
$$ \mu(x)\,\frac{dy}{dx} + \mu'(x)\,y = \mu(x)\,\frac{dy}{dx} + \mu(x)\,P(x)\,y $$
が成り立つ必要があります。両辺から $\mu(x)\,\dfrac{dy}{dx}$ を消去すると
$$ \mu'(x) = \mu(x)\,P(x) $$
という $\mu(x)$ に関する条件式が得られます。これが積分因子の満たすべき微分方程式です。
積分因子 $\mu(x)$ の導出
$\mu(x)$ が満たすべき方程式
$$ \frac{d\mu}{dx} = \mu(x)\,P(x) $$
は、実は $\mu$ に関する 変数分離形 の微分方程式です。前提記事で学んだ変数分離法を使って解きましょう。
ステップ1: 変数を分離します。
$$ \frac{d\mu}{\mu} = P(x)\,dx $$
ステップ2: 両辺を積分します。
$$ \int \frac{d\mu}{\mu} = \int P(x)\,dx $$
ステップ3: 左辺を計算します。
$$ \ln|\mu| = \int P(x)\,dx $$
ここで積分定数は、後で見るように最終結果に影響しないため省略できます（理由は後述）。
ステップ4: 両辺の指数関数を取ります。
$$ |\mu| = \exp\!\left(\int P(x)\,dx\right) $$
$\mu(x)$ は方程式の両辺に掛ける因子なので、正の値を取れば十分です。したがって
$$ \boxed{\mu(x) = \exp\!\left(\int P(x)\,dx\right)} $$
これが 積分因子 の公式です。
積分定数を省略できる理由: 仮に積分定数 $C_0$ を含めると $\mu(x) = e^{C_0}\exp\!\left(\int P(x)\,dx\right)$ となりますが、この $e^{C_0}$ は方程式の両辺に共通して掛かるだけなので、最終的な解には影響しません。したがって $C_0 = 0$ として問題ありません。
一般解の公式の導出
積分因子 $\mu(x) = \exp\!\left(\int P(x)\,dx\right)$ を用いて、1階線形ODEの一般解を導出します。
ステップ1: 元の方程式の両辺に $\mu(x)$ を掛けます。
$$ \mu(x)\,\frac{dy}{dx} + \mu(x)\,P(x)\,y = \mu(x)\,Q(x) $$
ステップ2: 左辺は、積分因子の定義から積の微分に書き換えられます。
$$ \frac{d}{dx}\bigl[\mu(x)\,y\bigr] = \mu(x)\,Q(x) $$
ステップ3: 両辺を $x$ で積分します。
$$ \mu(x)\,y = \int \mu(x)\,Q(x)\,dx + C $$
ここで $C$ は積分定数（任意定数）です。
ステップ4: 両辺を $\mu(x)$ で割ります。$\mu(x) = e^{\int P(x)\,dx} > 0$ なので常に割れます。
$$ \boxed{y = \frac{1}{\mu(x)}\left[\int \mu(x)\,Q(x)\,dx + C\right]} $$
これが1階線形常微分方程式の 一般解 です。$\mu(x) = \exp\!\left(\int P(x)\,dx\right)$ を代入すれば
$$ y = e^{-\int P(x)\,dx}\left[\int e^{\int P(x)\,dx}\,Q(x)\,dx + C\right] $$
と書くこともできます。1階のODEなので任意定数は1つ（$C$）であり、初期条件 $y(x_0) = y_0$ を1つ与えれば特殊解が確定します。
具体例1: 定数係数の場合
次の初期値問題を解きます。
$$ \frac{dy}{dx} + 2y = 6, \quad y(0) = 1 $$
ステップ1: 標準形との対応を確認します。
$$ P(x) = 2, \quad Q(x) = 6 $$
ステップ2: 積分因子を計算します。
$$ \mu(x) = \exp\!\left(\int 2\,dx\right) = \exp(2x) = e^{2x} $$
ステップ3: 両辺に $\mu(x) = e^{2x}$ を掛けます。
$$ e^{2x}\,\frac{dy}{dx} + 2e^{2x}\,y = 6e^{2x} $$
ステップ4: 左辺を積の微分にまとめます。
$$ \frac{d}{dx}\bigl[e^{2x}\,y\bigr] = 6e^{2x} $$
ステップ5: 両辺を積分します。
$$ e^{2x}\,y = \int 6e^{2x}\,dx = 6 \cdot \frac{e^{2x}}{2} + C = 3e^{2x} + C $$
ステップ6: 両辺を $e^{2x}$ で割ります。
$$ y = 3 + Ce^{-2x} $$
ステップ7: 初期条件 $y(0) = 1$ を適用します。
$$ 1 = 3 + C \cdot e^{0} = 3 + C $$
$$ C = -2 $$
したがって、特殊解は
$$ \boxed{y(x) = 3 – 2e^{-2x}} $$
です。$x \to \infty$ で $y \to 3$ に収束することが分かります。これは 定常解（斉次でない部分の特殊解 $y_p = 3$）に向かって指数関数的に近づく過渡応答です。
具体例2: 変数係数の場合
次の方程式を解きます。
$$ \frac{dy}{dx} + \frac{1}{x}\,y = x^2, \quad (x > 0) $$
ステップ1: 標準形との対応を確認します。
$$ P(x) = \frac{1}{x}, \quad Q(x) = x^2 $$
ステップ2: 積分因子を計算します。
$$ \mu(x) = \exp\!\left(\int \frac{1}{x}\,dx\right) = \exp(\ln x) = x $$
ここで $x > 0$ を仮定しているため、$\ln|x| = \ln x$ です。
ステップ3: 両辺に $\mu(x) = x$ を掛けます。
$$ x\,\frac{dy}{dx} + y = x^3 $$
ステップ4: 左辺を積の微分にまとめます。
$$ \frac{d}{dx}[x\,y] = x^3 $$
実際に左辺を確認すると $\dfrac{d}{dx}[xy] = x\dfrac{dy}{dx} + y$ であり、確かに一致しています。
ステップ5: 両辺を積分します。
$$ x\,y = \int x^3\,dx = \frac{x^4}{4} + C $$
ステップ6: 両辺を $x$ で割ります。
$$ \boxed{y = \frac{x^3}{4} + \frac{C}{x}} $$
これが一般解です。$C/x$ の項は $x \to 0^+$ で発散するため、初期条件の値によって解の振る舞いが大きく変わります。
物理への応用1: RC回路の過渡応答
抵抗 $R$ とキャパシタ $C_{\text{cap}}$ が直列に接続された回路に、時刻 $t = 0$ で一定電圧 $V_0$ を印加する場合を考えます。キルヒホッフの電圧則（KVL）から
$$ R\,i(t) + \frac{1}{C_{\text{cap}}}\int_0^t i(\tau)\,d\tau = V_0 $$
が成り立ちます。キャパシタ両端の電圧 $v(t) = \dfrac{1}{C_{\text{cap}}}\displaystyle\int_0^t i(\tau)\,d\tau$ と置き、$i(t) = C_{\text{cap}}\,\dfrac{dv}{dt}$ を代入すると
$$ RC_{\text{cap}}\,\frac{dv}{dt} + v = V_0 $$
が得られます。時定数 $\tau = RC_{\text{cap}}$ を用いて書き直すと
$$ \frac{dv}{dt} + \frac{1}{\tau}\,v = \frac{V_0}{\tau} $$
これは標準形であり、$P(t) = 1/\tau$、$Q(t) = V_0/\tau$ です。
積分因子の計算:
$$ \mu(t) = \exp\!\left(\int \frac{1}{\tau}\,dt\right) = e^{t/\tau} $$
方程式の両辺に $\mu(t)$ を掛ける:
$$ \frac{d}{dt}\bigl[e^{t/\tau}\,v\bigr] = \frac{V_0}{\tau}\,e^{t/\tau} $$
両辺を積分:
$$ e^{t/\tau}\,v = \frac{V_0}{\tau} \cdot \tau\,e^{t/\tau} + C = V_0\,e^{t/\tau} + C $$
$e^{t/\tau}$ で割る:
$$ v(t) = V_0 + Ce^{-t/\tau} $$
初期条件 $v(0) = 0$（キャパシタが初期放電状態）を適用:
$$ 0 = V_0 + C \implies C = -V_0 $$
したがって
$$ \boxed{v(t) = V_0\!\left(1 – e^{-t/\tau}\right)} $$
$t = \tau$ のとき $v(\tau) = V_0(1 – e^{-1}) \approx 0.632\,V_0$ であり、時定数 $\tau$ は電圧が最終値の約63.2%に達するまでの時間を表します。
物理への応用2: ニュートンの冷却法則
温度 $T(t)$ の物体が周囲温度 $T_{\text{env}}$ の環境に置かれたとき、ニュートンの冷却法則は
$$ \frac{dT}{dt} = -k(T – T_{\text{env}}) $$
と表されます。ここで $k > 0$ は冷却定数（物体の材質・形状・周囲の流体に依存）です。これを展開すると
$$ \frac{dT}{dt} + kT = kT_{\text{env}} $$
これも1階線形ODEの標準形であり、$P(t) = k$、$Q(t) = kT_{\text{env}}$ です。
積分因子の計算:
$$ \mu(t) = e^{kt} $$
一般解の導出:
$$ \frac{d}{dt}\bigl[e^{kt}\,T\bigr] = kT_{\text{env}}\,e^{kt} $$
$$ e^{kt}\,T = kT_{\text{env}} \cdot \frac{e^{kt}}{k} + C = T_{\text{env}}\,e^{kt} + C $$
$$ T(t) = T_{\text{env}} + Ce^{-kt} $$
初期条件 $T(0) = T_0$ を適用:
$$ T_0 = T_{\text{env}} + C \implies C = T_0 – T_{\text{env}} $$
したがって
$$ \boxed{T(t) = T_{\text{env}} + (T_0 – T_{\text{env}})\,e^{-kt}} $$
物体の温度は、初期温度 $T_0$ から周囲温度 $T_{\text{env}}$ に向かって指数関数的に近づきます。$T_0 > T_{\text{env}}$（物体が熱い）でも $T_0 < T_{\text{env}}$（物体が冷たい）でも同じ公式が使えます。
Python実装: RC回路の解析解と数値解の比較
RC回路の過渡応答を題材に、積分因子法で求めた解析解と scipy.integrate.solve_ivp による数値解を比較します。あわせて、ニュートンの冷却法則も可視化します。
import numpy as np
import matplotlib.pyplot as plt
from scipy.integrate import solve_ivp
# ============================================================
# 1. RC回路の過渡応答
# ============================================================
# パラメータ
R = 1e3 # 抵抗 [Ω]
C_cap = 1e-3 # キャパシタ [F]
V0 = 5.0 # 印加電圧 [V]
tau = R * C_cap # 時定数 [s]
# 微分方程式: dv/dt = (V0 - v) / tau
def rc_circuit(t, v):
return (V0 - v[0]) / tau
# 時間の設定
t_span = (0, 5 * tau) # 5τまで計算
t_eval = np.linspace(0, 5 * tau, 500)
# 数値解（solve_ivp）
sol_rc = solve_ivp(
rc_circuit,
t_span,
[0.0], # 初期条件 v(0) = 0
t_eval=t_eval,
method="RK45"
)
# 解析解
v_exact = V0 * (1 - np.exp(-t_eval / tau))
# ============================================================
# 2. ニュートンの冷却法則
# ============================================================
# パラメータ
T0 = 90.0 # 初期温度 [°C]（熱いコーヒー）
T_env = 22.0 # 周囲温度 [°C]
k_cool = 0.1 # 冷却定数 [1/min]
# 微分方程式: dT/dt = -k(T - T_env)
def cooling_law(t, T):
return -k_cool * (T[0] - T_env)
# 時間の設定
t_cool_span = (0, 60) # 60分
t_cool_eval = np.linspace(0, 60, 500)
# 数値解
sol_cool = solve_ivp(
cooling_law,
t_cool_span,
[T0],
t_eval=t_cool_eval,
method="RK45"
)
# 解析解
T_exact = T_env + (T0 - T_env) * np.exp(-k_cool * t_cool_eval)
# ============================================================
# 3. 可視化
# ============================================================
fig, axes = plt.subplots(1, 2, figsize=(14, 5))
# --- RC回路 ---
ax1 = axes[0]
ax1.plot(t_eval * 1000, v_exact, "-", linewidth=2.5, label="Analytical solution", color="royalblue")
ax1.plot(sol_rc.t * 1000, sol_rc.y[0], "o", markersize=3, label="Numerical (solve_ivp)", color="tomato", alpha=0.7)
ax1.axhline(y=V0, color="gray", linestyle="--", alpha=0.5, label=f"$V_0 = {V0}$ V")
ax1.axvline(x=tau * 1000, color="green", linestyle="--", alpha=0.5, label=f"$\\tau = {tau*1000:.0f}$ ms")
ax1.set_xlabel("Time [ms]", fontsize=12)
ax1.set_ylabel("Voltage $v(t)$ [V]", fontsize=12)
ax1.set_title("RC Circuit Transient Response", fontsize=13)
ax1.legend(fontsize=10)
ax1.grid(True, alpha=0.3)
ax1.set_ylim(-0.3, V0 * 1.15)
# --- ニュートンの冷却法則 ---
ax2 = axes[1]
ax2.plot(t_cool_eval, T_exact, "-", linewidth=2.5, label="Analytical solution", color="royalblue")
ax2.plot(sol_cool.t, sol_cool.y[0], "o", markersize=3, label="Numerical (solve_ivp)", color="tomato", alpha=0.7)
ax2.axhline(y=T_env, color="gray", linestyle="--", alpha=0.5, label=f"$T_{{\\mathrm{{env}}}} = {T_env}$°C")
ax2.set_xlabel("Time [min]", fontsize=12)
ax2.set_ylabel("Temperature $T(t)$ [°C]", fontsize=12)
ax2.set_title("Newton's Law of Cooling", fontsize=13)
ax2.legend(fontsize=10)
ax2.grid(True, alpha=0.3)
plt.tight_layout()
plt.savefig("first_order_linear_ode.png", dpi=150, bbox_inches="tight")
plt.show()
# 誤差の確認
print("=== 数値解と解析解の最大誤差 ===")
print(f"RC回路: {np.max(np.abs(sol_rc.y[0] - v_exact)):.2e} V")
print(f"ニュートン冷却: {np.max(np.abs(sol_cool.y[0] - T_exact)):.2e} °C")
このコードを実行すると、左のグラフにRC回路のキャパシタ電圧の時間変化が、右のグラフにコーヒーの温度変化が描画されます。どちらのグラフでも、解析解（実線）と solve_ivp による数値解（マーカー）がほぼ完全に一致することが確認できます。
RC回路では時定数 $\tau = RC = 1\,\text{s}$ の約5倍の時間でほぼ定常状態 $V_0 = 5\,\text{V}$ に到達します。ニュートンの冷却法則では、90°Cのコーヒーが約60分で周囲温度22°Cに近づいていく様子が見て取れます。
まとめ
本記事では、1階線形常微分方程式の解法である 積分因子法 を解説しました。
- 1階線形ODEの標準形は $\dfrac{dy}{dx} + P(x)\,y = Q(x)$ です。
- 積分因子 $\mu(x) = \exp\!\left(\int P(x)\,dx\right)$ を掛けることで、左辺を $\dfrac{d}{dx}[\mu y]$ という積の微分の形にまとめられます。
- 一般解は $y = \dfrac{1}{\mu(x)}\left[\displaystyle\int \mu(x)\,Q(x)\,dx + C\right]$ で与えられます。
- 積分因子法は、RC回路の過渡応答やニュートンの冷却法則など、物理・工学の幅広い問題に適用できます。
- 解析解は指数関数的な過渡応答の形を取り、定常値に向かって収束するパターンが共通して現れます。
積分因子法をマスターすれば、1階線形ODEは体系的に解くことができます。次のステップとして、以下のトピックに進んでいくとよいでしょう。
- 変数分離形 — 非線形だが分離可能な微分方程式の解法
- 2階定数係数線形ODE — 特性方程式を用いた解法（振動・減衰系の記述）
- ラプラス変換 — 初期値問題を代数的に解くための強力な手法