どうも、木村（@kimu3_slime）です。
今回は、Julia（SymPy）で重積分を計算する方法を紹介します。
準備
SymPyを使うので、持っていなければインストールしておきましょう。
1 2 | using Pkg Pkg.add("SymPy") |
準備として、以下のコードを実行しておきます。
1 | using SymPy |
重積分を計算する方法
今回の記事において、記号として利用する変数を用意しておきます。
1 2 | @vars x y z r θ real=true @vars R positive=true |
「integrate(関数,範囲1,範囲2)」で重積分（逐次積分）が計算できます。
1 2 | f1 = 1- x^2 -y^2 integrate(f1,(x, 0, 1),(y,0,1)) |
\[ \begin{aligned}- x^{2} – y^{2} + 1\end{aligned} \]
\[ \begin{aligned}\frac{1}{3}\end{aligned} \]
関数\(f\)の\(x\)方向の重心は
\[ \begin{aligned}\bar{x}:= \frac{\int _D xf(x,y) dx dy}{\int _D f(x,y) dx dy}\end{aligned} \]
と定義されます。これを関数として定義し、計算してみましょう。
1 2 3 4 5 | function center(f) return integrate(x*f,(x, 0, 1),(y,0,1))/integrate(f,(x, 0, 1),(y,0,1)) end center(f1) |
\[ \begin{aligned}\frac{1}{4}\end{aligned} \]
積分する領域は正方形でなくても、範囲を数式として指定できれば計算できます。三角形領域\(D=\{(x,y) \mid 0 \leq y \leq x, 0 \leq x \leq 1\}\)における重積分を試してみましょう。
1 | integrate(f1,(y,0,x),(x, 0, 1)) |
\[ \begin{aligned}\frac{1}{6}\end{aligned} \]
変数\(x\)について先に積分しても、結果は同じです（フビニの定理）。積分範囲の指定に気をつけましょう。
1 | integrate(f1,(x,y,1),(y, 0, 1)) |
\[ \begin{aligned}\frac{1}{6}\end{aligned} \]
半径\(R\)の円盤の面積を
\[ \begin{aligned}\int _{D_R} 1 dx dy \end{aligned} \]
を計算することによって求めてみましょう。
1 2 | integrate(Sym(1), (y, -sqrt(R^2-x^2), sqrt(R^2-x^2)) ) integrate(Sym(1), (y, -sqrt(R^2-x^2), sqrt(R^2-x^2)),(x, -R, R) ) |
\[ \begin{aligned}2 \sqrt{R^{2} – x^{2}}\end{aligned} \]
\[ \begin{aligned}\pi R^{2}\end{aligned} \]
きちんと円の面積が求められています。
ガウス積分の2乗を
\[ \begin{aligned}\int_{-\infty} ^\infty \int_{-\infty} ^\infty e^{-(x^2+y^2)}dxdy\end{aligned} \]
として求めてみましょう。積分範囲が無限（広義積分）であっても、2つのオー「oo」で無限大を指定できます。
1 | integrate(exp(-x^2 - y^2), (x, -oo, oo), (y, -oo, oo)) |
\[ \begin{aligned}\pi\end{aligned} \]
円盤領域で直交座標によって計算すると、適切な結果を表示してくれません。
1 | integrate(exp(-x^2 - y^2), (y, -sqrt(R^2-x^2), sqrt(R^2-x^2)) ,(x, -R, R) ) |
\[ \begin{aligned}\sqrt{\pi} \int_{- R}^{R} e^{- x^{2}} \operatorname{erf}{\left(\sqrt{R^{2} – x^{2}} \right)}\, dx\end{aligned} \]
\(\operatorname{erf}\)は誤差関数（ガウス積分の一部分）であり、積分の結果が単純になっていません。
これは極座標変換して計算することで解決します。
\[ \begin{aligned}\int_D f(x,y)dxdy = \int_I f(r\cos\theta, r \sin \theta) r dr d\theta\end{aligned} \]
被積分関数の変数を「数式.subs(代入する変数と値)」によって置き換えましょう。
1 2 | f2 =exp(-x^2 - y^2) f3 = simplify(f2.subs([(x, r* cos(θ)),(y,r* sin(θ))])) |
\[ \begin{aligned}e^{- x^{2} – y^{2}}\end{aligned} \]
\[ \begin{aligned}e^{- r^{2}}\end{aligned} \]
そして積分すれば
1 | integrate(f3 *r , (r, 0, R) ,(θ, 0, 2*PI) ) |
\[ \begin{aligned}2 \pi \left(\frac{1}{2} – \frac{e^{- R^{2}}}{2}\right)\end{aligned} \]
で、\(R\to \infty\)の極限を取れば
1 | limit(integrate(f3 *r , (r, 0, R) ,(θ, 0, 2*PI) ) , R, oo) |
\[ \begin{aligned}\pi \end{aligned} \]
さきほどと同じ結果が得られました。
ちなみに、ガウス積分の値自体は、1変数の広義積分としてそのまま計算できます。
参考：Julia（SymPy）で1変数関数を積分する方法（多項式、指数対数、三角関数）
以上、Julia（SymPy）で重積分を計算する方法を紹介してきました。
1変数関数だけでなく多変数の積分も、積分範囲に気をつければ、簡単に計算できますね。
木村すらいむ（@kimu3_slime）でした。ではでは。
コロナ社 (2020-03-26T00:00:01Z)
￥7,353 (コレクター商品)
こちらもおすすめ
Julia（SymPy）で1変数関数を積分する方法（多項式、指数対数、三角関数）
