Kepler.html

Kepler Problem

We demonstrate the numerical solution of Newton's equation for the Kepler problem in Cartesian coordinates.

A unit system is chosen to make the equations dimensionless, and we assume that the Sun is infinitely heavy, i.e., we just consider the motion of a single body around the center of attraction.

This worksheet contains several sections:

1) In the first section we introduce a system of coordinates that render the solar system dimensionless. This is achieved by choosing three scales:

the AU (average Sun-Earth distance) as a distance measure;

the solar mass M_s as a mass measure;

the year (time it takes for Earth to go around the Sun) as a time measure.

We solve the Kepler problem numerically for this problem, and demonstrate a few features, such as the virial theorem (relation between kinetic and potential energies).

2) We then explore the problem of radial motion, including a discussion of stability of circular orbits as a function the force law;

3) We then conclude with a brief discussion of the scattering problem.

We begin with section 1:

For the solar system we have the following parameters:

Solar mass: 2*10^{30} kg

Gravitational constant: 6.67*10^{-11} SI (N*m^2/kg^2)

Earth mass: 6*10^{24} kg

Sun-Earth distance 1.5 * 10^{11} m

Earth's orbit time: 3.15*10^{7}s

> restart;

> G:=6.67*10^(-11);

> M_s:=1.99*10^(30);

> M_e:=5.97*10^(24);

> G*M_s*M_e;

> AU:=1.5*10^(11);

> G*M_s*M_e/AU^2;

Earth's attraction towards the Sun is rather large when expressed in Newtons, but the acceleration in m/s^2 is not very large:

> G*M_s/AU^2;

Note that this gravitational acceleration is balanced at the center of the Earth by the centrifugal force due to the rotation about the Sun. The only net effect on an object at the surface of the Earth is in the form of a tidal force.

The period is given as:

> T_e:=evalf(365.26*24*60*60);

> v_e:=evalf(2*Pi*AU)/T_e;

Earth moves with about 30 km/s around the Sun.

Now we introduce a unit system in which the Earth's orbit radius is of order unity, the masses are measured in multiples of Sun's mass, and times are measured in terms of orbit times. We can pick three dimensionful quantities and set the scales as follows [remember: SI=MKS(A)]:

M = M' * M_s

T = T' * T_e

X = X' * AU

and with the substitution we express equations in terms of the new variables.

