The mathematical pendulum Based on a MAPLE worksheet by H.J. L\303\274dde making use of N.J. Giordano, Computational Physics, Chapter 3 Modified by J. Reinhardt, 2007, 2009 In this worksheet we study one of the simplest dynamical systems, the mathematical pendulum. A mass m is attached to a massless string (to be more precise: to a rigid stick) of length l and is subject to the force of gravitation. The dynamical state this system is describe by two variables, the angle of deflections \316\270 and the corresponding angular velocity \317\211= d\316\270/dt.The restoring force is proportional to sin \316\270; and thus is nonlinear. In addition there can be a frictional force proportional to the angular velocity. The trajectories of the free, undisturbed pendulum are very simple. However, if the pendulum is subject to a time-varying external perturbation the situation changes drastically. For certain parameters the driven pendulum exhibits irregular nonperiodic solutions and is an exemplary system for studying the properties of chaos. In this worksheet we study 1) the free mathematical pendulum, 2) the trajectories of a driven pendulum, 3) an analysis of chaotic motion.

restart: with(plots): We compare the exact restoring force sin \316\270 with the power series approximation (Taylor expansion) to various orders: s:=taylor(sin(theta),theta=0,8); sp:=array(1..4): for i from 1 to 4 do s:=taylor(sin(theta),theta=0,2*i): sp[i]:=convert(s,polynom): od: plot([seq(sp[i](theta),i=1..4),sin(theta)],theta=0..2*Pi,y=-1..1); For small angles \316\270<0.4 the linear approximation is reliable. LUklbXJvd0c2Iy9JK21vZHVsZW5hbWVHNiJJLFR5cGVzZXR0aW5nR0koX3N5c2xpYkdGJzZALUkjbWlHRiQ2JVEjV2VGJy8lJ2l0YWxpY0dRJmZhbHNlRicvJSxtYXRodmFyaWFudEdRJ25vcm1hbEYnLUkjbW9HRiQ2LVEifkYnRjIvJSZmZW5jZUdGMS8lKnNlcGFyYXRvckdGMS8lKXN0cmV0Y2h5R0YxLyUqc3ltbWV0cmljR0YxLyUobGFyZ2VvcEdGMS8lLm1vdmFibGVsaW1pdHNHRjEvJSdhY2NlbnRHRjEvJSdsc3BhY2VHUSYwLjBlbUYnLyUncnNwYWNlR0ZJLUYsNiVRJnNvbHZlRidGL0YyRjUtRiw2JVEkdGhlRidGL0YyRjUtRiw2JVEkdHdvRidGL0YyRjUtRiw2JVEoc3lzdGVtc0YnRi9GMkY1LUYsNiVRI29mRidGL0YyRjUtRiw2JVElT0RFc0YnRi9GMkY1LUYsNiVRLG51bWVyaWNhbGx5RidGL0YyLUY2Ni1RIixGJ0YyRjkvRjxRJXRydWVGJ0Y9Rj9GQUZDRkVGRy9GS1EsMC4zMzMzMzMzZW1GJ0Y1LUYsNiVRJ3Rha2luZ0YnRi9GMkY1LUYsNiVRImdGJ0YvRjItRjY2LVEiL0YnRjJGOUY7L0Y+Rl9vRj9GQUZDRkUvRkhRLDAuMTY2NjY2N2VtRicvRktGXXAtRiw2JVEibEYnRi9GMi1GNjYtUSI9RidGMkY5RjtGPUY/RkFGQ0ZFL0ZIUSwwLjI3Nzc3NzhlbUYnL0ZLRmZwLUkjbW5HRiQ2JFEiMUYnRjItRjY2LVEiOkYnRjJGOUY7Rj1GP0ZBRkNGRUZlcEZncEY1LUknbXNwYWNlR0YkNiYvJSdoZWlnaHRHUSYwLjBleEYnLyUmd2lkdGhHUSYwLjBlbUYnLyUmZGVwdGhHRmRxLyUqbGluZWJyZWFrR1ElYXV0b0YnLUYsNiNRIUYnLyUrZXhlY3V0YWJsZUdGMUYy ps1:=diff(theta(t),t)=omega(t); ps2:=diff(omega(t),t)=-om2*sin(theta(t)); ph2:=diff(omega(t),t)=-om2*theta(t); om2:=1.: om0:=0.: th0:=120.*(Pi/180.): sols:=dsolve({ps1,ps2,theta(0)=th0,omega(0)=om0},{theta(t),omega(t)},type=numeric): solx:=dsolve({ps1,ph2,theta(0)=th0,omega(0)=om0},{theta(t),omega(t)},type=numeric): fs:=odeplot(sols,[t,theta(t)],0..100,numpoints=500, color=red): fx:=odeplot(solx,[t,theta(t)],0..100,numpoints=500, color=blue): display(fs,fx,view=[0..20,-Pi..Pi],axes=boxed,labels=[t,theta],title=`Mathematical Pendulum (red) vs harmonic approximation (blue)`); For large amplitudes a clear difference is seen. The pendulum no longer undergoes harmonic oscillation and its period grows to infinity as \316\270LUklbXJvd0c2Iy9JK21vZHVsZW5hbWVHNiJJLFR5cGVzZXR0aW5nR0koX3N5c2xpYkdGJzYkLUklbXN1YkdGJDYlLUkjbWlHRiQ2I1EhRictRiM2Jy1JI21uR0YkNiRRIjBGJy8lLG1hdGh2YXJpYW50R1Enbm9ybWFsRicvJSdpdGFsaWNHUSV0cnVlRicvJStmb3JlZ3JvdW5kR1EsWzIwMCwwLDIwMF1GJy8lLHBsYWNlaG9sZGVyR0Y9L0Y5USdpdGFsaWNGJy8lL3N1YnNjcmlwdHNoaWZ0R1EiMEYnRjg= approaches 180 degrees. We also can inspect the phase space plot. fs:=odeplot(sols,[theta(t),omega(t)],0..100,numpoints=500, color=red): fx:=odeplot(solx,[theta(t),omega(t)],0..100,numpoints=500, color=blue): display(fs,fx,view=[-Pi..Pi,-Pi..Pi],axes=boxed,labels=[theta,omega], title=`Mathematical Pendulum (red) vs harmonic approximation (blue)`); Try out various initial conditions! The trajectories of a dynamical system in a two-dimensional phase space (x,y) can be understood if one plots the derivatives dy/dx = (dy/dt)/(dx/dt) in the form of little arrows on a grid. These arrows are tangent vectors to the trajectories going through the respective grid points. MAPLE offers a command which generates this type of tangent field plots automatically. restart: with(plots): with(DEtools): om2:=1.: ps1:=diff(theta(t),t)=omega(t); ps2:=diff(omega(t),t)=-om2*sin(theta(t)); dfieldplot([ps1,ps2],[theta,omega],t=0..1,theta=-2*Pi..+2*Pi,omega=-6..+6,dirgrid=[30,30],color=blue,arrows=large,axes=boxed,title=`Mathematical pendulum: Tangent field plot`);

