[SciPy-user] Usage of integrate.quadrature

[Anne Archibald]

On 14/06/07, Nils Wagner <nwagner at iam.uni-stuttgart.de> wrote:

> How can I use integrate.quadrature to compute
>
> \int\limits_0^1 h(xi) d xi
>
> or
>
> \int\limits_0^1 h(xi) h^T(xi) d xi,
>
> where h(x) is a v e c t o r-valued function, e.g.

The short answer is, you can't. Also, you should be using
integrate.quad if you can, but it won't do this either.

> def h(xi,le):
>     """ Shape functions """
>     tmp = zeros(4,float)
>     tmp[0] = 1-3*xi**2+2*xi**3
>     tmp[1] = (xi-2*xi**2+xi**3)*le
>     tmp[2] = 3*xi**2-2*xi**3
>     tmp[3] = (-xi**2+xi**3)*le
>     return tmp
>
> integrate.quadrature(?,0,1,args=(le))
>
> Any pointer would be appreciated.

Well, vector-valued integrals are by definition computed
componentwise, so if your function really does look like this you
could just integrate each piece separately (using
scipy.integrate.quadrature or scipy.integrate.quad, which is faster
and safer, though these are such simple functions...). In fact if this
is your function, you can find the integrals analytically, which will
be vastly faster and more robust (if your calculus is rusty, use MAPLE
or the online integrator or whatever).

If you have a function that unavoidably computes a vector every time -
perhaps it's based on a Fourier transform (though remember the FT is
linear!) or some other process - I don't think scipy has anything that
will do the integral for you. But it will not be hard to be as smart
as scipy.integrate.quadrature: all it does is Gaussian integration to
higher and higher order until it converges or gives up. To do Gaussian
quadrature, pick the orthogonal polynomials that correspond to your
weight function (probably the Legendre), evaluate your vector function
at the roots given by the orthogonal polynomial object, add the
results up weighted by the weights given by the orthogonal polynomial
object, and you've got the answer. To make it adaptive, just keep
increasing the order until you think it's converged as much as it's
going to.

Anne

P.S. Be careful with high-order orthogonal polynomials - the
recurrence relations scipy uses to compute with them begin to have
terrible roundoff error when you get to high orders (above maybe
fifty). I had to implement the Chebyshev polynomials myself as cos(n
arccos x) to get enough accuracy for my problem. -A