Plotting and animating movement of particle(s) inside velocity field in Python

Question

I have some code which correctly draws a vector field I want. I now want to plot and eventually animate the movement of one(or several) particles in that vector field. Now, I know I need to integrate with odeint to get the positions of a particle I place into the grid, but any tutorial or piece of code I find assumes I want to draw a parameter in relation to time... Now, I guess I could calculate for x and y individually and plot them, but there has to be a more efficient way? Do i calculate a vector product(u*v) and draw in relation to that? I guess not. Actually, I am struggling with the required parameters for odeint. So, let's say i want to draw the movement of a particle which has an initial position of X = 0.5 and Y = 0.5, in time intervals of 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()

I tried to massage the code using various data to give me at least something, but the ususal error I had was in the odeint line regarding data so I just left x0 and y0 blank as I suspect the error is there. Feel free to correct the remaining code if there is another error.

Also, how would I go about doing drawing the path of say, 5 particles, set 5 different initial conditions as a touple, a matrix, just type them out?

Thank you for your time in advance!

xdze2 · Accepted Answer · 2018-09-01 22:18:59Z

Here is part of the code with some modifications to make odeint works.

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, the initial state, have to be a vector (i.e. a 1d array), and the y argument of the dy/dt = velocity_field function too. Therefore, x and y have to be packed together, and unpacked in the function.

For the multiple solutions with different initial condition, a simple solution is to mix the solving part and the plot: (it could be better if the computation is long to separate the two, for instance by storing the solution in a list, and using another loop for the plot)

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');

which gives:

Collectives™ on Stack Overflow

Plotting and animating movement of particle(s) inside velocity field in Python

1 Answer 1

Comments

Your Answer

Hot Network Questions

Collectives™ on Stack Overflow

1 Answer 1

Comments

Your Answer

Sign up or log in

Post as a guest

Related