Let us consider which types of trajectories are possible. Starting with an initial angular velocity LUklbXJvd0c2Iy9JK21vZHVsZW5hbWVHNiJJLFR5cGVzZXR0aW5nR0koX3N5c2xpYkdGJzYkLUklbXN1YkdGJDYlLUkjbWlHRiQ2JlEnJiM5Njk7RicvJSdpdGFsaWNHUSZmYWxzZUYnLyU2c2VsZWN0aW9uLXBsYWNlaG9sZGVyR1EldHJ1ZUYnLyUsbWF0aHZhcmlhbnRHUSdub3JtYWxGJy1GIzYnLUkjbW5HRiQ2JFEiMEYnRjgvRjNGNy8lK2ZvcmVncm91bmRHUSxbMjAwLDAsMjAwXUYnLyUscGxhY2Vob2xkZXJHRjcvRjlRJ2l0YWxpY0YnLyUvc3Vic2NyaXB0c2hpZnRHUSIwRidGOA===0 we find oscillatory motion (librations). If the pendulum is give a sufficiently strong initial kick it will undergo rotational motion. Since the system is conservative, which of the two types of motion is realized depends on the value of the total energy E= T + V. The phase space portrait of the pendulum shows the possible trajectories in dependence of the total energy. restart: with(plots): We first define the kinetic and potential energy of the pendulum: T:=1/2*m*l^2*omega^2; V:=m*g*l*(1-cos(theta)); For various values of the total energy E=T+U the contour plot shows the possible trajectories (lines of constant energy). The limiting trajectory with E=2gl corresponds to a pendulum which starts from the vertical position, undergoes one rotation and comes to rest again in the vertical position. For E>2gl the pendulum rotates and for E<2gl it oscillates. g:=9.81: l:=9.81: m:=1.: # thus g/l=1 contourplot(T+V,theta=-3/2*Pi..3/2*Pi,omega=-5..5,contours=[20,50,100,2*g*l,400,600,800,1000,1200],numpoints=1000,axes=boxed,filled=true, coloring=[green,black],title=`Phase space portrait of the mathematical pendulum`);

