1

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: "IntegrationWarning: The integral is probably divergent, or slowly convergent." (it is known to be convergent) or "IntegrationWarning: The maximum number of subdivisions (50) has been achieved."

The equation is from: http://dx.doi.org/10.1029/WR003i001p00241

It is also available here: http://www.aqtesolv.com/papadopu.htm

Thank you for any assistance you have in numerical integration of fussy functions in python...

CODE EXAMPLE

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()
tomatessechees
  • 45
  • 2
  • 6
  • Try breaking up the integral into pieces. For example, you could sum the integrals over, say, [0.0, 0.1], [0.1, 1], [1, 10], and [10, inf]; experiment with the intervals to find one that works. Also, when you get the warning about the maximum number of subdivisions, try increasing that limit by using the `limit` argument, e.g. `limit=200`. – Warren Weckesser Jan 19 '18 at 18:02
  • The provided example is not minimal. I see `quad` used in `s()`, and `s()` is then plotted; no idea what's happening there. Try picking out _one_ function that you want to integrate and maybe show a plot. – Nico Schlömer Feb 15 '18 at 23:07
  • @Nico, It is a very specific function I am trying to integrate (`int_f`, see the link in the comment) for a variety of values of `t` and `r`. There is only _one_ function being integrated, so I do not see why you take issue with the example. This example shows how `s` jumps around with changing `r` and `t` values due to `quad` having problems with this integrand. – tomatessechees Feb 17 '18 at 14:24
  • I think one could simplify the example by writing something like: "Here's the function `int_f`, it looks so and so, and I want to integrate it. However, I'm getting the error such and such when running the above code." Anything more than that (like the parameters `t` and `r`) is distracting, imho. – Nico Schlömer Feb 17 '18 at 14:28

1 Answers1

6

I have bad news for you: Oscillatory integrands must be dealt with by special methods, for example Filon's method. Reading through the scipy docs, I see no way to do this. However, it appears MPMath implements precisely what you are looking for, and even gives an example using Bessel functions.

user14717
  • 4,757
  • 2
  • 44
  • 68
  • Thank you for this. I have not used that package before and will give it a go. I agree that scipy looks to be somewhat limited with respect to the implementation of various numerical methods adapted to a variety of "special" cases. – tomatessechees Jan 22 '18 at 13:13