How to integrate a discontinuous function by tplquad or nquad

[Safinazsalem]

I try to triple integrate the following function but it gives many errors:

IntegrationWarning: The integral is probably divergent, or slowly convergent.

IntegrationWarning: The maximum number of subdivisions (50) has been achieved.
  If increasing the limit yields no improvement it is advised to analyze 
  the integrand in order to determine the difficulties. 

IntegrationWarning: The algorithm does not converge.  Roundoff error is detected
  in the extrapolation table.  It is assumed that the requested tolerance
  cannot be achieved, and that the returned result (if full_output = 1) is 
  the best which can be obtained.

The function is f(t1,t2,m)

(-(50*t1**2*(-16*(16 + 20*t1 + 6*t1**2 + 4*t1**3 + 3*t1**4) + (1 + t1)**2*(448 + 1568*t1 + 2400*t1**2 + 1920*t1**3 + 864*t1**4 + 240*t1**5 + 72*t1**6 + 42*t1**7 + 18*t1**8 + 3*t1**9)*math.exp((4*(1 + t1))/(2 + t1))) + m*(1 + t1)*(-128 + 1536*t1**2*math.exp((4*(1 + t1))/(2 + t1)) + 2112*t1**4*math.exp((4*(1 + t1))/(2 + t1)) + 1008*t1**5*math.exp((4*(1 + t1))/(2 + t1)) + 312*t1**6*math.exp((4*(1 + t1))/(2 + t1)) + 114*t1**7*math.exp((4*(1 + t1))/(2 + t1)) + 60*t1**8*math.exp((4*(1 + t1))/(2 + t1)) + 21*t1**9*math.exp((4*(1 + t1))/(2 + t1)) + 3*t1**10*math.exp((4*(1 + t1))/(2 + t1)) + 64*t1*(-1 + 6*math.exp((4*(1 + t1))/(2 + t1))) + 32*t1**3*(-1 + 78*math.exp((4*(1 + t1))/(2 + t1)))))/(4.*(1 + t1)*(2 + t1)*(-48 - 96*t1 + 8*(-3 + 2*m**2)*t1**2 + 8*(3 + 2*m**2)*t1**3 + (6 + 4*m**2)*t1**4 + 6*t1**5 + 3*t1**6))) * ((50 * t1**2) +5* (-(50*t1**2*(-16*(16 + 20*t1 + 6*t1**2 + 4*t1**3 + 3*t1**4) + (1 + t1)**2*(448 + 1568*t1 + 2400*t1**2 + 1920*t1**3 + 864*t1**4 + 240*t1**5 + 72*t1**6 + 42*t1**7 + 18*t1**8 + 3*t1**9)*math.exp((4*(1 + t1))/(2 + t1))) + (k - m)*(1 + t1)*(-128 + 1536*t1**2*math.exp((4*(1 + t1))/(2 + t1)) + 2112*t1**4*math.exp((4*(1 + t1))/(2 + t1)) + 1008*t1**5*math.exp((4*(1 + t1))/(2 + t1)) + 312*t1**6*math.exp((4*(1 + t1))/(2 + t1)) + 114*t1**7*math.exp((4*(1 + t1))/(2 + t1)) + 60*t1**8*math.exp((4*(1 + t1))/(2 + t1)) + 21*t1**9*math.exp((4*(1 + t1))/(2 + t1)) + 3*t1**10*math.exp((4*(1 + t1))/(2 + t1)) + 64*t1*(-1 + 6*math.exp((4*(1 + t1))/(2 + t1))) + 32*t1**3*(-1 + 78*math.exp((4*(1 + t1))/(2 + t1)))))/(4.*(1 + t1)*(2 + t1)*(-48 - 96*t1 + 8*(-3 + 2*(k - m)**2)*t1**2 + 8*(3 + 2*(k - m)**2)*t1**3 + (6 + 4*(k - m)**2)*t1**4 + 6*t1**5 + 3*t1**6)))) 

I wonder is there a way to make python integrate this function smoothly.

I mean in quad one can pass various parameters to integrate.quad to change the absolute tolerance by epsrel or epsabs. Or one can discard the singularity points by points=. But in triple integral how to do so by tplquad or nquad . Also may be using if: else will be helpful with this function, but I wonder how ?

