Algorithm to calculate trajectories from a vector field

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?

Answer 1

[Jon]

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.

[Ashfaq Ahmed]

Hi, thank you for the comment it helps!

And yes, my vector field is not a function, rather it's a data set (matrix size 71X301) and I will have to approach with it while including some interpolation. Could you kindly tell me how can I approach for my case?

[Jon]

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.

[Ashfaq Ahmed]

Hi Jon, I could formulated the interpolated grid!

clear all; close all; clc;

[X,Y] = meshgrid(0:6,0:6);
U = 0.25*X;
V = 0.5*Y;
[Xq,Yq] = meshgrid(-6:0.25:6);
Uq = interp2(X,Y,U,Xq,Yq);
Vq = interp2(X,Y,V,Xq,Yq);
figure(1),
subplot(1,2,1)
quiver(X,Y,U,V,0)
title('Plotting the Vectors on a Coarse Grid','fontweight','bold', 'fontsize',20);
grid on
subplot(1,2,2)
quiver(Xq,Yq,Uq,Vq);
grid minor
title('Linear Interpolation Using Finer Grid','fontweight','bold', 'fontsize',20);

Since the interpolation is done, now I have more values for u and v in more points. Could you kindly explain me this sentence you wrote? "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. "

I think I am very close to the answer!

[Jon]

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

[Ashfaq Ahmed]

Hi, thank you so, sooo much! Now I get to see the trajectories at any point I want on this grid.

[Jon]

Yes, but as I said, there still seem to be some numerical solution issues, as you would hope that if you made the final time longer you would see the trajectory run all the way to the edge of the grid. Instead you get some values for the trajectory that look reasonable and then all NaN (not a number). Making the final time longer actually seems to make things worse (maybe somehow it increases the step size). Hopefully you can figure out how to resolve this problem, or as I said maybe it isn't an issue with your real data.

I did try using different solvers, ode15s ode23, etc but that didn't seem to help. I also tried adjusting the error tolerances and the value that the interp2 extrapolates to when points go outside of the defined grid. I didn't find any real improvement with any of these.

[Ashfaq Ahmed]

Hi @Jon, I think I have actually solved the issue about what you were mentioning. You are right - the NaN value appeared because of the characteristics of the vector field itself, since it was diverging to infinity as the x and y increase I tried it with different vector field and the code worked like a charm! For example,

[X,Y] = meshgrid(-3*pi:pi/4:3*pi,-3*pi:pi/4:3*pi);
U = sin(Y);
V = cos(X);
[Xq,Yq] = meshgrid(-3*pi:pi/4:3*pi,0.5,-3*pi:pi/4:3*pi);

And for t = 3 sec, the trajectory line was -

When I increased the t, the trajectory forms a circle (which makes sense.)

This whole code is giving me a fascinating idea! Suppose, instead of 1 particle, if I want to place 5 particles at the same time in the different location (of course, within the grid), what changes should I bring to the code? For instance, I want 5 particles initally starting from (x,y) = (1,3), (4,1), (3,5), (2,8), and (6,2)

Could you please give me an idea on this? Thank you so much for your extraordinary support so far! :D

[Jon]

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]

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]

Can you give me an idea on which function should I use to make an animation of the particle trajectories along with the time?

[Jon]

While looking for how to animate, I stumbled across this https://www.mathworks.com/help/matlab/ref/streamline.html

Looks like there is a built in function to plot the stream lines! So maybe you don't need to write all of the code yourself after all. Lesson is it is probably better to spend more time searching MATLAB documentation and less time jumping into coding. Oh well, I guess we learned something about computing streamlines :)

Regarding the animation, please see this https://www.mathworks.com/help/matlab/animation-1.html

In general for searching for MATLAB functions and how to do things I find it most effective to just use Google search with Matlab included in the key words, so for example: search the terms matlab streamline. This often seems to turn up better results then searching within the MATLAB help using their searchbar.

[Ashfaq Ahmed]

Thanks Jon. I don't know how can I properly expresss how grateful I am for your constant support!

[Jon]

You're very welcome. Good luck with your project

[Ashfaq Ahmed]

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]

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.