非斉次方程式は、外力や入力が加わったシステムの応答を記述します。例えば、強制振動(外力付きバネ系)、外部電圧源を持つ回路、制御入力を持つプラントなど、実際の工学問題の大半は非斉次です。
本記事では、非斉次線形ODEの一般解の構造を理解し、未定係数法と定数変化法の2つの系統的な解法を習得します。
本記事の内容
- 非斉次方程式の一般解 = 斉次一般解 + 特殊解
- 未定係数法(入力が多項式・指数・三角関数の場合)
- 定数変化法(一般的な手法)
- ロンスキアン
- Pythonでの実装
前提知識
一般解の構造
非斉次方程式
$$ ay” + by’ + cy = g(x) $$
の一般解は、次の2つの和です。
$$ \boxed{y = y_h + y_p} $$
- $y_h$: 斉次方程式 $ay” + by’ + cy = 0$ の一般解(2つの任意定数を含む)
- $y_p$: 非斉次方程式の特殊解(任意定数を含まない、1つの解)
未定係数法
$g(x)$ が多項式・指数関数・三角関数またはそれらの積の場合に有効な手法です。
ルール
| $g(x)$ の形 | $y_p$ の仮定形 |
|---|---|
| $P_n(x)$($n$ 次多項式) | $A_n x^n + A_{n-1} x^{n-1} + \cdots + A_0$ |
| $e^{rx}$ | $Ae^{rx}$ |
| $\cos\omega x$ or $\sin\omega x$ | $A\cos\omega x + B\sin\omega x$ |
| $P_n(x)e^{rx}$ | $(A_n x^n + \cdots + A_0)e^{rx}$ |
| $e^{rx}\cos\omega x$ | $e^{rx}(A\cos\omega x + B\sin\omega x)$ |
修正ルール: 仮定した $y_p$ が斉次解 $y_h$ と重複する場合、$x^s$ を掛けます($s$ は重複の重複度)。
例: $y” + 3y’ + 2y = 10e^{3x}$
斉次解: $y_h = C_1 e^{-x} + C_2 e^{-2x}$
$y_p = Ae^{3x}$ と仮定($\lambda = 3$ は特性根でないので修正不要)。
$$ \begin{align} y_p’ &= 3Ae^{3x}, \quad y_p” = 9Ae^{3x} \\ 9Ae^{3x} + 9Ae^{3x} + 2Ae^{3x} &= 10e^{3x} \\ 20A &= 10 \\ A &= \frac{1}{2} \end{align} $$
$$ y = C_1 e^{-x} + C_2 e^{-2x} + \frac{1}{2}e^{3x} $$
例: $y” + y = \sin x$(共鳴の例)
斉次解: $y_h = C_1\cos x + C_2\sin x$
$\sin x$ は斉次解に含まれるので、$y_p = x(A\cos x + B\sin x)$ と修正。
代入して係数を比較すると $A = -1/2, B = 0$。
$$ y_p = -\frac{x}{2}\cos x $$
定数変化法(パラメータ変化法)
$g(x)$ がどのような関数でも適用できる一般的な手法です。
斉次解を $y_h = C_1 y_1 + C_2 y_2$ として、$y_p = u_1(x) y_1 + u_2(x) y_2$ と仮定します。
$u_1, u_2$ は次の連立方程式を解いて求めます。
$$ \begin{cases} u_1′ y_1 + u_2′ y_2 = 0 \\ u_1′ y_1′ + u_2′ y_2′ = \frac{g(x)}{a} \end{cases} $$
ロンスキアン $W$ を用いた公式:
$$ W = \begin{vmatrix} y_1 & y_2 \\ y_1′ & y_2′ \end{vmatrix} = y_1 y_2′ – y_2 y_1′ $$
$$ u_1′ = -\frac{y_2 g(x)}{aW}, \quad u_2′ = \frac{y_1 g(x)}{aW} $$
積分して $u_1, u_2$ を求め、$y_p = u_1 y_1 + u_2 y_2$ とします。
例: $y” + y = \sec x$
$y_1 = \cos x, y_2 = \sin x$, $W = \cos^2 x + \sin^2 x = 1$
$$ \begin{align} u_1′ &= -\sin x \cdot \sec x = -\tan x \\ u_2′ &= \cos x \cdot \sec x = 1 \end{align} $$
$$ \begin{align} u_1 &= \int -\tan x \, dx = \ln|\cos x| \\ u_2 &= \int 1 \, dx = x \end{align} $$
$$ y_p = \cos x \ln|\cos x| + x\sin x $$
Pythonでの実装
import numpy as np
import matplotlib.pyplot as plt
from scipy.integrate import solve_ivp
x = np.linspace(0, 10, 500)
# 例1: y'' + 3y' + 2y = 10e^{3x}, y(0)=0, y'(0)=0
def ode1(t, y):
return [y[1], -3*y[1] - 2*y[0] + 10*np.exp(3*t)]
sol1 = solve_ivp(ode1, [0, 2], [0, 0], t_eval=np.linspace(0, 2, 300))
# 解析解: y = C1*e^{-x} + C2*e^{-2x} + 0.5*e^{3x}
# y(0)=0: C1+C2+0.5=0, y'(0)=0: -C1-2C2+1.5=0 → C1=-0.5, C2=0.5... wait
# Actually let me just compute: C1+C2=-0.5, -C1-2C2=-1.5 → C1=0.5, C2=-1
t1 = sol1.t
y_exact1 = 0.5*np.exp(-t1) - np.exp(-2*t1) + 0.5*np.exp(3*t1)
# 例2: y'' + y = sin(x) (共鳴), y(0)=0, y'(0)=0
def ode2(t, y):
return [y[1], -y[0] + np.sin(t)]
sol2 = solve_ivp(ode2, [0, 30], [0, 0], t_eval=np.linspace(0, 30, 500))
# 例3: y'' + y = sec(x), y(0)=1, y'(0)=0
def ode3(t, y):
return [y[1], -y[0] + 1/np.cos(t)]
sol3 = solve_ivp(ode3, [0, 1.4], [1, 0], t_eval=np.linspace(0, 1.4, 300))
fig, axes = plt.subplots(1, 3, figsize=(15, 4))
axes[0].plot(sol1.t, sol1.y[0], 'b-', linewidth=2, label='数値解')
axes[0].plot(t1, y_exact1, 'r--', label='解析解')
axes[0].set_title("$y''+3y'+2y=10e^{3x}$\n(未定係数法)")
axes[0].set_xlabel('x')
axes[0].set_ylabel('y')
axes[0].legend()
axes[0].grid(True, alpha=0.3)
axes[1].plot(sol2.t, sol2.y[0], 'b-', linewidth=2)
axes[1].set_title("$y''+y=\\sin x$\n(共鳴: 振幅増大)")
axes[1].set_xlabel('x')
axes[1].set_ylabel('y')
axes[1].grid(True, alpha=0.3)
axes[2].plot(sol3.t, sol3.y[0], 'b-', linewidth=2)
axes[2].set_title("$y''+y=\\sec x$\n(定数変化法)")
axes[2].set_xlabel('x')
axes[2].set_ylabel('y')
axes[2].grid(True, alpha=0.3)
plt.tight_layout()
plt.savefig('nonhomogeneous_ode.png', dpi=150, bbox_inches='tight')
plt.show()
共鳴の例($y” + y = \sin x$)では、入力の周波数が系の固有振動数と一致するため、振幅が時間とともに線形に増大していく様子が確認できます。
まとめ
- 非斉次ODEの一般解は $y = y_h + y_p$(斉次一般解 + 特殊解)
- 未定係数法: $g(x)$ が多項式・指数・三角関数のとき有効、仮定形を代入して係数決定
- 斉次解と重複する場合は $x^s$ を掛けて修正
- 定数変化法: 任意の $g(x)$ に適用可能、ロンスキアンを用いて公式的に計算
次のステップとして、以下の記事も参考にしてください。