double gyre - compute absolute and relative dispersion

Hello everyone! I wrote a routine to integrate numerical particle trajectories for a double gyre flow and now I am supposed to compute the absolute and relative dispersion based on the computed trajectories. Can anyobdy help me in building the for cycle? I believe the cycle has three indices ( i = length of the x-domain, j = length of the y-domain and, k = timesteps) but i have no idea how to build it. you can find attached the scripts that I am using. If you run" integrateDGtemp.m" you can quickly understand what it is about.

doublegyreVEC.m contains the double gyre vector field

rk4singlestep.m contains the Runge Kutta 4 scheme for solving the ODE

integrateDGtemp.m is the script used to plot the initial particle positions and then the double gyre field.

Also, attached here you can see the formulas that I am supposed to use for computing absolute and relative dispersion.

Thank you a lot if you can help.

Answer 1

[Jiri Hajek]

Hi, seems like you are almost there... But you have to make up your mind on the meaning of your dispersion formulae. In statistics, dispersion is not a concrete function, but rather a quality of a data set, which may be expressed by any of several statistical parameters. (I guess that the i, j, k indices all describe the same set of M particles and consequently the expressions in brackes have the meaning of eulerian distances between pairs of particles.)

To evaluate your formulae, you'll need just a few basic MATLAB functions:

• diff to compute the differences of positions (the bracketed terms)
• .* dot product of same-size vectors
• sum of a vector

Using these, you will be able to calculate single element of the A2 and R2 matricex. Then, to complete these matrices at a time instant, you would need two more for loops.

R2

[Tiziano Bagnasco]

Thank you so much for your answer... However, could you pleace help me in setting up a for cycle? I am very confused.. Mostly because I have not saved the x and y particle positions for every timestep in a proper format. Thanks

[Jiri Hajek]

Yes, your script would benefit from some enhancement and clear data structure. But you see, the scope and detail of your code is entirely up to you; here we can just discuss issues that cause confusion. To your question:

It seems not necessary to store particle positions for every timestep, to calculate dispersion. You may need storage for plotting or reporting, but let us focus on the algorithm, where you need just the current particle coordinates, say

currentX1  % 1xM array
currentX2  % 1xM array

and initial coordinates

initialX1 % 1xM array
initialX2 % 1xM array

So, your domain has two dimensions, x1 and x2. Then, your dispersion is a 2x2 matrix. You need to calculate each element of the dispersion matrix separately. Now, within your timestep, you need to add a loop to go through the M particles.

A2 = zeros(2);  % initialize the dispersion matrix
for m = 1:M
    A2(1,1) = A2(1,1) + (currentX1(m)-initialX1(m))^2/M;
    A2(2,2) = A2(2,2) + (currentX2(m)-initialX2(m))^2/M;
    A2(1,2) = A2(1,2) + (currentX1(m)-initialX1(m))*(currentX2(m)-initialX2(m))/M;
    A2(2,1) = A2(1,2);  % symmetry
end
% Here, you have the instantaneous absolute dispersion in A2

Similarly for the relative dispersion.