https://atatat.hatenablog.com/entry/sympy23

　今回はSymPyを用いて偏微分方程式を解くPDEツールについて紹介したい。基本的には以前紹介した微分方程式を解く方法と似ているが、現状1次のPDEなど制約が多いようだ。

atatat.hatenablog.com

1. 偏微分方程式ソルバーpdsolve

　微分方程式ではdsolveを用いていたが。偏微分方程式ではpdsolveを用いる。手順は偏微分方程式eqを定義して、pdsolve(eq)とすれば解が得られる。

$\displaystyle 2f(x,y)+3\frac{\partial f(x,y)}{\partial x}+2\frac{\partial f(x,y)}{\partial y}=0$

を解く。

まず今回利用するモジュール、関数を以下のようにロードする。

from sympy import Function, symbols, exp
from sympy.solvers.pde import checkpdesol, pdsolve, classify_pde

次にシンボリック変数、関数を定義し、上の偏微分方程式をeqとして定義する。

x,y=symbols('x y')
f=Function('f')
eq=2*f(x,y)+3*f(x,y).diff(x)+4*f(x,y).diff(y)
eq
>2*f(x, y) + 3*Derivative(f(x, y), x) + 4*Derivative(f(x, y), y)

pdsolve(eq)とすれば解が得られる。

sol=pdsolve(eq)
sol
>Eq(f(x, y), F(4*x - 3*y)*exp(-6*x/25 - 8*y/25))

偏微分方程式の解は

$\displaystyle f(x,y)=F(4x-3y)\exp \left( -\frac{6x}{25}-\frac{8y}{25} \right)$

ということがわかる。

2. 解の確認

基本的にpdsolveで得られた解は偏微分方程式を満たすが、偏微分の性質から微分方程式よりもさまざまな解が候補になりえる。そこで解が偏微分方程式を満たしているかどうかchkpdesol(eq,sol)で確認することができる。満たしていればTrue,満たしていなければFalseが返される。

まず、pdsolveで得られた解でチェックする。

checkpdesol(eq,sol)
>(True, 0)

解のうちF(4*x-3*y)は任意の関数を指定できるので、例えばF(x,y)=(4*x - 3*y)が解になっているか確認するには

 sol1= (4*x - 3*y)*exp(-6*x/25 - 8*y/25)
checkpdesol(eq,sol1)
>(True, 0)

同様にF(x,y)=(4*x - 3*y)**2とした場合も

sol2= (4*x - 3*y)**2*exp(-6*x/25 - 8*y/25)
checkpdesol(eq,sol2)
>(True, 0)

さらに線形関数以外を選択する場合、例えばF(x,y)=exp(4*x - 3*y)

sol3= exp(4*x - 3*y)*exp(-6*x/25 - 8*y/25)
checkpdesol(eq,sol3)
>(True, 0)

任意の関数の引数が4*x-3*y以外だと偏微分方程式は満たさない。例えばF(x,y)=exp(4*x-3)とすると

sol4= exp(4*x-3)*exp(-6*x/25 - 8*y/25)
checkpdesol(eq,sol4)
>(False, 12*exp(94*x/25 - 8*y/25 - 3))

3. 任意定数の偏微分方程式

上で解いた偏微分方程式の係数を任意の定数にした場合も同様に解が得られる。

$\displaystyle a\frac{\partial f(x,y)}{\partial x}+b\frac{\partial f(x,y)}{\partial y}+cf(x,y)=0$

a,b,cをシンボリック変数として定義して、上と同様にpdsolve()を適用する。

a,b,c=symbols('a b c')
eq_gen=a*f(x,y).diff(x)+b*f(x,y).diff(y)+c*f(x,y)
sol_gen=pdsolve(eq_gen)
sol_gen
>Eq(f(x, y), F(-a*y + b*x)*exp(-c*(a*x + b*y)/(a**2 + b**2)))

得られる解は

$\displaystyle f(x,y)=F(bx-ay)\exp \left( -\frac{c ( ax+by) }{a^{2}+b^{2}} \right)$

上の場合はa=3,b=4,c=2だったので一般解になっている。

最後にclassify_pdeにより偏微分方程式の分類を表示させることができる。今回解いた方程式は1st_linear_constant_coeff_homogeneousとなっている。

classify_pde(eq_gen)
>('1st_linear_constant_coeff_homogeneous',)

4. まとめ

　SimPyのpdfsolveにより条件付きながら偏微分方程式の解を得る方法について紹介した。