A similar procedure is used in modern physics to scale the atomic orbital problem; one notable exception is that one uses the electron mass as a scale. In the solar Kepler problem this is not done, because the planet's mass drops out of the equation of motion. In the quantum problem it does not, because the force of attraction is electric, i.e., not proportional to the lighter particle's (i.e., the electron's) mass.

What is our unit of velocity in these natural units for the Earth-around-the-Sun motion problem?

> unitLength:=AU; # introduced unit #1

> unitTime:=T_e; # introduced unit #2

> unitVel:=unitLength/unitTime; # a derived unit!

The velocities V' will come out as multiples of about 4.7 km/s, which means that Earth's cruising speed [about the Sun] becomes 2 Pi!

> v_e/unitVel;

What is the unit of force?

The natural unit should be the force experienced by Earth due to gravity.

The unit for acceleration is clearly the velocity unit divided by a time unit:

> unitAcc:=unitVel/unitTime; # a derived unit!

> unitMass:=M_s; # introduced unit #3 !

> unitForce:=unitMass*unitAcc; # a derived unit!

Given this unit for the acceleration, it is of interest to express the acceleration that Earth experiences from the gravitational attraction towards the Sun in this unit:

> (G*M_s/AU^2)/unitAcc;

The force experienced by Earth is small in these units:

> (G*M_s*M_e/AU^2)/unitForce;

We can find the unit for the gravitational constant from its dimensions:

> unitG:=unitForce*(unitLength)^2/unitMass^2;

Now the centrifugal force in SI units: we assume a point at which the velocity and position are at right angles, so that the magnitude of angular momentum is given by the product of magnitudes of position and velocity times mass:

> L0:=(M_e*AU*v_e);

> L0/(unitMass*unitVel*unitLength); # angular momentum in our units:

> subs(r=AU,-diff(L0^2/(2*M_e*r^2),r));

The gravitational attraction in SI units:

> -G*M_s*M_e/AU^2;

The balance between gravitational attraction and the centrifugal force is not perfect, it is at the level of the number of significant digits carried, which measn that it implies that Earth's orbit is not a perfect circle, and the AU averages over an ellipse with eccentricity close to 0. Earth's eccentricity is listed as 0.017. Mars has an eccentricity of 0.092, Mercury has 0.2

Now let us figure out the balancing in the unit system introduced. We calculate the centrifugal force from differentiating the centrifugal potential expressed as

L^2/(2*m*r^2):

> subs(r=AU,diff(L0^2/(2*M_e*r^2),r))*(unitMass*unitLength^3/(unitMass*unitLength*unitVel)^2);

> (-G*M_s*M_e/AU^2)*(unitLength^2/(unitMass^2*unitG));

> G/unitG;

Now, this establishes what value one should use for G*M_s when working in this unit system:

We take a fresh start, and observe how the Earth (an object of mass m) moves about a fixed gravitational center using a Cartesian coordinate system.

> restart; GM:=39.17: m:=6*10^(24)/(2*10^(30)):

The inverse power for the potential (here we can play if we wish):

> n:=10/10;

> V:=-GM*m/r^n;

From the conservation of the orientation of the angular momentum vector we know that the motion is confined to a plane, i.e., the problem is two-dimensional. We label the coordinates as x and y .

> V:=subs(r=sqrt(x^2+y^2),V);

The force is calculated as follows:

> Fx:=-diff(V,x);

> Fy:=-diff(V,y);

Note that both components of the force depend on both x and y . We make this more obvious by unapplying the two expressions into functions:

> FX:=unapply(Fx,x,y): FY:=unapply(Fy,x,y):

We can now evaluate the force vector at some location. Let us demonstrate that the force is central, i.e., points towards the origin of the coordinate system:

> F0:=[FX(1,2),FY(1,2)];

The y-component of the force vector at this location given by r = [1,2] is twice as large as the x-component, and this is sufficient for the force vector to be proportional to the position vector.

The Newton equation leads to a pair of coupled second-order equations:

> NE1:=m*diff(x(t),t$2)=FX(x(t),y(t));

> NE2:=m*diff(y(t),t$2)=FY(x(t),y(t));

Let us define some intial condition: we pick a point on the positive x -axis, with some velocity in the y-direction, but no velocity component in the x -direction, i.e., the initial position and velocity vectors are at right angles:

> x0:=1: y0:=0: vx0:=0: vy0:=6.28: # these conditions illustrate Earth's orbit

> x0:=1: y0:=0: vx0:=0: vy0:=4.: # these conditions give substantial ellipticity

> IC:=x(0)=x0,D(x)(0)=vx0,y(0)=y0,D(y)(0)=vy0;

> sol:=dsolve({NE1,NE2,IC},{x(t),y(t)},numeric,output=listprocedure):

> X:=subs(sol,x(t)): Y:=subs(sol,y(t)):

> [X(1),Y(1)];

This confirms that a trajectory is being calculated. So let us produce a graph. We need to put X(t) and Y(t) in quotes, because they evaluate to numbers only for numerical values of t .

> plot(['X(t)','Y(t)'],t=0..3,title="x(t) and y(t)",color=[blue,green],thickness=3);

OK, now we plot this parametrically:

> plot(['X(t)','Y(t)',t=0..3],title="parametric plot",color=red,thickness=3,scaling=constrained);

Exercise 1:

Change the intial conditions for the equations and produce the graphs. Watch how the period of the motion changes with the four parameters that determine the initial condition.

Exercise 2:

Change the potential function to a slightly different distance dependence from 1/ r . Observe what happens to the motion.

Let us visualize better what is going on by showing the vectorial quantities along the trajectory.

> Np:=30; T:=0.5;

> Lr:=[seq([X(T/Np*i),Y(T/Np*i)],i=1..Np)]:

> with(plots):

Warning, the name changecoords has been redefined

> PLT:=plot(Lr,style=point): display(PLT,scaling=constrained);

> VX:=subs(sol,diff(x(t),t)): VY:=subs(sol,diff(y(t),t)):

> with(plottools):

> scale_v:=0.1: scale_F:=300: # trial and error are required to adjust these

> for i from 1 to Np do: t_i:=i*T/Np:

> PLv[i] := arrow([X(t_i),Y(t_i)], evalm(scale_v*vector([VX(t_i),VY(t_i)])), .05, .2, .1, color=green):
od:

> for i from 1 to Np do: t_i:=i*T/Np:

> PLF[i] := arrow([X(t_i),Y(t_i)], evalm(scale_F*vector([FX(X(t_i),Y(t_i)),FY(X(t_i),Y(t_i))])), .05, .2, .1, color=red):
od:

> for i from 1 to Np do: t_i:=i*T/Np:

> PLr[i] := arrow([0,0],vector([X(t_i),Y(t_i)]), .05, .2, .1, color=blue):
od:

> display(seq(display(PLT,PLv[i],PLF[i],PLr[i]),i=1..Np),insequence=true,scaling=constrained,view=[-1..1,-1..1]);

Let us demonstrate what happens to the potential and kinetic energies as a function of time:

> T_kin:=t->m/2*(VX(t)^2+VY(t)^2);

> T_kin(0),T_kin(0.25),T_kin(0.5);

> V_pot:=t->subs(x=X(t),y=Y(t),V);

> E_tot:=t->T_kin(t)+V_pot(t);

> E_tot(0),E_tot(0.25),E_tot(0.5);

> plot(['T_kin(t)','V_pot(t)','E_tot(t)'],t=0..1,color=[red,blue,green],thickness=3,axes=boxed);

We observe that the total energy is conserved, and that the kinetic and potential energies as a function of time have to be related in some way to achieve this.

When the planet is close to the center of attraction it moves faster (has higher kinetic energy) on account of the fact that the potential energy has been lowered.

One can establish a relationship between the average values over a period (this is called the virial theorem in classical mechanics, it's details depend on the force law, and it has a counterpart in quantum mechanics).

We lift an approximate value for the period from the graph:

> Period:=0.5;

> T_avg:=evalf(Int(T_kin,0..Period,5,_NCrule))/Period;

> E_avg:=evalf(Int(E_tot,0..Period,5,_NCrule))/Period;

> V_avg:=evalf(Int(V_pot,0..Period,5,_NCrule))/Period;

> V_avg/E_avg;

Exercise 3:

Carry out these calculations for different initial conditions and verify that the relationship between these average values is a feature of Newton's law of gravitation

(inverse-squared force law). Ensure a correct choice of period by looking at the graph. Explore what happens for slightly different power laws,

such as n =12/10, or n =8/10, etc.

>

2) The radial motion problem

The radial problem for r ( t ) is derived from angular momentum conservation. It allows us to determine the trajectory from a single differential equation.

First we calculate the angular momentum:

> with(linalg):

Warning, the protected names norm and trace have been redefined and unprotected

> L0:=evalm(m*crossprod([x0,y0,0],[vx0,vy0,0]));

> t0:=2;

> evalm(m*crossprod([X(t0),Y(t0),0],[VX(t0),VY(t0),0]));

We observe that the calculated trajectory does observe the conservation of angular momentum approximately. This is a test if the numerical accuracy of the algorithm. We also see that the angular momentum vector is perpendicular to the x - y plane.

The Euclidean norm is called the Frobenius norm in the linalg package (square-root of sum of elements squared).

> L_mag:=norm(L0,frobenius);

The one-dimensional problem of solving for the distance from the center of attraction as a function of time:

> Vr:=-GM*m/r^n;

> Veff:=Vr+L_mag^2/(2*m*r^2);

> Feff:=unapply(-diff(Veff,r),r);

> DE:=m*diff(r(t),t$2)=Feff(r(t));

> r0:=sqrt(x0^2+y0^2);

> vr0:=0; # this is for our original choice of perpendicular position and velocity vectors.

> IC_r:=r(0)=r0,D(r)(0)=vr0:

> #sol_r:=dsolve({DE,IC_r},r(t)); # doesn't work well

> sol_r:=dsolve(DE,r(t));

This is not very useful!

> sol_r:=dsolve({DE,IC_r},r(t),numeric,output=listprocedure):

> R:=subs(sol_r,r(t)):

> plot('R(t)',t=0..2,thickness=3);

We have almost balanced the centrifugal against the gravitational forces in the case of Earth's orbit.

For the other initial condition a substantial variation in radial distance results.

We proceed with a discussion of the stability of circular orbits for different power-law potentials.

We choose simpler units by picking an attractive potential of form -k/r^n, and use a mass of unit 1 to move in this force field

>

> restart;

> n:=3/2; # try this with n=3 and observe; then try n=2.0001; also observe what happens for n=1.999x...

The stability condition: the centripetal acceleration needed for uniform circular motion is provided by the attractive force.

> STC:=0=Lsq/(m*r0^3)-n*k/r0^(n+1);

> Lsq:=solve(STC,Lsq);

> V_eff:=Lsq/(2*m*r^2)-k/r^n;

> F_eff:=unapply(diff(-V_eff,r),r);

> NE:=m*diff(r(t),t$2)=F_eff(r(t));

> m:=1; k:=1;

> r0:=1;

> IC:=r(0)=r0,D(r)(0)=-1/100;

> sol:=dsolve({NE,IC},r(t),numeric,output=listprocedure):

> R:=subs(sol,r(t)):

> plot('R(t)',t=0..100,thickness=3);

A harmonic solution which drives the particle back towards oscillations about the circle with radius r0, when circular orbits are stable. Otherwise we get runaway solutions.

There is a simple graphical answer to the question as to why the problem of circular orbits becomes unstable:

> P1:=plot(1/(2*r^2)-1/r,r=0..10,color=red):

> P2:=plot(1/(2*r^2)-1/r^3,r=0..10,color=blue):

> plots[display](P1,P2,view=[0..5,-2..1],thickness=3,title="Effective Potential");

The blue curve corresponds to the unstable case. We are also perhaps interested in a graph of the radial force.

> P3:=plot(-diff(1/(2*r^2)-1/r,r),r=0..10,color=red,numpoints=300):

> P4:=plot(-diff(1/(2*r^2)-1/r^3,r),r=0..10,color=blue,numpoints=300):

> plots[display](P3,P4,view=[0..5,-1..0.5],thickness=3,title="Effective Force");

The red curve represents a restoring force with equilibrium at r=1.

Another situation that works, i.e., which has stable circular orbits is the case of confining (linear or higher-power) potentials in three dimensions:

> P1:=plot(1/(2*r^2)+r,r=0..10,color=red):

> P2:=plot(1/(2*r^2)+r^2,r=0..10,color=blue):

> plots[display](P1,P2,view=[0..5,0..10],thickness=3,title="Effective Potential");

> P3:=plot(-diff(1/(2*r^2)+r,r),r=0..10,color=red,numpoints=300):

> P4:=plot(-diff(1/(2*r^2)+r^2,r),r=0..10,color=blue,numpoints=300):

> plots[display](P3,P4,view=[0..5,-10..10],thickness=3,title="Effective Force");

3) The scattering problem

Let us return to the Kepler problem, and consider positive-energy solutions. These correspond to orbits that pass by the center of attraction once only (true comets). These solutions are also relevant for charged particle scattering (Rutherford scattering).

Suppose a particle is very far away from the scattering center so that we can ignore the attraction. The 'free' particle is characterized by:

the asymptotic velocity vector, and associated kinetic energy;

position vector, which together with the velocity vector and the particle mass determines the angular momentum vector.

To illustrate the idea we use Maple to generate some numerical solutions.

> restart; with(linalg):

Warning, the protected names norm and trace have been redefined and unprotected

> m:=1: r0:=[-100,5,0]: v0:=[1,0,0]:

> L0:=evalm(m*crossprod(r0,v0));

> Lsq:=dotprod(L0,L0);

> n:=1;

> V:=-1/r^n;

> V:=subs(r=sqrt(x^2+y^2),V);

> T0:=m/2*dotprod(v0,v0);

> E0:=evalf(T0+subs(x=r0[1],y=r0[2],V));

We observe that the particle is not entirely 'free' at a distance of 100 units. In fact, the inverse-distance-squared force law does not fall quickly enough at large distances to really allow the particle to be ever truly free at finite separations.

> F_x:=unapply(diff(-V,x),x,y);