使用 scipy.integrate.nquad 实现数值积分

Question

我有这个带有相关极限的二维积分。该函数可以在Python中定义为

def func(gamma, u2, u3):
    return (1-1/(1+gamma-u3-u2))*(1/(1+u2)**2)*(1/(1+u3)**2)

其中u3的范围是0到gamma（正实数），u2的范围是0到gamma-u3。

如何使用 scipy.integrate.nquad 实现此功能？我尝试阅读文档，但并不容易理解，尤其是我对 Python.

比较陌生

扩展：我想为任意值 K 实现数值积分，在这种情况下，被积函数由 (1-1/(1+gamma-uk-....-u2))*(1/(1+uK)**2)*...*(1/(1+u2)**2) 给出。我编写了采用动态参数数量的函数，如下所示：

def integrand(gamma, *args):
    '''
    inputs:
     - gamma
     - *args = (uK, ..., u2)

    Output:
     - (1-1/(1+gamma-uk-....-u2))*(1/(1+uK)**2)*...*(1/(1+u2)**2)
    '''
    L = len(args)
    for ll in range(0, L):
        gamma -= args[ll]
    func = 1-1/(1+gamma)
    for ll in range(0, L):
        func *= 1/((1+args[ll])**2)
    return func

但是，我不确定如何对范围执行相同的操作，我将为范围设置一个函数，其中 uK 的范围从 0 到 gamma，u_{K-1}的范围从0到gamma-uK，....，u2的范围从0到gamma-uK-...-u2。

Answer 1

这是一个更简单的方法，使用 scipy.integrate.dblquad 而不是 nquad：

Return the double (definite) integral of func(y, x) from x = a..b and y = gfun(x)..hfun(x).

from  scipy.integrate import dblquad

def func(u2, u3, gamma):
    return (1-1/(1+gamma-u3-u2))*(1/(1+u2)**2)*(1/(1+u3)**2)


gamma = 10

def gfun(u3):
    return 0

def hfun(u3):
    return gamma-u3

dblquad(func, 0, gamma, gfun, hfun, args=(gamma,))

似乎gfun和hfun不接受额外的参数，所以gamma必须是一个全局变量。

使用nquad，经过多次试验和错误：

from  scipy.integrate import nquad

def func(u2, u3, gamma):
    return (1-1/(1+gamma-u3-u2))*(1/(1+u2)**2)*(1/(1+u3)**2)

def range_u3(gamma):
    return (0, gamma)

def range_u2(u3, gamma):
    return (0, gamma-u3)

gamma = 10
nquad(func, [range_u2, range_u3], args=(gamma,) )

tplquad源代码的有用引用：

# nquad will hand (y, x, t0, ...) to ranges0
# nquad will hand (x, t0, ...) to ranges1

并且从 nquad 文档中，变量的顺序是相反的（dblquad 相同）：

Integration is carried out in order. That is, integration over x0 is the innermost integral, and xn is the outermost

具有 k 嵌套集成的一般案例：

from  scipy.integrate import nquad
import numpy as np

def func(*args):
    gamma = args[-1]
    var = np.array(args[:-1])

    return (1-1/(1+gamma-np.sum(var)))*np.prod(((1+var)**-2))

def range_func(*args):
    gamma = args[-1]
    return (0, gamma-sum(args[:-1]))

gamma, k = 10, 2
nquad(func, [range_func]*k, args=(gamma,) )

使用 scipy.integrate.nquad 实现数值积分

Implementing numerical integration using scipy.integrate.nquad

python

scipy

numerical-integration