We study a pendulum with velocity-dependent friction and an external driving force which oscillates periodically in time with angular velocity \316\251. The resulting forms of motion depend sensitively on strength of the driving force LUklbXJvd0c2Iy9JK21vZHVsZW5hbWVHNiJJLFR5cGVzZXR0aW5nR0koX3N5c2xpYkdGJzYkLUklbXN1YkdGJDYlLUkjbWlHRiQ2JVEiRkYnLyUnaXRhbGljR1EldHJ1ZUYnLyUsbWF0aHZhcmlhbnRHUSdpdGFsaWNGJy1GIzYnLUYvNiVRI2V4RidGMkY1RjIvJStmb3JlZ3JvdW5kR1EsWzIwMCwwLDIwMF1GJy8lLHBsYWNlaG9sZGVyR0Y0RjUvJS9zdWJzY3JpcHRzaGlmdEdRIjBGJy9GNlEnbm9ybWFsRic=, the period \316\251 and the friction constant \316\261. We study the three values LUklbXJvd0c2Iy9JK21vZHVsZW5hbWVHNiJJLFR5cGVzZXR0aW5nR0koX3N5c2xpYkdGJzYlLUklbXN1YkdGJDYlLUkjbWlHRiQ2JVEiRkYnLyUnaXRhbGljR1EldHJ1ZUYnLyUsbWF0aHZhcmlhbnRHUSdpdGFsaWNGJy1GIzYoLUYvNiVRI2V4RidGMkY1RjIvJStmb3JlZ3JvdW5kR1EsWzIwMCwwLDIwMF1GJy8lLHBsYWNlaG9sZGVyR0Y0LyUwZm9udF9zdHlsZV9uYW1lR1ElVGV4dEYnRjUvJS9zdWJzY3JpcHRzaGlmdEdRIjBGJ0ZCL0Y2USdub3JtYWxGJw===0, 0.5, 1.2. restart: with(plots): pf1:=diff(theta(t),t)=omega(t); pf2:=diff(omega(t),t)=-om2*sin(theta(t))-alpha*omega(t)+Fex*sin(Omega*t); om2:=1.: alpha:=0.5: Omega:=2./3.: # system parameters th0:=Pi/2: om0:=0.: # initial conditions F:=array(1..3): # strength of driving force F[1]:=0.: F[2]:=0.5: F[3]:=1.2: mycolor:=['red','green','blue']: f:=array(1..3): Now the differential equations are solved three times: for i from 1 to 3 do Fex:=F[i]: sol:=dsolve({pf1,pf2,theta(0)=th0,omega(0)=om0},{theta(t),omega(t)},type=numeric,output=listprocedure): f[i]:=odeplot(sol,[t,theta(t)],0..60,numpoints=500,color=mycolor[i]); od: display(f[1],f[2],f[3],view=[0..60,-4..+4],axes=boxed,labels=[t,theta], title= `Driven Pendulum for Fex=0,0.5,1.2`); f1:=array(1..3): for i from 1 to 3 do Fex:=F[i]: sol:=dsolve({pf1,pf2,theta(0)=th0,omega(0)=om0},{theta(t),omega(t)},type=numeric,output=listprocedure): f1[i]:=odeplot(sol,[theta(t),omega(t)],0..60,numpoints=500,color=mycolor[i]): od: display(f1[1],f1[2],f1[3],view=[-4..4,-Pi..Pi],axes=boxed,labels=[theta,omega], title= `Driven Pendulum for Fex=0,0.5,1.2: Phase space plot for t=0..60`); Three different types of trajectories are seen. For LUklbXJvd0c2Iy9JK21vZHVsZW5hbWVHNiJJLFR5cGVzZXR0aW5nR0koX3N5c2xpYkdGJzYlLUklbXN1YkdGJDYlLUkjbWlHRiQ2JVEiRkYnLyUnaXRhbGljR1EldHJ1ZUYnLyUsbWF0aHZhcmlhbnRHUSdpdGFsaWNGJy1GIzYoLUYvNiVRI2V4RidGMkY1RjIvJStmb3JlZ3JvdW5kR1EsWzIwMCwwLDIwMF1GJy8lLHBsYWNlaG9sZGVyR0Y0LyUwZm9udF9zdHlsZV9uYW1lR1ElVGV4dEYnRjUvJS9zdWJzY3JpcHRzaGlmdEdRIjBGJ0ZCL0Y2USdub3JtYWxGJw===0 the trajectory approaches a stable fixed point attractor. The trajectory for LUklbXJvd0c2Iy9JK21vZHVsZW5hbWVHNiJJLFR5cGVzZXR0aW5nR0koX3N5c2xpYkdGJzYlLUklbXN1YkdGJDYlLUkjbWlHRiQ2JVEiRkYnLyUnaXRhbGljR1EldHJ1ZUYnLyUsbWF0aHZhcmlhbnRHUSdpdGFsaWNGJy1GIzYoLUYvNiVRI2V4RidGMkY1RjIvJStmb3JlZ3JvdW5kR1EsWzIwMCwwLDIwMF1GJy8lLHBsYWNlaG9sZGVyR0Y0LyUwZm9udF9zdHlsZV9uYW1lR1ElVGV4dEYnRjUvJS9zdWJzY3JpcHRzaGlmdEdRIjBGJ0ZCL0Y2USdub3JtYWxGJw===0.5 approaches a limit cycle. For LUklbXJvd0c2Iy9JK21vZHVsZW5hbWVHNiJJLFR5cGVzZXR0aW5nR0koX3N5c2xpYkdGJzYlLUklbXN1YkdGJDYlLUkjbWlHRiQ2JVEiRkYnLyUnaXRhbGljR1EldHJ1ZUYnLyUsbWF0aHZhcmlhbnRHUSdpdGFsaWNGJy1GIzYoLUYvNiVRI2V4RidGMkY1RjIvJStmb3JlZ3JvdW5kR1EsWzIwMCwwLDIwMF1GJy8lLHBsYWNlaG9sZGVyR0Y0LyUwZm9udF9zdHlsZV9uYW1lR1ElVGV4dEYnRjUvJS9zdWJzY3JpcHRzaGlmdEdRIjBGJ0ZCL0Y2USdub3JtYWxGJw===1.2 the trajectory shows irregular behaviour. Will it stabilize at later times? Let us follow the trajectories for a time interval ten times longer. f2:=array(1..3): for i from 1 to 3 do Fex:=F[i]: sol:=dsolve({pf1,pf2,theta(0)=th0,omega(0)=om0},{theta(t),omega(t)},type=numeric,output=listprocedure): f2[i]:=odeplot(sol,[theta(t),omega(t)],0..600,numpoints=2000,color=mycolor[i]): od: display(f2[1],f2[2],f2[3],view=[-5..70,-Pi..Pi],axes=boxed,labels=[theta,omega], title= `Driven Pendulum for Fex=0,0.5,1.2: Phase space plot for t=0..600`); Obviously the third trajectory remains irregular, even if we follow it for a long time. We are in the region of chaos. Try out various initial conditions and system parameters.

