title: Nesting numerical integrals in python
body: <p>I'm currently working on computing a numerical integral in python which is really four nested numerical integrals. Each integral is a function of a different variables, but the problem I'm encountering is that the limits are dependent on the integral its nested within (with the outermost integral having proper float ranges).</p>
  
  <p>For example, say I have a function f(a,b,c,d) with limits amin, amax, bmin, bmax, cmin, cmax, dmin, dmax. amin and amax are float. bmin and max are functions of a, cmin and cmax are functions of b, and dmin and dmax are functions of c. If scipy.integrate.quad() was just a for loop, then for each step, it would be possible to pass the value of a (or b, c, or d) into the limits such that they become floats.</p>
  
  <p>Is there a way to do this with scipy.integrate.quad? So far, I've tried simple nesting:</p>
  
  <pre><code>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
  </code></pre>
  
  <p>In this case, I still get the same error I was encountering before <em>ValueError: Invalid limits given</em> which seems to arise because I have symbols which are not being integrated over. I also tried using nquad, but I encountered the same errors.</p>
  
  <p>This is my current attempt at a process to integrate iteratively outward, only doing the numerics at the end:</p>
  
  <pre><code>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
  </code></pre>
  
  <p>The difficulty is that min1, max1 are functions of var2 and the function I'm actually working with appears to not have an analytic solution.</p>
  
  <p>Would it also help if I instead integrated from the outermost in as opposed to the innermost out?</p>
  
  <p>Thanks to Kazemakase, whose answer was helpful in sorting out my issue! I solved it using code slightly modified from what he wrote, so I've included that here as a reference to anyone else with a similar issue in the future.</p>
  
  <pre><code>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
  </code></pre>