nano_exit

基礎的なことこそ、簡単な例が必要だと思うのです。

sympyで(平面の)Greenの定理を確認

sympyの練習を兼ねて、平面に対するGreenの定理で、問題を解いてみる。
平面のグリーンの定理 [物理のかぎしっぽ]

以下の積分を求めてみる。

\displaystyle
\int \! \int_D ( 1 - ( x^2 + y^2) ) dx dy
\\
\displaystyle
D = \left\{ (x,y) | x^2 + y^2 \le 1 \right\}
原点で最大値を取り、等方向的で、境界でゼロを持つような、何かしらの密度を積分する、と思えば解りやすいだろうか。

普通に極座標に変換して解くと、面積素片

\displaystyle
dx \, dy
=
\left|
\begin{array}{cc}
  \frac{ \partial x }{ \partial r } & \frac{ \partial x }{ \partial \theta } \\ 
  \frac{ \partial y }{ \partial r } & \frac{ \partial y }{ \partial \theta } 
\end{array}
\right|
dr \, d\theta
=
\left|
\begin{array}{cc}
  \cos \theta & - r \sin \theta \\ 
  \sin \theta &  r \cos \theta
\end{array}
\right|
dr \, d\theta
=
r\, dr \, d\theta
であるから(ヤコビアンの中身は転置してても大丈夫)、

\displaystyle
\int \! \int_D ( 1 - ( x^2 + y^2 ) ) dx dy
  = \int^1_0 r dr \int^{2\pi}_0 \, d\theta \, ( 1 - r^2 )
\\
\displaystyle
  = 2 \pi \left[ \frac{ r^2 }{ 2 } - \frac{ r^4 }{ 4 } \right]^1_0
  = \frac{ \pi }{ 2 }

この積分に対する一連の操作をsympyでやると、

from sympy import *
init_session()

r, theta, c, s = symbols( "r θ c s" )
c = cos( theta )
s = sin( theta )

xrt = r * c
yrt = r * s

f0 = 1 - ( x**2 + y**2 )
f1 = f0.subs( [ ( x, xrt ), ( y, yrt) ] )

Jacob = Matrix([
  [ Derivative( xrt, r ), Derivative( xrt, theta ) ],
  [ Derivative( yrt, r ), Derivative( yrt, theta ) ]
])
Jacob = simplify( Jacob.det().doit() )

I = Integral( f1 * Jacob, ( r, 0, 1 ), ( theta, 0, 2 * pi ) )
I.doit()
# pi / 2

積分区間の決定は自分で判断する必要があるが、大方、機械的に処理できるのではないだろうか。
途中式を吐き出すようにすれば、もう少し分かったような気がするようになると思う。

注意点として、*.subs()で関係式を代入して式変形する際、Derivativeにそれで代入しようとしても、「積分した後の変数に代入しようとする」ため、例えば以下のようにしようとするとゼロ行列になる。

Jacob = Matrix([
  [ Derivative( x, r ), Derivative( x, theta ) ],
  [ Derivative( y, r ), Derivative( y, theta ) ]
])
Jacob = Jacob.subs([ ( x, xrt ), ( y, yrt ) ]).doit()
# zero matrix

Integral()も同様に、*.subs()を使っても積分範囲に代入されて、被積分関数を弄ることが出来ない。

この積分を、Greenの定理で線積分に直して計算してみる。
Greenの定理は以下の式で与えられる。

\displaystyle
\int \! \int_D \left( \frac{ \partial Q }{ \partial x } - \frac{ \partial P }{ \partial y } \right) dx dy
  =  \oint_C \left( P \, dx + Q \, dy \right)

したがって、例えば以下のように P, Qを定めると、

\displaystyle
\frac{ \partial Q }{ \partial x } = \frac{ 1 }{ 2 } - x^2
\rightarrow Q = \frac{ x }{ 2 } - \frac{ x^3 }{ 3 }
\\
\displaystyle
  - \frac{ \partial P }{ \partial y } = \frac{ 1 }{ 2 } - y^2
\rightarrow P = - \left( \frac{ y }{ 2 } - \frac{ y^3 }{ 3 } \right)
\\
\displaystyle
\therefore
\oint_C \left( P \, dx + Q \, dy \right) = 
  - \oint_C \left( \frac{ y }{ 2 } - \frac{ y^3 }{ 3 } \right) dx
  + \oint_C \left( \frac{ x }{ 2 } - \frac{ x^3 }{ 3 } \right) dy
\\
\displaystyle
= \int^{2\pi}_0 \left( \frac{ \sin \theta }{ 2 } - \frac{ \sin^3 \theta }{ 3 } \right) \sin \theta d\theta
  + \oint_C \left( \frac{ \cos \theta }{ 2 } - \frac{ \cos^3 \theta }{ 3 } \right) \cos \theta d\theta

これを計算するとちゃんと \pi/2になることを手で示しても良いが、計算ミスだらけで効率が悪いので、sympyにやってもらうことにする。

P = - integrate( 1/2 - y**2, y )
integrand = P.subs([ ( x, xrt ), ( y, yrt ) ]) * diff( xrt, theta )
integrand = integrand.subs([ ( r, 1 ) ])
int_range = ( theta, 0, 2 * pi )
P_int = integrate( integrand, int_range )
# pi / 4

Q = integrate( 1/2 - x**2, x )
integrand = Q.subs([ ( x, xrt ), ( y, yrt ) ]) * diff( yrt, theta )
integrand = integrand.subs([ ( r, 1 ) ])
int_range = ( theta, 0, 2 * pi )
Q_int = integrate( integrand, int_range )
# pi / 4

I = P_int + Q_int
# pi / 2

 O_x, P_yによる被積分関数の分割を、上の例と逆にしても、同じ答えが得られる。

P = - integrate( 1/2 - x**2, y )
integrand = P.subs([ ( x, xrt ), ( y, yrt ) ]) * diff( xrt, theta )
integrand = integrand.subs([ ( r, 1 ) ])
int_range = ( theta, 0, 2 * pi )
P_int = integrate( integrand, int_range )
# pi / 4

Q = integrate( 1/2 - y**2, x )
integrand = Q.subs([ ( x, xrt ), ( y, yrt ) ]) * diff( yrt, theta )
integrand = integrand.subs([ ( r, 1 ) ])
int_range = ( theta, 0, 2 * pi )
Q_int = integrate( integrand, int_range )
# pi / 4

I = P_int + Q_int
# pi / 2

もっと言えば、 Q_x被積分関数を全部押し付けてもちゃんと求まる。

P_int = 0

Q = integrate( 1 - ( x**2 + y**2 ), x )
integrand = Q.subs([ ( x, xrt ), ( y, yrt ) ]) * diff( yrt, theta )
integrand = integrand.subs([ ( r, 1 ) ])
int_range = ( theta, 0, 2 * pi )
Q_int = integrate( integrand, int_range )
# pi / 2

I = P_int + Q_int
# pi / 2

今はsympyで数値計算的なチェックになってしまっているが、代入を駆使して途中式のチェックとかが出来るとかなり便利だとは思う。
(ただ、手で書くのとどっちが速いか、という問題はあるが。。。)