Let us study what happens if we slightly vary the inital conditions. For this we integrate the ODEs twice, with initial conditions which differ by a tiny amount. We plot the absolute value of the difference between the values of \316\270(t) for the two calculations. restart: with(plots): pf1:=diff(theta(t),t)=omega(t); pf2:=diff(omega(t),t)=-om2*sin(theta(t))-alpha*omega(t)+Fex*sin(Omega*t); om2:=1: alpha:=0.5: Omega:=2/3: # system parameters th0:=Pi/2: om0:=0.: # inital conditions F:=array(1..3): F[1]:=0.: # strength of driving force F[2]:=0.5: F[3]:=1.2: mycolor:=['red','green','blue']: f3:=array(1..3): for i from 1 to 3 do Fex:=F[i]: sol:=dsolve({pf1,pf2,theta(0)=th0,omega(0)=om0},{theta(t),omega(t)},type=numeric,output=listprocedure): th:=subs(sol,theta(t)): sol1:=dsolve({pf1,pf2,theta(0)=1.000001*th0,omega(0)=1.000001*om0},{theta(t),omega(t)},type=numeric,output=listprocedure): th1:= subs(sol1,theta(t)): deltheta := th1-th: f3[i]:=logplot(abs(deltheta(t)),t=0..100,numpoints=500,color=mycolor[i]): od: display(f3[1],f3[2],f3[3],view=[0..50,1e-10..0.1],axes=boxed,labels=[t,deltatheta], title= `Driven Pendulum for Fex=0,0.5,1.2: Init. values differ by 1ppm`); We note a striking difference between the cases with regular motion (red: LUklbXJvd0c2Iy9JK21vZHVsZW5hbWVHNiJJLFR5cGVzZXR0aW5nR0koX3N5c2xpYkdGJzYlLUklbXN1YkdGJDYlLUkjbWlHRiQ2JVEiRkYnLyUnaXRhbGljR1EldHJ1ZUYnLyUsbWF0aHZhcmlhbnRHUSdpdGFsaWNGJy1GIzYoLUYvNiVRI2V4RidGMkY1RjIvJStmb3JlZ3JvdW5kR1EsWzIwMCwwLDIwMF1GJy8lLHBsYWNlaG9sZGVyR0Y0LyUwZm9udF9zdHlsZV9uYW1lR1ElVGV4dEYnRjUvJS9zdWJzY3JpcHRzaGlmdEdRIjBGJ0ZCL0Y2USdub3JtYWxGJw===0, stable fixed point attractor, green: LUklbXJvd0c2Iy9JK21vZHVsZW5hbWVHNiJJLFR5cGVzZXR0aW5nR0koX3N5c2xpYkdGJzYlLUklbXN1YkdGJDYlLUkjbWlHRiQ2JVEiRkYnLyUnaXRhbGljR1EldHJ1ZUYnLyUsbWF0aHZhcmlhbnRHUSdpdGFsaWNGJy1GIzYoLUYvNiVRI2V4RidGMkY1RjIvJStmb3JlZ3JvdW5kR1EsWzIwMCwwLDIwMF1GJy8lLHBsYWNlaG9sZGVyR0Y0LyUwZm9udF9zdHlsZV9uYW1lR1ElVGV4dEYnRjUvJS9zdWJzY3JpcHRzaGlmdEdRIjBGJ0ZCL0Y2USdub3JtYWxGJw===0.5, limit cycle attractor) and the chaotic case (blue, LUklbXJvd0c2Iy9JK21vZHVsZW5hbWVHNiJJLFR5cGVzZXR0aW5nR0koX3N5c2xpYkdGJzYlLUklbXN1YkdGJDYlLUkjbWlHRiQ2JVEiRkYnLyUnaXRhbGljR1EldHJ1ZUYnLyUsbWF0aHZhcmlhbnRHUSdpdGFsaWNGJy1GIzYoLUYvNiVRI2V4RidGMkY1RjIvJStmb3JlZ3JvdW5kR1EsWzIwMCwwLDIwMF1GJy8lLHBsYWNlaG9sZGVyR0Y0LyUwZm9udF9zdHlsZV9uYW1lR1ElVGV4dEYnRjUvJS9zdWJzY3JpcHRzaGlmdEdRIjBGJ0ZCL0Y2USdub3JtYWxGJw===1.2). In the first two cases the difference \316\264\316\270(t) drops to zero as the trajectories approach the attractor. In the latter case, however, the difference on the average increases exponentially (note the logarthmic scale). This is a quantitative criterion for chaotic motion: Neighboring trajectories with slightly different initial conditions are moving away from each other with exponentially increasing distance. The slope of this increase is called the Lyapunov exponent. If it is positive the system is chaotic. Try out different initial conditions in sol1. Study the dependence of the Lyapunov exponent on on the driving force LUklbXJvd0c2Iy9JK21vZHVsZW5hbWVHNiJJLFR5cGVzZXR0aW5nR0koX3N5c2xpYkdGJzYmLUklbXN1YkdGJDYlLUkjbWlHRiQ2JVEiRkYnLyUnaXRhbGljR1EldHJ1ZUYnLyUsbWF0aHZhcmlhbnRHUSdpdGFsaWNGJy1GIzYoLUYvNiVRI2V4RidGMkY1RjIvJStmb3JlZ3JvdW5kR1EsWzIwMCwwLDIwMF1GJy8lLHBsYWNlaG9sZGVyR0Y0LyUwZm9udF9zdHlsZV9uYW1lR1ElVGV4dEYnRjUvJS9zdWJzY3JpcHRzaGlmdEdRIjBGJ0YyRkJGNQ==.

