Calculating perihelia using Maple

Question

I have the differential equation

DE:=$\frac{d^2}{d\phi^2}(\frac{1}{r(\phi)})+\frac{1}{r(\phi)}=1+\frac{3}{64r(\phi)^2}$,

with initial conditions

ICS:=$r(0)=\frac{2}{3}, (r)'(0)=0$.

Solving for $r(\phi)$ and plotting the orbit with polar coordinates

sol:=dsolve({DE,ICs},numeric,output=listprocedure); 
r_sol:=rhs(sol[2]);
polarplot(r_sol(phi),phi=0..10*Pi;

I observe an orbit about the origin that has precessing perihelion. I want to calculate, using Maple, how much the perhelion precesses by per revolution. How do I find the value of $\phi$ for each of the perihelia?

Answer 1

An efficient way to do this in Maple is to use the events option of dsolve(...,numeric) to recognize when r(phi) is at a local minimum.

A less efficient way (in general, though perhaps not so noticeable here) would be to use the Optimization:-Minimize command on a procedure returned by dsolve(...,numeric) for evaluating r(phi).

restart;

# We augment the system of differential equations with
# a new function ddr(phi) which equals the second derivative
# of r(phi) w.r.t. phi.

DE := ddr(phi) = diff(r(phi),phi,phi),
      diff(1/r(phi),phi,phi) + 1/r(phi) = 1 + 3/(64*r(phi)^2):
ICs := ddr(0)=0, r(0)=2/3, D(r)(0)=0:

# Use an event, which halts the solver whenever both these
# conditions attain:
#
# 1) diff(r(phi),phi) = 0  ie. r(phi) is at a local min/max
#
# 2) ddr(phi)>0            ie. second derivative test
#
# Hence the solving will halt when r(phi) is at a local minimum,
# which characterizes perihelion.

sol := dsolve( {DE,ICs}, numeric, output=listprocedure,
               events=[ [ [ diff(r(phi),phi), ddr(phi)>0 ], halt ] ] ):

# The following is a more robust way to pick off the procedures
# returned by `dsolve`, as it doesn't depend on using
# positions in the above list (as did your `op` approach).
# Relying on position is fragile when variable names change or
# new outputs (like `ddr(phi)` here) are introduced.

phi_sol:=eval(phi,sol):
r_sol:=eval(r(phi),sol):
dr_sol:=eval(diff(r(phi),phi),sol):
ddr_sol:=eval(ddr(phi),sol):

# Now we can see how the event works.
# Any of the procedures will stop at this value of phi.

phi_sol(1000.0);
    Warning, cannot evaluate the solution further right
         of 6.6194338, event #1 triggered a halt

                          6.61943388363658

phi_sol(last);

                          6.61943388363658

r_sol(last);

                          0.666666666609673

dr_sol(last);    # first derivative

                                          -18
                       9.67481032353862 10   

ddr_sol(last);   # second derivative

                          0.175347222241220

# If we clear the event then we can compute further (say,
# up to the next halt). See below, for that done in a loop.

# But first we'll go around a few revolutions, obtaining
# a polar plot of r(phi) versus phi.

n := 11;  # How many revolutions we'll do.

                               n := 11

# We temporarily disable the halting event, so that we can
# plot r(phi). Then we re-enable it.

r_sol(eventdisable={1});
P:=plots:-polarplot(r_sol(phi),phi=0..n*2*Pi):
r_sol(eventenable={1});

P;

# Now, with the event re-enabled, we'll go around for `n` revolutions
# once again, starting with the ICs. But this time it will halt,
# whenever the event conditions are satisfied.
#
# Each time it halts we'll store the phi value, clear the event,
# and proceed again.

oldwarnlevel:=interface(warnlevel=0): # Temporarily suppress the warning.

old:=0:
T:='T':  # Name of a table to store the perihelion points.
for i from 1 to n do
  new:=phi_sol(old+5*Pi/2); # Attempt to go more a more than 2*Pi around
  T[i]:=new;
  try
    phi_sol(eventclear);
  catch:
  end try;
  print(sprintf("Halted at %4.2f instead of %4.2f + %a = %4.2f",
                new, old, 2*Pi, evalf(old+2*Pi)));
  old:=new;
end do:

interface(warnlevel=oldwarnlevel):  # Re-instate the original warning level.

           "Halted at 6.62 instead of 0.00 + 2*Pi = 6.28"
          "Halted at 13.24 instead of 6.62 + 2*Pi = 12.90"
          "Halted at 19.86 instead of 13.24 + 2*Pi = 19.52"
          "Halted at 26.48 instead of 19.86 + 2*Pi = 26.14"
          "Halted at 33.10 instead of 26.48 + 2*Pi = 32.76"
          "Halted at 39.72 instead of 33.10 + 2*Pi = 39.38"
          "Halted at 46.34 instead of 39.72 + 2*Pi = 46.00"
          "Halted at 52.96 instead of 46.34 + 2*Pi = 52.62"
          "Halted at 59.57 instead of 52.96 + 2*Pi = 59.24"
          "Halted at 66.19 instead of 59.57 + 2*Pi = 65.86"
          "Halted at 72.81 instead of 66.19 + 2*Pi = 72.48"

# Compute r(phi) for each stored phi value.

phiL:=convert(T,list):
r_sol(eventdisable={1});
rL:=map(r_sol,phiL):
r_sol(eventenable={1});

evalf[4](phiL); # A quick look at them.

       [6.619, 13.24, 19.86, 26.48, 33.10, 39.72,
        46.34, 52.96, 59.57, 66.19, 72.81]

# Make a list of lists, for point-plotting.

pts:=[seq( [rL[i], phiL[i]], i=1..nops(phiL) )]:

evalf[4](pts); # A quick look at them.

   [[0.6667, 6.619], [0.6667, 13.24], [0.6667, 19.86],
    [0.6667, 26.48],[0.6667, 33.10], [0.6667, 39.72],
    [0.6667, 46.34], [0.6667, 52.96], [0.6667, 59.57],
    [0.6667, 66.19], [0.6667, 72.81]]

# Now we'll make a point-plot from each pair.

# Just use  color=blue  if your Maple version lacks ColorTools:
colorlist:=[seq( ColorTools:-Color([1-i,i,1.0]),
                 i=0.1..0.9, (0.9-0.1)/(nops(pts)-1) )]:

Pts := seq( plot([pts[i]], color=colorlist[i],
                 axiscoordinates=polar, coords=polar,
                 style=point,
                 symbolsize=20, symbol=solidcircle),
            i=1..nops(pts) ):
#Pts;

plots:-display(P, Pts);

# We could make several kinds of animation. Here's one:

r_sol(eventdisable={1});
A := plots:-animate( plots:-display,
                     [ 'plots:-polarplot'(r_sol(phi),phi=0..floor(ii)*2*Pi,
                                      color=black), Pts[floor(ii)] ],
                     ii=1..n, frames=3*n, paraminfo=false ):
r_sol(eventenable={1});

plots:-display(A, insequence=true, coordinateview=[0.0..1.7, 0..2*Pi]);

# A less efficient way (in general, though perhaps not so
# noticeable here) would be to apply the Optimization:-Minimize
# command to r(phi).

r_sol(eventdisable={1});
 old_peri := 0.0:
 for i from 1 to n do
   optsol := Optimization:-Minimize( r_sol, old_peri+Pi/2 .. old_peri+Pi/2+2*Pi );
   new_peri := optsol[2][1];
   print(sprintf("Minimum for r(phi) found at phi=%4.2f instead of %4.2f + %a = %4.2f",
                new_peri, old_peri, 2*Pi, evalf(old_peri+2*Pi)));
   old_peri := new_peri:
 end do:
r_sol(eventenable={1});

    "Minimum for r(phi) found at phi=6.62 instead of 0.00 + 2*Pi = 6.28"
   "Minimum for r(phi) found at phi=13.24 instead of 6.62 + 2*Pi = 12.90"
   "Minimum for r(phi) found at phi=19.86 instead of 13.24 + 2*Pi = 19.52"
   "Minimum for r(phi) found at phi=26.48 instead of 19.86 + 2*Pi = 26.14"
   "Minimum for r(phi) found at phi=33.10 instead of 26.48 + 2*Pi = 32.76"
   "Minimum for r(phi) found at phi=39.72 instead of 33.10 + 2*Pi = 39.38"
   "Minimum for r(phi) found at phi=46.34 instead of 39.72 + 2*Pi = 46.00"
   "Minimum for r(phi) found at phi=52.96 instead of 46.34 + 2*Pi = 52.62"
   "Minimum for r(phi) found at phi=59.57 instead of 52.96 + 2*Pi = 59.24"
   "Minimum for r(phi) found at phi=66.19 instead of 59.57 + 2*Pi = 65.86"
   "Minimum for r(phi) found at phi=72.81 instead of 66.19 + 2*Pi = 72.48"