SymPyの使い方9 ～微分方程式の解 - つれづれなる備忘録

1. 基本的な微分方程式の定義、解法
2. 2階微分方程式の解法
3. 一般的な微分方程式の解法
4. まとめ

　今回はSymPyを用いて微分方程式の解を求める方法について紹介する。

1. 基本的な微分方程式の定義、解法

　定義が未知の関数、例えばf(x),g(x)を定義するにはsymbols(関数名,cls=Function)を用いて

f,g=symbols('f g',cls=Function)

とする。導関数は

f(x).diff(x)
>Derivative(f(x), x)

とすれば、定義が未知の導関数を表現することができる。

微分方程式を定義するには、前回の方程式の解を求めたときと同様にEqを利用することができる。

atatat.hatenablog.com

ここで最も基本的な微分方程式

$\displaystyle \frac{d f(x)}{dx} = a f(x)$

を定義するには以下のようにする。

diffeq1=Eq(f(x).diff(x),a*f(x))
deffeq1

微分方程式を解くにはdsolve('方程式',求める関数)を利用する。上の微分方程式を解くには

dsolve(diffeq1,f(x))
>Eq(f(x), C1*exp(a*x))

Eq(f(x),C1*exp(a*x))はf(x)=C1exp(ax)という意味になるので、正しい解が得られていることがわかる。

2. 2階微分方程式の解法

2階微分も上と同様の手順で解くことができる。例えば

$\displaystyle \frac{d^{2} f(x)}{dx^{2}}=k^{2}f(x), k>0$

を解くには、まず、k>0という指定をする。

k=symbols('k',positive=True)

次に微分方程式を以下のように定義して解く。

diffeq2=Eq(f(x).diff(x,x),k**2*f(x))
dsolve(diffeq2,f(x))
>Eq(f(x), C1*exp(-k*x) + C2*exp(k*x))

また、振動解が得られる

$\displaystyle \frac{d^{2} f(x)}{dx^{2}}=-k^{2}f(x), k>0$

は以下のように解く。

diffeq3=Eq(f(x).diff(x,x),-k**2*f(x))
dsolve(diffeq3,f(x))
>Eq(f(x), C1*sin(k*x) + C2*cos(k*x))

3. 一般的な微分方程式の解法

　一般的な微分方程式は、解けるものに限りがあるが、例えば

$\displaystyle \frac{d f(x)}{dx}+p(x)f(x)=0$

p(x)は任意の関数とする。上の微分方程式を以下のように定義し解くことができる。

p=symbols('p',cls=Function)
diff1st_gen=Eq(f(x).diff(x)+p(x)*f(x),0)
>Eq(f(x), C1*exp(-Integral(p(x), x)))

右辺に任意の関数がある非斉次型の微分方程式

$\displaystyle \frac{d f(x)}{dx}+p(x)f(x)=q(x)$

は同様の手順で解くことはできるが、f(x)=という陽の形式にならない。

q=symbols('q',cls=Function)
diff1st_gen=Eq(f(x).diff(x)+p(x)*f(x),q(x))
>Eq((exp(Integral(p(x), x)) - Integral(p(x)*exp(Integral(p(x), x)), x))*f(x) + Integral((f(x)*p(x) - q(x))*exp(Integral(p(x), x)), x), C1)

なお知られている一般解は

$\displaystyle f(x)=\left( \int q(x') \exp \left( \int p(x'') dx'' \right) dx' + C \right) \exp \left(-\int p(x') dx' \right)$

また解は知られているが、dsolveで解けない微分方程式もある。例えば

$\displaystyle \frac{d^{2}f(x)}{dx^{2}} -xp(x)\frac{df(x)}{dx}+p(x)f(x)=0$

を解くために、以下のように実行しようとするとエラーになる。

diff2nd_gen=Eq(f(x).diff(x,x)-x*p(x)*f(x).diff(x)+p(x)*f(x),0)
dsolve(diff2nd_gen,f(x))
>NotImplementedError: solve: Cannot solve -x*p(x)*Derivative(f(x), x) + f(x)*p(x) + Derivative(f(x), x, x)

4. まとめ

　今回はSymPyを用いて微分方程式の解を求める方法について紹介した。