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]
|
