[SciPy-User] Quadrature of a vector-valued function

[Rob Falck]

Thanks for the quick reply.  The odeint solution is worth considering,
but I suspect it uses more intervals than a quadrature typically
would.  Also, the point about the adaptive one possibly being faster
when not vectorized is a good point, although in my case I don't think
any of the equations is substantially "faster" than the others.  I was
able to kludge fixed_quad to get the results I want.  This should be
extensible to the variable-step quadrature method, since it just calls
fixed-quad repeatedly.  Right now this doesn't utilize the "vectorize"
utility in scipy.integrate, but it gets the job done.

import numpy as N
from scipy.special.orthogonal import p_roots

def fixed_quadv(func,a,b,args=(),n=5,vec_func=False):
    """
Compute a definite integral using fixed-order Gaussian quadrature.

Integrate `func` from a to b using Gaussian quadrature of order n.

Parameters
----------
func : callable
A Python function or method to integrate.  Must have a form compatible
with func(x,*args).  If x may be a vector, then vec_func option should
be
True.  If func is vector-valued, then it should return an array of
shape
(size(x),2).
a : float
Lower limit of integration.
b : float
Upper limit of integration.
args : tuple, optional
Extra arguments to pass to function, if any.
n : int, optional
Order of quadrature integration. Default is 5.
vec_func : boolean, optional
True if func accepts vector input for x.

Returns
-------
val : float
Gaussian quadrature approximation to the integral

See Also
--------
fixed_quad: fixed-order Gaussian quadrature
quad : adaptive quadrature using QUADPACK
dblquad, tplquad : double and triple integrals
romberg : adaptive Romberg quadrature
quadrature : adaptive Gaussian quadrature
romb, simps, trapz : integrators for sampled data
cumtrapz : cumulative integration for sampled data
ode, odeint - ODE integrators

"""
    [x,w] = p_roots(n)
    x = N.real(x)
    ainf, binf = map(N.isinf,(a,b))
    if ainf or binf:
        raise ValueError("Gaussian quadrature is only available for "
                "finite limits.")
    if vec_func:
        fx = (((b-a)*x)+(a+b))/2.0
        accum = N.dot(w,func(fx,*args))
        y = accum*(b-a)/2.0
    else:
        accum = 0
        for i in range(n):
            accum += w[i] * func((((b-a)*x[i])+(a+b))/2.0, *args )
        y = accum*(b-a)/2.0
    return y, None

def vector_func(x):
    result = N.zeros([N.size(x),2])
    a = 3*x**2+2*x+1
    result[:,0] = a
    b = 2*x+3
    result[:,1] = b
    return result

print fixed_quadv(vector_func,0,2,vec_func=True)


On May 25, 3:29 pm, Warren Weckesser <warren.weckes... at enthought.com>
wrote:
> Hi Rob,
>
> On Wed, May 25, 2011 at 2:14 PM, Rob Falck <robfa... at gmail.com> wrote:
> > Hello,
>
> > I'm trying to use a quadrature to integrate several equations, and all
> > have similar up-front calculations.  However, quad as implemented in
> > scipy.integrate only seems capable of handling scalar functions,
> > meaning I have to break my vector-valued function up into several
> > scalar functions.  Also, since quad is adaptive, I can't necessarily
> > store the common calculations because theres no guarantee that quad
> > will use the same grid when integrating each function.
>
> For some use cases, splitting up the function is preferable, especially if
> the evaluating the functions is computationally expense.  If one component
> is changing rapidly while another isn't, extra time would be spent making
> unnecessary calls to the nice function.
>
>
>
>
>
> > As an example, in the following code I would like the integral of
> > vector_func to return N.array([14,10]), but instead it errors out :
>
> > import numpy as N
> > from scipy.integrate import quad
>
> > def scalar_func(x):
> >    return 3*x**2+2*x+1
>
> > def vector_func(x):
> >    a = 3*x**2+2*x+1
> >    b = 2*x+3
> >    return N.array([a,b])
>
> > print quad(scalar_func,0,2)
> > print quad(vector_func,0,2)
>
> > The output generated is:
>
> > (14.0, 1.5543122344752192e-13)
> > Traceback (most recent call last):
> >  File "quadv.py", line 15, in <module>
> >    print quad(vector_func,0,2)
> >  File "/Library/Frameworks/Python.framework/Versions/2.7/lib/
> > python2.7/site-packages/scipy/integrate/quadpack.py", line 245, in
> > quad
> >    retval =
> > _quad(func,a,b,args,full_output,epsabs,epsrel,limit,points)
> >  File "/Library/Frameworks/Python.framework/Versions/2.7/lib/
> > python2.7/site-packages/scipy/integrate/quadpack.py", line 309, in
> > _quad
> >    return
> > _quadpack._qagse(func,a,b,args,full_output,epsabs,epsrel,limit)
> > quadpack.error: Supplied function does not return a valid float.
>
> > Any ideas how to best approach this?  The other integration routines
> > (quadrature, romberg, etc) don't seem to be able to do this either.
>
> Not sure if this is how to "best" approach this, but if you really want to
> do it, you could use odeint:
>
> In [369]: from scipy.integrate import odeint
>
> In [370]: result = odeint(lambda v, t: vector_func(t), [0,0], [0,2])
>
> In [371]: result[1]
> Out[371]: array([ 14.00000002,  10.        ])
>
> Warren
>
> > _______________________________________________
> > SciPy-User mailing list
> > SciPy-U... at scipy.org
> >http://mail.scipy.org/mailman/listinfo/scipy-user
>
>
>
> _______________________________________________
> SciPy-User mailing list
> SciPy-U... at scipy.orghttp://mail.scipy.org/mailman/listinfo/scipy-user