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

[Warren Weckesser]

Hi Rob,

On Wed, May 25, 2011 at 2:14 PM, Rob Falck <robfalck 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-User at scipy.org
> http://mail.scipy.org/mailman/listinfo/scipy-user
>
-------------- next part --------------
An HTML attachment was scrubbed...
URL: <http://mail.scipy.org/pipermail/scipy-user/attachments/20110525/650d69c0/attachment.html>