Here is the code I use:

from scipy import integrate 
import math
import scipy

k=3;a=0.1;nab=10**-9;n=0.94;kz= 0.002;
# Define the function to be integrated
def f(m,t1,t2):
    return (-(50*t1**2*(-16*(16 + 20*t1 + 6*t1**2 + 4*t1**3 + 3*t1**4) + (1 + t1)**2*(448 + 1568*t1 + 2400*t1**2 + 1920*t1**3 + 864*t1**4 + 240*t1**5 + 72*t1**6 + 42*t1**7 + 18*t1**8 + 3*t1**9)*math.exp((4*(1 + t1))/(2 + t1))) + m*(1 + t1)*(-128 + 1536*t1**2*math.exp((4*(1 + t1))/(2 + t1)) + 2112*t1**4*math.exp((4*(1 + t1))/(2 + t1)) + 1008*t1**5*math.exp((4*(1 + t1))/(2 + t1)) + 312*t1**6*math.exp((4*(1 + t1))/(2 + t1)) + 114*t1**7*math.exp((4*(1 + t1))/(2 + t1)) + 60*t1**8*math.exp((4*(1 + t1))/(2 + t1)) + 21*t1**9*math.exp((4*(1 + t1))/(2 + t1)) + 3*t1**10*math.exp((4*(1 + t1))/(2 + t1)) + 64*t1*(-1 + 6*math.exp((4*(1 + t1))/(2 + t1))) + 32*t1**3*(-1 + 78*math.exp((4*(1 + t1))/(2 + t1)))))/(4.*(1 + t1)*(2 + t1)*(-48 - 96*t1 + 8*(-3 + 2*m**2)*t1**2 + 8*(3 + 2*m**2)*t1**3 + (6 + 4*m**2)*t1**4 + 6*t1**5 + 3*t1**6))) * ((50 * t1**2) +5* (-(50*t1**2*(-16*(16 + 20*t1 + 6*t1**2 + 4*t1**3 + 3*t1**4) + (1 + t1)**2*(448 + 1568*t1 + 2400*t1**2 + 1920*t1**3 + 864*t1**4 + 240*t1**5 + 72*t1**6 + 42*t1**7 + 18*t1**8 + 3*t1**9)*math.exp((4*(1 + t1))/(2 + t1))) + (k - m)*(1 + t1)*(-128 + 1536*t1**2*math.exp((4*(1 + t1))/(2 + t1)) + 2112*t1**4*math.exp((4*(1 + t1))/(2 + t1)) + 1008*t1**5*math.exp((4*(1 + t1))/(2 + t1)) + 312*t1**6*math.exp((4*(1 + t1))/(2 + t1)) + 114*t1**7*math.exp((4*(1 + t1))/(2 + t1)) + 60*t1**8*math.exp((4*(1 + t1))/(2 + t1)) + 21*t1**9*math.exp((4*(1 + t1))/(2 + t1)) + 3*t1**10*math.exp((4*(1 + t1))/(2 + t1)) + 64*t1*(-1 + 6*math.exp((4*(1 + t1))/(2 + t1))) + 32*t1**3*(-1 + 78*math.exp((4*(1 + t1))/(2 + t1)))))/(4.*(1 + t1)*(2 + t1)*(-48 - 96*t1 + 8*(-3 + 2*(k - m)**2)*t1**2 + 8*(3 + 2*(k - m)**2)*t1**3 + (6 + 4*(k - m)**2)*t1**4 + 6*t1**5 + 3*t1**6)))) 


# Define the limits of integration
a = 0.1
b = 5


# Perform the integration
result, error = integrate.tplquad(f, 0.1, 2, lambda x: 0.1, lambda x: 2, lambda x,y: 0.1, lambda x,y: 2 )

print("result:", result)

Any help is appreciated!

[Gribouillis]

I see that the function contains denominators. Does the function have poles in the integration interval?

[Safinazsalem]

The function is given by:

$ f(t1,t2,m) = \psi(k,t) \phi(k,t) + 5 \psi(k,t) \psi(k-m ,t)$

where $\psi(k,t)$ is given in the attachments, and $\phi(t,k)= 50 * t**2 +k$.


You notice that $\psi(k,t)$ may has many divergences in the dominator , so I wonder how to solve that