Coding Crowell's Method for Orbiting Bodies

2 views (last 30 days)
I'm trying to make a function that takes some initial conditions and models of planetary bodies and it spits out a plot of what happens over time.
To do this I'm using the following differential equation,
Which is the formula given on Wikipedia (Orbit modeling - Wikipedia). The are vectors that represent the location of the bodies, and n is the number of these bodies, G is the gravitational constant 6.672e-11 . Now I'm using the substitutions, , and , to turn this into into a system of first order ODEs. First, before defining the ODE function I'll explain how the inputs work, I use this function here (the plotting part is unfinished, I'm just testing if the caluclations actually work),
function [T,V] = OrbModel(m,x0,v0,t,evalnum)
% INPUT EXPLANATIONS
% OUTPUT EXPLANATIONS
dim = length(x0(:,1));
p = length(m);
T = linspace(t(1),t(2),evalnum);
y0 = [x0(:);v0(:)];
[~,V]=ode45(@(t,y)ODEFunc(y,m),T,y0);
switch dim
case 2
plotres=zeros(1,p);
hold on
for i=1:p:dim*p
plot(V(:,i),V(:,i+1));
end
hold off
case 3
end
m is a vector that contains the masses of each planet. is a n by p matrix, where n is the dimension I'm working in, e.g. in 2 dimensions their positions are described by 2 dimensional vectors and p are the number of planets. Therefore, the first column of is the initial coordinates of the first planet, etc. is the exact same as but initial velocities instead. Now ideally I'd want to pass the matrix into the ode45( . ) function but I can't so I combined and into a single vector, in the code above. Now ODEFunc( . ) is the following,
function dy=ODEFunc(y,m)
G = 6.672e-11;
p = length(m);
dim = length(y)/(2*p);
dimc=dim-1;
dy = zeros(dim*p*2,1);
for i=1:p
dy(i*dim-dimc:i*dim) = y(i*dim-dimc+dim*p:i*dim+dim*p);
for k=[1:i-1, i+1:p]
dy(i*dim-dimc+p*dim:i*dim+p*dim) = dy(i*dim-dimc+p*dim:i*dim+p*dim)...
+G*m(k)*(y(k*dim-dimc:k*dim)-y(i*dim:i*dim))/norm(y(k*dim-dimc:k*dim)-y(i*dim-dimc:i*dim))^3;
end
end
end
which (hopefully) resembles the function used in Crowell's Method with the substitutions above. Now with these two functions I'm testing it by plotting an Earth-Sun system, the inputs I'm using are m=[1.989e30, 5.972e24], x0=[0 0; 1.4756e11 0], v0=[0 0; 0 30000], t and evalnum are unimportant. These inputs I got from a quick search on google, more specifically I'm assuming the sun to be at the center and the Earth , and the velocity of the sun to be 0 and the earth has an upwards velocity of m/s. With these inputs I get the following plot,
The orange line is Earth, and it starts at the bottom right. Now it seems to have relatively appropriate behavior of planetary bodies, arcing back towards the sun at (0,0), but then the earth just leaves, off to the right to infinity (at least it seems like that). Finally, my question is, since it seems to exhibit decent behavior it's really hard to tell if the code itself is incorrect or just my initial conditions are bad. As such, I'm wondering which one it is. If it's the initial conditions how would I go about finding the correct ones?

Accepted Answer

James Tursa
James Tursa on 29 Nov 2020
Edited: James Tursa on 29 Nov 2020
I don't have the patience to go through all of your dim and dimc logic. It would seem to be much easier to code and debug if you reshaped things first. E.g.,
function dy=ODEFunc(y,m)
G = 6.672e-11;
p = length(m);
y = reshape(y,[],2*p); % columns of positions and velocities
x = y(:,1:p); % pick off positions
v = y(:,p+1:end); % pick off velocities
dv = zeros(size(v)); % allocate velocity derivative matrix
for i=1:p
for k=[1:i-1, i+1:p]
dv(:,i) = dv(:,i) + _____; % insert code here to calculate acceleration of i'th planet
end
end
dy = [v(:);dv(:)];
end
You just need to fill in the blank, using x(:,i) and x(:,k) as the positions of the i'th and k'th planets respectively.
  2 Comments
Anthony P
Anthony P on 29 Nov 2020
Edited: Anthony P on 29 Nov 2020
Alright so I changed the function section to this,
function dy=ODEFunc(y,m)
G = 6.672e-11;
p = length(m);
y = reshape(y,[],2*p);
dim = length(y(:,1));
dy = zeros(dim,2*p);
for i=1:p
dy(:,i) = y(:,i+p);
for k=[1:i-1, i+1:p]
dy(:,i+p) = dy(:,i+p) + G*m(k)*(y(:,k)-y(:,i))/norm(y(:,k)-y(:,i))^3;
end
end
dy=dy(:);
end
Which is arguably a lot more readable (funnily enough this was what my code originally was and I only changed it because ode couldn't take matrices, and I forgot reshape existed). In any case, now my orbits are behaving a lot more nicely having begin circular around the sun, it seems like my previous code was incorrect.
Thanks and sorry for wasting your time (and possibly others too).
James Tursa
James Tursa on 29 Nov 2020
Edited: James Tursa on 29 Nov 2020
If your problem was solved then it wasn't a waste!
Note that you really don't need dim. You can just do:
dy = zeros(size(y));

Sign in to comment.

More Answers (0)

Categories

Find more on Numerical Integration and Differential Equations in Help Center and File Exchange

Community Treasure Hunt

Find the treasures in MATLAB Central and discover how the community can help you!

Start Hunting!