Standard Map on a Torus 

Frank Wang, fwang@lagcc.cuny.edu 

Department of Mathematics, LaGuardia Community College 

The City University of New York 

Abstract:  Maple worksheet for producing a standard map on a torus is developed.  A brief introduction of Hamiltonian mechanics using a simple pendulum system as an example is provided, followed by a discussion of a kicked rotor system which consists of a simple pendulum with the potential energy turned on in delta-function pulses.  The integration of the kicked rotor problem yields the standard map, which has characteristics of a large class of systems.  Plotting the standard map on a torus facilitates a three-dimensional visualization the Kolmogorov-Arnold-Moser (KAM) theory of chaos in a Hamiltonian system.  An animation demonstrating chaos created by homoclinic tangle near a hyperbolic point is contained in this worksheet.    

Introduction 

Classically a Hamiltonian function for a particle is the sum of the kinetic and potential energies, 

. 

The dynamics of the system are governed by Hamiltonian's equations, [1,2] 

Simple Pendulum 

For a simple pendulum, the potential energy is .  The Hamiltonian equations give 

>

restart:

 

>	Eq1 := diff(p(t),t)=-sin(q(t));

 

>	Eq2 := diff(q(t),t)=p(t);

 

>	DEtools[DEplot]([Eq1, Eq2], [q(t), p(t)], t=-20..20, q(t)=-1.2*Pi..1.2*Pi, p(t)=-2.2..2.2, [[q(0)=Pi/10, p(0)=0], [q(0)=Pi/6, p(0)=0], [q(0)=Pi/4, p(0)=0], [q(0)=Pi/2, p(0)=0], [q(0)=Pi-.01, p(0)=0], [q(0)=-Pi+.01, p(0)=0], [q(0)=0, p(0)=2.15], [q(0)=0, p(0)=-2.15]], linecolor=black, stepsize=0.1);

 

The Hamiltonian for the simple pendulum system is time-independent, thus the energy is conserved.  The trajectories in the phase space are contours of constant energy.  The curve emanating from the unstable equilibrium point , a hyperbolic point, is called a separatrix, which separates the "libration" (or oscillation) motion, appearing as a set of elliptical contours, and the "rotation" motion, appearing as a set of curves that do not cross the axis.  Because in coordinate space is equivalent , we form a cylinder by gluing together opposite sides of the phase space.  Three representative curves are plotted on this phase space cylinder. 

>	Soln1 := dsolve({Eq1, Eq2, q(0)=Pi/3, p(0)=0}, {q(t), p(t)}, numeric, output=listprocedure);

 

>	Soln2 := dsolve({Eq1, Eq2, q(0)=Pi-0.01, p(0)=0}, {q(t), p(t)}, numeric, output=listprocedure);

 

>	Soln3 := dsolve({Eq1, Eq2, q(0)=0, p(0)=-2.5}, {q(t), p(t)}, numeric, output=listprocedure);

 

>	with(plots):

 

Warning, the name changecoords has been redefined 

>	pt := plot3d([cos(u), sin(u), v], u=0..2*Pi, v=-3..3, style=hidden):

 

>	po := spacecurve([cos(rhs(Soln1(t)[3])), sin(rhs(Soln1(t)[3])), (rhs(Soln1(t)[2]))], t=0..8, color=blue, thickness=2):

 

>	ps := spacecurve([cos(rhs(Soln2(t)[3])), sin(rhs(Soln2(t)[3])), (rhs(Soln2(t)[2]))], t=-12..12, color=black, thickness=2, numpoints=400):

 

>	pr := spacecurve([cos(rhs(Soln3(t)[3])), sin(rhs(Soln3(t)[3])), (rhs(Soln3(t)[2]))], t=-4..4, color=red, thickness=2, numpoints=400):

 

>	display(pt, po, ps, pr);

 

>

 

On the cylinder it is apparent that libration or oscillation (blue curve) and rotation (red curve) are closed orbits of two different topological types; the separatrix (black curve) divides the phase space into three different regions.  Every point in the phase space must follow an orbit, either libration or rotation motion, thus there is no chaos in an integrable system.    

Kicked Rotor and Standard Map 

Under perturbation of the potential function, the separatrix might be broken down into a disordered region of chaos.  The mechanism to create chaos through homoclinic tangles is demonstrated in an animation in the last section.   

Consider a Hamiltonian with the potential energy turned on in delta-function pulses,   

. 

We may consider this "kicked rotor" system as a simple pendulum with potential energy turned on in pulses or "kicks" every unit time.  Integrating Hamiltonian's equations for this potential, one obtains[3, 4] 

These two equations known as the standard map.  One can apply periodic boundary conditions by taking the fractional part of the values of and , and the map becomes a square; we also use symbols and for the variables in our worksheet.  If we paste the opposite sides together, it forms a torus.     

