# Calculating perihelia using Maple

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);
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?

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(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(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;
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"