Algorithm to calculate trajectories from a vector field

Question

Ashfaq Ahmed 2022-1-10

0
链接

此问题的直接链接

https://ww2.mathworks.cn/matlabcentral/answers/1626260-algorithm-to-calculate-trajectories-from-a-vector-field

评论： Jon 2022-1-26

Dear altruists,

Suppose, I have a two-dimensional vector field, i.e., for each point (x, y) I have a vector (u, v), whereas u and v are functions of x and y.

This vector field canonically defines a set of trajectories, i.e. a set of paths a particle would take if it follows along the vector field. In the following image, the vector field is depicted in blue, and there are four trajectories (which are my expected outcome), depicted in dark red:

I am looking for an algorithm which will give me a trajectory of a virtual particle I imagine from anypoint of that vector field. The trajectories must satisfy some kind of minimum denseness in the plane (for every point in the plane we must have a 'nearby' trajectory), or some other condition to get a reasonable set of trajectories.

I could not find anything useful on Google on this, so I posted it here because I always get responses from this community :)

Before I start devising such an algorithm by myself: Are there any known algorithms for this problem? What is their name, for which keywords do I have to search?

0 个评论
显示 -2更早的评论隐藏 -2更早的评论

请先登录，再进行评论。

请先登录，再回答此问题。

Answer 1

Jon 2022-1-10

1
链接

此回答的直接链接

https://ww2.mathworks.cn/matlabcentral/answers/1626260-algorithm-to-calculate-trajectories-from-a-vector-field#answer_872130

编辑：Jon 2022-1-10

The trajectories you are seeking are solutions to a set of ordinary differential equations with initial conditions given by the starting point of the particle. You can solve ode's using MATLAB ode solvers see https://www.mathworks.com/help/matlab/ordinary-differential-equations.html

These functions will solve a system of ode's dx/dt = f(x). Since these functions only a require that you create a function that gives the derivatives (dx/dt) as a function of the system state (x). The same fuction that you are using to generate your vector fields can be used to provide the needed derivatives at a given state values. If your vector fields are from experimental data, you could probably use the data you have and some interpolation to calculate the derivatives needed for the ode solver.

15 个评论
显示 13更早的评论隐藏 13更早的评论

Jon 2022-1-10

Lets consider the velocities at each point as having components u = dx/dt and v = dy/dt, these are the two values that you must supply to your ODE solver for any given value of x and y.

Suppose from your experiments you have a two dimensional array U of horizontal velocity components u where the rows are y locations and the columns are x locations. Similarly you have an array V of the vertical velocity components. You can then use MATLAB's interp2 to interpolate to get u and v values at query points xq, yq, within the array. You will also probably need to make use of meshgrid, to supply suitable arrays X and Y of x and y locations to the interp2 function.

To get details on interp2 and meshgrid, type respectively doc interp2 and doc meshgrid, on your command line, or otherwise search your MATLAB help.

Once you work out how to interpolate a value of u and v given query points xq,yq you will need to put this functionality into your own MATLAB function that outputs a two element vector xdot, whose first element is u and second element is v. You then supply this function to the MATLAB ODE solver that you choose, e.g. ode45

Hope this gets you started.

Jon 2022-1-11

编辑：Jon 2022-1-11

在 MATLAB Online 中打开

odefun.m

I have attached a function that provides the derivatives (velocity components), and the script below, modified from your example, will use this function to overlay the velocity of a particle.

% define example velocity field
% in actual application this would be from experimental data
[X,Y] = meshgrid(0:6,0:6);
U = 0.25*X;
V = 0.5*Y;
[Xq,Yq] = meshgrid(-6:0.25:6);
% display the velocity field
figure(1),
quiver(X,Y,U,V,0)
title('Vector Field with Particle Trajectory ');
grid on
hold on
% assign initial position for particle to be tracked
p0(1) = 1;
p0(2) = 3;
% set final time for simulation
tF = 2;
% use ode solver to obtain trajectories
% need to pass the function that computes derivatives as a 
% function handle using anonymous function in order to pass
% additional parameters U,V,X,Y
extrap = max(max(U(:),max(V(:)))); % velocity value to use when outside of grid
f = @(t,p) odefun(t,p,U,V,X,Y,extrap);
tspan = [0,tF];
% assign options for solver
options = odeset('RelTol',1e-10,'AbsTol',1e-10);
[t,p] = ode45(f,tspan,p0);
% overlay the particle track
plot(p(:,1),p(:,2))
hold off

It runs, but there still seem to be some issues with the ode45 solver, as the trajectory values start getting to be computed as NaN's before the final time is reached. So you get an initial part of the trajectory, but then no further points. I think that the underlying problem is that your vector field represents an unstable system (the velocities get bigger the futher you are from the origin). I think this gives the solver some problems. Maybe if you're actual vector field does not behave this way you will not have problems.

Maybe another MATLAB answers reader has some ideas here ?

Anyhow this should at least get you started

Jon 2022-1-13

在 MATLAB Online 中打开