One needs to make a distinction between the trajectory in a phase space and points in a map.  In a time-dependent system, there is an extra dimension.  The map is a reduced phase plot, or a Poincare surface of section, representing the intersection of the trajectory with the planes .     

For a standard map on a torus, in addition to libration and rotation motions similar to a simple pendulum system, see plots below, a new type of motion is a region of disorganized scatter of points.  The reader should vary the parameter K and initial condition p[1] and q[1] in the worksheet to discover the dynamical behavior of the standard map.  The libration motion is generated by choosing K:=0.9, p[1]:=0.2 and q[1]:=0; the rotation motion K:=0.9, p[1]:=0.7 and q[1]:=0.5.    

>

K := 0.9;

 

>	imax := 10000; p := Vector(imax): q := Vector(imax):

 

>	p[1] := 0.1; q[1] := 0.5;

 

>	for i from 1 to imax-1 do p[i+1] := (evalf(p[i]-K/(2*Pi)*sin(2*Pi*q[i]))):  q[i+1] := (q[i] + p[i+1]): end do:

 

>	plot([seq([frac(q[n]), frac(p[n])], n=1..imax)], style=point, symbol=point, color=black, axes=boxed, view=[-1..1, -1..1]);

 

>	x := (u,v) -> (c+a*cos(v))*cos(u);

 

>	y := (u,v) -> (c+a*cos(v))*sin(u);

 

>	z := (u,v) -> a*sin(v);

 

>	a := 2: c := 5:

 

>	pd := plot3d([x(u,v), y(u,v), z(u,v)], u=0..2*Pi, v=0..2*Pi, scaling=constrained, style=hidden):

 

>	pm := pointplot3d({seq([x(2*Pi*p[i], 2*Pi*q[i]), y(2*Pi*p[i], 2*Pi*q[i]), z(2*Pi*p[i], 2*Pi*q[i])], i=1..imax)}, color=black):

 

>	display([pd, pm]);

 

>

 

From the above figure, chaotic motion is bounded in the phase space.  To further investigate the standard map, we choose 40 initial points uniformly distributed in the phase space, and map each point forward for 400 times.  Again the reader is encouraged to vary the K value to explore this map.        

>	K := 0.8; imax := 400; nr := 20; theta1 := 0; theta2 := 0.5;

 

>	r := Matrix(2*nr,imax): theta := Matrix(2*nr,imax):

 

>	for j from 1 to nr do  r[j,1] := j/nr:  theta[j,1] := theta1: end do:

 

>	for j from 1 to nr do for i from 1 to imax-1 do r[j,i+1] := (evalf(r[j,i]-K/(2*Pi)*sin(2*Pi*theta[j,i]))):  theta[j,i+1] := (theta[j,i] + r[j,i+1]): end do: end do:

 

>	for j from 1 to nr do  r[j+nr,1] := j/nr:  theta[j+nr,1] := theta2: end do:

 

>	for j from 1 to nr do for i from 1 to imax-1 do r[j+nr,i+1] := (evalf(r[j+nr,i]-K/(2*Pi)*sin(2*Pi*theta[j+nr,i]))):  theta[j+nr,i+1] := (theta[j+nr,i] + r[j+nr,i+1]): end do: end do:

 

>	plot([seq(seq([frac(theta[j,n]), frac(r[j,n])], n=1..imax), j=1..2*nr)], style=point, symbol=point, color=black, axes=boxed, view=[-1..1, -1..1]);

 

>	pd2 := plot3d([x(u,v), y(u,v), z(u,v)], u=0..2*Pi, v=0..2*Pi, scaling=constrained, style=patchnogrid, lightmodel=light1):

 

>	pm := pointplot3d({seq(seq([x(2*Pi*r[j,i], 2*Pi*theta[j,i]), y(2*Pi*r[j,i], 2*Pi*theta[j,i]), z(2*Pi*r[j,i], 2*Pi*theta[j,i])], i=1..imax), j=1..2*nr)}, color=black):

 

>	display([pd2, pm]);

 

>

 

In this map, the elliptical curves analogous to libration in a simple pendulum system is a structure of "nonlinear resonance" or "resonance island," and the circles winding around the torus analogous to rotation in a simple pendulum system is a structure of "invariant torus" or "KAM torus" (a circle is a 1-torus, not to confuse with the phase space 2-torus), after Kolmogorov, Arnold and Moser.  We show two maps below, for and .         

Without perturbation (), is constant and we should see circles (invariant 1-tori) winding around the phase space 2-torus.  The KAM theorem states that under a small perturbation most of the invariant 1-tori remain.  This fact is shown in the map for .  As increases, invariant tori becomes distorted and the resonance islands are formed.  When is near 0.8, chaotic motion appears prominently.  However, the invariant tori confine chaotic motion in a bounded region.  When is near 1, all the tori are broken and chaotic motion can wander further.  With an even higher value of all the resonance islands are distorted and the motion is entirely chaotic.        

Homoclinic Tangle 

