[SciPy-User] Warnings with QUAD
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
<etrade.griffiths@dsl.pipex.com> wrote:
> 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.
In some cases the integrate.quad might still give correct answers in
spite of the warning and singularity.
Josef
> _______________________________________________
> SciPy-User mailing list
> SciPy-User@scipy.org
> http://mail.scipy.org/mailman/listinfo/scipy-user
>
More information about the SciPy-User
mailing list