A useful tool for studying dynamical systems are stroboscopic maps. Instead of looking at the full time-dependent trajectories, only snapshots at discrete moments in time are taken, each of them represented by a single point in phase space. For a system with a periodic external driving force it is natural to choose equidistant sampling times corresponding to the driving time period. For the driven pendulum we will take sampling times defined by > \316\251 t = n \317\200;. The stroboscopic map is a special case of a Poincar\303\251 section used universally for studying nonlinear dynamical systems. Below we study the driven pendulum problem by integrating three coupled first-order ODEs and plotting points in the phase space at times LUklbXJvd0c2Iy9JK21vZHVsZW5hbWVHNiJJLFR5cGVzZXR0aW5nR0koX3N5c2xpYkdGJzYkLUklbXN1YkdGJDYlLUkjbWlHRiQ2JVEidEYnLyUnaXRhbGljR1EldHJ1ZUYnLyUsbWF0aHZhcmlhbnRHUSdpdGFsaWNGJy1GIzYlLUYvNiVRIm5GJ0YyRjVGMkY1LyUvc3Vic2NyaXB0c2hpZnRHUSIwRicvRjZRJ25vcm1hbEYn = n \317\200/ \316\251. For practical reasons we split the whole trajectory into jj=200 blocks with npt=50 sampling points each. When plotting the phase space diagram, the periodic angle variable \316\270 is projected into the interval [-\317\200 .. \317\200]; using the function 'normangle' defined below. restart: with(plots): ode1:=diff(theta(t),t)=omega(t); ode2:=diff(omega(t),t)=-om2*sin(theta(t))-alpha*omega(t)+Fex*sin(phi(t)); ode3:=diff(phi(t),t)=Omega; normangle:= x->x-2*Pi*trunc((x/Pi+sign(x))/2): # normalize angles to interval -Pi..+Pi om2:=1:alpha:=0.5: Fex:=1.20: Omega:=2./3.: # Parameters for chaotic dynamics th0:=Pi/2: om0:=0: ph0:=0.: jj:=200: # jj blocks with 50 pts each jjtrans:=10: # The first 500 points are discarded initcond:=theta(0)=th0, omega(0)=om0, phi(0)=ph0: # initial conditions The following calculation is time consuming. for j from 1 to jj do # loop over jj vars:=theta(t), omega(t), phi(t): sol:=dsolve({ode1,ode2,ode3,initcond},{vars},type=numeric,output=listprocedure); assign(sol); theta:=theta(t); omega:=omega(t); phi:=phi(t); npt:=50; # number of stroboscopic points in one block pts:=seq([normangle(theta(Pi*i/Omega)),omega(Pi*i/Omega)],i=1..npt); # (i) only points at 'stroboscopic times' are taken; (ii) theta is reduced to the interval [-Pi..Pi] plt(j):=plot([pts],view=[-Pi..Pi,-2..+2],color=blue); th0:=theta(Pi*npt/Omega); om0:=omega(Pi*npt/Omega): ph0:=phi(Pi*npt/Omega): unassign(t,theta,omega,phi); initcond:=theta(0)=th0,omega(0)=om0,phi(0)=ph0; # initial conditions for the next block od: display([seq(plt(j),j=jjtrans..jj)],symbol=POINT,style=point,view=[-Pi..Pi,-2..2],axes=boxed,labels=[theta,omega], title= `Poincare Section for Fex=1.2`); The Poincar\303\251 section shows a distorted but quite regular and localized geometrical object. Closer inspection reveals that the apparent lines consist of bundles of lines which again consist of narrower bundles, and so on ad infinitum. This is the typical property of a fractal object. The picture shows a strange attractor characteristic for a chaotic system. Show that the shape of the strange attractor does not depend on the initial condition. Have a look at the 'finestructure' of the strange attactor by 'zooming in'. This can be done by changing the viewing window interactively or in the display command. Study the change of the Poincar\303\251 section when going from regular to chaotic motion. Suggestion: Compare LUklbXN1Ykc2Iy9JK21vZHVsZW5hbWVHNiJJLFR5cGVzZXR0aW5nR0koX3N5c2xpYkdGJzYlLUkjbWlHRiQ2JVEiRkYnLyUnaXRhbGljR1EldHJ1ZUYnLyUsbWF0aHZhcmlhbnRHUSdpdGFsaWNGJy1JJW1yb3dHRiQ2JC1GLDYlUSNleEYnRi9GMi9GM1Enbm9ybWFsRicvJS9zdWJzY3JpcHRzaGlmdEdRIjBGJw===1.14 and LUklbXN1Ykc2Iy9JK21vZHVsZW5hbWVHNiJJLFR5cGVzZXR0aW5nR0koX3N5c2xpYkdGJzYlLUkjbWlHRiQ2JVEiRkYnLyUnaXRhbGljR1EldHJ1ZUYnLyUsbWF0aHZhcmlhbnRHUSdpdGFsaWNGJy1JJW1yb3dHRiQ2JC1GLDYlUSNleEYnRi9GMi9GM1Enbm9ybWFsRicvJS9zdWJzY3JpcHRzaGlmdEdRIjBGJw===1.15.

