[SciPy-user] Gaussian quadrature error
Anne Archibald
peridot.faceted@gmail....
Wed Sep 17 11:00:09 CDT 2008
2008/9/14 Dan Murphy <chiefmurph@comcast.net>:
> I am trying out the integrate.quadrature function on the function f(x)=e**x
> to the left of the y-axis. If the lower bound in not too negative, I get a
> reasonable answer, but if the lower bound is too negative, I get 0.0 as the
> value of the integral. Here is the code:
>
>
>
> from scipy import *
>
>
>
> def f(x):
>
> return e**x
>
>
>
> integrate.quadrature(f,-10.0,0.0) # answer is (0.999954600065,
> 3.14148596026e-010)
>
>
>
> but
>
>
>
> integrate.quadrature(f,-1000.0,0.0) # yields (8.35116510531e-090,
> 8.35116510531e-090)
>
>
>
> Note that 'val' and 'err' are equal. Is this a bug in quadrature?
No, unfortunately. It is a limitation of numerical quadrature in
general. Specifically, no matter how adaptive the algorithm is, it can
only base its result on a finite number of sampled points of the
function. If these points are all zero to numerical accuracy, then the
answer must be zero. So if you imagine those samples are taken at the
midpoints of 10 intervals evenly spaced between -1000 and 0, then the
rightmost one returns a value of e**(-50), which is as close to zero
as makes no nevermind. You might be all right if this were an adaptive
scheme and if it used the endpoints, since one endpoint is guaranteed
to give you one. But not using the endpoints is a design feature of
some numerical integration schemes.
The take-home lesson is that you can't just use numerical quadrature
systems blindly; you have to know the features and limitations of the
particular one you're using. Gaussian quadrature can be very accurate
for smooth functions, but it has a very specific domain of
applicability. scipy.optimize.quad is a little more general-purpose by
intent (and necessarily a little less efficient when Gaussian
quadrature will do) but it can be tricked too.
A more specific take-home lesson is to try to normalize your problem
as much as possible, so that all quantities you feed your integrator
are of order unity. Yes, it's a pain to have to handle scale factors
yourself, particularly in the normal case when you're solving a family
of related problems. But you'll get much more reliable performance.
Anne
More information about the SciPy-user
mailing list