I want to use vectorized computation (with python and numpy) to do some calculation involving bessel functions.

As is often the case, I need to calculate the sums of bessel functions of (infinit) differnt order. I need to do that for a whole vector field, and would like to use vectorized computation. Actually, both, the bessel orders and the dimensions of the field would be handy to do vectorized, but they intersect. Here is a minimal, dummed down demo of the issue:

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)


What would need to happen here is, that of each value of R the bessel functions with the order-vector would need to be calculated. It would need to create an additional dimension in the vector, so to speak.

Instead I get the following error message which I understand, but want to avoid:

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

How can I do that? One way would be to iterate through the space dimensions manually (not vectorized) and use the vectorisation for the inner loop where the sums over the bessel functions are calculated. But to me it would be preferable if there was a way to indicate to numpy that at this point an additional dimention in the involved vectors is needed and collapse this dimension again when the infinite sum is calculated.

Background, the math:

In reality I want to calulate a magnetic field with three vector components in every point in space of my meshgrid. In every point I need to calculate an infinit sum of besselfunctions of infinit orders:

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

with $r$, $z$ and $\theta$ as input parameters (space coordinates), $a',p, q$ as constant arguments and $\eta_m=|mq|$ as an abreviation.

  • $\begingroup$ It might help if you can write out what you are aiming to do using Mathjax, as you seem to know what you want to do mathematically, but just aren't sure how to express it with Numpy. As is, it's not clear to me what you are trying to do: m is a 20-element array and R is a (200,400,200)-element array, so they can't be multiplied because the lengths don't coincide along any dimension. $\endgroup$ – Tyberius Dec 9 '20 at 22:24
  • $\begingroup$ @Tyberius added math background. $\endgroup$ – Andreas Schuldei Dec 10 '20 at 8:22

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service, privacy policy and cookie policy

Browse other questions tagged or ask your own question.