sympyの練習を兼ねて、平面に対するGreenの定理で、問題を解いてみる。
平面のグリーンの定理 [物理のかぎしっぽ]
以下の積分を求めてみる。
原点で最大値を取り、等方向的で、境界でゼロを持つような、何かしらの密度を積分する、と思えば解りやすいだろうか。
普通に極座標に変換して解くと、面積素片は
であるから(ヤコビアンの中身は転置してても大丈夫)、
この積分に対する一連の操作を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の定理は以下の式で与えられる。
したがって、例えば以下のようにを定めると、
これを計算するとちゃんとになることを手で示しても良いが、計算ミスだらけで効率が悪いので、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
による被積分関数の分割を、上の例と逆にしても、同じ答えが得られる。
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
もっと言えば、に被積分関数を全部押し付けてもちゃんと求まる。
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で数値計算的なチェックになってしまっているが、代入を駆使して途中式のチェックとかが出来るとかなり便利だとは思う。
(ただ、手で書くのとどっちが速いか、という問題はあるが。。。)