2
$\begingroup$

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?

$\endgroup$
4
$\begingroup$

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;

enter image description here

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

enter image description here

# 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]);

enter image description here

# 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"
$\endgroup$

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service, privacy policy and cookie policy

Not the answer you're looking for? Browse other questions tagged or ask your own question.