Блог

Мне нужно построить графики контуров и колебаний скалярных и векторных полей, определенных на неравномерной сетке в координатах (r, theta).

В качестве минимального примера проблемы, с которой я столкнулся, рассмотрим график контура потоковой функции для магнитного диполя, контуры такой функции являются линиями тока соответствующего векторного поля (в данном случае магнитного поля).

Приведенный ниже код берет неравномерную сетку в координатах (r, theta), отображает ее на декартову плоскость и строит контурный график функции потока.

 import numpy as np
import matplotlib.pyplot as plt

r = np.logspace(0,1,200)
theta = np.linspace(0,np.pi/2,100)

N_r = len(r)
N_theta = len(theta)

# Polar to cartesian coordinates
theta_matrix, r_matrix = np.meshgrid(theta, r)
x = r_matrix * np.cos(theta_matrix)
y = r_matrix * np.sin(theta_matrix)

m = 5
psi = np.zeros((N_r, N_theta))

# Stream function for a magnetic dipole
psi = m * np.sin(theta_matrix)**2 / r_matrix

contour_levels = m * np.sin(np.linspace(0, np.pi/2,40))**2.

fig, ax = plt.subplots()
# ax.plot(x,y,'b.')  # plot grid points
ax.set_aspect('equal')
ax.contour(x, y, psi, 100, colors='black',levels=contour_levels)
plt.show()
  

Однако по какой-то причине график, который я получаю, выглядит неправильно:

Если я поменяю местами x и y в вызове функции contour, я получу желаемый результат:

То же самое происходит, когда я пытаюсь создать график колчана векторного поля, определенного на той же сетке и отображенного на плоскость xy, за исключением того, что замена x и y в вызове функции больше не работает.

Похоже, я где-то допустил глупую ошибку, но я не могу понять, что это такое.

If psi = m * np.sin(theta_matrix)**2 / r_matrix then psi увеличивается при переходе от 0 к pi / 2, а psi уменьшается при увеличении r.

Таким образом, контурная линия для psi должна увеличиваться в r по мере увеличения theta. В результате получается кривая, которая идет против часовой стрелки, поскольку она расходится от центра. Это согласуется с первым опубликованным вами графиком и результатом, возвращенным первой версией вашего кода с

 ax.contour(x, y, psi, 100, colors='black',levels=contour_levels)
  

Альтернативный способ подтвердить достоверность результата — посмотреть на график поверхности psi :

 import numpy as np
import matplotlib.pyplot as plt
import mpl_toolkits.mplot3d.axes3d as axes3d

r = np.logspace(0,1,200)
theta = np.linspace(0,np.pi/2,100)

N_r = len(r)
N_theta = len(theta)

# Polar to cartesian coordinates
theta_matrix, r_matrix = np.meshgrid(theta, r)
x = r_matrix * np.cos(theta_matrix)
y = r_matrix * np.sin(theta_matrix)

m = 5

# Stream function for a magnetic dipole
psi = m * np.sin(theta_matrix)**2 / r_matrix

contour_levels = m * np.sin(np.linspace(0, np.pi/2,40))**2.

fig = plt.figure()
ax = fig.add_subplot(1, 1, 1, projection='3d')
ax.set_aspect('equal')

ax.plot_surface(x, y, psi, rstride=8, cstride=8, alpha=0.3)
ax.contour(x, y, psi, colors='black',levels=contour_levels)
plt.show()