もう結構古い記事なのですが、以前SciPyで定積分をやる方法を記事にしたことがあります。
参考: scipyで定積分
最近、2変数関数の積分をやる機会があったのでこの機会に重積分をSciPyで行う方法を紹介します。SciPyのintegrateモジュールには、重積分用の関数が複数あります。
– dblquad (2変数関数の定積分)
– tplquad (3変数関数の定積分)
– nquad (一般のn変数の定積分)
dblquadの使い方
順番に説明していきます。まずは2重積分のdblquadです。関数の定義は次のようになっています。
scipy.integrate.dblquad(func, a, b, gfun, hfun, args=(), epsabs=1.49e-08, epsrel=1.49e-08)
必須なのは、積分対象のfunc, 外側の積分区間のa, b、そして内側の積分区間を示す、gfun, hfunです。
gfunとhfunは名前からわかる通り、定数ではなく関数です。これにより内側の積分の積分区間を変数にすることができます。つまり以下のような積分区間の積分ができます。
$$\int_0^1\int_0^y xy \,dxdy$$
例に挙げたのでこれを実装してみましょう。ちなみに解は$1/8=0.125$です。funcの定義は、内側の関数を第1引数にする必要があるので注意してください。
from scipy.integrate import dblquad
def f1(x, y):
return x*y
def x0(y):
return 0
def x1(y):
return y
print(dblquad(f1, 0, 1, x0, x1))
# (0.125, 5.515032205777789e-15)想定通りですね。積分結果と推定誤差が返ってくるのは1変数の積分と同様です。
内側の積分区間も定数から定数までだよ、要するに長方形領域で積分したいよ、って場合はgfun, hfun に定数を返す関数を返してください。
tplquadの使い方
続いて、3変数向けのtplquadです。これもdblquadとかなり似てる感じで使えます。積分変数が一個増えているので上限加減の指定がもう一個ある感じです。
scipy.integrate.tplquad(func, a, b, gfun, hfun, qfun, rfun, args=(), epsabs=1.49e-08, epsrel=1.49e-08)
たとえば次の積分をやってみましょう。
$$\int_0^1\int_0^z\int_0^{y+z} xyz \,dxdydz.$$
ちなみに答えは$17/144=0.11805555…$となるはずです。
引数の順番に注意が必要なので慎重にコーディングしてください。
from scipy.integrate import tplquad
def f2(x, y, z):
return x*y*z
def x0(y, z):
return 0
def x1(y, z):
return y + z
def y0(z):
return 0
def y1(z):
return z
print(tplquad(f2, 0, 1, y0, y1, x0, x1))
# (0.11805555555555557, 2.1916761217856673e-14)バッチリですね。
nquadの使い方
最後に一般のn変数を積分できるnquadの使い方を紹介します。
引数の形式が先ほどの二つと少し違います。
scipy.integrate.nquad(func, ranges, args=None, opts=None, full_output=False)
funcにn変数関数を渡して、rangesに積分区間を渡すことになります。rangesは配列で、1変数目から順番に区間の下限上限の2値の配列を格納しておけば良いです。また、ここにも一応関数を使うことはできます。
これはシンプルな例で、定数関数1を超立方体区間で積分してみました。
from scipy.integrate import nquad
def f3(w, x, y, z):
return 1
print(nquad(f3, [[-2, 2], [-2, 2], [-2, 2], [-2, 2]]))
# (256.0, 2.8421709430404007e-12)
$4^4=256$になりましたね。
ここで急にシンプルな例を出したのには事情がありまして、変数の数が多くなるとやはり積分は困難なようで、ちょっと複雑な例になると規定の反復回数をこなしても必要な精度に届かずWarningが出たりするケースが多々あります。
どうしても計算したい場合は limit パラメーター等をいじっての対応になりますのでドキュメントを参照しながら調整してみてください。
(僕も実運用で必要になったら改めて調査して紹介しようと思います。)