Help regarding 2 body planetary simulation
Show older comments
Hi,
I've been tasked to create a simple 2D planetary simulation involving a stationary sun and 3 orbiting planets. I started out trying to create a simulation involving the sun and 1 planet only and simplified their masses and the constant G. I believe the acceleration should be constant which means once I find the acceleration I will reuse it to find the new velocities and positions.
I'm quite new to MatLab and as such the code I want should be simple. We are tasked to solve this using Newton's gravitational laws. Furthermore the planets should only be affected by the gravitational force from the sun (and the sun shouldn't be affected by any gravitational force from the planets).
Here is my attempt so far:
%%Constants
global G m_0 m_1
m_0=1000; % Sun mass
m_1=1; % Earth mass
G=10^(-3); %Gravitational constant
n=1000; % Steps
dt=0.1; % Stepsize
%%Equations
% R_se = sqrt(x^2 + y^2); % Distance between earth and sun
% F_se = G*m_0*m_1/R_se^2; % The force the sun affects the earth with
% a_se =G*m_0/R_se^2; % a_es = F_es/m_1 Found with Newtons 2. law f=m*a
% xdotdot= -G*m_0*x/R_se^(3/2) % The acc. in x
% ydotdot= -G*m_0*y/R_se^(3/2) % The acc. in y
%%Start pos. and velocity for sun and earth
%Sun - stationary
x_sun=0;
y_sun=0;
vx_sun=0;
vy_sun=0;
%%Earth - these values have been randomly chosen but I wish to use Ephimerides (?) from Nasa with actual values once I get my code working
x=-10;
y=11;
vx=9;
vy=-5;
Calculations
%%Finding start values for radius and acc. (Because im assuming acc. to be constant - which im afraid may not actually be true)
R_se = sqrt(x^2 + y^2);
ax= -G*m_0*x/R_se^(3/2);
ay= -G*m_0*y/R_se^(3/2);
for n=1:n-1
%%Finding new values for radius, velocity and position.
vx_new=vx+ax*dt;
vy_new=vy+ay*dt;
x_new=x+vx_new*dt;
y_new=y+vy_new*dt;
R_se_new=sqrt(x_new^+y_new^2);
vx=vx_new;
vy=vy_new;
x=x_new;
y=y_new;
R_se=R_se_new;
hold on
plot(x_new,y_new),'.');
axis([-100 100 -100 100]);
drawnow;
end
So my thinking here is that first I find the acceleration which will be constant. Then I find the new velocity in the x and y plane by using the old velocity. From that I find the new positions using the old position and the new velocity times the stepsize. And then i find a new radius between the earth and sun by using the new positions.
Once this is over I say that the old values are the new ones so that when it loops back the old values will be the new ones and the newer values will be determined by the new ones. (If you know what I mean).
The problem:
When I run this I do get a figure but depending on the start values i choose for earth i get 2 moving lines of dots only going in the x plane. I have tried changing the values so many times with no luck.
First of all the sun shouldn't be moving and should only appear as one stationary dot and second of all the earth should be moving in a circle around the sun.
Sorry for the big wall of text!
3 Comments
Image Analyst
on 17 Dec 2016
REA
on 17 Dec 2016
Image Analyst
on 17 Dec 2016
Yes that's better, though you might have gotten rid of the double spacing which you probably put in when you didn't know how to format it properly. Be sure to read John's answer below.
Accepted Answer
John D'Errico
on 17 Dec 2016
Edited: John D'Errico
on 17 Dec 2016
1 vote
First of all, why do you need global variables? You are not calling anything anyway!
Learn to pass variables in as arguments to functions that will need them. Globals are a novice programming crutch that you never really need, but will cause you problems in the future, by making your code buggy and far more difficult to debug.
Next, you keep on talking about acceleration being a constant. Nothing could be further from the truth! While the gravitational constant is a constant, acceleration will vary with time! If the acceleration on your hypothetical earth really was a constant, then your planet would move off in a straight line, at faster and faster velocities, eventually exceeding the speed of light. I hope you are planning on adding relativistic motion to this simulation. Otherwise, when the earth exceeds the speed of light, things might get interesting.
Looking at your code, in fact, you set ax and ay outside of the loop, not changing them ever. So that is exactly what you are doing. Spaceship earth, reporting for duty! This planet is headed outside of the solar system, and soon, the galaxy.
The point is, in order for a planet to maintain an orbit, the accelerations in the x and y directions will indeed vary. In fact, the acceleration of your planet will always be pointed inwards, towards the sun. But the direction of that acceleration will vary, depending where the planet is in its orbit around the sun.
You can think of this in a polar coordinates form. So the acceleration is always radially inwards. This may be what you are thinking, that in polar coordinates, as long as the radius stays constant, then acceleration is indeed constant. But the direction of the radial vector is constantly changing as you move around a circle. The velocity will be orthogonal to the acceleration, if the planet is in orbit. So your planet stays constantly accelerating towards the sun, but stays in orbit. Yes, it seems vaguely paradoxical, but this is how orbital motion works.
So, yes, you need to make G a constant. They call it the gravitational constant for a good reason. But you cannot make ax and ay constant!
I'd also point out that the position of the sun is itself not constant. So a better model has the sun reacting to the forces on it, moving in tiny circles itself. Essentially, the sun and the earth move in concert, each acting on the forces they feel.
Finally, one thing you would realize is that over a long period of time, the predicted positions of those bodies would accumulate errors. You are using Euler's method to integrate what is essentially a system of differential equations. It would be a better choice to use a better method for this problem. So you might formulate the problem as a set of differential equations, then feed that to an ODE solver like ODE45.
I'm not saying ODE45 is what you truly need to use for your problem, but it would be a better way to solve the general problem. The good thing is then you allow the solver to do all the work. Good programming practice uses tools written by experts to solve your problems. Don't get the idea that you need to write these solvers yourself. In fact, too many people seem to think they do. They learned how to solve the problem that way, in some basic class on programming or mathematics, thus ODE solvers, optimization, whatever. Then when it is time for them to do some serious work, they think they need to code up Newton's method themselves. WRONG.
4 Comments
John D'Errico
on 17 Dec 2016
Edited: John D'Errico
on 17 Dec 2016
Better now, but you have lots of questions, and still a few problems. So in random order...
Yes, I am sure you are expected to use Euler's method. It is adequate for the purpose, which here is more to get you coding, learning the skills of programming, as well as teaching you something about numerical methods.
So, yes, it would be better to use an ODE solver. You can formulate everything in terms of differential equations, passing them into a tool like ODE45. But that also does not teach you those programming skills I mentioned before. So, you are stuck with Euler's method. Hey, it was good enough for Euler, and he was no slouch. :)
As far as the sun is concerned, think of it like the earth equations of motion. The earth "feels" a force on it based on the distance to the sun, and the gravitational constant. Thus that force will be proportional to
F = (m1*m2)/d^2
where d is the distance between the bodies.
So what does the sun feel? THE SAME FORCE! But what are the accelerations on the two bodies? Here we can use
F = m*a
to learn that while the sun and the earth both feel the same force,the relative accelerations on the two bodies are widely different, because the two masses are different by a factor of roughly 333000 (according to Google.) Thus since the mass of the sun is so much larger than the earth, the acceleration the sun takes on is so much smaller, by the same proportion.
The point is, to a decent approximation we can decide to treat the sun as motionless, but in fact, we can view both bodies as orbiting around a point at their common center of mass. This is one of the important skills of mathematical modeling - to learn which terms you can ignore as unimportant, and which ones are absolutely crucial to include in your model.
Ok, back to your code. Does MATLAB understand what you are doing? HELL NO! Your computer understands nothing at all. The blasted things are dumb as a fence post. :) So the question is, did you write the correct code? Closer than you were before perhaps, but not yet there. A quick glance shows me these lines:
ax_new= -G*m_0*x/R_se^(3/2); % Find the acc_x for the earth
ay_new= -G*m_0*y/R_se^(3/2); % Find the acc_y for the earth
They are inside the loop now, but x and y stay constant. While you are updating x_new and y_new, you never used them in the computation of those accelerations.
So as the earth moves, you also need to update the values of x and y in your formula to obtain new acceleration components. (I've not checked the formula itself to see if it is correct, beyond not being updated properly.)
John D'Errico
on 18 Dec 2016
I agree that making the sun motionless is a valid approximation. I even explained why that is a valid approximation, and how to know why this would be so. My point was that it IS an approximation, and it might be one of the things one would add into a more complete model of a solar system.
It looks like you are getting closer to a viable model, now that the accelerations are being computed with each iteration.
To ensure it is working properly, I would first VERY carefully check the parameters in the model, since you are using a somewhat bastardized unit system.
m_0 and m_1 I can accept, since those have units of kg.
Next, I would look at G, and carefully verify it is correct for the unit system in terms of Au as a distance unit, with a time unit of days.
Just a programming style point here, I would NEVER use m_0 and m_1 here as variable names. Instead, m_sun and m_earth are just as easy to type, and make your code far easier to verify and debug. Numbered variables like that can be confusing, and will lead to headaches in the future.
So as long as you have just two objects, m_sun and m_earth work well enough. If you were going to add planets though, then you might use a vector (of length 9) for the planet masses. Store the positions, velocities, accelerations each in their own vectors or arrays. Even better would be to store the planet positions in a 9x2 array, thus column 1 would be the x-coordinate, column 2 is y. Then the earth will be found at xy(3,:).
Again, this is what you want to start thinking of for the future, if you would expand this model.
You are getting much closer now though.
More Answers (0)
Categories
Find more on Earth and Planetary Science 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!
