Message Boards Message Boards

A period function for anharmonic oscillations

Brad Klee, Wolfram Research
Posted 8 years ago

I just updated to Mathematica 11 and celebrating with a short note on A276816, which I am busy editing to fix my overly-liberal use of "Coefficient". Here's a Mma11-working copy of my best discovery yet, a function for determining the period of a one dimensional anharmonic oscillation with arbitrary-valued potential expansion coefficients, roughly following the method of Plane Pendulum and Beyond.

  1. Define Ansatz for polar coordinate Phase Space Trajectory

    RExp[n_] := Expand[b Plus[R[0], Total[b^# R[#] & /@ Range[n]]]]

  2. Expansion function for calculating R[n]

    RCalc[n_] := With[{basis = Subtract[Tally[Join[Range[n + 2], #]][[All, 2]], 
   Table[1, {n + 2}]] & /@ IntegerPartitions[n + 2][[3 ;; -1]]}, 
       Total@ReplaceAll[ Times[-2, Multinomial @@ #, v[Total[#]],   
            Times @@ Power[RSet[# - 1] & /@ Range[n + 2], #]] & /@ 
                  basis, {Q^2 -> 1, v[2] -> 1/4}]]

  3. Calculate differential time dependence

    dt[n_] := With[{exp = Normal[Series[-1/(1 + x) /. 
       x -> Total[(2 # v[#] RExp[n - 1]^(# - 2) & /@  Range[3, n + 2])], {b, 0, n}]]}, 
          Expand@ReplaceAll[Coefficient[exp, b, #] & /@ Range[n], R -> RSet]]

  4. List ring generators

    RingGens[n_] := Times @@ (v /@ #) & /@ (IntegerPartitions[n] /. x_Integer :> x + 2)

  5. Construct coefficient triangle

    tri[m_] :=  MapThread[ Function[{a, b}, Times[-# /. v[n_] :> Q^n /. 
       Q^n_ :> Binomial[n, n/2], (1/2) Coefficient[a, #]] & /@ b],
         {dt[2 m][[2 #]] & /@ Range[m], RingGens[2 #] & /@ Range[m]}]

  6. Expand period from integer triangle

    PeriodExpansion[tri_, n_] := ReplaceAll[ 1 + Dot[MapThread[
     Dot, {tri, 2 RingGens[2 #] & /@ Range[n]}], (2 h)^( Range[n])], 
                {v[m_] :> (v[m]*(1/2)^m)}]

  7. Preload function values

    RSet[0] = 1; Set[RSet[#], Expand@RCalc[#]] & /@ Range[2*10];

  8. Transform potential to set of v[_] coefficients

    vRep[potential_, n_] := With[{coeffs = Coefficient[Series[potential, {x, 0, n}], x, #] & /@  Range[2, n]},
  v[#] -> coeffs[[# - 1]]/2/(coeffs[[1]] 2) & /@ Range[2, n]]

Now let's look at the period function up to order 5 in dimensionless energy

    PeriodExpansion[tri[5], 5]
    Out[]= 1 + 2 h (15 v[3]^2 - 3 v[4]) + 
     4 h^2 ((3465 v[3]^4)/4 - 945/2 v[3]^2 v[4] + (105 v[4]^2)/4 + 
        105/2 v[3] v[5] - (15 v[6])/4) + 
     8 h^3 ((255255 v[3]^6)/4 - 225225/4 v[3]^4 v[4] + 
        45045/4 v[3]^2 v[4]^2 - (1155 v[4]^3)/4 + 15015/2 v[3]^3 v[5] - 
        3465/2 v[3] v[4] v[5] + (315 v[5]^2)/8 - 3465/4 v[3]^2 v[6] + 
        315/4 v[4] v[6] + 315/4 v[3] v[7] - (35 v[8])/8) + 
     16 h^4 ((334639305 v[3]^8)/64 - 101846745/16 v[3]^6 v[4] + 
        72747675/32 v[3]^4 v[4]^2 - 3828825/16 v[3]^2 v[4]^3 + (
        225225 v[4]^4)/64 + 14549535/16 v[3]^5 v[5] - 
        3828825/8 v[3]^3 v[4] v[5] + 675675/16 v[3] v[4]^2 v[5] + 
        675675/32 v[3]^2 v[5]^2 - 45045/32 v[4] v[5]^2 - 
        3828825/32 v[3]^4 v[6] + 675675/16 v[3]^2 v[4] v[6] - 
        45045/32 v[4]^2 v[6] - 45045/16 v[3] v[5] v[6] + (3465 v[6]^2)/
        64 + 225225/16 v[3]^3 v[7] - 45045/16 v[3] v[4] v[7] + 
        3465/32 v[5] v[7] - 45045/32 v[3]^2 v[8] + 3465/32 v[4] v[8] + 
        3465/32 v[3] v[9] - (315 v[10])/64) + 
     32 h^5 ((29113619535 v[3]^10)/64 - 45176306175/64 v[3]^8 v[4] + 
        11712375675/32 v[3]^6 v[4]^2 - 2342475135/32 v[3]^4 v[4]^3 + 
        305540235/64 v[3]^2 v[4]^4 - (2909907 v[4]^5)/64 + 
        1673196525/16 v[3]^7 v[5] - 1405485081/16 v[3]^5 v[4] v[5] + 
        305540235/16 v[3]^3 v[4]^2 v[5] - 14549535/16 v[3] v[4]^3 v[5] + 
        305540235/64 v[3]^4 v[5]^2 - 43648605/32 v[3]^2 v[4] v[5]^2 + 
        2297295/64 v[4]^2 v[5]^2 + 765765/32 v[3] v[5]^3 - 
        468495027/32 v[3]^6 v[6] + 305540235/32 v[3]^4 v[4] v[6] - 
        43648605/32 v[3]^2 v[4]^2 v[6] + 765765/32 v[4]^3 v[6] - 
        14549535/16 v[3]^3 v[5] v[6] + 2297295/16 v[3] v[4] v[5] v[6] - 
        135135/64 v[5]^2 v[6] + 2297295/64 v[3]^2 v[6]^2 - 
        135135/64 v[4] v[6]^2 + 61108047/32 v[3]^5 v[7] - 
        14549535/16 v[3]^3 v[4] v[7] + 2297295/32 v[3] v[4]^2 v[7] + 
        2297295/32 v[3]^2 v[5] v[7] - 135135/32 v[4] v[5] v[7] - 
        135135/32 v[3] v[6] v[7] + (9009 v[7]^2)/128 - 
        14549535/64 v[3]^4 v[8] + 2297295/32 v[3]^2 v[4] v[8] - 
        135135/64 v[4]^2 v[8] - 135135/32 v[3] v[5] v[8] + 
        9009/64 v[6] v[8] + 765765/32 v[3]^3 v[9] - 
        135135/32 v[3] v[4] v[9] + 9009/64 v[5] v[9] - 
        135135/64 v[3]^2 v[10] + 9009/64 v[4] v[10] + 
        9009/64 v[3] v[11] - (693 v[12])/128)

With so many terms, each a possible hiding place for error, it's worthwhile to validate by testing on well-know potential functions. Let's try two with minima at x = 0, of pendulum and Coulomb type:

PotentialList = {1 - Cos[x], 1/(2 (x + 1)^2) - 1/(x + 1)}
(PeriodExpansion[tri[5], 5] /. vRep[#, 12]) & /@ PotentialList
Out[]={ 1 + (1/8)*h + (9/256)h^2 + (25/2048)h^3 + (1225/262144)h^4 + (3969/2097152)h^5 , 
         1 + 3 h + (15/2)h^2 + (35/2)h^3 + (315/8)h^4 + (693/8)h^5   }

The first expansion contains the leading terms of EllipticK(m), under the transformation h->2m. Two or three of these expansion coefficients can easily be measured in a simple experiment using a USB computer mouse .

The second expansion eventually converges like 

1 + 3 h + (15/2)h^2 + (35/2)h^3 + (315/8)h^4 + (693/8)h^5 + ... = (1 - 2 h)^(-3/2)

The converged function can be used to determine the exact period of a Kepler orbit. The oscillation of planets are well known in astronomy, but what's more, the method here develops and applies to precession as in this planetary precession demonstration.

The period of a pendulum and a planet in just one function... Working more on this as I finish up my PhD dissertation in the next one or two years.

Best Regards.
POSTED BY: Brad Klee
4 Replies
Sort By:

Hi Moderation Team,

Thanks for the plaudit! But if we're going to make this a featured post, let's go ahead and add an animation.

Performing standard numerical time evolution along the pendulum C.E.S. allows us to compute a "phase space race":

Phase Space Race

This image shows trajectories for six energies, with ten approximations at each energy level, using the same coloring scheme where precision increases toward blue.

Especially for the outer higher-energy trajectories, the poorly-approximated green curves race ahead of schedule, finishing before their due time. The blue curves linger a while longer around the critical points, and arrive at the finish line in more reasonable time.

Since the period is defined as the time for one complete cycle around phase space, the energy dependence of the period is very obvious in the animation.

POSTED BY: Brad Klee

enter image description here - you have earned "Featured Contributor" badge, congratulations !

This is a great post and it has been selected for the curated Staff Picks group. Your profile is now distinguished by a "Featured Contributor" badge and displayed on the "Featured Contributor" board.
POSTED BY: EDITORIAL BOARD

Hi Todd,

Same to you... I was reading through your ABC book, cool idea.

I agree that the code is difficult to follow without pictures, but remind you that the pictures are already out there in numerous references. And in the future we are likely to see more references, with pictures, better explanations, etc.

Actually the real difficulty is with the optimization, but this is also the main bragging point. Whereas only 10 orders of magnitude are given in the Wolfram Demonstration "Exact and Approximate Relativistic Corrections to the Orbital Precession of Mercury", in the above code we can easily compute 20 orders of magnitude.

If you want to see some pictures, we might as well make some.

  • Add phase angle dependence, and list trajectories through 20 orders of magnitude

    TrajectoryList = (RExp[#] /. R -> RSet /. v[n_] :> Q^n v[n]) & /@ 
Range[20];

  • User Coulomb potential expansion coefficients for v[ _ ]

    KeplerPhaseSpace =  MapThread[#1 /. #2 &, {TrajectoryList, 
vRep[PotentialList[[2]], 2 + #] & /@ Range[20]}];

  • List plot functions for 10 Conserved Energy Surface (CES) contours in phase space

    PlotFList = Flatten[KeplerPhaseSpace /. {Q -> Cos[\[Phi]], b -> #/15} & /@ Range[10]];

  • Plot the CES contours. Precision increases from green to blue.

    Show[MapThread[  PolarPlot[#1, {\[Phi], 0, 2 Pi}, PlotStyle -> #2] &, {PlotFList, 
Flatten[Table[  RGBColor[0, (20 - #)/20, #/20] & /@ Range[20], {10}]]}], 
PlotRange -> All, ImageSize -> 800]

CES Contours in Phase Space

  • CES Coordinates at order-20 precision

    f20 = {# Cos[\[Phi]], # Sin[\[Phi]], 1/2 b^2} &@(TrajectoryList[[20]] /. 
vRep[PotentialList[[2]], 22] /. Q -> Cos[\[Phi]]);

  • Place potential function into phase space with zero momentum

    Potential = {x, 0, 1/2 + 1/(2 (1 + x)^2) - 1/(1 + x)}

  • Plot CES over phase space with potential function as a tangent curve

    Show[  ParametricPlot3D[f20, {\[Phi], 0, 2 \[Pi]}, {b, 0, 10/15}, 
ImageSize -> 600, Axes -> False, Boxed -> False, PlotRange -> All],
ParametricPlot3D[Potential, {x, -1/2, 3}], ImageSize -> 1000 ]

CES Over phase space

These are just plots related to the "Rexp" function. With some additional development of the "dt" function, we could even animate motion of the oscillator along a CES.

Happy Holidays !
POSTED BY: Brad Klee

Brad,

Glad to see you are still producing good stuff. I remember your work on aperiodic tilings. Have to admit though it is hard to follow this without pictures. Unlike OEIS it is easy to attach images on Community.
POSTED BY: Todd Rowland

Group Abstract Group Abstract