[Python] Why does this code takes a long time to be excuted then it gives IntegrationWarning

The following python code takes long time to be executed, then it gives the following error, and never gives an output:

/usr/local/lib/python3.11/dist-packages/scipy/integrate/_quadpack_py.py:1260: IntegrationWarning: The integral is probably divergent, or slowly convergent.
quad_r = quad(f, low, high, args=args, full_output=self.full_output,
/usr/local/lib/python3.11/dist-packages/scipy/integrate/_quadpack_py.py:1260: 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. If the position of a
local difficulty can be determined (singularity, discontinuity) one will
probably gain from splitting up the interval and calling the integrator
on the subranges. Perhaps a special-purpose integrator should be used.
quad_r = quad(f, low, high, args=args, full_output=self.full_output,


So any help is appreciated to optimize the code and make it executed probably.

Here is the code:

from functools import partial
import numpy as np
import scipy.special as special
from scipy import integrate
from scipy.special import kn, hyp1f1 # Import hyp1f1
import matplotlib.pyplot as plt
import math
import time, sys
from scipy.integrate import nquad
from functools import partial


a=0.01;nab=10**-9;kz=0.002;H=10**-5;start= 100 *H; end=10**9 *H;step=10**8 *H;
ti =10 *H;
tf= 100 *H;
########################################################


def psi1(m,t1):
return 0.1 + t1 * (0.1 + t1 * (-6.4576 + m**2 * (-0.625 - 21.66203 * t1) - 100.282 * t1 - 1.25 * m**4 * t1))

def phi1(m,t1,k):
return 0.1 + t1 * (10.30555 + 6.1 * (k-m)**6 * t1**2 + (k-m)**4 * t1 * (3. + 148.022222 * t1) + t1 * (224.953549 + 4034.491 * t1) + (k-m)**2 * (1. + t1 * (50.45555 + 1284.838 * t1)))

def psik1(m,t1,k):
return 0.1 + t1 * (0.1 + t1 * (-6.4576 + (k-m)**2 * (-0.625 - 21.66203 * t1) - 100.282 * t1 - 1.25 * (k-m)**4 * t1))


# f1(m,t1,k) = - psi1(m,t1) phi1(k-m,t1) - 5 psi1(m,t1) psi1(k-m,t1)

def f1(m,t1,k):
return - psi1(m, t1) * phi1(m, t1, k) - 5.0 * psi1(m, t1) * psik1(m,t1,k)

#############################################

def psi2(m,t2):
return 0.1 + t2 * (0.1 + t2 * (-6.4576 + m**2 * (-0.625 - 21.66203 * t2) - 100.282 * t2 - 1.25 * m**4 * t2))

def phi2(m,t2,k):
return 0.1 + t2 * (10.30555 + 6.1 * (k-m)**6 * t2**2 + (k-m)**4 * t2 * (3. + 148.022222 * t2) + t2 * (224.953549 + 4034.491 * t2) + (k-m)**2 * (1. + t2 * (50.45555 + 1284.838 * t2)))

def psik2(m,t2,k):
return 0.1 + t2 * (0.1 + t2 * (-6.4576 + (k-m)**2 * (-0.625 - 21.66203 * t2) - 100.282 * t2 - 1.25 * (k-m)**4 * t2))


# f2(m,t2,k) = - psi2(m,t2) phi2(k-m,t2) - 5 psi2(m,t2) psi2(k-m,t2)

def f2(m,t2,k):
return - psi2(m, t2) * phi2(m, t2, k) - 5.0 * psi2(m, t2) * psik2(m,t2,k)

#############################################


def psim2(m,t2,k):
return 0.1 + t2 * (0.1 + t2 * (-6.4576 + (k-m)**2 * (-0.625 - 21.66203 * t2) - 100.282 * t2 - 1.25 * (k-m)**4 * t2))

def phim2(m,t2):
return 0.1 + t2 * (10.30555 + 6.1 * m**6 * t2**2 + m**4 * t2 * (3. + 148.022222 * t2) + t2 * (224.953549 + 4034.491 * t2) + m**2 * (1. + t2 * (50.45555 + 1284.838 * t2)))

def psikm2(m,t2):
return 0.1 + t2 * (0.1 + t2 * (-6.4576 + m**2 * (-0.625 - 21.66203 * t2) - 100.282 * t2 - 1.25 * m**4 * t2))


# f(k,(k-m),t2) = - psim2(m,t2) phim2(m,t2) - 5 psi2(m,t2) psi2(k-m,t2)

def fm2(m,t2,k):
return - psim2(m, t2,k) * phim2(m, t2) - 5.0 * psim2(m, t2,k) * psikm2(m,t2)

#############################################

def ffmk(t1,t2,m,k):
return f1(m,t1,k) * (fm2(m,t2,k) + f1(m,t1,k))


###############################################

def parf(t1,t2,m,mu,a,k):
# Avoid division by zero by adding a small value (epsilon) to t1 and t2
epsilon = 1e-10 # A small value to avoid division by zero
return (H * a)**2 * k**3 * m**3 *(1/(t1*H)) * (1/(t2*H)) * (4/9)**2 * nab**2 * (1/kz)*0.2 * (1 - 2 * mu**2 + mu**4) * (1/(k**2 - 2 * k * mu * m + m**2)**(3/2))


def ptot(t1,t2,m,mu,a,k):
return parf(t1,t2,m,mu,a,k) * ffmk(t1,t2,m,k)


def rho(t1,t2,m,mu,a,k):
return k**2 * (H * a * 10*9)**2 * ptot(t1,t2,m,mu,a,k)

def omega(t1,t2,m,mu,a,k):
return k**2 * ptot(t1,t2,m,mu,a,k) * 10**18 / (H * 10**12 )**2

###############################################

# Define the integration bounds as a list of tuples
ranges = [[-H, a] ,[-H, a],[0.1, 1],[-1,1]]
rang= [[-H, a] ,[-H, a],[0.1, 1]]


X = np.arange(0.1* H,end,step)
X2 = np.arange(0.1,10,0.1)

####################################################
# Fixing the variable name from t_f to tf
delta_t1 = - H + a
delta_t2 = -ti + a

####################################################

def p21(a, k):
"""
Calculates the integral of p(t, y, s) with respect to s from 1 to 2.
"""
# Use quad to integrate p with respect to s
integral_result, error = nquad(rho,[[-H, a] ,[-H, a] ,[0.1, 1],[-1,1]],args=(a,k) )
return integral_result

def p22(a, k):
"""
Calculates the integral of p(t, y, s) with respect to s from 1 to 2.
"""
# Use quad to integrate p with respect to s
epsilon = 1e-10
integral_result1, error1 = nquad(rho,[[-ti + epsilon, a] ,[-ti + epsilon, a] ,[0.1, 1],[-1,1]],args=(a,k) )
return integral_result1

def F1(x):
# Initialize res with the same shape as x to store results
res = np.zeros_like(x, dtype=float)
epsilon = 1e-10
for i,val in enumerate(x):
y,err = nquad(p21, [[H, a]] ,args=(val/(2*np.pi),))
res= y / delta_t1
return res


def F2(x):
# Initialize res with the same shape as x to store results
res1 = np.zeros_like(x, dtype=float)
epsilon = 1e-10
for i,val in enumerate(x):
y,err = nquad(p22, [[ti + epsilon, a]] ,args=(val/(2*np.pi),))
res1= y / delta_t2
return res1




X2=X/H

# Convert X2 to a NumPy array if it's not already
X2 = np.asarray(X2)

# Replace NaN or Inf values in F(X), F2(X), and F3(X) with 0
F1X = np.nan_to_num(F1(X))
F2X = np.nan_to_num(F2(X))


##########################################
plt.plot(X,F1X, color='b', label='$\eta_a=-H$') # Plot using FX
plt.plot(X,F2X, color='red', label='$\eta_a=-tf$') # Plot using FX
plt.xlabel(" $f~ Hz$ ")
plt.ylabel(r"$\rho_h(k,\tau \rightarrow 0)$") # Use raw string and correct LaTeX syntax
plt.xscale('log')
plt.yscale('log')
plt.legend()
plt.savefig('Pi-k.png')

Continue reading...