title: Numerical integration of highly oscillating 1-D integrand (containing Bessel functions) in python
body: <p>I am attempting to numerically evaluate a real-valued integrand composed of several Bessel functions (of the first and second kind). The integrand is oscillating and decaying, and needs to be evaluated between 0 and +∞. So far, my attempts using the scipy.integrate subpackage (quad and fixed_quad) have been unsuccessful. The evaluated value jumps around when in reality it should be smooth. For certain sets of parameter values, I also receive warnings: <em>"IntegrationWarning: The integral is probably divergent, or slowly convergent."</em> (it is known to be convergent) or <em>"IntegrationWarning: The maximum number of subdivisions (50) has been achieved."</em></p>
  
  <p>The equation is from: <a href="http://dx.doi.org/10.1029/WR003i001p00241" rel="nofollow noreferrer">http://dx.doi.org/10.1029/WR003i001p00241</a></p>
  
  <p>It is also available here: <a href="http://www.aqtesolv.com/papadopu.htm" rel="nofollow noreferrer">http://www.aqtesolv.com/papadopu.htm</a></p>
  
  <p>Thank you for any assistance you have in numerical integration of fussy functions in python...</p>
  
  <p>CODE EXAMPLE</p>
  
  <pre><code>import matplotlib.pyplot as plt
  import numpy as np
  import pandas as pd
  from scipy import special as sps
  import scipy.integrate as integrate
  
  # define constants and variables (SI mks units):
  r_w = 0.15
  r_c = 0.16
  b = 10
  S_s = 1E-6
  Q = 0.001
  S = S_s*b
  K=1E-8
  T=K*b
  alpha = (r_w**2)*S/r_c**2
  def r_D(r):
      return r/r_w
  def u(r,t):
      return r**2*S/(4*T*t)
  def B(beta):
      return beta*sps.jv(0,beta) - 2*alpha*sps.jv(1,beta)
  def A(beta):
      return beta*sps.yn(0,beta) - 2*alpha*sps.yn(1,beta)
  def intd_f(beta,r,t): 
      TOP = (1-np.exp(-beta**2*r_D(r)**2)/(4*u(r,t)))*( sps.jv(0,beta*r_D(r))*A(beta) - sps.yn(0,beta*r_D(r))*B(beta) )
      BOT = (A(beta)**2 + B(beta)**2)*beta**2
      return TOP/BOT
  def s(r,t):
      banana = (2*alpha*Q)/(np.pi**2*K*b)
      apple = integrate.quad(intd_f, 0, np.inf, args=(r,t))[0]
      #apple = integrate.fixed_quad(intd_f, 0, 1E100, args=(r,t))[0] # another option I have tried
      return banana*apple
  
  #%% simple example usage...
  r=np.arange(1,10,.1)
  t=60*60*24*pd.Series([1/24,1,365,3650])
  plt.figure()
  for tt in t:
      print('time=',tt)
      snow=[]
      for rr in r:
          print('r=',rr)
          snow.append(s(rr,tt))
      plt.subplot(2,1,1)
      plt.plot(r,snow,label='t='+str(tt/(60*60*24))+'d')
      plt.subplot(2,1,2)
      plt.semilogy(r,np.abs(snow))    
  plt.subplot(2,1,1)
  plt.legend()
  </code></pre>