[SciPy-user] Gaussian quadrature error

Anne Archibald peridot.faceted@gmail....
Mon Oct 6 21:23:32 CDT 2008


2008/10/6 Dan Murphy <chiefmurph@hotmail.com>:
> Thanks for the clear explanation, Anne. If I understand the second take-home
> lesson -- feed integrators unity-order quantities -- then instead of calling
> quadrature as I did:
>    integrate.quadrature(f,-1000.0,0.0)
> I should scale my function f so that my call looks more like
>    integrate.quadrature(f,-1.0,0.0)
> Was that what you meant?

Yes, though with the particular function you used - exp(x) - you then
have the problem of exceedingly rapid changes within the region of
integration. If what you actually wanted was to integrate from
-infinity to zero, then you might do better to rescale your x
coordinate nonlinearly, say by x' = 1/(1-x) so that the integration
becomes one from zero to one without a sharp edge near one side.

A further take-home lesson is that Gaussian quadrature is extremely
efficient for functions that behave like high-order polynomials -
which exp(x) does on small intervals but not on large intervals.
However, you can use customized Gaussian quadrature based on a
different set of polynomials to integrate functions that look like a
polynomial times a weight factor.

> Also, when you say "scipy.optimize.quad is a little more general-purpose"
> did you mean "scipy.integrate.quad"? Indeed, the call
>    integrate(quad(f,-1000.0,0.0)
> worked great!

Oops. Yes. It's still possible to trick it, but scipy.integrate.quad
is based on old, well-tested code that is designed to be fairly robust
against peculiar functions. It will be a little slower than Gaussian
quadrature for functions that are extremely smooth, but then Gaussian
quadrature suffers very badly when handed functions like abs(x).

Anne


More information about the SciPy-user mailing list