System of differential equations

  • MATLAB
  • Thread starter dRic2
  • Start date
  • #1
dRic2
Gold Member
813
198
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:

Answers and Replies

  • #2
12,497
6,284
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?
 
  • #3
dRic2
Gold Member
813
198
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.
 
  • #4
12,497
6,284
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?
 
  • #6
12,497
6,284
  • Like
Likes dRic2
  • #7
DrClaude
Mentor
7,554
3,894
but It gives wrong results.
You'll have to elaborate.
 
  • #8
dRic2
Gold Member
813
198
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.
 
  • #9
kreil
Insights Author
Gold Member
668
67
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.
 
  • Like
Likes AVBs2Systems and dRic2
  • #10
dRic2
Gold Member
813
198
Thanks @Krell. I was just experimenting with Matlab to get used to it, but it is very useful to know!
 

Related Threads on System of differential equations

Replies
1
Views
1K
Replies
1
Views
2K
Replies
11
Views
33K
Top