Исследование траектории частицы в потенциальном поле¶
Задаем потенциал¶
Задаем 3-х мерную синусоидальную решетку в параболической потенциальной яме.
import numpy as np def V(x, y, z): return np.cos(10*x) + np.cos(10*y) + np.cos(10*z) + 2*(x**2 + y**2 + z**2)
Создаем прямоугольную 3-х мерную сетку для построения потенциала. Сетка содержит 100 интервалов и расположена в кубе по X, Y, Z от –2 до 2. Вычисляем потенциал для точек этой сетки.
X, Y, Z = np.mgrid[-2:2:100j, -2:2:100j, -2:2:100j] VV = V(X,Y,Z)
Используем модуль mlab (с.м. mlab: Python scripting for 3D plotting) для интерактивной визуализации данных. Визуализируем потенциал изоповерхностями:
from mayavi import mlab mlab.contour3d(X, Y, Z, VV) mlab.show()
После выполнения построение окно Mayavi обладает полнофункциональным интерактивным интерфейсом и позволяет поворачивать объекты, менять вид представления данных, параметры окна, цветовые схемы, плотность изолиний и т.д. Поменяем параметры изоповерхностей. Это можно сделать двойным нажатием на IsoSurface в дереве структур Mayavi (если вы запускаете скрипт из питона, нужно кликнуть на значке Mayavi в верхнем углу, чтобы появилось меню). Это открывает диалог, изменяющий параметры контуров изоповерхности. Хорошая картинка получается после выключения автоматического создания контуров и выбора –0,5 для начального контура и 15 для последнего. Можно поменять Colors and legends в дереве объектов и выбрать другую цветовую схему - LUT (Look Up Table).
Для лучшего представления потенциала, хотелось бы увидеть больше контуров, но новые контуры будут закрывать внутренние. Одно из решений использование плоскости сечения. Нажмите правой кнопкой на элементе IsoSurface и добавьте элемент ScalarCutPlane в одном из подменю раздела “Add module”. Можно перемещать плоскость сечения кликая на нее мышкой и удерживая кнопку нажатой.
Можно манипулировать с объектами Mayavi через открытое окно питон. Добавим легенду в наше окно.
mlab.colorbar(title='Potential', orientation='vertical')
Можно сделать то-же самое через меню LUT в интерактивном режиме.
Задаем силы¶
Для изучения движения частицы в заданном потенциале, необходимо получить соответствующие силы, как градиент потенциала. Создаем простейшую функцию градиента:
def gradient(f, x, y, z, d=0.01): fx = f(x+d, y, z) fx_ = f(x-d, y, z) fy = f(x, y+d, z) fy_ = f(x, y-d, z) fz = f(x, y, z+d) fz_ = f(x, y, z-d) return (fx-fx_)/(2*d), (fy-fy_)/(2*d), (fz-fz_)/(2*d)
Для проверки работы функции градиента, визуализируем векторное поле. Чтобы не перегружать картинку множеством векторов построим градиенты только в плоскости X=50, и каждой третьей точке по сетке.
Vx, Vy, Vz = gradient(V, X[50, ::3, ::3], Y[50, ::3, ::3], Z[50, ::3, ::3]) mlab.quiver3d(X[50, ::3, ::3], Y[50, ::3, ::3], Z[50, ::3, ::3], Vx, Vy, Vz, scale_factor=-0.2, color=(1, 1, 1)) mlab.show()
Траектории¶
Теперь можно использовать scipy для интегрирования траекторий. Сначала зададим динамический поток, функцию, которая возвращает изменение каждой переменной за один шаг по времени. Эта функция используется в численном решателе ОДУ и задает динамику системы. В нашем случае изменение координат определятся соответствующими скоростями, а изменение скоростей – силами. Сила состоит из градиента потенциала (V) с обратным знаком и силы трения, пропорциональной скорости.
def flow(r, t): """ The dynamical flow of the system """ x, y, z, vx, vy, vz = r fx, fy, fz = gradient(V, x-.2*np.sin(6*t), y-.2*np.sin(6*t+1), z-.2*np.sin(6*t+2)) return np.array((vx, vy, vz, -fx - 0.3*vx, -fy - 0.3*vy, -fz - 0.3*vz))
Создаем траекторию. Начальные данные - нулевые.
from scipy.integrate import odeint # Initial conditions R0 = (0, 0, 0, 0, 0, 0) # Times at which we want the integrator to return the positions: t = np.linspace(0, 50, 500) R = odeint(flow, R0, t)
Рисуем траектории после удаления секущей плоскости и векторного поля, щелкая на соответствующие элементы и выбирая «удалить». Рисуем траекторию через параметрический график со своей цветовой схемой.
x, y, z, vx, vy, vz = R.T trajectory = mlab.plot3d(x, y, z, t, colormap='hot', tube_radius=None) mlab.colorbar(trajectory, title='Time', orientation='vertical') mlab.show()