Eu tenho essa integral bidimensional com limites dependentes. A função pode ser definida em Python como

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

onde os limites deu3 é de 0 agamma (número real positivo) e os limites deu2 é de 0 agamma-u3.

Como posso implementar isso usandoscipy.integrate.nquad? Eu tentei ler a documentação, mas não foi fácil de seguir, especialmente eu sou relativamente novo no Python.

Extensão: gostaria de implementar uma integração numérica para uma arbitragemK, onde o integrando, neste caso, é fornecido por(1-1/(1+gamma-uk-....-u2))*(1/(1+uK)**2)*...*(1/(1+u2)**2). Eu escrevi a função que recebe um número dinâmico de argumentos da seguinte maneira:

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

No entanto, não sei como fazer o mesmo para os intervalos, onde terei uma função para os intervalos, em queuK varia de 0 agamma, u_{K-1} varia de 0 agamma-uK, ....,u2 varia de 0 agamma-uK-...-u2.