はじめに
物理シミュレーションや数理モデルを見ていると、よく微分方程式が出てきます。
たとえば、次のような形です。
\frac{dx}{dt}=f(x)
これは、現在の値 $x$ に対して、次にどちら向きにどれくらい変化するかを表しています。
ただ、式だけを見ても、実際にどのような時間発展になるのかは意外とイメージしにくいです。
実際の微分方程式では、複数の項が組み合わさっていることが多いです。
そのため、いきなり全体を理解しようとすると少し大変ですが、
それぞれの項が単独ではどのような変化をしたがるのか
を少し分解して眺めてみると、式の見え方が変わることがあります。
そこでこの記事では、かなり単純な例として、
\frac{dx}{dt}=0,\quad 1,\quad x,\quad x^2,\quad x^3
を考えます。
目的は、厳密な数学的議論というより、
右辺 $f(x)$ の形が変わると、時間発展の雰囲気がどう変わるのか
を Python で眺めてみることです。
実装コードは Google Colab からもお試しいただけます：
https://colab.research.google.com/drive/1Rc5PTj224bswVqSw2Vzy33qhEVJyNj6g?usp=sharing
今回眺める微分方程式
今回は、次の5つの微分方程式を扱います。
\frac{dx}{dt}=0
\frac{dx}{dt}=1
\frac{dx}{dt}=x
\frac{dx}{dt}=x^2
\frac{dx}{dt}=x^3
初期条件はすべて
x(0)=x_0
とします。
この並びはかなり単純ですが、
- 変化しない
- 一定速度で変化する
- 今の値に比例して変化する
- 非線形に変化する
という違いを見るには、ちょうどよい例です。
解析解
それぞれの解析解は次のようになります。
| 微分方程式 | 解析解 | 雰囲気 |
|---|---|---|
| $\frac{dx}{dt}=0$ | $x(t)=x_0$ | 変化しない |
| $\frac{dx}{dt}=1$ | $x(t)=x_0+t$ | 一定速度で増える |
| $\frac{dx}{dt}=x$ | $x(t)=x_0 e^t$ | 指数関数的に増える |
| $\frac{dx}{dt}=x^2$ | $x(t)=\frac{x_0}{1-x_0t}$ | $x_0>0$ では有限時間で発散する |
| $\frac{dx}{dt}=x^3$ | $x(t)=\frac{x_0}{\sqrt{1-2x_0^2t}}$ | $x_0\neq0$ では有限時間で発散する |
ここで、$\frac{dx}{dt}=x^2$ の場合、$x_0>0$ なら
t=\frac{1}{x_0}
で発散します。
一方で、$x_0<0$ の場合は、正の時間方向では発散せず、$0$ に近づいていきます。
また、$\frac{dx}{dt}=x^3$ の場合は、$x_0\neq0$ なら
t=\frac{1}{2x_0^2}
で発散します。
正の初期値から始めると $+\infty$ へ、負の初期値から始めると $-\infty$ へ向かいます。
つまり、非線形な項では、初期値の大きさだけでなく、符号によっても時間発展の様子が変わることがあります。
同じ初期値で時間発展を比べる
まずは、初期値を $x_0=0.2$ に固定して、それぞれの時間発展を同じ図に描いてみます。
import numpy as np
import matplotlib.pyplot as plt
def sol_zero(t, x0):
return x0 * np.ones_like(t)
def sol_one(t, x0):
return x0 + t
def sol_x(t, x0):
return x0 * np.exp(t)
def sol_x2(t, x0):
"""
dx/dt = x^2
x0 > 0 のとき t = 1/x0 で発散する。
発散後は通常の解として描かない。
"""
denom = 1 - x0 * t
y = x0 / denom
if x0 > 0:
y = np.where(denom > 0, y, np.nan)
return y
def sol_x3(t, x0):
"""
dx/dt = x^3
x0 != 0 のとき t = 1/(2*x0^2) で発散する。
発散後は通常の解として描かない。
"""
denom = 1 - 2 * x0**2 * t
y = np.full_like(t, np.nan, dtype=float)
valid = denom > 0
y[valid] = x0 / np.sqrt(denom[valid])
return y
models = {
r"$dx/dt=0$": sol_zero,
r"$dx/dt=1$": sol_one,
r"$dx/dt=x$": sol_x,
r"$dx/dt=x^2$": sol_x2,
r"$dx/dt=x^3$": sol_x3,
}
x0 = 0.2
t = np.linspace(0, 5, 1000)
fig, ax = plt.subplots(figsize=(8, 5))
for label, func in models.items():
with np.errstate(divide="ignore", invalid="ignore"):
y = func(t, x0)
# 見やすさのため、表示範囲外は消す
y = np.where(np.isfinite(y), y, np.nan)
y = np.where(np.abs(y) < 10, y, np.nan)
ax.plot(t, y, label=label)
ax.axhline(0, color="black", linewidth=0.8)
ax.set_xlabel("t")
ax.set_ylabel(r"$x(t)$")
ax.set_title(fr"Time evolution from $x_0={x0}$")
ax.set_ylim(-0.5, 8)
ax.grid(alpha=0.3)
ax.legend()
plt.tight_layout()
plt.show()
この図を見ると、右辺の形によって時間発展の雰囲気が変わることがわかります。
$\frac{dx}{dt}=0$ では、時間が経っても変化しません。
$\frac{dx}{dt}=1$ では、一定速度で直線的に増えていきます。
$\frac{dx}{dt}=x$ では、値が大きくなるほど増加速度も大きくなるため、指数関数的に増えていきます。
一方で、少し注意したいのが $\frac{dx}{dt}=x^3$ です。
一見すると、$x^3$ は非線形なので激しく変化しそうに見えます。
しかし、初期値が $x_0=0.2$ の場合、
x_0=0.2,\quad x_0^2=0.04,\quad x_0^3=0.008
です。
そのため、初期時刻では $x^3$ の項はかなり小さくなります。
つまり、$x^3$ は常に強い項というより、
$|x|$ が小さいときは弱く、$|x|$ が大きくなると急激に強くなる項
として見るとわかりやすいです。
初期値を変えてみる
次に、初期値を負から正まで変えてみます。
import numpy as np
import matplotlib.pyplot as plt
def sol_zero(t, x0):
return x0 * np.ones_like(t)
def sol_one(t, x0):
return x0 + t
def sol_x(t, x0):
return x0 * np.exp(t)
def sol_x2(t, x0):
"""
dx/dt = x^2
x0 > 0 のとき t = 1/x0 で発散する。
発散後は通常の解として描かない。
"""
denom = 1 - x0 * t
y = x0 / denom
if x0 > 0:
y = np.where(denom > 0, y, np.nan)
return y
def sol_x3(t, x0):
"""
dx/dt = x^3
x0 != 0 のとき t = 1/(2*x0^2) で発散する。
発散後は通常の解として描かない。
"""
denom = 1 - 2 * x0**2 * t
y = np.full_like(t, np.nan, dtype=float)
valid = denom > 0
y[valid] = x0 / np.sqrt(denom[valid])
return y
models = {
r"$dx/dt=0$": sol_zero,
r"$dx/dt=1$": sol_one,
r"$dx/dt=x$": sol_x,
r"$dx/dt=x^2$": sol_x2,
r"$dx/dt=x^3$": sol_x3,
}
x0_list = [-1.0, -0.5, -0.2, 0.2, 0.5, 1.0]
t = np.linspace(0, 5, 1000)
fig, axes = plt.subplots(2, 3, figsize=(15, 8), sharex=True, sharey=True)
axes = axes.ravel()
for ax, x0 in zip(axes, x0_list):
for label, func in models.items():
with np.errstate(divide="ignore", invalid="ignore"):
y = func(t, x0)
# 見やすさのため、表示範囲外は消す
y = np.where(np.isfinite(y), y, np.nan)
y = np.where(np.abs(y) < 10, y, np.nan)
ax.plot(t, y, label=label)
ax.axhline(0, color="black", linewidth=0.8)
ax.axhline(1, color="gray", linewidth=0.8, linestyle="--")
ax.axhline(-1, color="gray", linewidth=0.8, linestyle="--")
ax.set_title(fr"$x_0={x0}$")
ax.set_xlabel("t")
ax.set_ylabel(r"$x(t)$")
ax.set_ylim(-6, 6)
ax.grid(alpha=0.3)
axes[-1].legend(loc="upper right", fontsize=9)
plt.tight_layout()
plt.show()
図中の破線は $x=1$ と $x=-1$ を表しています。
これは、$x$, $x^2$, $x^3$ の大小関係が切り替わる目安です。
たとえば、
|x|<1
の範囲では、高次の項ほど小さくなります。
具体的には、$x=0.2$ のとき、
|x|=0.2,\quad x^2=0.04,\quad |x^3|=0.008
です。
そのため、小さい値の近くでは、$x^3$ のような高次の項はかなりおとなしく見えます。
一方で、
|x|>1
では、高次の項ほど大きくなります。
つまり、非線形項は、小さい値の近くでは静かですが、$|x|$ が $1$ を超えたあたりから急に効きやすくなります。
また、負の初期値を見ると、$x^2$ と $x^3$ の違いもわかりやすくなります。
\frac{dx}{dt}=x^2
では、右辺は常に $0$ 以上です。
そのため、$x<0$ から始めても、$x$ は正の方向へ進みます。
ただし、負の初期値から始めた場合、すぐに正の値へ飛ぶわけではありません。
$\frac{dx}{dt}=x^2$ の解は、負の値から始めると、時間とともに $0$ に近づいていきます。
一方で、
\frac{dx}{dt}=x^3
では、$x<0$ のとき右辺も負になります。
そのため、負の値から始めると、さらに負の方向へ進み、有限時間で $-\infty$ に発散します。
このように、同じ非線形項でも、偶数乗か奇数乗かによって時間発展の向きが変わります。
かなり単純な例ですが、
- $x^2$ は常に正方向へ押す
- $x^3$ は符号を保ったまま、正なら正へ、負なら負へ押す
という見方ができます。
発散後の枝について
ここで、可視化するときに少し注意が必要です。
たとえば、$\frac{dx}{dt}=x^2$ の解析解は
x(t)=\frac{x_0}{1-x_0t}
です。
$x_0>0$ の場合、$t=1/x_0$ で分母が $0$ になり、解は $+\infty$ に発散します。
この時刻を越えて解析式をそのまま評価すると、分母が負になるため、グラフ上では負側に別の枝が現れます。
たとえば $x_0=1$ なら、
x(t)=\frac{1}{1-t}
なので、$t=1$ を越えると式の値は負になります。
しかし、これは発散後の式を機械的に描いているだけです。
通常の意味での時間発展が、$+\infty$ に飛んだあとに $-\infty$ から戻ってくる、という意味ではありません。
そのため、この記事のプロットでは、発散時刻を越えた部分は表示しないようにしています。
このような処理を入れないと、有限時間発散を持つ解では、見かけ上の枝が出てしまうことがあります。
まとめ
この記事では、
\frac{dx}{dt}=0,\quad 1,\quad x,\quad x^2,\quad x^3
という単純な微分方程式を例に、解析解と時間発展を可視化しました。
今回の例からわかるように、右辺 $f(x)$ の形が変わるだけで、時間発展の雰囲気はかなり変わります。
特に、
- $\frac{dx}{dt}=0$ は変化しない
- $\frac{dx}{dt}=1$ は一定速度で変化する
- $\frac{dx}{dt}=x$ は指数関数的に変化する
- $\frac{dx}{dt}=x^2$ は常に正方向へ押す
- $\frac{dx}{dt}=x^3$ は符号を保ったまま外側へ押す
- 高次の項は $|x|<1$ では弱く、$|x|>1$ で急に効きやすい
- 有限時間で発散する解では、解析式をそのまま描くと、見かけ上の枝が現れることがある
という雰囲気が見えました。
実際の物理シミュレーションや数理モデルでは、もっと複雑な項が組み合わさって現れます。
しかし、複雑な方程式をいきなり全体として理解しようとする前に、
それぞれの項が、単独ではどのような変化を好むのか
を眺めておくと、式の読み方が少し楽になります。
今回はかなり素朴な例でしたが、微分方程式に馴染むための第一歩として、こういう単純な形を可視化して眺めてみました。
微分方程式の項の雰囲気を掴むときの参考になれば嬉しいです。