问题描述
x,y = np.meshgrid(np.linspace(-8,8,30),np.linspace(-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
这是我的速度场，这就是它的外观velocity field
我尝试了一些在与我的问题类似的其他问题上找到的提示，但是我得到了没有任何意义的空白图表或线条。我想这部分是由于图表中间的零。 我想要的是一种从具有初始方向的初始点发送无质量粒子并观察它如何在该场中移动的方法。
谢谢！
解决方法
这取决于你真正需要什么。我可以提供我写的东西，但可能有更好的 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,y,V):
i,j = indxs(p,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[:,i+1]
W = W + a * (1-b) * V[:,j+1,i] + a * b * V[:,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,V,t_step):
k1 = VF(p,V)
k2 = VF(p + t_step*k1/2,V)
k3 = VF(p + t_step*k2/2,V)
k4 = VF(p + t_step*k3,V)
return p + t_step * (k1 + 2*k2 + 2*k3 + k4) / 6
def VF_trajectory(p,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],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,n_iterations)
# this function actually allows you to generate a plot of the flow-trajectories of
# a mesh-grid vector field
plt.streamplot(x,V[0,:])
# 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,y_bottom,y_top])
plt.show()