El hamiltoniano en las variables p y x

H =  ((p^2)/(2m))+ ((m ((ω0)^2))/2) (x^2)+ ϵ (p^4)

H0 =  ((p^2)/(2m))+ ((m ((ω0)^2))/2) (x^2)

H1 =  (p^4)

Las transformaciónes AA para el oscilador lineal y el hamiltoniano H1 escrito en esas variables:

x[θ_, J_]:= (Sqrt[((2 J)/(m ω0))])Sin[θ];
p[θ_, J_]:=(Sqrt[2 J m ω0])Cos[θ];

H1[θ_, J_]:= ((p[θ, J])^4);

El retrato de fases:

ϵpr = -0.1;
ContourPlot[((p^2)/2)+((x^2)/2)+ ϵpr (p^4), {x, -2, 2}, {p, -2,2}, ContourShading->None, Contours->24]

Cálculo del valor medio de H1

Integrate[(1/(2π))H1[θ, i], {θ, 0, 2π}]

H1medio[i_]:= (3/2) (i^2) (m^2)  ((ω0)^2);

Cálculo de la corrección a la frecuencia:

D[H1medio[i], i]

ω[i_]:= ω0 + 3 i (m^2) ϵ ((ω0)^2);

Cálculo de la corrección a la transformación canónica  que relaciona  las variables AA de H0 a las de H1, llamadas aquí θ, J y ϕ, i, respectivamente

Integrate[(1/(ω0))(H1medio[i] - H1[θp, i]), {θp, 0, θ}]

y1[θ_, i_]:= -(1/8) (i^2) (m^2)  ω0 (8 Sin[2 θ]+Sin[4 θ]);

Las ecuaciones de transformación entre las nuevas variables de AA y las antiguas

Jota[ϕ_, i_]:= i + ϵ Derivative[1,0][y1][ϕ, i];
Theta[ϕ_, i_]:= ϕ - ϵ Derivative[0,1][y1][ϕ, i];

Expansión hasta orden ϵ de las variables originales  x y p como funciones del tiempo

(Sqrt[((2 i)/(m ω0))])Collect[(1/((Sqrt[2]) (Sqrt[(i/(m ω0))])))Normal[Series[x[Theta[ϕ, i], Jota[ϕ, i]], {ϵ, 0, 1}]], ϵ, Simplify]
(Sqrt[2 i m ω0])Collect[(1/(Sqrt[2 i m ω0]))Normal[Series[p[Theta[ϕ, i], Jota[ϕ, i]], {ϵ, 0, 1}]], ϵ, Simplify]

xϵ[ϕ_, i_, ϵ_]:= (Sqrt[2]) (Sqrt[(i/(m ω0))]) (Sin[ϕ]+(3/4) i (m^2) ϵ ω0 (3+2 Cos[2 ϕ]) Sin[ϕ]);
pϵ[ϕ_, i_, ϵ_]:= (Sqrt[2]) (Sqrt[i m ω0]) (Cos[ϕ]+(1/4) i (m^2) ϵ ω0 (-6 Cos[ϕ]+Cos[3 ϕ]));

Cálculo de las constantes i y ϕ0 a orden ϵ a partir de las condiciones iniciales (sería largo de explicar)

Collect[Normal[Series[((m ((ω0)^2))/2)((xϵ[ϕ0, i, ϵ])^2)+ (1/(2m))((pϵ[ϕ0, i, ϵ])^2), {ϵ, 0, 1}]], ϵ, Simplify]

PowerExpand[Collect[Normal[Series[PowerExpand[ArcTan[m ω0((xϵ[ϕ0, i, ϵ])/(pϵ[ϕ0, i, ϵ]))]], {ϵ, 0, 1}]], ϵ, FullSimplify]]

Phi0 = ArcTan[m ω0 (x0/p0)] + ϵ Phi1;
I0 =(1/(ω0))( ((m ((ω0)^2))/2)(x0^2)+ (1/(2m))(p0^2)) + ϵ I1;

exp1 = Collect[Normal[Series[I0 ω0-(1/2) (I0^2) (m^2) ϵ ((ω0)^2) (4 Cos[2 Phi0]+Cos[4 Phi0]), {ϵ, 0, 1}]], ϵ, FullSimplify]
exp2 = Collect[TrigExpand[Collect[Normal[Series[Phi0+I0 (m^2) ϵ ω0 Cos[Phi0] (4+Cos[2 Phi0]) Sin[Phi0], {ϵ, 0, 1}]], ϵ, FullSimplify]], ϵ, FullSimplify]

sol = FullSimplify[Solve[{Coefficient[exp1, ϵ] == 0, Coefficient[exp2, ϵ] == 0}, {I1, Phi1}]][[1]]