Up to now two typical properties of chaotic motion have been studied: (i) extreme sensitivity of trajectories on the initial conditions, (ii) characteristic geometric structure of the Poincar\303\251 sections. Now we investigate the transition from regular to chaotic motion in dependence of control parameters. We compare four solutions of the driven coupled pendulum for the force constants LUklbXJvd0c2Iy9JK21vZHVsZW5hbWVHNiJJLFR5cGVzZXR0aW5nR0koX3N5c2xpYkdGJzYzLUklbXN1YkdGJDYlLUkjbWlHRiQ2JVEiRkYnLyUnaXRhbGljR1EldHJ1ZUYnLyUsbWF0aHZhcmlhbnRHUSdpdGFsaWNGJy1GIzYkLUYvNiVRI2V4RidGMkY1L0Y2USdub3JtYWxGJy8lL3N1YnNjcmlwdHNoaWZ0R1EiMEYnLUkjbW9HRiQ2LVEiPUYnRj0vJSZmZW5jZUdRJmZhbHNlRicvJSpzZXBhcmF0b3JHRkgvJSlzdHJldGNoeUdGSC8lKnN5bW1ldHJpY0dGSC8lKGxhcmdlb3BHRkgvJS5tb3ZhYmxlbGltaXRzR0ZILyUnYWNjZW50R0ZILyUnbHNwYWNlR1EsMC4yNzc3Nzc4ZW1GJy8lJ3JzcGFjZUdGVy1JI21uR0YkNiRRJDEuNEYnRj0tRkM2LVEiLEYnRj1GRi9GSkY0RktGTUZPRlFGUy9GVlEmMC4wZW1GJy9GWVEsMC4zMzMzMzMzZW1GJy1GQzYtUSJ+RidGPUZGRklGS0ZNRk9GUUZTRlxvL0ZZRl1vLUZlbjYkUSUxLjQ1RidGPUZobkZgby1GZW42JFElMS40N0YnRj1GaG5GYG8tRmVuNiRRJDEuNUYnRj0tRkM2LVEiLkYnRj1GRkZJRktGTUZPRlFGU0Zcb0Zjby1JJ21zcGFjZUdGJDYmLyUnaGVpZ2h0R1EmMC4wZXhGJy8lJndpZHRoR1EmMC4wZW1GJy8lJmRlcHRoR0ZlcC8lKmxpbmVicmVha0dRKG5ld2xpbmVGJy1GYXA2JkZjcEZmcEZpcC9GXHFRJWF1dG9GJy1GLzYjUSFGJ0Y9restart: with(plots): pf1:=diff(theta(t),t)=omega(t):pf2:=diff(omega(t),t)=-om2*sin(theta(t))-alpha*omega(t)+Fex*sin(Omega*t): om2:=1.: alpha:=0.5: Omega:=2./3.: th0:=Pi/2: om0:=0.: F:=array(1..4): F[1]:=1.4: F[2]:=1.45: F[3]:=1.47: F[4]:=1.5: mycolor:=['red','green','blue','black']: normangle:= x->x-2*Pi*trunc((x/Pi+sign(x))/2): # normalize angles to interval -Pi..+Pi for j from 1 to 4 do Fex:=F[j]: sol:=dsolve({pf1,pf2,theta(0)=th0,omega(0)=om0},{theta(t),omega(t)},type=numeric, output=listprocedure): assign(sol); theta:=theta(t); npt:=600; # number of points pts:=seq([i/4,normangle(theta(i/4))-(j-1)*2*Pi],i=1..npt); # shift successive curves by 2Pi outplot(j):=plot([pts],color=mycolor[j]): unassign(t,theta,omega); od: display([seq(outplot(j),j=1..4)],view=[50..150,-7*Pi..Pi],axes=boxed,labels=[t,theta], title= `Driven pendulum for Fex=1.4, 1.45, 1.47, 1.5`); For the given parameters the pendulum undergoes a modulated rotational motion in synchrony with the driving frequency (one rotation per time period T=2\317\200/\316\251. Following a transient region, for LUklbXN1Ykc2Iy9JK21vZHVsZW5hbWVHNiJJLFR5cGVzZXR0aW5nR0koX3N5c2xpYkdGJzYlLUkjbWlHRiQ2JVEiRkYnLyUnaXRhbGljR1EldHJ1ZUYnLyUsbWF0aHZhcmlhbnRHUSdpdGFsaWNGJy1JJW1yb3dHRiQ2JC1GLDYlUSNleEYnRi9GMi9GM1Enbm9ybWFsRicvJS9zdWJzY3JpcHRzaGlmdEdRIjBGJw===1.4 the motion is periodic. With increasing driving force the system undergoes a bifurcation and at LUklbXN1Ykc2Iy9JK21vZHVsZW5hbWVHNiJJLFR5cGVzZXR0aW5nR0koX3N5c2xpYkdGJzYlLUkjbWlHRiQ2JVEiRkYnLyUnaXRhbGljR1EldHJ1ZUYnLyUsbWF0aHZhcmlhbnRHUSdpdGFsaWNGJy1JJW1yb3dHRiQ2JC1GLDYlUSNleEYnRi9GMi9GM1Enbm9ybWFsRicvJS9zdWJzY3JpcHRzaGlmdEdRIjBGJw===1.45 we find a period doubling. A second period doubling follows (to a period of 4T) as observed for LUklbXN1Ykc2Iy9JK21vZHVsZW5hbWVHNiJJLFR5cGVzZXR0aW5nR0koX3N5c2xpYkdGJzYlLUkjbWlHRiQ2JVEiRkYnLyUnaXRhbGljR1EldHJ1ZUYnLyUsbWF0aHZhcmlhbnRHUSdpdGFsaWNGJy1JJW1yb3dHRiQ2JC1GLDYlUSNleEYnRi9GMi9GM1Enbm9ybWFsRicvJS9zdWJzY3JpcHRzaGlmdEdRIjBGJw===1.47. Such a cascade of period doublings leads to nonperiodic chaotic motion.

