使用 scipy.integrate.nquad 实现数值积分 [英] Implementing numerical integration using scipy.integrate.nquad

问题描述

我有这个具有相关限制的二维积分.该函数可以在 Python 中定义为

I have this 2-dimensional integral with dependent limits. The function can be defined in Python as

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.

where the limits of u3 is from 0 to gamma (a positive real number), and the limits of u2 is from 0 to gamma-u3.

如何使用 scipy.integrate.nquad?我试图阅读文档，但并不容易理解，尤其是我对 Python 比较陌生.

How can I implement this using scipy.integrate.nquad? I tried to read the documentation, but it was not easy to follow, especially I am relatively new to Python.

扩展:我想对任意 K 实现数值积分，在这种情况下，被积函数由 (1-1/(1+gamma-uk-....-u2))*(1/(1+uK)**2)*...*(1/(1+u2)**2).我编写了一个接受动态数量参数的函数，如下所示:

Extension: I would like to implement a numerical integration for an arbiraty K, where the integrand in this case is given by (1-1/(1+gamma-uk-....-u2))*(1/(1+uK)**2)*...*(1/(1+u2)**2). I wrote the function that takes a dynamic number of arguments as follows:

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.

However, I am not sure how to do the same for the ranges, where I will have one function for the ranges, where uK ranges from 0 to gamma, u_{K-1} ranges from 0 to gamma-uK, ...., u2 ranges from 0 to gamma-uK-...-u2.

推荐答案

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

Here is a simpler method using scipy.integrate.dblquad instead of nquad:

从 x = a..b 返回 func(y, x) 的双(定)积分和y = gfun(x)..hfun(x).

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 必须是一个全局变量.

It seems that gfun and hfun do not accept the extra arguments, so gamma has to be a global variable.

使用 nquad，经过多次尝试和错误:

Using nquad, after many trial and error:

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源代码的有用引用:

Useful quote from the source code of tplquad:

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

并且从 nquad 文档来看，变量的顺序是相反的(对于 dblquad 也是如此):

And from the nquad documentation, the order of the variables is reversed (same for dblquad):

集成是按顺序进行的.即x0上的积分是最里面的积分，xn是最外面的

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

具有 k 嵌套集成的通用案例:

Generic case with k nested integrations:

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 实现数值积分的文章就介绍到这了，希望我们推荐的答案对大家有所帮助，也希望大家多多支持IT屋！