System of differential equations

  • MATLAB
  • Thread starter dRic2
  • Start date
  • #1
dRic2
Gold Member
875
233
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
13,578
7,566
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
875
233
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
13,578
7,566
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
13,578
7,566
  • #8
dRic2
Gold Member
875
233
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
68
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
875
233
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
R
Replies
1
Views
2K
Replies
11
Views
33K
  • Last Post
Replies
1
Views
4K
Replies
2
Views
21K
Replies
1
Views
2K
Replies
1
Views
898
  • Last Post
Replies
6
Views
2K
  • Last Post
Replies
3
Views
791
Top