0

I have this 2-dimensional integral with dependent limits. The function can be defined in Python as

def func(gamma, u2, u3):
    return (1-1/(1+gamma-u3-u2))*(1/(1+u2)**2)*(1/(1+u3)**2)

where the limits of u3 is from 0 to gamma (a positive real number), and the limits of u2 is from 0 to gamma-u3.

How can I implement this using scipy.integrate.nquad? I tried to read the documentation, but it was not easy to follow, especially I am relatively new to Python.

Extension: I would like to implement a numerical integration for an arbiraty K, where the integrand in this case is given by (1-1/(1+gamma-uk-....-u2))*(1/(1+uK)**2)*...*(1/(1+u2)**2). I wrote the function that takes a dynamic number of arguments as follows:

def integrand(gamma, *args):
    '''
    inputs:
     - gamma
     - *args = (uK, ..., u2)

    Output:
     - (1-1/(1+gamma-uk-....-u2))*(1/(1+uK)**2)*...*(1/(1+u2)**2)
    '''
    L = len(args)
    for ll in range(0, L):
        gamma -= args[ll]
    func = 1-1/(1+gamma)
    for ll in range(0, L):
        func *= 1/((1+args[ll])**2)
    return func

However, I am not sure how to do the same for the ranges, where I will have one function for the ranges, where uK ranges from 0 to gamma, u_{K-1} ranges from 0 to gamma-uK, ...., u2 ranges from 0 to gamma-uK-...-u2.

BlackMath
  • 932
  • 7
  • 24

1 Answers1

0

Here is a simpler method using scipy.integrate.dblquad instead of nquad:

Return the double (definite) integral of func(y, x) from x = a..b and y = gfun(x)..hfun(x).

from  scipy.integrate import dblquad

def func(u2, u3, gamma):
    return (1-1/(1+gamma-u3-u2))*(1/(1+u2)**2)*(1/(1+u3)**2)


gamma = 10

def gfun(u3):
    return 0

def hfun(u3):
    return gamma-u3

dblquad(func, 0, gamma, gfun, hfun, args=(gamma,))

It seems that gfun and hfun do not accept the extra arguments, so gamma has to be a global variable.

Using nquad, after many trial and error:

from  scipy.integrate import nquad

def func(u2, u3, gamma):
    return (1-1/(1+gamma-u3-u2))*(1/(1+u2)**2)*(1/(1+u3)**2)

def range_u3(gamma):
    return (0, gamma)

def range_u2(u3, gamma):
    return (0, gamma-u3)

gamma = 10
nquad(func, [range_u2, range_u3], args=(gamma,) )

Useful quote from the source code of tplquad:

# nquad will hand (y, x, t0, ...) to ranges0
# nquad will hand (x, t0, ...) to ranges1

And from the nquad documentation, the order of the variables is reversed (same for dblquad):

Integration is carried out in order. That is, integration over x0 is the innermost integral, and xn is the outermost

Generic case with k nested integrations:

from  scipy.integrate import nquad
import numpy as np

def func(*args):
    gamma = args[-1]
    var = np.array(args[:-1])

    return (1-1/(1+gamma-np.sum(var)))*np.prod(((1+var)**-2))

def range_func(*args):
    gamma = args[-1]
    return (0, gamma-sum(args[:-1]))

gamma, k = 10, 2
nquad(func, [range_func]*k, args=(gamma,) )
xdze2
  • 3,986
  • 2
  • 12
  • 29
  • My example of double integral is a special case. I am interested in the more general case of (K-1)-nested integral. What does the `args=(gamma,)` do? – BlackMath Aug 01 '18 at 04:54
  • Also, is there anyway to generalize the above code for `nquad` for an arbitrary integer `K`? In this case the function will be `def func(gamma, u2, u3, ..., uK): return (1-1/(1+gamma-uk-....-u2))*(1/(1+uK)**2)*...*(1/(1+u2)**2)`, and the limits for uK is from 0 to gamma, for u_{K-1} from 0 to gamma-uK, ...., and for u2 from 0 to gamma-uK-...-u3. – BlackMath Aug 01 '18 at 04:57
  • `gamma` is not a variable of the integration but a parameter. The `args` option is used to pass extra parameters to the function, so `gamma` – xdze2 Aug 01 '18 at 05:05
  • It should be possible to write it for `k` number of variables... by using a generic range function accepting arbitrary input arguments: `range(*args): ...` – xdze2 Aug 01 '18 at 05:12
  • I think `func` is easy to construct for a generic `K`, but I cannot say the same for the range function. I mean, how to have one range function? and how to call it from the `quad` method? – BlackMath Aug 01 '18 at 05:26
  • I updated the original post by adding an extension to the original problem, which is a generalization to the special case presented originally. I appreciate if you could give me some more hints on how to handle the ranges for the generic case. – BlackMath Aug 01 '18 at 05:57
  • Thanks for the update. I like how you vectorized the code using numpy. However, in your code k is actually (k-1). From your code, I understand that `gamma` is the last elemnet in `args`, right? I get how `func` works, but I still don't understand how `range_func` works? When we call the `quand` method, what is the sequence of steps that happen? – BlackMath Aug 01 '18 at 18:38
  • `nquad` start by calling the first range function only with `gamma` as argument, it obtains the first range (0, gamma). The first integration can start, however because the integrations are nested, to do so the other integrations have to be performed. It is like a recursive call of a function. [Here](https://stackoverflow.com/a/51305099/8069403) is an example of the integration over a volume using 3 nested `quad`, it may be easier to understand. Your right: `k-1` is `k` in my code, and `gamma` is the last element of `args` – xdze2 Aug 02 '18 at 08:07
  • OK, so, the first call is with gamma only, but how the range function knows that the next range is from `0` to `gamma-uK`, and so one? In other words, how `args` changes from `gamma` to `uK, gamma` ... `uK, ..., u2, gamma` in the code? It was more clear to me for the special case `K=3`. – BlackMath Aug 02 '18 at 14:20
  • Can we adapt the above code to the case presented here: https://stackoverflow.com/questions/51958138/nested-integral-in-python? – BlackMath Aug 23 '18 at 03:24