python - using SciPy to integrate a function that returns a matrix or array

Question

I have a symbolic array that can be expressed as:

from sympy import lambdify, Matrix

g_sympy = Matrix([[   x,  2*x,  3*x,  4*x,  5*x,  6*x,  7*x,  8*x,   9*x,  10*x],
                  [x**2, x**3, x**4, x**5, x**6, x**7, x**8, x**9, x**10, x**11]])

g = lambdify( (x), g_sympy )

So that for each x I get a different matrix:

g(1.) # matrix([[  1.,   2.,   3.,   4.,   5.,   6.,   7.,   8.,   9.,  10.],
      #         [  1.,   1.,   1.,   1.,   1.,   1.,   1.,   1.,   1.,   1.]])
g(2.) # matrix([[  2.00e+00,   4.00e+00,   6.00e+00,   8.00e+00,   1.00e+01, 1.20e+01,   1.40e+01,   1.60e+01,   1.80e+01,   2.00e+01],
      #         [  4.00e+00,   8.00e+00,   1.60e+01,   3.20e+01,   6.40e+01, 1.28e+02,   2.56e+02,   5.12e+02,   1.02e+03,   2.05e+03]])

and so on...


I need to numerically integrate g over x, say from 0. to 100. (in the real case the integral does not have an exact solution) and in my current approach I have to lambdify each element in g and integrate it individually. I am using quad to do an element-wise integration like:

ans = np.zeros( g_sympy.shape )
for (i,j), func_sympy in ndenumerate(g_sympy):
    func = lambdify( (x), func_sympy)
    ans[i,j] = quad( func, 0., 100. )

There are two problems here: 1) lambdify used many times and 2) for loop; and I believe the first one is the bottleneck, because the g_sympy matrix has at most 10000 terms (which is not a big deal to a for loop).

As shown above lambdify allows the evaluation of the whole matrix, so I thought: "Is there a way to integrate the whole matrix?"

scipy.integrate.quadrature has a parameter vec_func which gave me hope. I was expecting something like:

g_int = quadrature( g, x1, x2 )

to get the fully integrated matrix, but it gives the ValueError: matrix must be 2-dimensional

---

EDIT: What I am trying to do can apparently be done in Matlab using quadv and has already been discussed for SciPy

The real case has been made available here.

To run it you will need:

• numpy
• scipy
• matplotlib
• sympy

Just run: "python curved_beam_mrs.py".

You will see that the procedure is already slow, mainly because of the integration, indicated by the TODO in file curved_beam.py.

It will go much slower if you remove the comment indicated after the TODO in file curved_beam_mrs.py.

The matrix of functions which is integrated is showed in the print.txt file.

Thank you!

Answer 1

The first argument to either quad or quadrature must be a callable. The vec_func argument of the quadrature refers to whether the argument of this callable is a (possibly multidimensional) vector. Technically, you can vectorize the quad itself:

>>> from math import sin, cos, pi
>>> from scipy.integrate import quad
>>> from numpy import vectorize
>>> a = [sin, cos]
>>> vectorize(quad)(a, 0, pi)
(array([  2.00000000e+00,   4.92255263e-17]), array([  2.22044605e-14,   2.21022394e-14]))

But that's just equivalent to explicit looping over the elements of a. Specifically, it'll not give you any performance gains, if that's what you're after. So, all in all, the question is why and what exactly you are trying to achieve here.