在python中嵌套数值积分答案

【问题标题】：Nesting numerical integrals in python在python中嵌套数值积分
【发布时间】：2017-05-11 21:23:52
【问题描述】：

我目前正在计算 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）的值传递到限制中，使它们变为浮点数。

有没有办法用 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

【问题讨论】：

@Gerry 非常感谢 cmets，我已经尝试解决这个问题。如果需要更多信息，请告诉我！

标签： python loops scipy nested numerical-integration

【解决方案1】：

嵌套调用quad 是正确的方法。然而，简单是不够的。传递给 quad 的每个参数都必须是一个函数 - 而不是先前调用的结果。这是一个在菱形区域abs(x1) + abs(x2) <= 1 上集成f = x1**2 - x2**3 的示例。

我们需要能够在-/+(1 - abs(x1)) 之间为x1 之间的任何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)

【讨论】：

谢谢，多玩了一点，它成功了！我已经发布了上述问题中的代码以供将来参考。
@PhysicistAbroad 很高兴这对您有所帮助。嵌套函数实际上很好地反映了嵌套积分的概念。

【解决方案2】：

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

例如，要计算 f(x,y)=1 在面积 0

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，这是应该的。

【讨论】：