# A period function for anharmonic oscillations

Posted 1 year ago
2744 Views
|
4 Replies
|
12 Total Likes
|
 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. Define Ansatz for polar coordinate Phase Space Trajectory RExp[n_] := Expand[b Plus[R[0], Total[b^# R[#] & /@ Range[n]]]] 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}]] 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]] List ring generators RingGens[n_] := Times @@ (v /@ #) & /@ (IntegerPartitions[n] /. x_Integer :> x + 2) 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]}] 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)}] Preload function values RSet[0] = 1; Set[RSet[#], Expand@RCalc[#]] & /@ Range[2*10]; 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.
4 Replies
Sort By:
Posted 1 year ago
 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 1 year ago
 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 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 ] 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 1 year ago
 - 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 1 year ago
 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":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.
Community posts can be styled and formatted using the Markdown syntax.