MATLAB Planetary Motion: Solving ODE45 & Plotting Orbits

[Scatterbrains]

Hello,

Ive got a problem in MATLAB and have searched the internet a lot for help with this but nothing really addresses my issue specifically.\. I wish to use the ode45 solver to calculate and plot the orbits of all the planets in the solar system. I can manage just fine to get them to go around the sun but am struggling with even understanding how to get them to interact with each other in terms of their gravity. I've put what I've done for Earth and venus below, please forgive the sloppy work I am not very good at it yet.
function Gravity1 ()

clear;
format long eng;

StartTime = 0;
TimeStep = 24 * 3600 * 5;
EndTime = 24 * 3600 * 365;


%Venus initial conditions

VX0 = 1.08e11; %Position X
VY0 = 0; %Position Y
VVX0 = 0; %Velocity X
VVY0 = 35000; %Velocity Y

%Earth intial conditions

EX0 = 1.496e11; %Position X
EY0 = 0; %Position Y
EVX0 = 0; %Velocity X
EVY0 = 29500; %Velocity Y


TVec = StartTime:TimeStep:EndTime;

%Column Vectors for initial conditions
EarthVec = [VX0; VY0; VVX0; VVY0];
VenusVec = [EX0; EY0; EVX0; EVY0];

%Built in solver
[TimeVec, EarthMat] = ode45(@GravityDE, TVec, EarthVec);
[TimeVec, VenusMat] = ode45(@GravityDE, TVec, VenusVec);


EarthXVec = EarthMat (:,1);
EarthYVec = EarthMat (:,2);
VenusXVec = VenusMat (:,1);
VenusYVec = VenusMat (:,2);

hold on;
plot (EarthXVec,EarthYVec, '.');
plot (VenusXVec,VenusYVec,'.');
function dYVec = GravityDE (TVec, YVec)


MStar = 1.99e30;
GravConst = 6.67259e-11;

X = YVec(1);
Y = YVec(2);
VX = YVec(3);
VY = YVec(4);

RangeSq = X.^2 + Y.^2;
Theta = atan2(Y,X);
AccMag = (GravConst .* MStar ./ RangeSq);
AccX = -AccMag .* cos (Theta);
AccY = -AccMag .* sin (Theta);

dYVec = [VX; VY; AccX; AccY];
Sorry about the minimal comments, I am stressed out and have butchered things here and there, I hope you can follow this. Any help on improvements would also be appreciated. I really want to understand this and figure out how to deal with more than 2 bodies.

Thanks

 

[superkan619]

That sounds great, "certainly I'm going to learn MATLAB someday"! But I don't think its possible to make planets revolve while considering gravitational interaction b/w them. This is so because I've learned that nobody has ever been able to solve the 3-BODY problem.

 

[D H]

You need to

• Make a *big* vector that incorporates all of the planetary states in one structure.
• Consider the Sun to be a planet. It orbits, too.
• Use a better form for the gravitational acceleration (i.e., use vectors only. That atan2 stuff forces an oversimplistic form.)
• Compute the gravitational interactions among all bodies. (i.e., you need a loop).

If you want to be at all realistic you need to go to 3D. The planets do not orbit in the same plane.

 

[D H]

But I don't think its possible to make planets revolve while considering gravitational interaction b/w them. This is so because I've learned that nobody has ever been able to solve the 3-BODY problem.

That is just wrong. Of course the N-body problem is solvable. It just can't be expressed in terms of elementary functions. One has to resort to numerical integration techniques to find a solution.

 

[Scatterbrains]

You need to

• Make a *big* vector that incorporates all of the planetary states in one structure.
• Consider the Sun to be a planet. It orbits, too.
• Use a better form for the gravitational acceleration (i.e., use vectors only. That atan2 stuff forces an oversimplistic form.)
• Compute the gravitational interactions among all bodies. (i.e., you need a loop).

If you want to be at all realistic you need to go to 3D. The planets do not orbit in the same plane.

Thanks for the reply. I think I can manage most of that however what I don't know how to do is to calculate the gravitational attraction between a planet at a position and a planet at another position. Would I use a similar function to the one I had before? the mind boggles.

My vector knowledge isn't hugely advanced unfortunately.

 

[belliott4488]

In general, the gravitational force on a body of mass m₁ at position r₁ due to a second body, of mass m₂ at position r₂ is given by

F = Gm₁m₂ (r₂ - r₁)/|r₂ - r₁|³

The vector part is just the inverse square dependence on the separation distance, in the direction from body 1 to body 2.

If this level of vector math is unfamiliar to you, then I think you might find a detailed simulation of this problem to be kind of difficult. Vectors are used precisely because they make problems like this easier.

 

[Scatterbrains]

In general, the gravitational force on a body of mass m₁ at position r₁ due to a second body, of mass m₂ at position r₂ is given by

F = Gm₁m₂ (r₂ - r₁)/|r₂ - r₁|³

