使用 SciPy 与可变被积函数成员集成并放入 meshgrid答案

【问题标题】：Integrate using SciPy with variable integrand members and put into meshgrid使用 SciPy 与可变被积函数成员集成并放入 meshgrid
【发布时间】：2015-02-15 23:50:02
【问题描述】：

我正在寻找一个矢量场，其中两个分量是通过数值积分确定的，并且积分是使用要绘制矢量的空间坐标计算的。

在下面的代码中，我定义了依赖于 r 和 z 的被积函数，然后在循环中定义它们的值并在每个 (r,z) 处计算两个积分。

我有两个问题：

1) 有没有更 Pythonic 的方法来使用变量 r 和 z 评估这些积分？

2) 存储积分值及其坐标以从meshgrids 生成quiver 图的最佳方法是什么？

import scipy.integrate as spi
import scipy.special as sps
import numpy as np

R = 2
V = 10

def Efield_r_integrand(k):
    return np.exp(-k*z)*k*R*V*sps.jv(1,k*r)*sps.jv(1,k*R)

def Efield_z_integrand(k):
    return np.exp(-k*z)*k*R*V*sps.jv(0,k*r)*sps.jv(1,k*R)

x_max = 3.0
z_max = 3.0
n_pts = 20

for i in xrange(n_pts):
    r = float(i)/float(n_pts)*r_max
    for j in xrange(n_pts):
        z = float(j)/float(n_pts)*z_max
        current_Efield_r = spi.quad(Efield_r_integrand,0,np.inf)[]
        current_Efield_z = spi.quad(Efield_z_integrand,0,np.inf)[]

【问题讨论】：

标签： python numpy matplotlib scipy

【解决方案1】：

你的代码已经相当pythonic了，除了结构

for i in xrange(n_pts):
    r = float(i)/float(n_pts)*r_max

这让人想起其他一些编程语言。在 Python 中，写起来会更 Pythonic

for r in np.arange(0, rmax, rmax/n_pts):

因为您不需要中间变量i。

话虽如此，评估在网格上定义的函数的积分是我不再编写双 for 循环的事情，而是让便利函数 np.vectorize 来处理：

import scipy.integrate as spi
import scipy.special as sps
import numpy as np
import matplotlib.pyplot as plt

def Efield_r_integrand(k, z, r, R, V):
    return np.exp(-k*z)*k*R*V*sps.jv(1,k*r)*sps.jv(1,k*R)

def Efield_z_integrand(k, z, r, R, V):
    return np.exp(-k*z)*k*R*V*sps.jv(0,k*r)*sps.jv(1,k*R)

x_max, z_max, n_pts = 3.0, 3.0, 20
R, V = 2, 10

Z, X = np.mgrid[0:z_max:n_pts*1j, 0:x_max:n_pts*1j] # note that this creates a grid with the boundaries (0 and 3.0) included! 

def integrate_on_grid(func, lo, hi, *args):
    """Returns a callable that can be evaluated on a grid."""
    return np.vectorize(lambda n,m: spi.quad(func, lo, hi, (n,m)+args)[0])

Efield_r, Efield_z = [integrate_on_grid(func, 0, np.inf, R, V)(Z,X)
    for func in (Efield_r_integrand, Efield_z_integrand)]

plt.quiver(X, Z, Efield_r, Efield_z)

最后一行展示了如何轻松使用迄今为止获得的结果来生成quiver 图。仍有一些重复可以删除：Efield_r_integrand 最多与Efield_z_integrand 相同，理想情况下可以通过让quad 了解如果调用被积函数的返回值是一个数组来解决这个问题，这意味着对数组的每个元素进行积分。然而，这不是它的工作方式，但我离题了。您可以再创建一个函数来获取这些方程中的公因子并调用它，但这取决于个人喜好。

【讨论】：