Using sympy.lambdify together with numpy arrays
has a severe issue when the lambdified function
is a constant. Let me show an example,
first I show what I would like to have
in any case:

x = sympy.Symbol("x")
f = sympy.lambdify(x, x**2, "numpy")

Now assume we want to evaluate f on a grid u:

u = numpy.linspace(-2,2,10)

where u.shape is (10,) and the result
f(u).shape has shape (10,) too. This is what I
expect when doing numerics. Now the example where
it goes wrong:

g = sympy.lambdify(x, sympy.diff(x**2,x,2), "numpy")

So the function is essentially a constant independent of x.
And this is the cause for the issue:

g(u) return a single float 2 and not a numpy array anymore!

This is not what I expect in numerics, although it's obvious
from a programmers point of view. (And the same issue is present
in plain python lambda expressions.)

Possible solutions:

• use numpy.vectorize:

  h = numpy.vectorize(sympy.lambdify(x, sympy.diff(x**2,x,2), "numpy"))

h(u) then gives
array([2, 2, 2, 2, 2, 2, 2, 2, 2, 2])

• Drawbacks: vectorize is SLOW, see also http://www.mail-archive.com/numpy-discussion@scipy.org/msg00587.html We should really only use it when absolutly necessary.

• use an additional pseudo-dependence:

  k = sympy.lambdify(x, sympy.diff(x**2,x,2) + 0*x, "numpy")

• This does not work because sympy simplifies the 0*x away

• A possible starting point for a patch:

def lambdify(symbols, expression):
    """Specialized lambdify routine avoiding the unused argument issue.
    A common sympy/python lambda expression 'f = lambda x: 2' has the
    problem that f(x) returns a float 2 even if x was a numpy array.
    Thus we loose all information on the shape of x.
    """
    f = sp.lambdify(symbols, expression, "numpy")

    if not all([ i in expression for i in symbols ]):
        f = np.vectorize(f)

    return f

Original issue for #5642: http://code.google.com/p/sympy/issues/detail?id=2543
Original author: https://code.google.com/u/107490137238222069432/