A homoclinic tangle is a fundamental mechanism by which chaos is created near the hyperbolic point .  A hyperbolic point has stable and unstable directions, determined by the Jacobian.  For the simple pendulum system, the stable and unstable separatrices connect smoothly.  For the standard map, at a large value of stable and unstable separatrices oscillate wildly about the hyperbolic point.  The unstable separatrix is constructed by taking 1000 random points randomly distributed close to the hyperbolic point and in the unstable direction, then map each point forward.  The stable separatrix is constructed similarly but in stable direction and map backward.                            

>

restart:

 

>	K := 1.5; fpr := 0; fpt := 0.5;

 

>	with(LinearAlgebra):

 

>	Jac := <<1+K| 1>, <K| 1>>;

 

>	(lam, vec) := Eigenvectors(Jac);

 

>	imax := 12; nr := 1000;

 

>	r1 := Matrix(nr,imax): theta1 := Matrix(nr,imax):

 

>	r2 := Matrix(nr,imax): theta2 := Matrix(nr,imax):

 

>	Rr := (RandomTools[Generate](list(float(range=0..0.002, method=uniform), nr))):

 

>	Rt := (RandomTools[Generate](list(float(range=0..0.002, method=uniform), nr))):

 

>	for j from 1 to nr do  r1[j,1] := fpr + abs(vec[1,1])*(Rr[j]-0.001):  theta1[j,1] := fpt + abs(vec[2,1])*(Rt[j]-0.001): end do:

 

>	for j from 1 to nr do for i from 1 to imax-1 do r1[j,i+1] := (evalf(r1[j,i]-K/(2*Pi)*sin(2*Pi*theta1[j,i]))):  theta1[j,i+1] := (theta1[j,i] + r1[j,i+1]): end do: end do:

 

>	for j from 1 to nr do  r2[j,1] := fpr + abs(vec[1,2])*(Rr[j]-0.001):  theta2[j,1] := fpt + abs(vec[2,2])*(Rt[j]-0.001): end do:

 

>	for j from 1 to nr do for i from 1 to imax-1 do theta2[j,i+1] := (theta2[j,i] - r2[j,i]): r2[j,i+1] := (evalf(r2[j,i]+K/(2*Pi)*sin(2*Pi*theta2[j,i+1]))): end do: end do:

 

>	with(plots):

 

Warning, the name changecoords has been redefined 

>

for i from 1 to imax do pf[i] := plot([seq([frac(theta1[j,i]), frac(r1[j,i])], j=1..nr)], style=point, symbol=point, color=blue, axes=boxed, view=[-1..1, -1..1]): pb[i] := plot([seq([frac(theta2[j,i]), frac(r2[j,i])], j=1..nr)], style=point, symbol=point, color=red, axes=boxed, view=[-1..1, -1..1]): pt[i] := display([pf[i], pb[i]]): end do:

 

>	display([seq(pt[i], i=1..imax)], insequence=true);

 

>	x := (u,v) -> (c+a*cos(v))*cos(u);

 

>	y := (u,v) -> (c+a*cos(v))*sin(u);

 

>	z := (u,v) -> a*sin(v);

 

>	a := 2: c := 5:

 

>	pd := plot3d([x(u,v), y(u,v), z(u,v)], u=0..2*Pi, v=0..2*Pi, style=hidden):

 

>

for i from 1 to imax do pf[i] := pointplot3d({seq([x(2*Pi*r1[j,i], 2*Pi*theta1[j,i]), y(2*Pi*r1[j,i], 2*Pi*theta1[j,i]), z(2*Pi*r1[j,i], 2*Pi*theta1[j,i])], j=1..nr)}, color=blue): pb[i] := pointplot3d({seq([x(2*Pi*r2[j,i], 2*Pi*theta2[j,i]), y(2*Pi*r2[j,i], 2*Pi*theta2[j,i]), z(2*Pi*r2[j,i], 2*Pi*theta2[j,i])], j=1..nr)}, color=red): ptd[i] := display([pf[i], pb[i], pd]): end do:

 

>	display([seq(ptd[i], i=1..imax)], insequence=true, scaling=constrained, orientation=[180, 45]);

 

>

 

Further theoretical background for this worksheet can be found in the references listed below, particularly [3, 4].   

References 

1. H. Goldstein, C. P. Poole and J. L. Safko, Classical Mechanics, 3rd ed., Addison-Wesley, 2002. 

2. V. I. Arnold, Mathematical Methods of Classical Mechanics, 2nd ed., Springer-Verlag, 1989. 

3. Todd Timberlake, "A computational approach to teaching conservative chaos," American Journal of Physics 72, 1002--1007 (2004), and references therein.  http://fsweb.berry.edu/academic/mans/ttimberlake/cc    

4. Linda E. Reichl, The Transition to Chaos: Conservative Classical systems and Quantum Manifestations, 2nd ed., Spring-Verlag, 2004.