[SciPy-user] Integration with precalculated values

[Anne Archibald]

2008/6/16 Mico Filós <elmico.filos at gmail.com>:
> Hi,
>
> I need to evaluate this function for an array of ys
>
>    f(y) = integral( exp(x*x)  * psi(x), A -y, B -y) # 2nd and 3rd
> arguments: lower and upper limits
>
> where
>
>    psi(x) = integral( exp(t*t) * ( 1 + erf(t) )**2 , -Inf, x)
>
> Apart from the fact that the integrand blows up easily, I would gain
> some speed by
> precalculating psi for the ranges I need for the computation of f(y).
> Is there any proper way to do that?
>
> I would really appreciate any hint or suggestion on how to evaluate
> the integral :)

Hmm.

I would start with an implementation you know to be reliable; it may
be fast enough, and if not, it can serve as a good test case. For this
I'd use scipy.integrate.quad, more or less the way you wrote it there.

If you want to accelerate things by precalculating psi, how you handle
it depends on the y values and the accuracy you need. The most direct
way assumes you have no control of the y values; then you can use a
fairly standard trick to replace psi(t) by a function that is faster
to calculate. (At least, I often use it.) Just compute the function at
a reasonable set of values (I'm afraid you'll have to figure out what
qualifies as "reasonable"), and use scipy.interpolate.splrep to fit a
spline through them. Then instead of evaluating psi(x) directly, you
can evaluate the spline (which will be very quick). Since your
function is always positive and probably varies a lot in magnitude,
I'd fit a spline to the log and then exponentiate on evaluation. Also,
you can drastically accelerate the evaluation of psi(t) at a sequence
of values by using an ODE integrator (say scipy.integrate.odeint)
once. (Since integration needs to walk across t values anyway, you
might as well use an integrator that can keep track of all the values.
Also, if it's psi that likes to blow up, odeint can cope with that.
See its docstring.) Given a fast way to evaluate psi(t), quad should
allow you to evaluate the function you actually want rather quickly.

If this is too slow, you can probably do something involving more
carefully discretizing your two functions, but it will probably
require you to restrict the range of y values.

Looking at your function, I strongly suspect things can be simplified
by rewriting it as a multidimensional integral and rethinking how you
do the integral. In particular, I'd look to see if anything can be
done with rewriting the erf as an integral, since its exp(-u**2) might
help keep the other exp(x**2) and exp(t**2) under control. But as
always, it's a tradeoff between work for you and work for the machine.

Anne