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

Question

我有一些代码可以正确绘制我想要的矢量场。我现在想要绘制并最终动画化该矢量场中一个（或多个）粒子的运动。现在，我知道我需要与 odeint 集成以获取我放置到网格中的粒子的位置，但是我发现的任何教程或代码片段都假设我想绘制与时间相关的参数......现在，我猜我可以单独计算 x 和 y 并绘制它们，但必须有更有效的方法吗？我是否计算向量积（u*v）并与之相关？我猜不会。实际上，我正在努力解决 odeint 所需的参数。因此，假设我想在 dt = 0.5 的时间间隔内绘制初始位置为 X = 0.5 和 Y = 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 个不同的初始条件设置为一个 touple、一个矩阵，只需输入它们？

提前感谢您的宝贵时间！

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

给出：