This worksheet is designed to create plots of the trajectories of test masses

orbiting in precessing elliptical orbits in the Schwarzschild metric

>

restart:

>	with(plots):

Warning, the name changecoords has been redefined

Conservation of angular momentum:

>	Leq:=diff(phi(tau),tau)=B/r(tau)^2;

Conservation of energy:

>	Eeq:=diff(t(tau),tau) = A/(1-2*m/r(tau));

>	H:=1-2*m/r(tau);

>	Hprime:=diff(1-2*m/r,r);

>	Hprimeeq:=subs(r=r(tau),Hprime);

Those are the equations coming from the metric and the conservation laws, i.e., time independence and spherical symmetry.

    Of course they are using spherical symmetry to draw everything in the equatorial plane.

 The equation below is the geodesic equation for r

>	r2eq:=diff(r(tau),tau$2)-B^2*H/r(tau)^3+(1/2)*(Hprimeeq/H)*(A^2-diff(r(tau),tau)^2);

 Then I want to study the energy equations as well, which come from the normalizaton of the 4-velocity.

>	Veff:=-(m/r(tau))+B^2/(2*r(tau)^2)-m*B^2/r(tau)^3;

>	r1eq:=diff(r(tau),tau)^2/2 + Veff=(A^2-1)/2;

>	ellipse:=subs({diff(r(tau),tau)=0},r1eq);

Now let me just get rid of the mass by setting it to the mass of the sun, as measured in kilometers:

>

m:=1.477;

  Next we state the values of the perihelion and aphelion distances of Mercury from the sun, again in kilometers:

>	rpm:=46*10^6; ram:=69.82*10^6;

We introduce the speed of light here, in km/sec, but it may well not be of any use.

>	c:=2.99792458*10^5;

We may then insert those values in the equation for the effective potential, as turn around points,

to determine the associated values of the constants.

>	ellperi:=subs(r(tau)=rpm,ellipse);

>	ellap:=subs(r(tau)=ram,ellipse);

>	pars:=[solve({ellperi,ellap},{A,B})];

>	Merc:=[op(pars[4])];

>	ellipsewab:=subs(Merc,ellipse);

>	testsol:=[solve(ellipsewab,{r(tau)})];

Notice that these numbers ought to have been our original ones, namely .46 x 10^8 and .6982 x 10^8.  They differ, already in the third place?

   Let's see if we can do any better?

>	Digits:=25;

>	ellperib:=subs(r(tau)=rpm,ellipse);ellapb:=subs(r(tau)=ram,ellipse);

>	parsb:=[solve({ellperib,ellapb},{A,B})];

>	Mercb:=parsb[4];

>	ellipsewabb:=subs(Mercb,ellipse);

>	testsb:=solve(ellipsewabb,{r(tau)});

AAAh Ha!  Now, it's actually doing the reasonable thing.

The first solution is then the correct values for the constants of the motion for Mercury.   We will insert them back into the equations:

  However, first I'm going to try some numbers for the constants that are somewhat simpler, which correspond to a much closer orbit.  

These will give us an orbit that "precesses" much more strongly.

>	B:=3.5695; m:=1; A:=0.9509661236;

>	Leq, r2eq;

   The three values below, then, also from the same book, tell me I'm going to start things up from rest at a turning point, where r(0)=11m.

>	ini:=r(0)=11, D(r)(0)=0, phi(0)=0;

>	Eq9:=dsolve({Leq,r2eq,ini},{r(tau),phi(tau)},numeric,output=listprocedure);

   The command below, where I'm printing out the results of a plot command in full, are mostly to show us how to create such a plot, and to

   interpret them.

>

>

>	PP1:=polarplot([rhs(Eq9(tau)[3]),rhs(Eq9(tau)[2]),tau=0..1],scaling=constrained);

>	nops(PP1), nops(op(1,PP1)), nops(op(1,op(1,PP1)));

We interpret all the numbers just above by saying that the dataset labeled as PP1 has 3 parts:  the first part is the complete PLOT routine, the

second part is the CONSTRAINED statement, which means that the vertical and horizontal axes should use the same scaling, and the third

is a statement about where the plot should be displayed.

   Then the plot statement itself has two parts:  they are the description of the CURVE to be plotted, and then the colours to use.

 Lastly the CURVE has some 49 points, which Maple thought was sufficient to generate the desired more-or-less smooth curve.

    Note that each of the 49 points is an x-coord, and a y-coord, which the "polarplot" command generated from the values of r and theta

           which were passed to it.

   I also show below what such a plot actually looks like.

>	polarplot([rhs(Eq9(tau)[3]),rhs(Eq9(tau)[2]),tau=0..300],scaling=constrained);

Notice the orbit is rather unexpected in its structure, since the "precession" is huge.  It took more than two turns around, in phi, in order to get back to its aphelion.  In order to see better how that works; let's make some other plots.

>

Eq9(100);

>	plot([rhs(Eq9(tau)[1]),rhs(Eq9(tau)[3]),tau=0..300]);

