Poincare Sections for the Walker-Ford 2-2 resonance Hamiltonian

 

>    ham:=-E+J1+J2-J1^2-3*J1*J2+J2^2+alpha*J1*J2*C;

ham := -E+J1+J2-J1^2-3*J1*J2+J2^2+alpha*J1*J2*C

>   

Here C stands for cos2(phi1-phi2) ; below I have inserted the value of that in terms of qi's and pi's:

>    hama:=J1+J2-J1^2-3*J1*J2+J2^2+alpha*(-J1*J2+(q1*q2+p1*p2)^2/2);

hama := J1+J2-J1^2-3*J1*J2+J2^2+alpha*(-J1*J2+1/2*(q1*q2+p1*p2)^2)

and below here I insert the values for the Ji's, to get it altogether in terms of the 4 coord's in phasespace

>    ha:=map(factor,collect(subs(J1=(p1^2+q1^2)/2,J2=(p2^2+q2^2)/2,hama),[p1,p2,q1,q2]));

ha := -1/4*p1^4-1/4*(3*p2^2-p2^2*alpha+2*q1^2-2+3*q2^2+alpha*q2^2)*p1^2+alpha*q1*q2*p2*p1+1/4*p2^4-1/4*(-2+3*q1^2+alpha*q1^2-2*q2^2)*p2^2-1/4*q1^4+1/4*(-3*q2^2+2+alpha*q2^2)*q1^2+1/2*q2^2+1/4*q2^4
ha := -1/4*p1^4-1/4*(3*p2^2-p2^2*alpha+2*q1^2-2+3*q2^2+alpha*q2^2)*p1^2+alpha*q1*q2*p2*p1+1/4*p2^4-1/4*(-2+3*q1^2+alpha*q1^2-2*q2^2)*p2^2-1/4*q1^4+1/4*(-3*q2^2+2+alpha*q2^2)*q1^2+1/2*q2^2+1/4*q2^4

>    H1:=map(factor,collect(subs(J1=(p1^2+q1^2)/2,J2=(p2^2+q2^2)/2,alpha=1,hama),[p1,p2,q1,q2]));

H1 := -1/4*p1^4-1/2*(-1+p2^2+2*q2^2+q1^2)*p1^2+q1*q2*p2*p1+1/4*p2^4-1/2*(-1-q2^2+2*q1^2)*p2^2-1/4*q1^4-1/2*(q2-1)*(q2+1)*q1^2+1/2*q2^2+1/4*q2^4
H1 := -1/4*p1^4-1/2*(-1+p2^2+2*q2^2+q1^2)*p1^2+q1*q2*p2*p1+1/4*p2^4-1/2*(-1-q2^2+2*q1^2)*p2^2-1/4*q1^4-1/2*(q2-1)*(q2+1)*q1^2+1/2*q2^2+1/4*q2^4

>    with(plots): with(DEtools):

>    hamilton_eqs(H1);

