How to properly solve a system of four first order ODEs

(Just for some background) Using the equations for the velocity of point vortices in a fluid, you obtain 2*N first order ODEs that describe the position. I'm trying to get the solution of the positions for just 2 point vortices to start off with, but can't see my error, but the result is in fact wrong right now. Any help would be awesome!

code is below, as well as the formula I am using

if true
%gamma is the strength of the vortices
gamma1 = 2;
gamma2 = 1;
%y1 and y2 become x(1) and x(2)
%x1 and x2 become x(3) and x(4)
f = @(t,x) [gamma2/2/pi*(-(x(1)-x(2)) / ((x(3)-x(4))^2 + (x(1)-x(2))^2));
  gamma1/2/pi*(-(x(2)-x(1)) / ((x(4)-x(3))^2 + (x(2)-x(1))^2));
  gamma2/2/pi*((x(3)-x(4)) / ((x(3)-x(4))^2 + (x(1)-x(2))^2));
  gamma1/2/pi*((x(4)-x(3)) / ((x(4)-x(3))^2 + (x(2)-x(1))^2))];
%Initial conditions are the initial positions
[t,x] = ode45(f, [0 10], [3 1 3 1]);
subplot(1,2,1)
plot(x(:,3),x(:,1))
xlabel("x1")
ylabel("y1")
subplot(1,2,2)
plot(x(:,4),x(:,2))
xlabel("x2")
ylabel("y2")

end

[Aarti Dwivedi]

Please format the code properly, and use the {}Code option for greater ease of reading.

[Star Strider]

‘... the result is in fact wrong right now.’

What is the correct result?

[2Lindsay2]

Solving these by hand is impossible, so I can't give you an exact answer to compare to. But in general, by looking at the graphs, the distance between the two points (the vortex locations) at any time is constant.

Furthermore, If both gamma1 and gamma2 are EQUAL and the vortex spins are in the SAME direction(designated by the sign of gamma), you get a circular path of constant R where each vortex moves around the midpoint of the imaginary line connecting the two.

For the numbers above (gamma1 > gamma2), It should still be circular path traveled by each, but the point on the imaginary line that they each travel around is shifted closer to vortex 1 of strength gamma1.

Based on the results I got, the vortices diverged away from each other. Which is wrong.

It's easier to see if you overlay the two plots using "hold on".

Could this be something to do with how ODE45 is solving the system?

[Star Strider]

We only have a segment of your code. Providing the complete function and the ode45 call as well as the way you plot the results would be helpful.

Fluid dynamics is not an area of my expertise, at least, so any additional information you can provide may also be helpful.

[2Lindsay2]

I've included the full code that I am using (hopefully it's all showing, if not let me know and I can just type it all in manually)(19 lines). To overlay the two plots I replaced everything after the ode45 call with

if true
  plot(x(:,3),x(:,1))
  hold on
  plot(x(:,4),x(:,2))
end

and get the following graph when the time span is take out to 100 seconds.

No worries on understanding the fluid mechanics here, I just can't figure out why I'm not getting a correct solution from the ode45 function.

[Star Strider]

Please copy and paste your entire ODE function code and your calling script with the ode45 call, including the initial conditions. Also include all constants you may use.

What should the figure look like?

[2Lindsay2]

Here is the m file for you to look at. This is all the code that I am using. I'm not using a function script, just the function handle in the one main script ODESolutions.m.

In the figure the center of curvature should be the same.

Answer 1

[David Goodmanson]

Hi 2L2,

I don't think you are doing yourself any favors with the notation

%y1 and y2 become x(1) and x(2)

%x1 and x2 become x(3) and x(4)

In your function, try swapping the first two rows with the second two rows, to obtain

f = @(t,x) [gamma2/2/pi*((x(3)-x(4)) / ((x(3)-x(4))^2 + (x(1)-x(2))^2));
            gamma1/2/pi*((x(4)-x(3)) / ((x(4)-x(3))^2 + (x(2)-x(1))^2));
            gamma2/2/pi*(-(x(1)-x(2)) / ((x(3)-x(4))^2 + (x(1)-x(2))^2));
            gamma1/2/pi*(-(x(2)-x(1)) / ((x(4)-x(3))^2 + (x(2)-x(1))^2));];
%Initial conditions are the initial positions
[t,x] = ode45(f, [0 100], [3 1 3 1]);

Running the final time out to 100 shows each pair going almost all the way around a circle.

[2Lindsay2]

Wow!!! I never would have guessed that switching rows would make a difference. Do you know why that is?

Thank you so much!

[David Goodmanson]

That's because your code did not correctly represent the formulas, and the corrected code happens to be the same as swapping the rows in the code you had.

The formulas don't say so, but u is a vector of x velocities v is a vector y velocities. Your basic setup has a four element vector [y1; y2; x1; x2] so the corresponding velocity vector is [vel_y1; vel_y2; vel_x1; vel_x2]; in terms of the formulas that is [v1; v2; u1; u2]. If you write new code with this starting point, it's the same as swapping the rows in what you had.