MATLAB Planetary Motion: Solving ODE45 & Plotting Orbits

In summary, the gravitational attraction between two bodies depends on their separation distance and the mass of the two bodies.
  • #1
Scatterbrains
8
0
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
 
Physics news on Phys.org
  • #2
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.
 
  • #3
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.
 
  • #4
superkan619 said:
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.
 
  • #5
D H said:
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.
 
  • #6
In general, the gravitational force on a body of mass m1 at position r1 due to a second body, of mass m2 at position r2 is given by

F = Gm1m2 (r2 - r1)/|r2 - r1|3

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.
 
  • #7
belliott4488 said:
In general, the gravitational force on a body of mass m1 at position r1 due to a second body, of mass m2 at position r2 is given by

F = Gm1m2 (r2 - r1)/|r2 - r1|3

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.
 
  • #8
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.
 
  • #9
Did you rewrite your GRAVITYDE function to work with vectors of length 8?
 
  • #10
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.
 
  • #11
D H said:
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.
 

FAQ: MATLAB Planetary Motion: Solving ODE45 & Plotting Orbits

1. How can I use MATLAB to solve ODE45 for planetary motion?

To solve ODE45 for planetary motion in MATLAB, you can use the ode45 function, which is specifically designed for solving ordinary differential equations. You will need to define the initial conditions and the function that represents the motion of the planet, and then use the ode45 function to solve the equation and plot the results.

2. Can I use MATLAB to plot the orbits of multiple planets?

Yes, you can use MATLAB to plot the orbits of multiple planets. You will need to define the initial conditions and the functions for each planet's motion, and then use the ode45 function to solve the equations and plot the orbits on the same graph.

3. How do I interpret the output of the ODE45 function in MATLAB?

The output of the ODE45 function in MATLAB is a vector containing the values of the dependent variables at each time step. These values can be used to plot the trajectory of the planet's motion, and can also be used to calculate other properties such as velocity and acceleration.

4. Is it possible to add other forces, such as gravitational attraction from other planets, to the ODE45 function in MATLAB?

Yes, it is possible to add other forces to the ODE45 function in MATLAB. You can add additional terms to the function that represents the motion of the planet, taking into account the forces from other bodies in the system. You will also need to define the initial conditions for these additional forces.

5. Can I use MATLAB to simulate the motion of a spacecraft in a planetary system?

Yes, you can use MATLAB to simulate the motion of a spacecraft in a planetary system. You will need to define the initial conditions and the functions for the motion of the spacecraft and the planets, and then use the ode45 function to solve the equations and plot the results. You can also add other forces, such as thruster firings or gravitational assists, to make the simulation more realistic.

Similar threads

Replies
2
Views
3K
Replies
2
Views
7K
Replies
2
Views
4K
Replies
1
Views
2K
Replies
1
Views
2K
Back
Top