Python-可変長引数を持つ関数を積分する

POINT

• 動的に変数の数が変わる関数を「定義」する方法．
• 動的に変数の数が変わる関数を「積分」する方法．

思いついた方法をメモしておきます．
他にもっと良い方法があるかもしれません．ぜひ教えて下さい！

やりたいこと

例えば

\begin{aligned}
f_1(x_{0}, x_{1}, \cdots , x_{N})
=\prod_{i=0}^{N} g_1(x_i)
\end{aligned}

という関数や

\begin{aligned}
f_2(x_{0}, x_{1}, \cdots , x_{N})
=\prod_{i=0}^{N-1} g_2(x_{i}, x_{i+1})
\end{aligned}

という関数を積分したい．

さらに，変数の数（関数$g$を掛ける回数）を色々変えて計算したい．このとき，変数の個数ごとに関数を定義するのではなく，1つの関数でやりくりしたい．

例1：1変数関数を掛け合わせる

上の$f_1$について考えましょう．

可変長引数をもつ関数

具体例として$g_1(x)=x$の場合を考えます．したがって，$f_1(x_{0}, x_{1}, \cdots , x_{N}) = x_{0} \times x_{1} \times \cdots \times x_{N}$となります．

関数は以下のように定義できます：

def g1(x):
    y = x
    return y

def f1(*args):
    y = 1
    for x in args:
        y = y * g1(x)
    return y

もちろん，$g_1$の定義を変えればこのタイプのどんな関数でも表現できます．

上の例では，この関数を利用することで

\begin{aligned}
f_1(1,2,3,\cdots ,N)=N!
\end{aligned}

が計算できます．確認してみましょう：

import numpy as np

def g1(x):
    y = x
    return y

def f1(*args):
    y = 1
    for x in args:
        y = y * g1(x)
    return y

n = 5  # integer
variable = np.arange(1, n+1)
# 5! = 120
print(str(n)+'! =', f1(*variable))

積分

$f_1$の積分は

\begin{aligned}
I
&=\int_{a_{0}}^{b_{0}}\cdots \int_{a_{N}}^{b_{N}}
f_1(x_{0}, x_{1}, \cdots , x_{N})
\,\mathrm{d}x_{0}\cdots \,\mathrm{d}x_{N} \\
&=\prod_{i=0}^{N} \biggl[ \int_{a_{i}}^{b_{i}} g_1(x_i) \,\mathrm{d}x_{i} \biggr]
\end{aligned}

となります．

ここで，また$g_1(x)=x$とすると

\begin{aligned}
I&=\prod_{i=0}^{N} \biggl[ \frac{(b_{i})^2 - (a_{i})^2}{2} \biggr]
\end{aligned}

と厳密に計算できます．そこで，次の具体例で「scipy.integrate.nquadによる積分結果」と「厳密解」を比較して一致することを確認しましょう：

\begin{aligned}
I
&=\int_{0}^{1}\,\mathrm{d}x_{0}
\int_{1}^{2} \,\mathrm{d}x_{1}
\int_{3}^{4}\,\mathrm{d}x_{2}\,
[x_{0} x_{1} x_{2}]
\end{aligned}

import numpy as np
from scipy import integrate

def g1(x):
    y = x
    return y

def f1(*args):
    y = 1
    for x in args:
        y = y * g1(x)
    return y

def exact(*args):
    y = 1
    for x in args:
        y = y * (x[1]**2 - x[0]**2)/2
    return y
    
interval_1 = [0, 1]
interval_2 = [1, 2]
interval_3 = [3, 4]
interval = [interval_1, interval_2, interval_3]

integral, err = integrate.nquad(f1, interval)

# Error= 0.0
print('Error=', integral - exact(*interval))

このように，「scipy.integrate.nquad」の被積分関数として「可変長引数を持つ関数」を指定すれば，渡す積分区間の個数を変えるだけで変数の数（$g_1$を掛ける回数）も自動で変えることができます．