+- Python Forum (https://python-forum.io)
+-- Forum: Python Coding (https://python-forum.io/forum-7.html)
+--- Forum: Homework (https://python-forum.io/forum-9.html)
+--- Thread: How to integrate a discontinuous function by tplquad or nquad (/thread-42201.html)
How to integrate a discontinuous function by tplquad or nquad - Safinazsalem - May-27-2024
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!
RE: How to integrate a discontinuous function by tplquad or nquad - Gribouillis - May-27-2024
I see that the function contains denominators. Does the function have poles in the integration interval?
RE: How to integrate a discontinuous function by tplquad or nquad - Safinazsalem - May-29-2024
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