# System of differential equations

• MATLAB
Gold Member
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

Code:
function out = sis(t, s)

% s(1) is x
% s(2) is y
% s(3) is v_x
% s(4) is v_y

global 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

Last edited:

jedishrfu
Mentor
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?

Gold Member
and is the v_z really v_y?

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
Mentor
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?

jedishrfu
Mentor
dRic2
DrClaude
Mentor
but It gives wrong results.
You'll have to elaborate.

Gold Member

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
Gold Member
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.

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