Nesting numerical integrals in python

Question

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).

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.

Is there a way to do this with scipy.integrate.quad? So far, I've tried simple nesting:

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

In this case, I still get the same error I was encountering before ValueError: Invalid limits given 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.

This is my current attempt at a process to integrate iteratively outward, only doing the numerics at the end:

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

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.

Would it also help if I instead integrated from the outermost in as opposed to the innermost out?

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.

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 thanks very much for the comments, I've tried to address it a little more. Please let me know if more is needed! — PhysicistAbroad
– PhysicistAbroad, Commented May 12, 2017 at 1:04

MB-F · Accepted Answer · 2017-05-12 07:53:05Z

1

Nesting calls to quad is the correct approach. However, simple is not enough. Each argument passed to quad needs to be a function - not a result of a previous call. Here is an example that integrates f = x1**2 - x2**3 over the diamond-shaped area abs(x1) + abs(x2) <= 1.

We need to be able to integrate x2 between -/+(1 - abs(x1)) for any x1 between +/-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

Instead of explicitly writing a function for each integral you can use nquad that does the wrapping for you. It takes an integration range for each variable that can also be a function depending on other variables integrated over:

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)

answered May 12, 2017 at 7:53

MB-F

23.8k5 gold badges71 silver badges127 bronze badges

Sign up to request clarification or add additional context in comments.

2 Comments

PhysicistAbroad Over a year ago

Thanks, played around with this a little more and it worked! I've posted the code that worked in my question above for future reference.

MB-F Over a year ago

@PhysicistAbroad I'm glad this helped you. Nested functions actually mirror the concept of nested integrals nicely.

user8153 · Accepted Answer · 2017-05-12 07:25:33Z

1

If I understand you correctly then nquad should be able to handle this integration. Note that its ranges argument allows you to specify functions that depend on the values of other integration variables.

For example, to calculate the integral of f(x,y)=1 over the area 0 <= y <= 1 and 0 <= x <= y (in other words, we calculate the area of the upper-left half triangle of the unit square):

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)))

This gives 0.5 for the value of the integral, as it should.

answered May 12, 2017 at 7:25

user8153

4,1251 gold badge11 silver badges18 bronze badges

Collectives™ on Stack Overflow

Nesting numerical integrals in python

2 Answers 2

2 Comments

Comments

Your Answer

Hot Network Questions

Collectives™ on Stack Overflow

2 Answers 2

2 Comments

Comments

Your Answer

Sign up or log in

Post as a guest

Related