CODEWITHSUNDEEP
pythonscipy
SuperCiocia
SuperCiocia

Reputation: 1971

Meshgrid with lambda functions

I want to plot a function in 3D, defined this way

from scipy.integrate import quad
from scipy.special import jn

integrand = lambda x, r: np.exp(-x**2) * jn(0, r*x) * x
intensity = lambda r: quad(lambda x: integrand(x, r), 0, 5)

such that

intensity(1)

gives me the value for r = 1.

I want to plot this in 3D as function of the radius in polar coordinates, so I define a meshgrid like this:

rho = np.linspace(0, 1.25, 50)
p = np.linspace(0, 2*np.pi, 50)
R, P = np.meshgrid(rho, p)
Z = intensity(R)

with the intention of then plotting it in 3D Cartesians by doing a change of coordinates:

X, Y = R*np.cos(P), R*np.sin(P)

surf = ax.plot_surface(X, Y, Z, cmap=plt.cm.YlGnBu_r)

However when I give intensity as argument which is not a single number, but an array, it complains that

quadpack.error: Supplied function does not return a valid float.

How can I combine a lambda function with a meshgrid?

Upvotes: 1

Views: 782

Answers (2)

Arend
Arend

Reputation: 19

I know i'm late to this topic but I have been running into a similar issue. I found a workaround, it is quick but not elegant. Therefore I'm also looking for a better solution. Basically I define this function in terms of a vectorised lambda function.

    from scipy.integrate import quad
    from scipy.special import jn
    import numpy as np

    integrand = np.vectorize(lambda x, r: np.exp(-x**2) * jn(0, r*x) * x)
    intensity = np.vectorize(lambda r: quad(lambda x: integrand(x, r), 0, 5))

The rest of the code will now work as intended because it is vectorised. Therefore it can be applied to a an array.

    rho = np.linspace(0, 1.25, 50)
    p = np.linspace(0, 2*np.pi, 50)
    R, P = np.meshgrid(rho, p)
    Z = intensity(R)

Upvotes: 1

Chris
Chris

Reputation: 531

Once you've made a meshgrid R is no longer a single value but an array.

Furthermore, scipy.integrate.quad returns a tuple of the value and the estimated error.

It's maybe a question of personal taste but I like to use lambda only for anonymous functions. Otherwise, I just get confused.

My solution is probably not the fastest but it's good enough given the current parameters I think.

import numpy as np
import matplotlib.pyplot as plt
from mpl_toolkits.mplot3d import Axes3D
from scipy.integrate import quad
from scipy.special import jn

def integrand(x, r): 
    return np.exp(-x**2) * jn(0, r*x) * x

def intensity(r):
    output = np.zeros_like(r)
    for i, ri in enumerate(r.flat):
        output.flat[i] = quad(lambda x: integrand(x, ri), 0, 5)[0]
    return output

rho = np.linspace(0, 1.25, 50)
p = np.linspace(0, 2*np.pi, 50)
R, P = np.meshgrid(rho, p)
Z = intensity(R)

fig = plt.figure()
ax = fig.add_subplot(111, projection='3d')

X, Y = R*np.cos(P), R*np.sin(P)

surf = ax.plot_surface(X, Y, Z, cmap=plt.cm.YlGnBu_r)

Upvotes: 2

Related Questions