Las condiciones iniciales para i y ϕ a partir de las condiciones iniciales para x y p, más la condición inicial para p si el dato es la velocidad v.

ϕ0[x0_, p0_, ϵ_] := ArcTan[m ω0 (x0/p0)] +  ϵ (1/2) (m^2) p0 x0 ω0 (-3-((2 (p0^2))/((p0^2)+(m^2) (x0^2) ((ω0)^2)))) ;
i0[x0_, p0_, ϵ_]:=(1/(ω0))( ((m ((ω0)^2))/2)(x0^2)+ (1/(2m))(p0^2) )+ϵ  ((5 (p0^4)-6 (m^2) (p0^2) (x0^2) ((ω0)^2)-3 (m^4) (x0^4) ((ω0)^4))/(8 ω0));
P0[v_, ϵ_]:= m v - 4 ϵ (m^4)(v^3);

Verifiaciones milagrosas:

PowerExpand[FullSimplify[Collect[Normal[Series[xϵ[ϕ0[x0, p0, ϵ], i0[x0, p0, ϵ], ϵ], {ϵ, 0,1}]], ϵ, TrigExpand], {x0 > 0, ω0> 0, p0>0, m> 0}]]

PowerExpand[FullSimplify[Collect[Normal[Series[pϵ[ϕ0[x0, p0, ϵ], i0[x0, p0, ϵ], ϵ], {ϵ, 0,1}]], ϵ, TrigExpand], {x0 > 0, ω0> 0, p0>0, m> 0}]]

A otra cosa

ω0 = m = 1;

Comparación de las soluciones a orden ϵ con el resultado de integrar las ecuaciones canónicas numéricamente y con la órbita no perturbada del oscilador lineal unidimensional.

x0 = 0;
v0 = 1;
ϵpr =Rationalize[-0.03]; (* ϵ es una cantidad negativa, para conservar el significado físico del hamiltoniano dada en el problema*)
ωpr = ω0 + 3 i0[x0, p0, 0] (m^2) ϵpr ((ω0)^2) ;

p0v = P0[v0, ϵpr];
p0 = 1;   (* se comparan trayectorias con la misma condición inciaial en el impulso. Sise elige p0 = p0v, se comparan trayectorias con la misma condición inicial en la velocidad*)

Phi0 = ϕ0[x0, p0, ϵpr];
I0 = i0[x0, p0, ϵpr];

Module[{x, p, t}, 
sol = NDSolve[{x'[t] == p[t] + 4 ϵpr ((p[t])^3), p'[t] ==- x[t],x[0] == x0, p[0] == p0}, {x, p}, {t, 0, 10 2(π/(ωpr))}][[1]]];
xsol = sol[[1]][[2]];
psol = sol[[2]][[2]];

g1 = ParametricPlot[{xϵ[ω0 t + ϕ0[x0, p0, 0], i0[x0, p0, 0], 0], pϵ[ω0 t + ϕ0[x0, p0, 0], i0[x0, p0, 0], 0]}, {t, 0, 2((2π)/(ω0))}, AspectRatio->Automatic, PlotStyle->Orange]; (* Naranja: órbita para el oscilador lineal sin perturbar*)

g2 = ParametricPlot[{xϵ[ωpr t + Phi0, I0, ϵpr], pϵ[ωpr t + Phi0, I0, ϵpr]}, {t, 0,2 ((2π)/(ωpr))}, AspectRatio->Automatic, PlotStyle->{Opacity[0.9],AbsoluteThickness[7], Blue}]; (*Azul: solución a orden ϵ*)

g3 = ParametricPlot[{xsol[t], psol[t]}, {t, 0,2 ((2π)/(ωpr))}, AspectRatio->Automatic, PlotStyle->{Yellow,Opacity[0.8], AbsoluteThickness[2]}]; (*Amarello: solución numérica de las ecs. completas*)

Show[g1, g2, g3, PlotRange-> All, AspectRatio-> Automatic, Background-> Black, ImageSize-> 600]

Comparación fina entre las solución numérica y la perturbativa: izq. las dos funciones en el mismo gráfico, der., su diferencia.

Posición:

Show[#, ImageSize-> 400]&/@ {Plot[{xsol[t],xϵ[ωpr t + Phi0, I0, ϵpr]} , {t, 0, 10π}],Plot[{xsol[t]-xϵ[ωpr t + Phi0, I0, ϵpr]} , {t, 0, 10π}]}

Impulso:

Show[#, ImageSize-> 400]&/@ {Plot[{psol[t],pϵ[ωpr t + Phi0, I0, ϵpr]} , {t, 0, 10 π}], Plot[{psol[t]-pϵ[ωpr t + Phi0, I0, ϵpr]} , {t, 0, 10 π}]}