MATLAB System of differential equations

[dRic2]

Hi, I was trying to solve the simplest problem of planetary motion (for one planet).

The equations should be:

$F_x = m \frac {d^2x} {dt^2} = G \frac {Mmx} {r^3}$
$F_y = m \frac {d^2y} {dt^2} = G \frac {Mmy} {r^3}$

where $r = \sqrt{x^2 + y^2}$

So I re-wrote the system like this:

$\frac {dx} {dt} = v_x$
$\frac {dy} {dt} = v_y$
$\frac {dv_x} {dt} = G \frac {Mmx} {r^3}$
$\frac {dv_y} {dt} = G \frac {Mmy} {r^3}$

I'm not sure how to implement this in matlab...

I tried this

function out = sis(t, s)

% s(1) is x
% s(2) is y
% s(3) is v_x
% s(4) is v_yglobal G Ms Mt
out = zeros(4,1);out(1) = s(3);

out(2) = s(4);r = sqrt(s(1)^2 + s(2)^2);out(3) = G * Ms * Mt * s(1) / (r^3);

out(4) = G * Ms * Mt * s(2) / (r^3);end

but It gives wrong results.

Thanks
Ric

 

[jedishrfu]

Why do you set out(1) = s(3) and out(2) = s(4) shouldn't they be something like out(1) = s(1) + t*s(3) and out(2) = s(2) + t*s(4) ?

and is the v_z really v_y?

 

[dRic2]

and is the v_z really v_y?

yes, my bad.

Why do you set out(1) = s(3) and out(2) = s(4) shouldn't they be something like out(1) = s(1) + t*s(3) and out(2) = s(2) + t*s(4) ?

I want to use ode45 instead of manually implement the algorithm and I though that is the way to do it.

 

[jedishrfu]

You misunderstood what I was asking. My understanding is that the s array is the state of the system and since s1 and s2 represent the position then you must update the position each time with a delta x and a delta y value.

The delta x is in fact $vx * \Delta t$ which in your case is t giving the $out1=s1+vx*t$ expression.

Do you see what I mean?

 

[dRic2]

Yes, but as far as I know it does not work that way... https://it.mathworks.com/help/matlab/math/choose-an-ode-solver.html https://it.mathworks.com/help/matlab/ref/ode45.html
Have I misunderstood the examples ?

 

[jedishrfu]

Maybe this stackexchange will help

https://stackoverflow.com/questions/28516339/basic-2-body-interaction-using-matlabs-ode45

They do the same scheme for out1 and out2 that you did so ignore mine. I used the Open Source Physics ode solvers and you had to provide a step function to modify the state vector at each time delta and it was done as I’ve shown it. The Matlab approach does it different so that is what you should follow.

 

[DrClaude]

but It gives wrong results.

You'll have to elaborate.

 

[dRic2]

Thanks @jedishrfu, very nice link!

First of all, I forgot the minus sign in front of $\frac GMmx/r^3$. And then, as concerned as I was with mathematics, I totally forgot about the physics of the problem! My initial conditions were wrong and physically inconsistent so the numerical integration kept failing.

 

[kreil]

Try ode113 instead of ode45, as it is a variable order solver that includes the results of several past time steps. It is better than ode45 for problems where high precision is required, such as orbital dynamics problems.

Example here:
https://www.mathworks.com/help/matlab/ref/ode113.html#bu6t62o-1_1

ode45 actually outputs more points than it computes - it uses interpolation to get a few extra points per time step, which results in nicer looking plots. That can be adjusted with the 'Refine' option.

 

[dRic2]

Thanks @Krell. I was just experimenting with Matlab to get used to it, but it is very useful to know!