具有 lambda 函数的网格网格

Question

我想在 3D 中绘制一个这样定义的函数

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)

这样

intensity(1) 

给我 r = 1 的值。

---

我想在 3D 中将其绘制为极坐标中半径的函数，因此我定义了一个这样的网格：

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

然后通过更改坐标以 3D 笛卡尔坐标绘制它：

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

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

但是，当我将intensity 作为不是单个数字而是数组的参数时，它会抱怨

quadpack.error: 提供的函数没有返回有效的浮点数。

如何将 lambda 函数与网格网格相结合？

Answer 1

一旦你创建了一个网格网格，R 就不再是一个单一的值，而是一个数组。

此外，scipy.integrate.quad 返回值和估计误差的元组。

这可能是个人品味的问题，但我喜欢仅将 lambda 用于匿名函数。否则，我只会感到困惑。

我的解决方案可能不是最快的，但考虑到我认为的当前参数，它已经足够了。

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)

Answer 2

我知道我迟到了这个话题，但我遇到了类似的问题。 我找到了一种解决方法，它很快但不优雅。因此，我也在寻找更好的解决方案。基本上我是根据向量化的 lambda 函数来定义这个函数的。

    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))

其余代码现在将按预期工作，因为它是矢量化的。因此它可以应用于数组。

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