Python 中速度场内粒子运动的绘图和动画

Question

我有一些代码可以正确绘制我想要的矢量场。我现在想绘制并最终为该矢量场中一个（或多个）粒子的运动制作动画。现在，我知道我需要与 odeint 集成以获得我放置在网格中的粒子的位置，但是我发现的任何教程或代码片段都假设我想绘制一个与时间相关的参数......现在，我猜我可以单独计算 x 和 y 并绘制它们，但必须有更有效的方法吗？我是否计算向量乘积 (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 个不同的初始条件设置为元组、矩阵，只需将它们键入即可？

提前感谢您抽出时间！

Answer 1

这是部分代码，经过一些修改后 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，初始状态，必须是向量（即一维数组），dy/dt = velocity_field 函数的 y 参数也是。因此，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');

给出：

Python 中速度场内粒子运动的绘图和动画

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

python

field

vector

scipy

odeint