我想使用矢量化计算（使用 python 和 numpy）来做一些涉及贝塞尔函数的计算。

通常情况下，我需要计算（无限）不同阶的贝塞尔函数的总和。我需要为整个向量场这样做，并且想使用向量化计算。实际上，贝塞尔阶数和场的维度都可以方便地进行矢量化，但它们相交。这是该问题的最小化演示：

import numpy as np
from scipy.special import kn, iv

order = 20
size = 200

r = np.arange(-size, size, 1.0)
theta = np.arange(0, 2 * np.pi, 2 * np.pi / size)
z = np.arange(0, 1, 1 / size)

R, T, Z = np.meshgrid(r, theta, z)

m = np.arange(order)

eta = m * R
bessel_term = iv(m, eta) * kn(m, eta)

print(sum(bessel_term))

这里需要发生的是，对于 R 的每个值，需要计算具有阶向量的贝塞尔函数。可以这么说，它需要在向量中创建一个额外的维度。

相反，我收到以下错误消息，我理解，但想避免：

Traceback (most recent call last):
    eta = m * R 
ValueError: operands could not be broadcast together with shapes (20,) (200,400,200) 

我怎样才能做到这一点？一种方法是手动迭代空间维度（未矢量化）并将矢量化用于计算贝塞尔函数总和的内部循环。但对我来说，如果有一种方法可以向 numpy 表明此时需要在所涉及的向量中增加一个维度，并在计算无限和时再次折叠这个维度，那将是更可取的。

背景，数学：

实际上，我想在我的网格空间的每个点中计算一个具有三个矢量分量的磁场。在每一点我都需要计算无穷阶的贝塞尔函数的无穷和：

$$B_r = j\frac{3 \mu_0 I_1}{4 \pi r} \sum_{m= ... -5, -2, 1, 4, 7...}^{\infty} 2 m I_m(\eta_m) K_m\left( \eta_m \frac{r}{a’} \right) + \frac{2 \pi r m q}{p} \left[ I_{m-1}(\eta_m) K_{m-1}\left(\eta_m \frac{r}{a’} \right) + I_{m+1}(\eta_m) K_{m+1}\left( \eta_m \frac{r}{a'} \right) \right] \exp[jm(\theta - 2 \pi z/p)]$$

和$r$,$z$和$\theta$作为输入参数（空间坐标），$a',p, q$作为常量参数和$\eta_m=|mq|$作为缩写。