The Sun-Earth System

Comp110 Course Module 6a

Our goal in this module is model the trajectory of the earth and the sun as they move due to Newton's law of gravity and his third law of motion.  Our approach will be to derive a system of differential equation from these laws, solve this system using Mathematica, and then visualize the resulting orbits.  One advantage of this computational approach is that we can then tinker with the law of gravity and see what happens as we modify this law.

While the derivation that I present may seem a little complicated at first, note that I will end up providing the basic differential equations that govern orbital motion.  However, you should still read my derivation closely, since you will need to modify these equations to add the moon to this system in Module 6b.

Deriving the equations of motion

Modeling the position of the Earth and Sun

Given that our task is to solve for the positions of the Earth and the Sun as they orbit each other due to the law of gravity, we must first settle on a coordinate system for representing the positions of the Earth and the Sun.  To this end, we will model the Earth and the Sun as points pe and ps whose positions vary as a function of time t.  If each of these points has an x-coordinate and a y-coordinate, we can view pe[t] and ps[t] as being time-varying functions that return the positions of the Earth and the Sun.

pe[t_]={xe[t],ye[t]};
ps[t_]={xs[t],ys[t]};

One nice result of treating the position of the Earth as a time-varying function is that the velocity of the Earth is simply the derivative of pe[t] with respect to t, i.e. the vector pe'[t].  Similarly, the acceleration of the Earth is simply the vector pe''[t], the second derivative of pe[t] with respect to t.

Newton's law of gravity

Newton's law of gravity states that the attractive force between two celestial bodies is proportional to the product of the masses of the two bodies divided by the square of the distance between the two bodies.  The constant of proportionality G is known as the gravitational constant. If the Sun has a mass ms and the Earth has a mass me , then the force of gravity is

((G ms me)/(r^2))

where r is the distance between the Sun and the Earth.  Newton's law of gravity is often referred to as an inverse square law, since gravity is proportional to (1/(r^2)).  In this module, we will experiment with modified laws of gravity where we replace (1/(r^2)) by (1/r) and (1/(r^3)).

Gravity has a direction as well as a magnitude.  Thus, the gravitational force fe[t] of the Sun on the Earth can be thought of as a vector whose length is the scalar ((G ms me)/(r^2)) and whose direction points from the Earth towards the Sun along the vector ps[t]-pe[t].

We can derive a more explicit formula for fe[t] in terms of ps[t] and pe[t] by building a unit-length vector from the Earth to the Sun and scaling this vector by ((G ms me)/(r^2)).  This unit vector has the form

((ps[t]-pe[t])/(leng[ps[t]-pe[t]])),

where the function leng[v] computes the length of a vector v={x,y} via

leng[{x_,y_}]:=(Sqrt[(x^2)+(y^2)])

Note that the expression leng[ps[t]-pe[t]] is simply the distance r in the expression ((G ms me)/(r^2)).  Therefore, the gravitational force fe[t] of the Sun on the Earth satisfies

fe[t]==((G ms me )/((leng[ps[t]-pe[t]])^2)) ((ps[t]-pe[t])/(leng[ps[t]-pe[t]]))==((G ms me (ps[t]-pe[t]))/((leng[ps[t]-pe[t]])^3))

Notice that the factor leng[ps[t]-pe[t]] in the denominator is cubed, even though gravity is an inverse square law.  The third factor in the denominator comes from multiplying the inverse square law by the unit vector
	((ps[t]-pe[t])/(leng[ps[t]-pe[t]])).
Notice too that the term ps[t] - pe[t] in the numerator does not cancel one of the factors of leng[ps[t] - pe[t]] in the denominator, since ps[t] - pe[t] is a vector whereas leng[ps[t] - pe[t]] is a scalar.

Newton's third law

To complete our derivation of the differential equations that govern motion, we must use Newton's third law of motion to tie the effect of gravitational force to the acceleration of the Earth.  According to Newton's third law,  a force applied to a body is equal to the body's mass times the body's acceleration.  Therefore, in our case, gravity causes the Earth to accelerate according to the relation

fe[t]=me pe''[t].

