Построение и анимация движения частиц внутри поля скоростей в Python - PullRequest
Новая
       14

Построение и анимация движения частиц внутри поля скоростей в Python

0 голосов
/ 01 сентября 2018

У меня есть код, который правильно рисует векторное поле, которое я хочу. Теперь я хочу построить и в конечном итоге анимировать движение одной (или нескольких) частиц в этом векторном поле. Теперь я знаю, что мне нужно интегрироваться с odeint, чтобы получить позиции частицы, которую я помещаю в сетку, но любой учебник или фрагмент кода, который я нахожу, предполагает, что я хочу нарисовать параметр относительно времени ... Теперь, я думаю, Я мог бы рассчитать для х и у индивидуально и построить их, но должен быть более эффективный способ? Я рассчитываю векторное произведение (u * v) и рисую по отношению к этому? Я думаю, нет. На самом деле, я борюсь с необходимыми параметрами для odeint. Итак, скажем, я хочу нарисовать движение частицы, которая имеет начальное положение X = 0,5 и Y = 0,5, во временных интервалах dt = 0,5. 

import numpy as np
import matplotlib.pyplot as plt
from matplotlib import animation
from scipy.integrate import odeint

def velocity_field(x, y, t):
    vx = -np.sin(2 * np.pi * x) * np.cos(2 * np.pi * y) - 0.1 * np.cos(2 * np.pi * t) * np.sin(
        2 * np.pi * (x - 0.25)) * np.cos(2 * np.pi * (y - 0.25))
    vy = np.cos(2 * np.pi * x) * np.sin(2 * np.pi * y) + 0.1 * np.cos(2 * np.pi * t) * np.cos(
        2 * np.pi * (x - 0.25)) * np.sin(
        2 * np.pi * (y - 0.25))
    return vx, vy

def entire_domain():
    xrange = (0, 1)
    yrange = (0, 1)
    mesh_sz_x = 50
    mesh_sz_y = 50
    dx = (xrange[1] - xrange[0]) / (mesh_sz_x - 1)
    dy = (yrange[1] - yrange[0]) / (mesh_sz_y - 1)

    x_mat, y_mat = np.mgrid[xrange[0]:xrange[1]:dx, yrange[0]:yrange[1]:dy]

    x_dot, y_dot = velocity_field(x=x_mat, y=y_mat, t=0)

    speed = np.sqrt(x_dot ** 2 + y_dot ** 2)

    u_n = x_dot / speed
    v_n = y_dot / speed

    plt.contourf(x_mat, y_mat, speed, 12, cmap=plt.get_cmap('viridis'),
                 interp='bicubic')

    plt.quiver(x_mat, y_mat, u_n, v_n  # data
               , color='black'
               , headlength=7
               , pivot='mid'
               ,
               )  # length of the arrows

    #This part is wrong
    '''
    x0 = ?????
    y0 = ?????
    t = np.arange(0, 100, 0.05)

    X = odeint(velocity_field, x0, y0, t)
    print(X)
    '''
    plt.show()

if __name__ == '__main__':
    entire_domain()

Я пытался помассировать код, используя различные данные, чтобы дать мне хоть что-то, но у меня была обычная ошибка в строке odeint относительно данных, поэтому я просто оставил пустыми x0 и y0, так как подозреваю, что ошибка есть. Не стесняйтесь исправлять оставшийся код, если есть другая ошибка.

Кроме того, как мне рисовать путь, скажем, 5 частиц, задать 5 различных начальных условий в виде матрицы, просто наберите их?

Спасибо, что уделили время заранее!

1 Ответ

0 голосов
/ 02 сентября 2018

Вот часть кода с некоторыми изменениями, чтобы odeint работал. 

import numpy as np
import matplotlib.pyplot as plt
from scipy.integrate import odeint

def velocity_field(xy, t):
    x, y = xy
    vx = -np.sin(2*np.pi * x) * np.cos(2*np.pi * y) \
        - 0.1*np.cos(2 * np.pi * t) * np.sin(2*np.pi*(x - 0.25)) * np.cos(2*np.pi*(y - 0.25))
    vy =  np.cos(2*np.pi * x) * np.sin(2*np.pi * y) \
       + 0.1*np.cos(2*np.pi * t) * np.cos(2*np.pi*(x - 0.25)) * np.sin(2*np.pi*(y - 0.25))
    return (vx, vy)


xy0 = (0, 0)
t_span = np.arange(0, 100, 0.05)

sol = odeint(velocity_field, xy0, t_span)

plt.plot(sol[:, 0], sol[:, 1]);
plt.axis('equal'); plt.xlabel('x'); plt.ylabel('y');

y0, начальное состояние, должно быть вектором (т.е. массивом 1d), а также аргументом y функции dy/dt = velocity_field. Следовательно, x и y должны быть упакованы вместе и распакованы в функции.

Для нескольких решений с различным начальным условием, простое решение состоит в том, чтобы смешать решающую часть и график: (может быть лучше, если вычисление будет длинным, чтобы разделить их, например, путем сохранения решения в списке, и используя другой цикл для сюжета) 

initial_conditons = [(1, .4), (1, 0), (4, .7)]
t_span = np.arange(0, 100, 0.05)

plt.figure(figsize=(6, 6))
for xy0 in initial_conditons:
    sol = odeint(velocity_field, xy0, t_span)
    plt.plot(sol[:, 0], sol[:, 1]);

plt.axis('equal'); plt.xlabel('x'); plt.ylabel('y');

, что дает:

nice example

...