This shows the radial oscillations well, as a function of tau.  We see it started at 11, went down to its

 minimum of about 5.5, at tau about 128 or so, and then back up to 11, at tau about 260 or so.

>	plot([rhs(Eq9(tau)[1]),rhs(Eq9(tau)[3]),tau=254.8..255.2],numpoints=200);

We can see that numerical inaccuracy is setting in, since the orbit should not get larger than 11; nonetheless, the

 maximum seems to be at approx. 254.9

>	Eq9(254.9);

  We see that at this FIRST return to the aphelion, the angle is already up to 15.4, i.e., almost 5 times Pi.

    Note that this is entirely consistent with our graph above, where that many turns were made.

            Now I would like to generate an animated plot, so that we can see it "happening":

>	for n from 1 to 20 do PP||n:=polarplot([rhs(Eq9(tau)[3]),rhs(Eq9(tau)[2]),tau=0..50*n], scaling=constrained); end do:

>	display([seq(PP||n,n=1..20)],insequence=true);

>	A:='A'; B:='B';

Now, let's go back and look at the one for Mercury:

>	val:=[op(Mercb)];

>	B:=rhs(val[2]); A:=rhs(val[1]);m:=1.477;

>

Leq,r2eq;

>	iniMerc:=r(0)=rpm, D(r)(0)=0, phi(0)=0;

>	EqMerc:=dsolve({Leq,r2eq,iniMerc},{r(tau),phi(tau)},numeric,output=listprocedure);

 We can now use this to create a plot of the orbit:

>	polarplot([rhs(EqMerc(tau)[3]),rhs(EqMerc(tau)[2]),tau=0..10^(12)],scaling=constrained);

Well, that seems not to have gotten us even half of all the way around, yet.  Let's push on furthr.

>	polarplot([rhs(EqMerc(tau)[3]),rhs(EqMerc(tau)[2]),tau=0..2*10^(12)],scaling=constrained);

That does begin to look like an orbit, and now we're most of the way there.

I wonder what is the appropriate value?

>	polarplot([rhs(EqMerc(tau)[3]),rhs(EqMerc(tau)[2]),tau=0..2.25*10^(12)],scaling=constrained);

Aah!  2.25 doesn't quite close, but it turns out that 2.3 does close.  We want to know what is phi when it closes back, i.e.,

 when r gets back to 46 MegaKm.

   Let's do some radial plots:

>	plot([rhs(EqMerc(tau)[1]),rhs(EqMerc(tau)[3]),tau=2.27*10^(12)..2.29*10^(12)]);

>	plot([rhs(EqMerc(tau)[1]),rhs(EqMerc(tau)[3]),tau=2.277*10^(12)..2.28*10^(12)]);

 So it appears that 2.2784 x 10^(12) is the most likely value I'm going to get.

We can also look at it in the polar mode.

>	polarplot([rhs(EqMerc(tau)[3]),rhs(EqMerc(tau)[2]),tau=2.27*10^(12)..2.29*10^(12)]);

There we see the perihelion point---notice that I took out the constraint on the scaling, so that now the scaling is

very different in the two directions!  It looks to me as if it's further along than 2.28

>	polarplot([rhs(EqMerc(tau)[3]),rhs(EqMerc(tau)[2]),tau=2.2781*10^(12)..2.2785*10^(12)],numpoints=500);

You can see that, here, we are beyond the level of accuracy of the numerical processes being used.  

Therefore I will decide that anywhere in there is close enough.

>	EqMerc(2.2782*10^(12)); EqMerc(2.27827*10^(12)); EqMerc(2.27829*10^(12));EqMerc(2.27834*10^(12)); EqMerc(2.278345*10^(12)); EqMerc(2.2783455*10^(12));EqMerc(2.278346*10^(12));

We see that the derivative in the next to the last one is  negative and really very small, and it turns back positive at the next place.

Therefore it seems reasonable to go ahead and take  tau = 2.2783455 x 10^(12) as the closest we're going to get, and then

we have phi at this point is 6.283179890518743544231355, for whatever all those digits are worth.

We will subtract that from 2Pi and see what we get:

>	delt:=evalf(6.283179890518743544231355-2*Pi);

Multiply this by 180/Pi to change to degrees and then by 3600 to change to seconds of arc:

>	deltsec:=delt*3600*180/evalf(Pi);

So this is a precession (negative) of -1.117 seconds per year.

I also decide to determine what the value of tau is, for the one orbit, in years, rather than kilometers.

>	secyr:=365.25*86400;

>	2.2783455*10^(12)/c/secyr;

>	polarplot([rhs(EqMerc(tau)[3]),rhs(EqMerc(tau)[2]),tau=0..2.25*10^(13)],scaling=constrained);

>	polarplot([rhs(EqMerc(tau)[3]),rhs(EqMerc(tau)[2]),tau=0..2.25*10^(14)],scaling=constrained);

Presumably the graph above, for quite a large number of orbits---about 100---does allow one to actually see some

             of the precession of  the orbit.

>