SymPyの使い方25 ～ PDE2 - つれづれなる備忘録

　今回も前回に続いてSymPyを用いて偏微分方程式を解くPDEツールについて紹介したい。

atatat.hatenablog.com

1. 偏微分方程式ソルバーpdsolveの制約
2. 右辺含む任意の1次偏微分方程式
3. 係数が関数の場合
4. まとめ

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

　前回触れたように現状はpdsolveは1次の偏微分のみ対応しており、2次以上は解くことができない。試しに以下のように2次の偏微分方程式を定義して解こうとするとNotImplementedError: psolve: Cannot solve Derivativeと表示されエラーになる。

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

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

sol=pdsolve(eq)
>NotImplementedError: psolve: Cannot solve Derivative(f(x, y), (x, 2)) + Derivative(f(x, y), (y, 2)

2. 右辺含む任意の1次偏微分方程式

前回のpde_1st_linear_constant_coeff_homogeneous以外で右辺に任意関数がある場合も解くことができる。

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

g=Function('g')
a,b,c=symbols('a b c')
gen1_eq=a*f(x,y).diff(x)+b*f(x,y).diff(y)+c*f(x,y)-g(x,y)

pdsolveで解くと、少し時間がかかって以下のように解が表示される。

sol_gen1=pdsolve(gen1_eq)
sol_gen1
>Eq(f(x, y), ((a**2 + b**2)*F(-a*y + b*x) + Integral(g(-a*b*y/(a**2 + b**2) + a*xi/(a**2 + b**2) + b**2*x/(a**2 + b**2), a**2*y/(a**2 + b**2) - a*b*x/(a**2 + b**2) + b*xi/(a**2 + b**2))*exp(c*xi/(a**2 + b**2)), (xi, a*x + b*y)))*exp(-c*(a*x + b*y)/(a**2 + b**2))/(a**2 + b**2))

上を少し整えると

$\displaystyle f(x,y)=\left( F(\eta)+\frac{\int^{\eta} g (\alpha,\beta) \exp ( \frac{c\xi}{a^{2}+b^{2}} ) d\xi }{a^{2}+b^{2}} \right) \exp \left( \frac{c\xi}{a^{2}+b^{2}} \right)$

$\displaystyle \xi=ax+by, \eta=-ay+bx$

$\displaystyle \alpha=\frac{a\xi+b\eta}{a^{2}+b^{2}},\beta=\frac{-a\eta+b\xi}{a^{2}+b^{2}}$

classify_pdeで方程式の型を表示すると1st_linear_constant_coeff_Integralが得られる。

classify_pde(gen1_eq)
>('1st_linear_constant_coeff', '1st_linear_constant_coeff_Integral')

dsolveのときと同様にhintオプションで1st_linear_constant_coeff_Integralを指定すると処理が早くなる。

sol_gen1=pdsolve(gen1_eq,hint='1st_linear_constant_coeff_Integral')

係数a,b,cやg(x,y)が任意でない場合は、hintオプションなしでも即座に解が計算される。

gen1_eq2= f(x,y).diff(x) + f(x,y).diff(y) + f(x,y) - exp(x + y)
sol2_gen1=pdsolve(gen1_eq2)
sol2_gen1
>Eq(f(x, y), F(x - y)*exp(-x/2 - y/2) + exp(x + y)/3

3. 係数が関数の場合

　係数a,b,cが以下のように関数の場合もpdsolveで解ける

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

ただし、上の偏微分方程式をそのまま解こうとすると非常に時間がかかるので、具体的にそれぞれ関数を指定したものを解いてみる。

gen2_eq=x*(f(x,y).diff(x)) - y*(f(x,y).diff(y)) + y**2*f(x,y) - y**2
classify_pde(gen2_eq)
>('1st_linear_variable_coeff',)

方程式の型は1st_linear_variable_coeffとなる。以下hintは指定しなくてもあまり時間は変わらないが、pdsolveで解が得られる。

 sol_gen2=pdsolve(gen2_eq,hint='1st_linear_variable_coeff')
 sol_gen2
>Eq(f(x, y), F(x*y)*exp(y**2/2) + 1)

4. まとめ

　今回はSimPyのpdfsolveを用いていろいろな1次偏微分方程式の解を得る方法について紹介した。