非斉次方程式の解法(未定係数法・定数変化法)

非斉次方程式は、外力や入力が加わったシステムの応答を記述します。例えば、強制振動(外力付きバネ系)、外部電圧源を持つ回路、制御入力を持つプラントなど、実際の工学問題の大半は非斉次です。

本記事では、非斉次線形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)$ に適用可能、ロンスキアンを用いて公式的に計算

次のステップとして、以下の記事も参考にしてください。