変数分離法

変数分離形とは，微分方程式の関数 \(f(x,y)\)が\(x\)と\(y\)のそれぞれの関数\(g(x),h(y)\)に分離できるときに利用できる解法である．すなわち，

\[ \begin{align} \frac{dy}{dx} = g(x)h(y) \end{align} \]

に式変形できるときに利用される方法である．変数分離形は常微分方程式の解法の中でも最も基本的で強力な解法の一つであるので，実際に計算し，身につけてほしい．

変数分離形の解法

与えられた一階常微分方程式が上記のような形状に式変形可能であるとき，次のステップを踏むことで未知関数 \(y=y(x)\) を計算できる．

まず，前述した分離形の微分方程式の両辺に\(dx\)をかけると以下のようになる

\[ \begin{align} dy = g(x)h(y)dx \end{align} \]

続いて，両辺を\(h(y)\)で割ると以下のようになる．

\[ \begin{align} \frac{1}{h(y)}dy = g(x)dx \end{align} \]

式を見てわかるように，\(x\)と\(y\)の項がそれぞれ両辺に移項できたことがわかる．最後に両辺を積分する．

\[ \begin{align} \int\frac{1}{h(y)}dy = \int g(x)dx \end{align} \]

この積分を実行することで微分方程式の解である未知関数 \(y(x)\)を求めることができる．

変数分離形の具体例

では，次の微分方程式を変数分離形で解いてみよう．以下の微分方程式は \(g(x)=1, h(y)=y\)のときの場合である．

\[ \begin{align} \frac{dy}{dx} = y \end{align} \]

解答はクリックで確認できる．

まず，左辺に \(y\) のみ，右辺に \(x\) の項のみとなるように整理すると，

\[ \begin{align} \frac{dy}{y} = dx \end{align} \]

となる．続いて，両辺を積分することを考える．

\[ \begin{align} \int \frac{dy}{y} = \int dx \end{align} \]

より，積分すると以下となる．

\[ \begin{align} \log |y| = x + C' \end{align} \]

ここで本来ならば，左辺の積分定数 \(C_1\) と右辺の積分定数 \(C_2\) が存在するがまとめて \(C'=C_1 + C_2\) としている．対数の定義より，\(\log |y| = x \rightarrow e^x = |y|\) であり，絶対値を外すと，\(y=\pm e^x\) となる．よって一般解は，

\[\begin{split} \begin{align} |y| &= e^{x + C'} \\ y &= \pm e^{C'} e^x \\ y & = C e^x \end{align} \end{split}\]

となる．ただし，\(C=\pm e^{C'}\) としている．以上より，与えられた微分方程式を変数分離形で解くことができた．

Pythonによる実装

では，上記の微分方程式をsympyで解く．単一のセルで作成したが，各行での動作は前回の講義のようにStep-by-stepで確認されたい．

from sympy import symbols, Eq, Derivative, Function, dsolve

# 変数と未知関数の定義
x = symbols('x')
y = Function('y')(x)

# dy/dxの定義
dy = Derivative(y, x)

# 微分方程式の定義
eq = Eq(dy, y)

# 微分方程式を解く
y_ = dsolve(eq, y)
y_

\[\displaystyle y{\left(x \right)} = C_{1} e^{x}\]

手計算と結果が一致することを確認されたい．最後に解集合を可視化する．

import numpy as np
from sympy import plotting
import matplotlib.pyplot as plt

N = 10

eqs = []
for c in np.linspace(-3, 3, N):
    eqs.append(y_.rhs.subs(symbols('C1'), c))
p = plotting.plot(*eqs, (x, -2, 2), show=False)

cm = plt.get_cmap('magma', N)
for i in range(cm.N):
    p[i].line_color = cm(i)
p.show()

置換による変数分離形への帰着

微分方程式によっては置換により変数分離形に帰着できる形がある．具体的に，

\[ \frac{dy}{dx} = f(y + ax + b) \]

のような微分方程式を考える．このような微分方程式では直接的に変数分離することはできない．

そこで，\(u = y + ax + b\) で置換する．このとき

\[ \frac{du}{dx} = \frac{dy}{dx} + a \]

となる．これを与えられた微分方程式に代入すると

\[\begin{split} \begin{align} \frac{du}{dx} - a &= f(u) \\ \frac{du}{dx} &= f(u) + a \\ \end{align} \end{split}\]

となり，変数分離形に帰着できる．

そして，

\[\begin{split} \begin{align} \frac{du}{dx} - a &= f(u) \\ \frac{1}{f(u) - a} du &= dx \\ \end{align} \end{split}\]

となり，両辺を積分すると

\[\begin{split} \begin{align} \int \frac{1}{f(u) - a} du &= \int dx \\ \int \frac{1}{f(u) - a} du &= x + C \\ \end{align} \end{split}\]

で一般解を導出できる．ただし，\(C\)は任意定数とする．

置換による変数分離形の具体例

次の微分方程式を解く．

\[ \frac{dy}{dx} = \frac{1}{y+x} - 1 \]

解答はクリックで確認できる．

ここで．\(u=y+x\)として一般解を求める．

このとき，両辺を \(x\) で微分したときの導関数は次のようになる．

\[ \frac{du}{dx} = \frac{dy}{dx} + 1 \]

これを与えられた微分方程式に代入し，変形すると．

\[\begin{split} \begin{align} \frac{du}{dx} - 1 &= \frac{1}{u} - 1 \\ \frac{du}{dx} &= \frac{1}{u} \end{align} \end{split}\]

これを変数分離形として解く．

\[\begin{split} \begin{align} \frac{du}{dx} - 1 &= \frac{1}{u} - 1 \\ u du &= dx \end{align} \end{split}\]

両辺を積分し，整理する．

\[\begin{split} \begin{align} \int u du &= \int dx \\ \frac{1}{2} u^2 &= x + C \\ \frac{1}{2} (y+x)^2 &= x + C \end{align} \end{split}\]

これを \(y\) について解く．

\[ y = -x \pm \sqrt{2x + 2C} \]

以上より一般解が得られた．

Pythonによる実装

では，上記の微分方程式をsympyで解く．

from sympy import symbols, Eq, Derivative, Function, dsolve

# 変数と未知関数の定義
x = symbols('x')
y = Function('y')(x)

# dy/dxの定義
dy = Derivative(y, x)

# 微分方程式の定義
eq = Eq(dy, 1 / (y + x) - 1)

# 微分方程式を解く
y_ = dsolve(eq, y)

このとき手計算した一般解からもわかるように二つの一般解を持つことに注意する．

y_[0]

\[\displaystyle y{\left(x \right)} = - x - \sqrt{C_{1} + 2 x}\]

y_[1]

\[\displaystyle y{\left(x \right)} = - x + \sqrt{C_{1} + 2 x}\]

import numpy as np
from sympy import plotting
import matplotlib.pyplot as plt

N = 10

eqs = []
for c in np.linspace(-3, 3, N):
    eqs.append(y_[0].rhs.subs(symbols('C1'), c))
    eqs.append(y_[1].rhs.subs(symbols('C1'), c))
p = plotting.plot(*eqs, (x, -2, 2), show=False)

cm = plt.get_cmap('magma', N*2)
for i in range(0, cm.N, 2):
    p[i].line_color = cm(i)
    p[i+1].line_color = cm(i)
p.show()