在 python 中嵌套数值积分

Question

我目前正在计算 python 中的数值积分，它实际上是四个嵌套的数值积分。每个积分都是不同变量的函数，但我遇到的问题是限制取决于它嵌套在其中的积分（最外面的积分具有适当的浮点范围）。

例如，假设我有一个函数 f(a,b,c,d)，其限制为 amin、amax、bmin、bmax、cmin、cmax、dmin、dmax。 amin 和 amax 是浮动的。 bmin和max是a的函数，cmin和cmax是b的函数，dmin和dmax是c的函数。如果 scipy.integrate.quad() 只是一个 for 循环，那么对于每一步，都可以将 a（或 b、c 或 d）的值传递到 limits 中，使它们成为浮点数。

有没有办法用 scipy.integrate.quad 做到这一点？到目前为止，我已经尝试过简单的嵌套：

def int4(func,min1,max1,min2,max2,min3,max3,min4,max4):
    finalfunc = scint.quad(scint.quad(scint.quad(scint.quad(func,min1,max1),min2,max2),min3,max3),min4,max4)
    return finalfunc

在这种情况下，我仍然遇到与 ValueError: Invalid limits given 之前遇到的相同错误，这似乎是因为我有未被整合的符号。我也尝试过使用nquad，但我遇到了同样的错误。

这是我目前尝试向外迭代迭代的过程，只在最后做数字：

def int4(func,var1,var2,var3,min1,max1,min2,max2,min3,max3,min4,max4):
    func2 = sym.integrate(func,(var1,min1,max1)
    func3 = sym.integrate(func2,(var2,min2,max2))
    func4 = sym.integrate(func3,(var3,min3,max3))
    finalfunc = scint.quad(func4,min4,max4) 
    return finalfunc

困难在于 min1、max1 是 var2 的函数，而我实际使用的函数似乎没有解析解。

如果我改为从最外层进入而不是从最内层向外集成，是否也有帮助？

感谢 Kazemakase，他的回答对解决我的问题很有帮助！我使用对他编写的代码稍作修改的代码解决了这个问题，因此我将其包含在此处作为将来遇到类似问题的其他人的参考。

import numpy as np
import scipy.integrate as si

def func(x1, x2, x3, x4):
    return x1**2 - x2**3+x3*x2 - x4*x3**3  

def int1():
    """integrates `int2` over x1"""
    a1, b1 = -1, 3
    def int2(x1):
        """integrates `func` over x2 at given x1.""" 
        #partial_func1 = lambda x2: func(x1, x2)
        b2 = 1 - np.abs(x1)
        a2 = -np.abs(x1**3)
        def int3(x2):
            a3 = x2
            b3 = -a3
            def int4(x3):
                partial_func = lambda x4: func(x1, x2, x3, x4)
                a4 = 1+np.abs(x3)
                b4 = - a4
                return si.quad(partial_func,a4,b4)[0]
            return si.quad(int4, a3, b3)[0]
        return si.quad(int3, a2, b2)[0]     
    return si.quad(int2, a1, b1)[0]

result = int1()  # -22576720.048151683

Answer 1

如果我理解正确，那么 nquad 应该能够处理此集成。请注意，其 ranges 参数允许您指定依赖于其他积分变量值的函数。

例如，要计算 f(x,y)=1 在面积 0 <= y <= 1 和 0 <= x <= y 上的积分（换句话说，我们计算面积单位正方形的左上半三角形):

from scipy.integrate import nquad
def f(x,y): return 1.
def x_integration_boundary(y): return (0,y)
nquad(f,(x_integration_boundary,(0,1)))

这给出了 0.5 的积分值，这是应该的。

Answer 2

嵌套调用 quad 是正确的方法。然而，简单是不够的。传递给 quad 的每个参数都需要是一个函数——而不是先前调用的结果。这是一个在菱形区域 abs(x1) + abs(x2) <= 1.

上积分 f = x1**2 - x2**3 的示例

对于 +/-1.

之间的任何 x1，我们需要能够在 -/+(1 - abs(x1)) 之间集成 x2

import numpy as np
import scipy.integrate as si


def func(x1, x2):
    return x1**2 - x2**3    

def int2(x1):
    """integrates `func` over x2 at given x1.""" 
    partial_func = lambda x2: func(x1, x2)
    b2 = 1 - np.abs(x1)
    a2 = -b2
    return si.quad(partial_func, a2, b2)[0]    

def int1():
    """integrates `int2` over x1"""
    a1, b1 = -1, 1
    return si.quad(int2, a1, b1)[0]

result = int1()  # 0.33333333333333337

您可以使用 nquad 为您进行包装，而不是为每个积分显式编写一个函数。它为每个变量取一个积分范围，它也可以是一个函数，取决于积分的其他变量：

def range1(x2):
    b1 = 1 - np.abs(x2)
    a1 = -b1
    return a1, b1    

result, err = si.nquad(func, [range1, (-1, 1)])
# (0.33333333333333337, 7.401486830834376e-15)

在 python 中嵌套数值积分

Nesting numerical integrals in python

python

loops

scipy

nested

numerical-integration