Help regarding 2 body planetary simulation

Question

REA 2016-12-17

0
链接

此问题的直接链接

https://ww2.mathworks.cn/matlabcentral/answers/317275-help-regarding-2-body-planetary-simulation

评论： John D'Errico 2016-12-18

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 个评论
显示 1更早的评论隐藏 1更早的评论

REA 2016-12-17

Like this?

Image Analyst 2016-12-17

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.

请先登录，再进行评论。

请先登录，再回答此问题。

Answer 1

John D'Errico 2016-12-17

1
链接

此回答的直接链接

https://ww2.mathworks.cn/matlabcentral/answers/317275-help-regarding-2-body-planetary-simulation#answer_247618

编辑：John D'Errico 2016-12-17

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 个评论
显示 2更早的评论隐藏 2更早的评论

REA 2016-12-17

编辑：REA 2016-12-17

在 MATLAB Online 中打开

Thanks so much for the answer!

I've removed the global part (honestly I didn't even know what it did)and i've edited my code a bit but I'm having problems and I'm not sure if MatLab understands what I'm trying to do. I believe the problem is the radius not updating and thus the acceleration doesn't update either etc.

for step=1:n-1
      R_se=sqrt(x^2 + y^2); % Find the radius between sun and earth
      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
      vx_new=vx+ax_new*dt; %Take "old" x_velocity and add new x_acc*stepsize
      vy_new=vy+ay_new*dt; %Same but with y
      x_new=x+vx_new*dt; %Take "old" x_pos. add new x_vel*stepsize
      y_new=y+vy_new*dt; %Same but with y
      x_new=x; %Rename the newly found x_pos to x
      y_new=y; %Same but with y
      ax_new=ax; %Same
      ay_new=ay; %Same
      vx_new=vx; %Same
      vy_new=vy; %Same
    hold on
    plot(x_new,y_new,'o');
    axis([-1000 1000 -1000 1000]);
    drawnow;
  end

I do get the figure but depending on the axis I choice I either can't see anything or just a very small ring (no idea whats going on). What I wanna know is, does MatLab understand what I'm trying to do (and if so, my code for calculating the new position over time is correct?).

And where does the sun come into play in all of this? I've assigned its start values but I don't use them anywhere and I'm not sure how to introduce it!

One last thing, I'm sure that there's a waaaay better way of creating a piece of code which you can copy-paste into multiple bodies instead of what i've written. Is this true and if so what would be the big change? Would it be to use an ODE solver?

PS. regarding the usage of Euler's method I'm quite sure this is how they expect us to solve this task. I don't know if we're allowed to use any other method (and if it would make sense) because of the level of MatLab experience we have!

John D'Errico 2016-12-17

编辑：John D'Errico 2016-12-17

在 MATLAB Online 中打开

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.)

REA 2016-12-18

编辑：REA 2016-12-18

在 MATLAB Online 中打开

Edit: I believe I have found a major flaw in my calculations. Since the acc. isn't constant I cannot use the usual kinematic equations for position and velocity! I do however believe that the equation for acc. is however still correct (or I could go back to the equation a_earth=(G*m_sun)/(sqrt(x_earth^2+y_earth^2))^2.

Thank you again for the great answer! I think you might've misunderstood me with regards to the sun. I do know that it is indeed affected by the gravity field of other objects but our assignment wants us to ignore this very small effect it has and thus I'd like to know if I should make correlating motion equations for the sun in which the sum would always be 0 (because of the start values all being 0). Or if I should just get MatLab to draw a yellow dot in Origo that does nothing at all.

Anyway, I looked at my code again and did a few tweaks (I added real values for mass of sun and earth and earth pos. taken from JPL Nasa and also the G constant) and it now runs (which is probably because I changed the order of some stuff around) but the orbit it takes is very odd. Also I have no idea how the axis should be (I'm assuming it's in AU and if so, then it should be around 10 so that when I add other planets I'd actually be able to see them (if their distance to the sun is further than 1 AU)

Heres my new code:

%%Constants
m_0=1.9891E30; % Sun mass in kg
m_1=5.972E24; % Earth mass in kg
G=6.674*10^(-11); % Gravitational constant %AU^3/(kg*day^2)
n=3650; %1 year is 36.5 because 1 step is 10 days (According to JPL Horizon data) so this is 100 years.
dt=1; %Stepsize is by the DATA from JPL aut. set to 10 days (dt=1 aka 10 days)
x=-1.667607819720000E-01; %AU
y=9.690699368667393E-01; %AU
vx=-1.723469549881712E-02; %AU/Day
vy=-2.976030157857586E-03; %AU/Day
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);
%%So I first find the radius and acc. of my given start pos. and velocities and then I find the new values through my loop!
    for step=1:n-1
        x_new=x+vx*dt
        y_new=y+vy*dt; 
        vx_new=vx+ax*dt
        vy_new=vy+ay*dt;
        R_se_new=sqrt(x_new^2 + y_new^2)
        ax_new= -G*m_0*x_new/R_se_new^(3/2)
        ay_new= -G*m_0*y_new/R_se_new^(3/2);   
        x=x_new;
        y=y_new;
        vx=vx_new;
        vy=vy_new;
        hold on
        plot(x_new,y_new,'.');
        axis([-1 1 -1 1]); % I think this is in AU and if so this should have the orbit go all around in the periphery of my figure in a perfect circle. If this is true I need it to be higher if I'm going to add other planets which are further away from the sun (and also to create a clearer view).
        drawnow;
    end

Last thing, I saw a video of some guy modelling gravity and he says when getting the values from JPL Nasa that G should be 1.4879*10^(-34) which units AU^3/(kg*day^2) instead of the value because the pos. and velocity data for earth is in AU for pos and AU / Day for velocity. And I just wanna make sure this is correct because the trajectory of my planet is closely tied to my G constant.

Thanks again for all the help!

John D'Errico 2016-12-18

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.

请先登录，再进行评论。

Help regarding 2 body planetary simulation

3 个评论
显示 1更早的评论隐藏 1更早的评论

采纳的回答

4 个评论
显示 2更早的评论隐藏 2更早的评论

更多回答（0 个）

另请参阅

类别

标签

Community Treasure Hunt

Help regarding 2 body planetary simulation

3 个评论 显示 1更早的评论隐藏 1更早的评论

采纳的回答

4 个评论 显示 2更早的评论隐藏 2更早的评论

更多回答（0 个）

另请参阅

类别

标签

Community Treasure Hunt

3 个评论
显示 1更早的评论隐藏 1更早的评论

4 个评论
显示 2更早的评论隐藏 2更早的评论