使用 SciPy 与可变被积函数成员进行积分并放入 meshgrid

Question

我想绘制一个矢量场，其中两个分量由数值积分确定，积分是使用要绘制矢量的 space 中的坐标计算的。

在下面的代码中，我将被积函数定义为依赖于 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)[]

Answer 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 理解如果 return调用被积函数的值是一个数组，这意味着对数组的每个元素进行积分。但是，这不是它的工作方式，但我离题了。您可以再创建一个函数来获取这些方程式中的公因数并调用它，但这取决于个人喜好。

使用 SciPy 与可变被积函数成员进行积分并放入 meshgrid

Integrate using SciPy with variable integrand members and put into meshgrid

python

numpy

matplotlib

scipy