Glad it is working well on a more realistic velocity field. The simplest way is to do multiple initial particle postions is just to make a loop around the ode solver feeding it new initial conditions. If you only want to see the trajectories displayed, and don't need to store the individual tracks you can just put the plotting inside the loop. Like this:

% display the velocity field
figure(1),
quiver(X,Y,U,V,0)
title('Vector Field with Particle Trajectory ');
grid on
hold on
% assign initial positions for particle to be tracked
% first rows are x and y coordinates, columns are different starting
% positions
p0 = [1 4 2;3 5 2];
% set final time for simulation
tF = 2;
% use ode solver to obtain trajectories
% need to pass the function that computes derivatives as a
% function handle using anonymous function in order to pass
% additional parameters U,V,X,Y
extrap = max(max(U(:),max(V(:)))); % velocity value to use when outside of grid
f = @(t,p) odefun(t,p,U,V,X,Y,extrap);
tspan = [0,tF];
% assign options for solver
options = odeset('RelTol',1e-10,'AbsTol',1e-10);
% loop through initial particle positions
numInit = size(p0,2); % number of initial positions is number of columns in p0
for k = 1:numInit
    [t,p] = ode45(f,tspan,p0(:,k));
    % overlay the particle track
    plot(p(:,1),p(:,2))
end
hold off

It is also possible to solve for multiple initial conditions all at once https://www.mathworks.com/help/matlab/math/solve-system-of-odes-with-multiple-initial-conditions.html which might be more efficient. The programming gets more complicated though, so unless you are having peformance/speed problems it probably isn't worth it.

If you need to save the individual tracks, you could put them into a cell array. Might be able to put them into a 3-d array with a page for each particle, or two 2-d array one for each velocity component, but I'm not sure if you can ensure that the number of time steps for each particle will be the same. Could make it a little more complicated to deal with.

Ashfaq Ahmed 2022-1-13

编辑：Ashfaq Ahmed 2022-1-13

在 MATLAB Online 中打开

Wonderful!! I was actually trying it this morning and I did somethinig very close to what you proposed (except, mine one was very basic).

figure(1),
quiver(X,Y,U,V,0,'autoscale','on'); 
title('Vector Field with Particle Trajectory starting from point (x,y)', 'fontsize',20)
axis equal; grid on; hold on
%% Assign the initial position for the particle to be tracked
p0= [1,3; -5,2; 8,1; -4,4; 6,3];
%% Set simulation time for the simulation
tF = 5;
%% Solving with ODE
extrap = max(max(U(:),max(V(:)))); % velocity value to use when outside of grid
f = @(t,p) odefun(t,p,U,V,X,Y,extrap);
tspan = [0,tF];
for i=1:5
options = odeset('RelTol',1e-10,'AbsTol',1e-10);
[t,p] = ode45(f,tspan,p0(i,:));
% overlay the particle track
plot(p(:,1),p(:,2),'LineWidth',4,'MarkerSize',10)
hold on
   
end
%% Function
function dpdt = odefun(~,p,U,V,X,Y,extrap)
xq = p(1);
yq = p(2);
% interpolate to get velocity components at query points
u = interp2(X,Y,U,xq,yq,"cubic");
v = interp2(X,Y,V,xq,yq,"cubic");
% assign derivative vector
dpdt = [u;v];
end

MATLAB is awesome!!

Ashfaq Ahmed 2022-1-26

编辑：Ashfaq Ahmed 2022-1-26

在 MATLAB Online 中打开

Hi @Jon, I have submitted the project with some changes and it went really good! I have, however, one confusion that I could never solve. Could you please tell me why did this line always make 41x2 matrix for p?

[t,p] = ode45(f,tspan,p0);

This was the comment you wrote on 11 January -

https://www.mathworks.com/matlabcentral/answers/1626260-algorithm-to-calculate-trajectories-from-a-vector-field#comment_1930410

Jon 2022-1-26

The rows in p are the time steps, the columns are the x and y positions. So I guess you understand why it had two columns. Regarding the number of rows, being 41, I this corresponds to the number of time steps the algorithm takes to go from 0 to tspan. The ode45 solver is a variable step size solver, so you don't have direct control of how many steps it takes. It depends upon the vector field and also the options for the error tolerances.

请先登录，再进行评论。

Algorithm to calculate trajectories from a vector field

0 个评论
显示 -2更早的评论隐藏 -2更早的评论

采纳的回答

15 个评论
显示 13更早的评论隐藏 13更早的评论

更多回答（0 个）

另请参阅

类别

标签

Community Treasure Hunt

Algorithm to calculate trajectories from a vector field

0 个评论 显示 -2更早的评论隐藏 -2更早的评论

采纳的回答

15 个评论 显示 13更早的评论隐藏 13更早的评论

更多回答（0 个）

另请参阅

类别

标签

Community Treasure Hunt

0 个评论
显示 -2更早的评论隐藏 -2更早的评论

15 个评论
显示 13更早的评论隐藏 13更早的评论