[diff(p1(t),t) = q1(t)*p1(t)^2-q2(t)*p2(t)*p1(t)+2*q1(t)*p2(t)^2+q1(t)^3+(q2(t)-1)*(q2(t)+1)*q1(t), diff(p2(t),t) = 2*q2(t)*p1(t)^2-q1(t)*p2(t)*p1(t)-q2(t)*p2(t)^2+1/2*(q2(t)+1)*q1(t)^2+1/2*(q2(t)-1)*q...
[diff(p1(t),t) = q1(t)*p1(t)^2-q2(t)*p2(t)*p1(t)+2*q1(t)*p2(t)^2+q1(t)^3+(q2(t)-1)*(q2(t)+1)*q1(t), diff(p2(t),t) = 2*q2(t)*p1(t)^2-q1(t)*p2(t)*p1(t)-q2(t)*p2(t)^2+1/2*(q2(t)+1)*q1(t)^2+1/2*(q2(t)-1)*q...
[diff(p1(t),t) = q1(t)*p1(t)^2-q2(t)*p2(t)*p1(t)+2*q1(t)*p2(t)^2+q1(t)^3+(q2(t)-1)*(q2(t)+1)*q1(t), diff(p2(t),t) = 2*q2(t)*p1(t)^2-q1(t)*p2(t)*p1(t)-q2(t)*p2(t)^2+1/2*(q2(t)+1)*q1(t)^2+1/2*(q2(t)-1)*q...
[diff(p1(t),t) = q1(t)*p1(t)^2-q2(t)*p2(t)*p1(t)+2*q1(t)*p2(t)^2+q1(t)^3+(q2(t)-1)*(q2(t)+1)*q1(t), diff(p2(t),t) = 2*q2(t)*p1(t)^2-q1(t)*p2(t)*p1(t)-q2(t)*p2(t)^2+1/2*(q2(t)+1)*q1(t)^2+1/2*(q2(t)-1)*q...

>    ic:=generate_ic(H1,{t=0,q2=0,q1=0,p2=0.816,energy=0.1},1);

ic := {[0., 1.246378248, .816, 0., 0.]}

>    P:=poincare(H1,t=0..500,ic,stepsize=0.05,iterations=3,scene=[q1=-2..2,p1=-2..2,p2=-2..2],3):

_____________________________________________________

`H = .99999999e-1`, `  Initial conditions:`, t = 0., p1 = 1.246378248, p2 = .816, q1 = 0., q2 = 0.

`Maximum H deviation : .4956100000e-3 %`

_____________________________________________________

`Time consumed: 123 seconds`

>    display(P,orientation=[-98,70],tickmarks=[3,3,3]);

[Maple Plot]

>    poincare(H1,t=0..1000,ic,stepsize=0.05,iterations=3,scene=[q2,p2]);

_____________________________________________________

`H = .99999999e-1`, `  Initial conditions:`, t = 0., p1 = 1.246378248, p2 = .816, q1 = 0., q2 = 0.

`Number of points found crossing the (q2,p2) plane: 532`

`Maximum H deviation : .9408000000e-3 %`

_____________________________________________________

`Time consumed: 25 seconds`

[Maple Plot]

>    ic:=generate_ic(H1,{t=0,q2=0,q1=0,p2=0.816,energy=0.2},1);

ic := {[0., 1.173288828, .816, 0., 0.]}

>    P:=poincare(H1,t=0..500,ic,stepsize=0.05,iterations=3,scene=[q1=-2..2,p1=-2..2,p2=-2..2],3):

_____________________________________________________

`H = .20000000`, `  Initial conditions:`, t = 0., p1 = 1.173288828, p2 = .816, q1 = 0., q2 = 0.

`Maximum H deviation : .2023200000e-3 %`

_____________________________________________________

`Time consumed: 135 seconds`

>    display(P,orientation=[-98,70],tickmarks=[3,3,3]);

[Maple Plot]

>    poincare(H1,t=0..1000,ic,stepsize=0.05,iterations=3,scene=[q1,p1]);

_____________________________________________________

`H = .20000000`, `  Initial conditions:`, t = 0., p1 = 1.173288828, p2 = .816, q1 = 0., q2 = 0.

`Number of points found crossing the (q1,p1) plane: 119`

`Maximum H deviation : .3907200000e-3 %`

_____________________________________________________

`Time consumed: 20 seconds`

[Maple Plot]

>    poincare(H1,t=0..1000,ic,stepsize=0.05,iterations=3,scene=[q1,p2]);

_____________________________________________________

`H = .20000000`, `  Initial conditions:`, t = 0., p1 = 1.173288828, p2 = .816, q1 = 0., q2 = 0.

`Number of points found crossing the (q1,p2) plane: 119`

`Maximum H deviation : .3907200000e-3 %`

_____________________________________________________

`Time consumed: 20 seconds`

[Maple Plot]

>    poincare(H1,t=0..1000,ic,stepsize=0.05,iterations=3,scene=[p1,p2]);

_____________________________________________________

`H = .20000000`, `  Initial conditions:`, t = 0., p1 = 1.173288828, p2 = .816, q1 = 0., q2 = 0.

`Number of points found crossing the (p1,p2) plane: 471`

`Maximum H deviation : .3925800000e-3 %`

_____________________________________________________

`Time consumed: 21 seconds`

[Maple Plot]

Let's try harder to find a truly periodic one.  I hunt with alpha=0.1

>    ha;

-1/4*p1^4-1/4*(3*p2^2-p2^2*alpha+2*q1^2-2+3*q2^2+alpha*q2^2)*p1^2+alpha*q1*q2*p2*p1+1/4*p2^4-1/4*(-2+3*q1^2+alpha*q1^2-2*q2^2)*p2^2-1/4*q1^4+1/4*(-3*q2^2+2+alpha*q2^2)*q1^2+1/2*q2^2+1/4*q2^4
-1/4*p1^4-1/4*(3*p2^2-p2^2*alpha+2*q1^2-2+3*q2^2+alpha*q2^2)*p1^2+alpha*q1*q2*p2*p1+1/4*p2^4-1/4*(-2+3*q1^2+alpha*q1^2-2*q2^2)*p2^2-1/4*q1^4+1/4*(-3*q2^2+2+alpha*q2^2)*q1^2+1/2*q2^2+1/4*q2^4

>    hap1:=subs(alpha=0.1,ha):

Now I want e=.1, alpha=.1,p10=0,q20=0,q10=1.144885896 and p20=.5317149904, numbers which were hopefully

acquired so as to do this.  They correspond to ii=.7967422731, and were calculated off in some other sheet, from the full,

perturbed Hamiltonian.

>    icp1:=generate_ic(hap1,{t=0,q1=1.144885896,q2=0,p2=.5317149909,energy=0.1},1);

icp1 := {[0., .4407416253e-4, .5317149909, 1.144885896, 0.]}

>    P:=poincare(hap1,t=0..500,icp1,stepsize=0.05,iterations=3,scene=[q1=-2..2,p1=-2..2,p2=-2..2],3):

_____________________________________________________

`H = .10000000`, `  Initial conditions:`, t = 0., p1 = .4407416253e-4, p2 = .5317149909, q1 = 1.144885896, q2 = 0.

`Maximum H deviation : .6790000000e-5 %`

_____________________________________________________

`Time consumed: 111 seconds`

>    display(P,orientation=[-98,70],tickmarks=[3,3,3]);

[Maple Plot]

I'm ecstatic; it actually worked!

    Let's look at some actual sections.

>    poincare(hap1,t=0..1000,icp1,stepsize=0.05,iterations=3,scene=[q1,p1]);

_____________________________________________________

`H = .10000000`, `  Initial conditions:`, t = 0., p1 = .4407416253e-4, p2 = .5317149909, q1 = 1.144885896, q2 = 0.

`Number of points found crossing the (q1,p1) plane: 239`

`Maximum H deviation : .1352000000e-4 %`

_____________________________________________________

`Time consumed: 16 seconds`

[Maple Plot]

>    poincare(hap1,t=0..1000,icp1,stepsize=0.05,iterations=3,scene=[p2,p1]);

_____________________________________________________

`H = .10000000`, `  Initial conditions:`, t = 0., p1 = .4407416253e-4, p2 = .5317149909, q1 = 1.144885896, q2 = 0.

`Number of points found crossing the (p2,p1) plane: 238`

`Maximum H deviation : .1349000000e-4 %`

_____________________________________________________

`Time consumed: 16 seconds`

[Maple Plot]

>    poincare(hap1,t=0..1000,icp1,stepsize=0.05,iterations=3,scene=[q1,p2]);

_____________________________________________________

`H = .10000000`, `  Initial conditions:`, t = 0., p1 = .4407416253e-4, p2 = .5317149909, q1 = 1.144885896, q2 = 0.

`Number of points found crossing the (q1,p2) plane: 239`

`Maximum H deviation : .1352000000e-4 %`

_____________________________________________________

`Time consumed: 16 seconds`

[Maple Plot]

>