If we substitute the definition of fe[t] into this equation and cancel out the factor of me on both sides, we arrive at the fundamental differential equation that describes the motion of the Earth based on the positions of the Earth and the Sun,

pe''[t]== ((G ms)/((leng[ps[t]-pe[t]])^3))(ps[t]-pe[t]).

Solving the equation of motion

Luckily for us, Mathematica is very good at solving differential equations of this type.  To set up and solve this equation, we must first choose correct values for the constants ms, me and G.

G= 4.981×(10^(-10));
ms=1.989×(10^30);
me=5.977×(10^24);

Here, the gravitational constant (G) is in units of kilograms-(kilometers^3)/(day^2).  Given that the mass of the Sun is so large in comparison to the mass of the Earth, we will make a simplifying assumption that the Sun simply stays fixed in one spot.  To make our life easier, let's assume that the Sun stays at the origin.

ps[t_]:={0,0}

At this point, we can let Mathematica derive the actual differential equations in terms of the unknown position {xe[t],ye[t]}.  The Mathematica variable sunEarthDiffEqs will contain these differential equations,

SetAttributes[Equal,Listable]

sunEarthDiffEqs=pe''[t]== ((G ms)/((leng[ps[t]-pe[t]])^3))(ps[t]-pe[t])

Note that the resulting equations are fairly simple (although they do contain a square root).  These equations involve two unknown functions xe[t] and ye[t] that we will solve for by invoking NDSolve.  To help NDSolve find a unique solution to these equations, we need to specify an initial position and an initial velocity for the Earth.

sunEarthInitConds = 
{xe[0]==1.495979×(10^8),
ye[0]==0,
xe'[0]==0,
ye'[0]==2.56721625×(10^6)}

Putting all of these ingredients together in the Mathematica function NDSolve, we can compute the orbit of the Earth for 1000 days.  For future use, we will store this orbit in the function earthTraj.

earthTraj[t_]=pe[t]/.
NDSolve[{sunEarthDiffEqs,sunEarthInitConds},{xe,ye},{t,0,1000}][[1]]

Observe that Mathematica produces a vector of functions depending on t that describe the Earth's orbit.

Visualizing the orbits

To complete our solution, we need to visualize the orbit of the Earth.  We will visualize this orbit in two ways:  as a static trajectory and as a dynamic graphic using Manipulate.  The function ParametricPlot does a nice job of plotting the circular trajectory of the Earth.

earthOrbit=ParametricPlot[earthTraj[t],{t,0,1000},AspectRatio->Automatic]

We can also create a dynamic graphics of the Earth orbiting the Sun using Manipulate.  The yellow disk is the Sun and the blue disk is the Earth.  Note that holding down Alt while you move the slider cause the Earth to move very slowly.

radE=6400 500;
radS=10 700000;

Manipulate[Show[
earthOrbit,
Graphics[{
RGBColor[1,1,0],
Disk[{0,0},radS],
RGBColor[0,0,1],
Disk[earthTraj[t],radE]},
PlotRange->1.5*(10^8),Axes->True]],
{t,0,1000}]

Orbital motion under a (1/(r^2)) law

Experimenting with different laws of gravity

In his famous book, The Principia, Isaac Newton introduced the inverse square law of universal gravitation.  In the same book, Newton also studied other possible laws of gravity, where the magnitude of the force depends on the distance in ways other than inversely with the square of the distance.  For example, in addition to studying a (1/(r^2)) law, Newton also analyzed laws depending on (1/(r^3)) and (1/r).  In the next part of this module, you are going to see what the planetary paths would look like in another universe, where the laws of gravity are slightly different from our own.  This ability to simulate other worlds is one of the great powers of the computer.  We will now perform experiments on the computer that are impossible to perform in the real world.

Orbital motion under a (1/(r^3)) law

sunEarthInitConds = 
{xe[0]==1.495979×(10^6),
ye[0]==0,
xe'[0]==0,
ye'[0]==2.104*(10^4)};

Orbital motion under a (1/r) law

sunEarthInitConds = 
{xe[0]==1.495979×(10^10),
ye[0]==0,
xe'[0]==0,
ye'[0]==3.14865*(10^10)};