python嵌套數值積分

Question

我目前正在計算python中的數值積分，它實際上是四個嵌套數值積分。每個積分是不同變量的函數，但我遇到的問題是限制取決於其嵌套內積分（最外面的積分具有適當的浮點範圍）。例如，假設我有一個具有限制amin，amax，bmin，bmax，cmin，cmax，dmin，dmax的函數f（a，b，c，d）。 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異常遇到了同樣的錯誤：給予無效的界限，這似乎發生，因爲我有沒有被符號整合在一起。我也試過使用nquad，但是我遇到了同樣的錯誤。

這是我在處理當前試圖反覆向外整合，只是做了NUMERICS底：

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

嵌套調用quad是正確的方法。但是，簡單是不夠的。傳遞給quad的每個參數都需要是一個函數 - 不是先前調用的結果。下面是一個將f = x1**2 - x2**3整合到菱形區域abs(x1) + abs(x2) <= 1上的示例。

我們需要能夠在-/+(1 - abs(x1))之間整合之間的任何x1+/-1之間。

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) 

Answer 2

如果我理解正確，那麼nquad應該能夠處理這種整合。請注意，其參數ranges允許您指定取決於其他集成變量值的函數。

例如，計算F（X，Y）= 1上的積分區域0 < = Y < = 1和0 < = X < = Y（換言之，我們計算上的區域單位正方形的-left半三角形）：

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，因爲它應該。