x,y = np.meshgrid(np.linspace(-8,8,30),np.linspace(-8,8,30))
q=3
w=10
freq=2
wavelenght=0.6
r=x**2+y**2
u=np.zeros((len(x),len(y)))
v=np.zeros((len(x),len(y)))
for i in range(0,len(x)):
for j in range (0,len(y)):
if (r[i,j]<=q**(3/4)):
x[i,j]=0
y[i,j]=0
if (r[i,j]>q**(3/4)):
u[i,j]=freq*wavelenght
这是我的速度场，这是速度场的样子
我在其他类似的问题上尝试了一些技巧，但我得到的是毫无意义的空白图表或线条。我想这部分是因为图中间的零。我希望有一种方法可以从一个具有初始方向的初始点发送无质量粒子，并观察它在这个场中是如何移动的。
谢谢！
这取决于您真正需要什么。我可以提供我写的东西，但可能有更好的python包或函数可以做得更好。由于我没有真正的向量场，我只是用一个网格网格向量场来测试代码，这个向量场来自于对网格顶点的多项式向量场的限制。随意使用您自己的网格矢量场
import numpy as np
import matplotlib.pyplot as plt
# left and right ends of the interval on the x-axis:
x_left = -8
x_right = 8
# bottom and top ends of the interval of the y-axis:
y_bottom = -8
y_top = 8
# number of points to be included in the mesh-grid along the x-axis
N_points_x = 200
# number of points to be included in the mesh-grid along the y-axis
N_points_y = 200
# generate the mesh-gird
x, y = np.meshgrid(np.linspace(x_left, x_right, N_points_x),np.linspace(y_bottom, y_top, N_points_y))
# a function that calculates the indices of the
# closest point from the mesh-grid to a given arbitrary point p
def indxs(p, x, y):
i = np.int(x.shape[1]*(p[0] - x[0,0])/(x[0,-1] - x[0,0]))
j = np.int(y.shape[0]*(p[1] - y[0,0])/(y[-1,0] - y[0,0]))
return i, j
# a simple quadratic interpolation of the mesh-grid vector filed
# to a quadratically interpolated vector field at a point p inside
# mesh-grid the square in which p is located
def VF(p, x, y, V):
i, j = indxs(p, x, y)
if 0 < i and i < x.shape[1]-1 and 0 < j and j < y.shape[0]-1:
a = ( p[0] - x[0, i] ) / (x[0, i+1] - x[0, i])
b = ( p[1] - y[j, 0] ) / (y[j+1, 0] - y[j, 0])
W = (1-a) * (1-b) * V[:, j, i] + (1-a) * b * V[:, j, i+1]
W = W + a * (1-b) * V[:, j+1, i] + a * b * V[:, j+1, i+1]
return W #/ np.linalg.norm(W) # you can also normalize the vector field to get only the trajecotry, without accounting for parametrization
else:
return np.array([0.0, 0.0])
# integrating the v#ector field one time step
# starting from a given point p
# uses Runge-Kutta 4 integrations, which
# allows you to sample the vector fields at four different near-by points and
# 'average' the result
def VF_flow(p, x, y, V, t_step):
k1 = VF(p, x, y, V)
k2 = VF(p + t_step*k1/2, x, y, V)
k3 = VF(p + t_step*k2/2, x, y, V)
k4 = VF(p + t_step*k3, x, y, V)
return p + t_step * (k1 + 2*k2 + 2*k3 + k4) / 6
def VF_trajectory(p, x, y, V, t_step, n_iter):
traj = np.empty( (2, n_iter), dtype = float)
traj[:, 0] = p
#m = 0
#while m < n_iter and (x[0,0] < traj[0,m] traj[1,m]):
for m in range(n_iter-1):
traj[:, m+1] = VF_flow(traj[:, m], x, y, V, t_step)
m = m+1
return traj
'''
This part is just to set up a test example, to run the code on.
When needed, generate your own vector field
'''
# initialize an empty array for the values of the vector filed at each point from the mesh-grid
V = np.empty((2, x.shape[1], y.shape[0]), dtype=float)
# a function that turns a planar polynomial vector field into a vector field on the mesh-grid
def U(x, y):
return -x + 0.3*x*y + 0.7*y**2, -0.5*y + 0.2*x*y - 0.1*x**2
# turns a planar polynomial vector field into a vector field on the mesh-grid
V[0,:,:], V[1,:,:] = U(x, y)
'''
End of the test example setup
'''
n_iterations = 700
# initial point
p0 = np.array([ -7.3, 3.1])
#p0 = np.array([-3, 5])
t_step = 0.005
trajectory = VF_trajectory(p0 , x, y, V, t_step, n_iterations)
# this function actually allows you to generate a plot of the flow-trajectories of
# a mesh-grid vector field
plt.streamplot(x, y, V[0,:,:], V[1,:,:])
# here I am plotting the trajectory generated by my functions on top of the plot of the
# flow-portrait of the mesh-grid vector field
plt.plot(trajectory[0, :], trajectory[1, :], '-r')
plt.axis('square')
plt.axis([x_left, x_right, y_bottom, y_top])
plt.show()