The vector part is just the inverse square dependence on the separation distance, in the direction from body 1 to body 2.

If this level of vector math is unfamiliar to you, then I think you might find a detailed simulation of this problem to be kind of difficult. Vectors are used precisely because they make problems like this easier.

That does look familiar, thanks. I could do all this stuff on paper just not terribly well in matlab. I'll get there eventually.

 

[Scatterbrains]

Ok now the next problem is that when i put the state vectors into a column as below:


%Venus initial conditions

VX0 = 1.08e11; %Position X
VY0 = 0; %Position Y
VVX0 = 0; %Velocity X
VVY0 = 35000; %Velocity Y

%Earth intial conditions

EX0 = 1.496e11; %Position X
EY0 = 0; %Position Y
EVX0 = 0; %Velocity X
EVY0 = 29500; %Velocity Y


InitialStateVec = [VX0; VY0; VVX0; VVY0; EX0; EY0; EVX0; EVY0]


I get an error in MATLAB when it tries to do the calculations with the GravityDE function:


? Error using ==> odearguments at 116
Solving GRAVITYDE requires an initial condition vector of length 8.


It works fine when I do just the 4 states for the 1 planet, but it doesn't like it for 2.


Once again, thanks for taking the time to read this.

 

[D H]

Did you rewrite your GRAVITYDE function to work with vectors of length 8?

 

[Scatterbrains]

Ok that was a problem, I thought I did rewrite it and saved it but MATLAB kept using an older version for some reason. Figured that bit out now.

I've tried to incorporate loops to calculate how each body affects the acceleration of the planet, in the case of 3 bodies there's a for loop with 2 nested for loops each of which calculated the acceleration due to the star and the other for the other planet. I am not sure if I did this correctly nor do I know how to set the for loop condition, which I assume is the time period over which this is all taking place...not sure!

I'll copy what I've changed so far below, this assuming the input for the function is the column containing the 4 states of the 2 planets. Any suggestions are welcome.



function dYVec = GravityDE (TVec, YVec)




%Calculate the individual accelerations due to bodies then add them
%together


%Calculations for Venus
for

X = YVec(1);
Y = YVec(2);
VX = YVec(3);
VY = YVec(4);
X2 = YVec(5);
Y2 = YVec(6);
VX2 = YVec(7);
VY2 = YVec(8);

%Acceleration due to star

for

MBody = 2.00e30;
[AccX,AccY] = Acceleration(X,X2,Y,Y2);
AccX1 = AccX;
AccY1=AccY;

end

%Acceleration due to Earth

for

Mbody = 5.98e24
[AccX,AccY] = Acceleration (X,X2,Y,Y2);
AccX2=AccX;
AccY2=AccY;

end

%Calculate sum of Acceleration Vectors
VenusAccX = AccX1 + AccX2;
VenusAccY = AccY1 + AccY2;
%Venus Velocity Vectors
VenusVX = VX;
VenusVY = VY;


end


%Calculations for Earth
for

X = YVec(5);
Y = YVec(6);
VX = YVec(7);
VY = YVec(8);
X2 = YVec(1);
Y2 = YVec(2);
VX2 = YVec(3);
VY2 = YVec(4);

%Acceleration due to star

for

MBody = 2.00e30;
[AccX,AccY] = Acceleration(X,X2,Y,Y2);
AccX1 = AccX;
AccY1=AccY;

end

%Acceleration due to Venus

for

Mbody = 5.98e24
[AccX,AccY] = Acceleration (X,X2,Y,Y2);
AccX2=AccX;
AccY2=AccY;

end

%Calculate sum of Acceleration Vectors

EarthAccX = AccX1 + AccX2;
EarthAccY = AccY1 + AccY2;
%Earth Velocity vectors

EarthVX = VX;
EarthVY = VY;


end



dYVec = [VenusVX; VenusVY; VenusAccX; VenusAccY; EarthVX; EarthVY; EarthAccX; EarthAccY];




function [AccX,AccY] = Acceleration(X,X2,Y,Y2)


GravConst = 6.67259e-11;
RangeSq = (X2-X).^2 + (Y2 -Y).^2;
Theta = atan2 ((X2-X),(Y2-Y));
AccMag = GravConst .* MBody ./ RangeSq;
PlntAccX = -AccMag .* cos (Theta);
PlntAccY = -AccMag .* sin (Theta);


I also thought in future I'd include the state of the star in the input matrix for this function so that it is easier to calculate. Though, I think I'm going to just make it so that the star is stationary in space as my head is already hurting from this alone. Thanks very much so far.
I'm still using the atan function. Once I figure this whole thing out I'll work on getting that sorted. Thanks very much so far.

 

[superkan619]

That is just wrong. Of course the N-body problem is solvable. It just can't be expressed in terms of elementary functions. One has to resort to numerical integration techniques to find a solution.

Thank You D H.