How to properly solve a system of four first order ODEs

Question

2Lindsay2 2018-7-26

0
链接

此问题的直接链接

https://ww2.mathworks.cn/matlabcentral/answers/412167-how-to-properly-solve-a-system-of-four-first-order-odes

评论： David Goodmanson 2018-7-29

(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

7 个评论
显示 5更早的评论隐藏 5更早的评论

Star Strider 2018-7-28

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 2018-7-28

ODESolutions.m

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 2018-7-28

0
链接

此回答的直接链接

https://ww2.mathworks.cn/matlabcentral/answers/412167-how-to-properly-solve-a-system-of-four-first-order-odes#answer_330776

在 MATLAB Online 中打开

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.

2 个评论
显示无隐藏无

2Lindsay2 2018-7-29

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 2018-7-29

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.

请先登录，再进行评论。