josef.pktd@gmai... josef.pktd@gmai...
Wed Apr 28 10:55:27 CDT 2010

```On Mon, Apr 26, 2010 at 12:47 PM, Etrade Griffiths
> Hi
>
> I am trying to integrate a function with an infinite upper limit but
> keep getting these warnings from scipy.integrate.quad:
>
> Warning: 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.
>
> Not quite sure what this means.  Increasing the number of
> sub-divisions had no effect and the components of the function seem
> reasonably well behaved.  As the following example shows, I get these
> messages over a wide range of values.  Would be grateful for any
> suggestions to stop these warnings - thanks in advance.
>
> import math
> import scipy.special
> import scipy.integrate
>
> # function to integrate
>
> def func(u, tDa):
>     num = 1.0 - math.exp(-u * u * tDa)
>     j0u = scipy.special.j0(u)
>     y0u = scipy.special.y0(u)
>     den = j0u * j0u + y0u * y0u
>     return num / den / u / u / u
>
> # main program
>
> print("    tD          wD")
>
> for tD in [1.0E-3,1.0E-2,1.0E-1,1.0E0,1.0E1,1.0E2,1.0E3,1.0E4,1.0E5]:
>     x = scipy.integrate.quad(func, 0.0, scipy.integrate.inf, args=(tD))
>     wD = 4.0 * x[0] / math.pi / math.pi
>     print("%10.3e  %10.3e" % (tD, wD) )

your function  a spike/singularity close to zero, e.g.
ud = np.linspace(1e-8,0.1)
fv = np.array([func(ui, 1.0E-3) for ui in ud])

import matplotlib.pyplot as plt
plt.plot(fv)

Someone offered a version of integrate.quad that can handle
singularities on the mailing list last year.

spite of the warning and singularity.

Josef

> _______________________________________________
> SciPy-User mailing list
> SciPy-User@scipy.org
> http://mail.scipy.org/mailman/listinfo/scipy-user
>
```