We determine a reference trajectory LUklbXJvd0c2Iy9JK21vZHVsZW5hbWVHNiJJLFR5cGVzZXR0aW5nR0koX3N5c2xpYkdGJzYlLUklbXN1YkdGJDYlLUkjbWlHRiQ2JVEieEYnLyUnaXRhbGljR1EldHJ1ZUYnLyUsbWF0aHZhcmlhbnRHUSdpdGFsaWNGJy1GIzYoLUYvNiVRInJGJ0YyRjVGMi8lK2ZvcmVncm91bmRHUSxbMjAwLDAsMjAwXUYnLyUscGxhY2Vob2xkZXJHRjQvJTBmb250X3N0eWxlX25hbWVHUSpIZWFkaW5nfjFGJ0Y1LyUvc3Vic2NyaXB0c2hpZnRHUSIwRidGQi9GNlEnbm9ybWFsRic=(t) and then solve the linearized problem 'around this trajectory'. restart: with(plots): deq1:=diff(xr(t),t)=yr(t); deq2:=diff(yr(t),t)=-om2*sin(xr(t))-alpha*yr(t)+Fex*sin(zr(t)); deq3:=diff(zr(t),t)=Omega; om2:=1.: alpha:=0.5: Omega:=2./3.: # system parameters Fex:=1.2: th0:=Pi/2: om0:=0.: # initial conditions The differential equations are solved for the reference trajectory sol:=dsolve({deq1,deq2,deq3,xr(0)=th0,yr(0)=om0,zr(0)=0.},{xr(t),yr(t),zr(t)},type=numeric,output=listprocedure): x1:=subs(sol,xr(t)): y1:=subs(sol,yr(t)): refplot:=odeplot(sol,[t,xr(t)],0..600,numpoints=2000): display(refplot,view=[0..600,-20..+20],axes=boxed,labels=[t,x], title=`Driven Pendulum for Fex=1.2`); The linearized differential equations. Since the third row of the Jacobian vanishes there are only two differential equations. Note that the reference trajectory (x(t),y(t)) enters through the time-dependent coefficients. ld1:=diff(eta1(t),t)=eta2(t); ld2:=diff(eta2(t),t)=-om2*cos(x1(t))*eta1(t)-alpha*eta2(t); Solution of the linear system for an arbitrary initial condition: sollin:=dsolve({ld1,ld2,eta1(0)=1.,eta2(0)=1.},{eta1(t),eta2(t)},type=numeric,output=listprocedure): e1:=subs(sollin,eta1(t)): e2:=subs(sollin,eta2(t)): ploteta:=logplot(sqrt(e1(t)*e1(t)+e2(t)*e2(t)),t=0..100): display(ploteta,view=[0..100,0.1..10^8],axes=boxed,labels=[t,`|eta|`], title= `Linear solution`); Approximate evaluation of the largest Lyapunov exponent \317\203 = LUklbXJvd0c2Iy9JK21vZHVsZW5hbWVHNiJJLFR5cGVzZXR0aW5nR0koX3N5c2xpYkdGJzYnLUknbXVuZGVyR0YkNiUtSSNtb0dGJDYtUSRsaW1GJy8lLG1hdGh2YXJpYW50R1Enbm9ybWFsRicvJSZmZW5jZUdRJnVuc2V0RicvJSpzZXBhcmF0b3JHRjcvJSlzdHJldGNoeUdGNy8lKnN5bW1ldHJpY0dGNy8lKGxhcmdlb3BHRjcvJS5tb3ZhYmxlbGltaXRzR1EldHJ1ZUYnLyUnYWNjZW50R0Y3LyUnbHNwYWNlR1EmMC4wZW1GJy8lJ3JzcGFjZUdRLDAuMTY2NjY2N2VtRictRiM2Ji1JI21pR0YkNiVRInRGJy8lJ2l0YWxpY0dGQi9GM1EnaXRhbGljRictRi82LVEtJnJpZ2h0YXJyb3c7RidGMi9GNlEmZmFsc2VGJy9GOUZZL0Y7RkIvRj1GWS9GP0ZZL0ZBRlkvRkRGWS9GRlEsMC4yNzc3Nzc4ZW1GJy9GSUZbby1GLzYtUSgmaW5maW47RidGMkZYRlovRjtGWUZmbkZnbkZobkZpbkZFL0ZJRkdGMi8lLGFjY2VudHVuZGVyR0ZZLUZONiNRIUYnLUYjNiVGZG8tSSZtZnJhY0dGJDYoLUYjNigtSSNtbkdGJDYkUSIxRidGMkZRLyUrZm9yZWdyb3VuZEdRK1swLDE2MCw4MF1GJy8lLHBsYWNlaG9sZGVyR0ZCLyU2c2VsZWN0aW9uLXBsYWNlaG9sZGVyR0ZCRlMtRiM2J0ZNRlEvRmNwUSxbMjAwLDAsMjAwXUYnRmVwRlMvJS5saW5ldGhpY2tuZXNzR1EiMUYnLyUrZGVub21hbGlnbkdRJ2NlbnRlckYnLyUpbnVtYWxpZ25HRmJxLyUpYmV2ZWxsZWRHRllGMkZkb0Yy LUklbXJvd0c2Iy9JK21vZHVsZW5hbWVHNiJJLFR5cGVzZXR0aW5nR0koX3N5c2xpYkdGJzYnLUkjbWlHRiQ2I1EhRictSSNtb0dGJDYtUSNsbkYnLyUsbWF0aHZhcmlhbnRHUSdub3JtYWxGJy8lJmZlbmNlR1EmdW5zZXRGJy8lKnNlcGFyYXRvckdGOC8lKXN0cmV0Y2h5R0Y4LyUqc3ltbWV0cmljR0Y4LyUobGFyZ2VvcEdGOC8lLm1vdmFibGVsaW1pdHNHRjgvJSdhY2NlbnRHRjgvJSdsc3BhY2VHUSYwLjBlbUYnLyUncnNwYWNlR0ZHLUYwNi1RMCZBcHBseUZ1bmN0aW9uO0YnRjMvRjdRJmZhbHNlRicvRjpGTi9GPEZOL0Y+Rk4vRkBGTi9GQkZOL0ZERk5GRUZILUkobWZlbmNlZEdGJDYkLUYjNiZGKy1GIzYlRistSSZtZnJhY0dGJDYoLUYjNiUtRlY2Ji1GIzYmLUYsNiVRJyYjOTUxO0YnLyUnaXRhbGljR0ZORjMtRlY2JC1GIzYlLUYsNiVRInRGJy9GY29RJXRydWVGJy9GNFEnaXRhbGljRidGW3BGXXBGM0ZbcEZdcEYzLyUlb3BlbkdRInxnckYnLyUmY2xvc2VHRmFwRltwRl1wLUYjNictRlY2Ji1GIzYoRl9vLUZWNiQtRiM2Jy1JI21uR0YkNiRRIjBGJ0YzRltwLyUrZm9yZWdyb3VuZEdRLFsyMDAsMCwyMDBdRicvJSxwbGFjZWhvbGRlckdGXHBGXXBGM0ZbcEZicUZlcUZdcEYzRl9wRmJwRltwRmJxRmVxRl1wLyUubGluZXRoaWNrbmVzc0dRIjFGJy8lK2Rlbm9tYWxpZ25HUSdjZW50ZXJGJy8lKW51bWFsaWduR0Zcci8lKWJldmVsbGVkR0ZORjNGK0YzRjNGMw== normeta:= t -> sqrt(e1(t)*e1(t)+e2(t)*e2(t)): lyapunov:= t-> log(normeta(t)/normeta(0))/t: lyapunov(100); lyapunov(200); lyapunov(500); LUklbXJvd0c2Iy9JK21vZHVsZW5hbWVHNiJJLFR5cGVzZXR0aW5nR0koX3N5c2xpYkdGJzYjLUkjbWlHRiQ2I1EhRic= It is difficult to reach high numerical accuracy.