python - 使用 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 比较陌生。

扩展：我想为 arbiraty 实现数值积分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。

标签： pythonscipynumerical-integration

Answer 1

这是一个使用scipy.integrate.dblquad代替的更简单的方法nquad：

从 x = a..b 和 y = gfun(x)..hfun(x) 返回 func(y, 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）：

整合有序进行。也就是说，x0 上的积分是最内层的积分，xn 是最外层的积分

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,) )