Spring 2005

# Mathematics 118r Dynamical Systems Spring 2005

Course Head: Oliver knill
Office: SciCtr 434
Email: knill@math.harvard.edu

 News Info Plan Time Assign Project Exam Show Script Lab Faq Lib Link

## Laboratory

 C programs: Xstar, an N-body solver The best mandelbrot zoomer is definitely xaos. On OS X, you can install it with "fink install xaos" (provided, you have Fink installed. Mike Creutz xtoys precompiled for linux Mike Creutz xtoys precompiled for OS X xtoys page chirikov, dynamics in the plane

 CAS: (Classdemo) Mathematica code snippets. Note that each example is independent of the others. If you should try to run all of them in the same notebook, there could be confusions between variables. Just deal with each case seperatly. Exterior billiards at a semicircle ```(* It is an open problem since 1960, whether there exists a *) (* convex exterior billiard table with an unbounded orbit *) (* the following lines compute an orbit of 20000 at the *) (* semi circle |arg(z)|<= pi/2, |z|<=1 *) (* starting at z0 = 2+i. The orbit seems to escape to infinity *) (* Serge Tabatchnikov had suggested that this exterior billiard *) (* table is unstable. The problem is open. *) (* Mathematica code: Oliver Knill, March 2005, Math118r Harvard *) T1[z_]:=Module[{}, r=Abs[z]; alpha=2 ArcCos[1/r]; z*Exp[I alpha]]; T2[z_]:=z+2(-I-z); T3[z_]:=z+2(I-z); T[z_]:=N[If[Re[z]<0 && Im[z]>-1,T2[z],If[Re[z]>=0 && Im[z]>=1, T3[z],T1[z]]]]; A1=Graphics[{RGBColor[1,1,0],Disk[{0,0},1,{-Pi/2,Pi/2}]}]; ss=NestList[T,2+I,20000]; coord[z_]:={Re[z],Im[z]}; A3=Graphics[{RGBColor[1,0,0],PointSize[0.00001],Map[Point,Map[coord,ss]]}]; S=Show[A1,A3,AspectRatio->1,PlotRange->All] ``` Plotting geodesics on Surfaces (Code from Barrett O'Neill). This example plots geodesics in the hyperbolic plane: ```chr[1,1,1][e_,f_,g_][u_,v_]:=Simplify[(D[e[u,v],v]*f[u,v]-2*D[f[u,v],u]*f[u,v]+D[e[u,v],u]*g[u,v])/(2*(e[u,v]*g[u,v]-f[u,v]^2))]; chr[1,1,2][e_,f_,g_][u_,v_]:=Simplify[(D[e[u,v],v]*g[u,v]-D[g[u,v],u]*f[u,v])/(2*(e[u,v]*g[u,v]-f[u,v]^2))]; chr[1,2,1][e_,f_,g_][u_,v_]:=chr[1,1,2][e,f,g][u,v]; chr[1,2,2][e_,f_,g_][u_,v_]:=Simplify[(-D[g[u,v],v]*f[u,v]+2*D[f[u,v],v]*g[u,v]-D[g[u,v],u]*g[u,v])/(2*(e[u,v]*g[u,v]-f[u,v]^2))]; chr[2,1,1][e_,f_,g_][u_,v_]:=Simplify[(-D[e[u,v],v]*e[u,v]+2*D[f[u,v],u]*e[u,v]-D[e[u,v],u]*f[u,v])/(2*(e[u,v]*g[u,v]-f[u,v]^2))]; chr[2,1,2][e_,f_,g_][u_,v_]:=Simplify[(-D[e[u,v],v]*f[u,v]+D[g[u,v],u]*e[u,v])/(2*(e[u,v]*g[u,v]-f[u,v]^2))]; chr[2,2,1][e_,f_,g_][u_,v_]:=chr[2,1,2][e,f,g][u,v]; chr[2,2,2][e_,f_,g_][u_,v_]:=Simplify[(D[g[u,v],v]*e[u,v]-2*D[f[u,v],v]*f[u,v]+D[g[u,v],u]*f[u,v])/(2*(e[u,v]*g[u,v]-f[u,v]^2))]; eh[u_,v_]:=1/v^2;fh[u_,v_]:=0;gh[u_,v_]:=1/v^2 cht[i_,j_,k_][t_]:=chr[i,j,k][eh,fh,gh][u,v]/.{u->u[t],v->v[t]} gdesolv[u0_,v0_,du0_,dv0_,t0_,tmin_,tmax_]:=NDSolve[{u'[t]==p[t],v'[t]==q[t], p'[t]+cht[1,1,1][t]*p[t]^2+2cht[1,1,2][t]*p[t]*q[t]+cht[1,2,2][t]*q[t]^2==0, q'[t]+cht[2,1,1][t]*p[t]^2+2cht[2,1,2][t]*p[t]*q[t]+cht[2,2,2][t]*q[t]^2==0, u[t0]==u0,v[t0]==v0,p[t0]==du0,q[t0]==dv0},{u,v,p,q},{t,tmin,tmax}] soln1=gdesolv[1,1,1,2,0,-3,3];S1=ParametricPlot[Evaluate[{u[t],v[t]}/.soln1],{t,-3,3}] soln2=gdesolv[1,1,2,0,0,-3,3];S2=ParametricPlot[Evaluate[{u[t],v[t]}/.soln2],{t,-3,3}] soln3=gdesolv[1,1,1.5,-0.3,0,-3,3];S3=ParametricPlot[Evaluate[{u[t],v[t]}/.soln3],{t,-3,3}] Show[{S1,S2,S3}] ``` The Pollard rho method to factor integers uses a quadratic map T(x) = x2+c in the finite set of numbers modulo n. It is an extremely simple but effective method. Much more effective then simple trial division. Here is the simplest implementation without any tuning ```FactorWithPollard[n_] := Module[{a=3, x=17, y=17, q=1}, While[q<2, Do[x=Mod[x*x-a,n]; y=Mod[y*y-a,n]; y=Mod[y*y-a,n]; q=Mod[q*(x-y),n],{i,20}]; q=Mod[GCD[n,x-y],n] ];q]; FermatNumber[n_]:=2^(2^n)+1; FactorWithPollard[FermatNumber[5]] ``` Status on factoring Fermat numbers An unsolved problem. One knows that for almost all numbers x=[a0;a1,a2,a3,...] the limit (a1 a2 a3 ... an)1/n is the Kinchin constant 2.685452001.... It looks as if this limit is also obtained for pi, but one does not know: ```Module[{n=10000},s=ContinuedFraction[Pi,n]; (Product[s[[k]],{k,n}])^(1/n)] ``` Beta expansion. Mathematica has built in the beta expansions to any rel base b>1. ```b=Sqrt[2]; S=RealDigits[Pi,b,1000]; A=S[[1]]; L=S[[2]]; N[Sum[A[[i]] b^(L-i),{i,1,Length[A]}]] ``` Continued fraction expansion. Good approximations of Pi can be obtained by doing a continued fraction expansion to some degree, then take that rational number. Example: ```FromContinuedFraction[ContinuedFraction[Pi,20]] ``` Symbolic dynamics of Baker map. Iterate the Baker map and show simultaneously the evolution of the corresponding code. ```LengthOfDisplayedShift=10; M=LengthOfDisplayedShift; T[{u_,v_}]:=If[u<1/2,{2u,v/2},{2u-1,(v+1)/2}]; Ti[{u_,v_}]:=If[v<1/2,{u/2,2v},{(u+1)/2,2v-1}]; f[{u_,v_}]:=If[u<1/2,0,1]; d[s_]:=Delete[s,1]; S[{u_,v_}]:={Reverse[d[Map[f,NestList[Ti,{u,v},M]]]],"|",Map[f,NestList[T,{u,v},M]]} Si[{x_,comma_,y_}]:={Sum[y[[k]] 2^(-k),{k,Length[y]}],Sum[x[[Length[x]-k]] 2^(-k-1),{k,0,Length[x]-1}]}; AnimateShift[{u_,v_},n_]:=Module[{a=u,b=v},Print["\n"];Do[Print[{a,b},"\t\t",Flatten[S[{a,b}]]]; {a,b}=T[{a,b}],{n}]]; AnimateShift[{7/16,0},8] ``` Symbolica. Frank Zizzas Symbolica ( local copy) allows to compute objects in symbolic dynamics. If you put put it into your Mathemaitca library of Standard packages like .../AddOns/StandardPackages/DiscreteMath/Symbolica.m you can run things like ```Needs["DiscreteMath`Symbolica`"]; G=MakeGraph[{{1,1,1},{1,0,1},{0,1,1}}] Show[G] NEntropy[G] ZetaFunction[G] ``` The Golden Ratio subshift ```Needs["DiscreteMath`Symbolica`"]; Show[MakeGraph[ForbiddenWords[{1,1}]]] ``` Mandelbrot set ```M=Compile[{x,y},Block[{z=x+I y,k=0},While[Abs[z]<2&&k<50,z=z^2+x+I y;++k];k]]; DensityPlot[50-M[x,y],{x,-2.,1.},{y,-1.5,1.5},PlotPoints->200,Mesh->False]; ``` Try this out where zk+c is replaced by zk+c for k>2. Julia set or rather the Prisoner set Kc. The parameter c=-0.12+0.74 i is the Douady rabbit. ```J=Compile[{u,v},Block[{z=u+I*v,k=0},While[Abs[z]<99&&k<50,z=z^2-.12+.74I;++k];k]]; DensityPlot[50-J[u,v],{u,-1,1},{v,-1,1},PlotPoints->200,Mesh->False]; ``` Make a Julia set GIF movie. We change the parameter c along the unit circle. ```F[c_]:=Module[{},J=Compile[{u,v},Block[{z=u+I*v,k=0},While[Abs[z]<99&&k<50,z=z^2+c;++k];k]]; DensityPlot[50-J[u,v],{u,-2,2},{v,-2,2},ColorFunction->Hue,Background->GrayLevel[0.0], PlotPoints->400,Frame->False,Mesh->False]]; Export["juliamovie.gif",Table[F[Exp[2Pi k I/60]], {k,60}],"GIF",ConversionOptions ->{"Loop"->True},ImageSize->400] ``` Spin model in porous medium suggested by Monika Schleier in HW 7: the alphabet is {-1,0,1}. A function of neighboring spin configurations determines the probability to have spin 1 in the next generation. ```M = 20; gens = 600; x=Table[Random[Integer,2]-1,{M},{M}]; m[k_]:=Mod[k-1,M]+1; y=x; Do[ Do[s=2-Sum[If[i != 0 || j != 0, x[[m[k + i],m[l + j]]]/(i^2+j^2)^(3/2),0],{i,-2,2},{j,-2,2}]; p = Exp[-s*x[[k,l]]]/(2 Cosh[s*x[[k,l]]]); If[x[[k,l]] == 0, y[[k,l]] = 0, y[[k,l]] = If[Random[Real, 2/p] < 1,-x[[k,l]],x[[k,l]]]],{k,M},{l,M}]; x=y; ListDensityPlot[(x+1)/2,Mesh->False,Frame->False,Axes->False,ColorFunctionScaling -> False],{gens}] ``` Life Note that the following Mathematica implementation of Life does not aim for speed but for clarity of what is done. When run in a Mathematica notebook, you can animate the sequence of graphics. ```SizeOfWorld=50; M=SizeOfWorld; x=Table[Random[Integer,1],{M},{M}]; m[k_]:=Mod[k-1,M]+1; Do[y=Table[x[[m[i+1],j]]+x[[m[i-1],j]]+x[[i,m[j+1]]]+x[[i,m[j-1]]]+ x[[m[i-1],m[j-1]]]+x[[m[i-1],m[j+1]]]+ x[[m[i+1],m[j-1]]]+x[[m[i+1],m[j+1]]],{i,M},{j,M}]; Do[Do[If[x[[i,j]]==1, If[Abs[y[[i,j]]-2.5]>1,x[[i,j]]=0], If[y[[i,j]]==3, x[[i,j]]=1]],{i,M}],{j,M}]; ListDensityPlot[1-x,Mesh->False,Axes->False],{k,10}] ``` Particle Cellular Automaton (see homework). ```M=300; n=3; x0=Table[Random[Integer,n],{M}]; T[x_]:=Mod[x*RotateRight[x]*RotateLeft[x]+1,n]; ListDensityPlot[Reverse[NestList[T,x0,M]], Frame->False, Mesh->False,ColorFunction->Hue,Axes->False] ``` One dimensional Cellular Automata ```M=100;V=CellularAutomaton[18,Table[Random[Integer,1],{k,M}],M]; ListDensityPlot[Reverse[V],Mesh->False,Frame->False,Axes->False] ``` Turing machines are dynamical systems which can do computations. They had been introduced by Alan Turing in 1936. Its not difficult to realize a Turing machine as a cellular automaton. A challenge is the construction of busy beavers: given n states, how many 1 can the machine write before it stops? Here is a machine of Heiner Marxen who has in 1989 constructed a 5 state machine which produces 4098 ones on the tape before it stops. The orbit of this dynamical system is severeal million timesteps long. The Mathematica implementation below is from 1994, when I taught this course at Caltech. It was inspired by an article by Roman Maeder in the Mathematica journal (Volume 3, Summer 1993), who claims there in the abstract that "every computer science student has at one point written a Turing Machine simulator". So, here is mine, written as a dynamical system: Warning: the busy beaver needs several hours to eat! If you do not have Mathematica, you find an other Turing machine implementations of mine here. ```F[{state_,symbol_,position_}]:=Module[{newsym,move,newstate}, If[state==1 && symbol==b, newsym=m ; move=-1; newstate=2]; If[state==1 && symbol==m, newsym=m ; move=-1; newstate=1]; If[state==2 && symbol==b, newsym=m ; move= 1; newstate=3]; If[state==2 && symbol==m, newsym=m ; move= 1; newstate=2]; If[state==3 && symbol==b, newsym=m ; move=-1; newstate=1]; If[state==3 && symbol==m, newsym=m ; move= 1; newstate=4]; If[state==4 && symbol==b, newsym=m ; move=-1; newstate=1]; If[state==4 && symbol==m, newsym=m ; move= 1; newstate=5]; If[state==5 && symbol==b, newsym=m ; move= 1; newstate=0]; If[state==5 && symbol==m, newsym=b ; move= 1; newstate=3]; {newstate,newsym,position+move}]; T[{state_,band_,position_}]:=Module[{newstate,newsym,newpos,newband=band}, {newstate,newsym,newpos}=F[{state,band[[position]],position}]; newband[[position]]=newsym; If[newpos>Length[band],newband=Append[newband,b]]; If[newpos<1,newband=Prepend[newband,b];newpos++]; {newstate,newband,newpos}]; Orbit[{state_,band_,position_}]:=Module[{z={state,band,position},k=0}, While[Not[state==0],z=T[z];k++;If[Mod[k,10]==0,Print[k]]];Print[z]]; Orbit[{1,{b},1}] ``` Condition for stability of periodic orbit of period 2. The following three lines verify that a period 2 orbit is stable if l is bigger then R2 or l is smaller then R1+R2, where l is the length of the orbit and R1,R2 are the radii of curvature at the two impact points. We just compute the trace of the Jacobean of T2 and see when its absolute value is smaller than 2 ```B1={{1,l},{0,1}}; B2={{1,0},{-2/R1,1}}; B3={{1,l},{0,1}}; B4={{1,0},{-2/R2,1}}; {Simplify[Tr[B4.B3.B2.B1] ==2], Simplify[Tr[B4.B3.B2.B1] ==-2]} ``` Proof of integrability of elliptic billiard: the next lines prove that the product between the distances of a ray to the two focal points is conserved in an ellipse. ```F1 = { Sqrt[a^2-b^2],0}; (* first focal point *) F2 = {-Sqrt[a^2-b^2],0}; (* second focal point *) s[t_]:={ a*Cos[t],b*Sin[t]}; (* ellipse *) n[t_]:={ b*Cos[t],a*Sin[t]}; (* normal *) v[t_]:={-a*Sin[t],b*Cos[t]}; (* tangent *) A0=s[t]; (* impact point *) B1=s[t]+n[t]+r*v[t]; (* previous impact point *) B2=s[t]+n[t]-r*v[t]; (* future impact point *) cross[A_,B_]:=A[[1]]*B[[2]]-A[[2]]*B[[1]]; DistancePointLine[P_,A_,B_]:=cross[(P-A),(B-A)]/Sqrt[(B-A).(B-A)]; d1 = DistancePointLine[F1,A0,B1]; d2 = DistancePointLine[F2,A0,B1]; e1 = DistancePointLine[F1,A0,B2]; e2 = DistancePointLine[F2,A0,B2]; Simplify[d1*d2-e1*e2] == 0 ``` Menger Sponge: (be careful by running this with n>3, the number of cubes grows exponentially). ```Menger[0,{i_,j_,k_}]:={Cuboid[{i,j,k}]}; Menger[n_,{i_,j_,k_}]:=Module[{s={}},Do[If[Not[u==v==2 || v==w==2 || u==w==2], s=Union[s,Menger[n-1,3^(n-1)*{u,v,w}+{i,j,k}]]],{u,3},{v,3},{w,3}];s]; Show[Graphics3D[Menger[2,{0,0,0}]],Boxed->False,Axes->False]; ``` The Lorentz attractor: ```Lorentz=NDSolve[{x'[t]==10(y[t]-x[t]),y'[t]==-x[t] z[t]+28x[t]-y[t], z'[t]==x[t]*y[t]-8z[t]/3,x[0]==z[0]==0,y[0]==.3},{x,y,z},{t,0,40},MaxSteps->5000]; ParametricPlot3D[Evaluate[{x[t],y[t],z[t]}/. Lorentz],{t,0,40}, PlotPoints->5000]; ``` Class demo: Drawing orbits of differential equations ```S=NDSolve[{x'[t]==y[t],y'[t]==-x[t]-(x[t]^2-1)*y[t], x[0]==0,y[0]==0.2},{x,y},{t,0,30}]; ParametricPlot[Evaluate[{x[t], y[t]}] /. S[[1]], {t,0,30}] ``` Class demo: The riddle of the first lecture ```T[x_]:=4 x (1-x); A=NestList[T,0.3,100]; S[x_]:=4 x -4x^2; B=NestList[S,0.3,100]; ListPlot[A-B] ``` Class demo: Drawing bifurcation diagram ```g=Compile[{c},({c,#} &) /@ Union[Drop[NestList[c*#(1-#) &,0.3,500],400]]]; s=Table[g[c],{c,2.5,4,0.0001}]; ListPlot[Flatten[s,1],PlotStyle->{AbsolutePointSize[.0001]}] ``` Class demo Cobweb construction (this three line version is a shortened version of a cobweb Mathematica code of Monika Schleier turned in with homework). ```T[x_]:=4 x(1-x); o={{0.3,T[0.3]}}; p={{0.3,0},{0.3,T[0.3]}}; Do[l=Last[Last[o]]; o=Append[o,{l,l}]; o=Append[o,{l,T[l]}],{8}]; Show[Plot[{T[x],x},{x,0,1}],Graphics[{Line[p],Line[o]}]]; ``` Class demo: This is the proof that the Tent map S and Ulam map T are conjugated. You have to run the 3 lines both with c=1 and c=-1. The proof makes use of the fact that two continuous functions f(x),g(x) on the interval [0,1] are the same, if their derivatives f',g' exist and are the same on (1/2,1] and on [0,1/2) and if f(1/2)=g(1/2). ```T[x_]:=4x(1-x); U[x_]:=(1/2)-ArcSin[1-2x]/Pi; c=1; S[x_]:=1-c 2 (x-1/2); f[x_]:=U[T[x]]; g[x_]:=S[U[x]]; A=f'[x]; B=g'[x]; Simplify[A^2==B^2] && f[0.5]==g[0.5] ``` Do you believe this "proof" done by a machine? If not, you can check the identity by hand, but it will take you longer. Zoom into the Henon attractor ```T[{x_,y_}]:=N[{1-1.4*x^2+y,0.3*x}]; S[p_]:=T[p];c=Last[NestList[S,{0,0},300]]; F=NestList[S,c,100000]; G[R_]:=ListPlot[F,PlotRange->R,PlotStyle->{PointSize[10^(-6)]}]; Do[G[{{c[[1]]-2^k,c[[1]]+2^k},{c[[2]]-2^k,c[[2]]+2^k}}],{k,-6,0}]; ``` Computation of stable and unstable manifolds in Henon map according to Francescini and Russo ```(* stable and unstable manifolds in Henon map a la Francescini and Russo *) (* mathematica implementation: Oliver knill *) a=1.4;b=0.3; T[{x_,y_}]:={1-a*x^2+y,b*x}; s=Solve[T[{x,y}]=={x,y},{x,y}]; {x0,y0}={x,y} /. s[[2]] H={{-2a*x0,1},{b,0}}; l=Eigenvalues[H][[1]] f[x_]:=Module[{n=Length[x]},Append[x,a*Sum[x[[j+1]]*x[[n-j+1]],{j,n-1}]/(-l^n+b l^(-n)-2a x0)]]; xx=Module[{y={x0,1}},Do[y=f[y],{8}];y]; x[t_]:=xx[[1]]+Sum[xx[[k+1]]*t^k,{k,1,Length[xx]-1}]; y[t_]:=b*x[t/l]; p[t_]:=Module[{n=7+Floor[Log[Abs[t]]/Log[Abs[l]]],q}, q={x[t/l^n],y[t/l^n]};Do[q=T[q],{n}];q]; S1=ParametricPlot[{p[t]},{t,-10,10},PlotPoints->1000, PlotRange->{{-1,1.5},{-0.5,0.5}}]; Ti[{x_,y_}]:={y/b,x-1+a y^2/b^2}; l=Eigenvalues[H][[2]]; f[x_]:=Module[{n=Length[x]},Append[x,a*Sum[x[[j+1]]*x[[n-j+1]],{j,n-1}]/(-l^n+b l^(-n)-2a x0)]]; xx=Module[{y={x0,1}},Do[y=f[y],{8}];y]; x[t_]:=xx[[1]]+Sum[xx[[k+1]]*t^k,{k,1,Length[xx]-1}]; y[t_]:=b*x[t/l]; p[t_]:=Module[{n=7+Floor[Log[Abs[t]]/Log[Abs[l]]],q}, q={x[t*l^n],y[t*l^n]};Do[q=Ti[q],{n}];q]; S2=ParametricPlot[{p[t]},{t,-10,10},PlotPoints->1000, PlotRange->{{-1,1.5},{-0.5,0.5}}]; Show[{S1,S2,Show[Graphics[{PointSize[0.03],Point[{x0,y0}]}]]},AspectRatio->1] ```
 WWW: Java, flash, javascript:
 URL: Software collection links:
 Math118r, Dynamical systems, Spring 2005, Oliver Knill, knill@math.harvard.edu. Department of Mathematics, Faculty of Art and Sciences, Harvard University, Background music credit: "Barocco", by "Rondo Veneziano" under the lead of Gian Piero Reverberi.