{VERSION 5 0 "IBM INTEL NT" "5.0" } {USTYLETAB {CSTYLE "Maple Input" -1 0 "Courier" 0 1 255 0 0 1 0 1 0 0 1 0 0 0 0 1 }{CSTYLE "2D Math" -1 2 "Times" 0 1 0 0 0 0 0 0 2 0 0 0 0 0 0 1 }{CSTYLE "2D Input" 2 19 "" 0 1 255 0 0 1 0 0 0 0 0 0 0 0 0 1 } {CSTYLE "" -1 256 "" 1 16 0 0 0 0 0 1 0 0 0 0 0 0 0 0 }{CSTYLE "" -1 257 "" 1 14 0 0 0 0 0 0 1 0 0 0 0 0 0 0 }{CSTYLE "" -1 258 "" 0 1 0 0 0 0 0 0 1 0 0 0 0 0 0 0 }{CSTYLE "" -1 259 "" 1 14 0 0 0 0 0 0 1 0 0 0 0 0 0 0 }{CSTYLE "" -1 260 "" 0 1 0 0 0 0 0 0 1 0 0 0 0 0 0 0 } {CSTYLE "" -1 261 "" 0 1 0 0 0 0 0 0 1 0 0 0 0 0 0 0 }{PSTYLE "Normal " -1 0 1 {CSTYLE "" -1 -1 "Times" 1 12 0 0 0 1 2 2 2 2 2 2 1 1 1 1 }1 1 0 0 0 0 1 0 1 0 2 2 0 1 }{PSTYLE "Normal" -1 256 1 {CSTYLE "" -1 -1 "Times" 1 16 0 0 0 1 2 1 2 2 2 2 1 1 1 1 }1 1 0 0 0 0 1 0 1 0 2 2 0 1 }{PSTYLE "" 0 257 1 {CSTYLE "" -1 -1 "" 0 1 0 0 0 0 0 0 1 0 0 0 0 0 0 0 }0 0 0 -1 -1 -1 0 0 0 0 0 0 -1 0 }} {SECT 0 {EXCHG {PARA 0 "" 0 "" {TEXT 256 17 "Fractal Dimension" }} {PARA 0 "" 0 "" {TEXT -1 0 "" }}{PARA 0 "" 0 "" {TEXT -1 79 "The analy sis of Poincare plots for the Duffing oscillator which is explored in \+ " }{TEXT 19 11 "Duffing.mws" }{TEXT -1 475 " indicates that for certai n parameter ranges the Poincare section represents an accumulation of \+ points. Our aim is to calculate the fractal dimension of this strange \+ attractor. We introduce a simple box counting algorithm in order to ca lculate the capacity dimension, and we test it on a fractal structure \+ whose dimension can be derived from the algorithm that prescribes it. \+ First, we have to construct an algorithm that generates this fractal, \+ which is called a Koch curve." }}{PARA 0 "" 0 "" {TEXT -1 0 "" }} {PARA 0 "" 0 "" {TEXT 257 10 "Koch curve" }}{PARA 0 "" 0 "" {TEXT -1 0 "" }}{PARA 0 "" 0 "" {TEXT -1 269 "A Koch curve can be generated by \+ an iterative algorithm. A basic step consists of taking a unit line, c utting it in three, and replacing the middle segment by an equilateral triangle without base, thereby extending the length of the original e ntire line segment by 4/3." }}{PARA 0 "" 0 "" {TEXT -1 0 "" }}{PARA 0 "" 0 "" {TEXT -1 95 "Suppose we define the curve by its vertices. We i mplement the algorithm by the following steps:" }}{PARA 0 "" 0 "" {TEXT -1 57 "1) find the points that cut the line into three segments; " }}{PARA 0 "" 0 "" {TEXT -1 112 "2) find the midpoint, and rotate it \+ to horizontal orientation, displace the y-coordinate by the required a mount;" }}{PARA 0 "" 0 "" {TEXT -1 64 "3) rotate the shifted midpoint \+ back to the original orientation;" }}{PARA 0 "" 0 "" {TEXT -1 55 "4) a ssemble a list that contains all required vertices." }}{PARA 0 "" 0 " " {TEXT -1 0 "" }}{PARA 0 "" 0 "" {TEXT -1 91 "This may not be the mos t efficient implementation, but it is starightforward to understand." }}{PARA 0 "" 0 "" {TEXT -1 0 "" }}}{EXCHG {PARA 0 "> " 0 "" {MPLTEXT 1 0 73 "restart; with(stats): # statistics package is used for least-s quares fit." }}}{EXCHG {PARA 0 "> " 0 "" {MPLTEXT 1 0 77 "Koch:=proc(L L) local xi,yi,xf,yf,x1,y1,x2,y2,xh,yh,phi,cf,sf,x3,y3,xhp,yhp,d;" }} {PARA 0 "> " 0 "" {MPLTEXT 1 0 55 "xi:=LL[1][1]: yi:=LL[1][2]: xf:=LL[ 2][1]: yf:=LL[2][2]:" }}{PARA 0 "> " 0 "" {MPLTEXT 1 0 49 "x1:=evalf(x i+(xf-xi)/3): y1:=evalf(yi+(yf-yi)/3):" }}{PARA 0 "> " 0 "" {MPLTEXT 1 0 53 "x2:=evalf(xi+2*(xf-xi)/3): y2:=evalf(yi+2*(yf-yi)/3):" }} {PARA 0 "> " 0 "" {MPLTEXT 1 0 43 "xh:=evalf((xf+xi)/2): yh:=evalf((yf +yi)/2):" }}{PARA 0 "> " 0 "" {MPLTEXT 1 0 32 "phi:=evalf(arctan(yf-yi ,xf-xi));" }}{PARA 0 "> " 0 "" {MPLTEXT 1 0 27 "cf:=cos(phi): sf:=sin( phi):" }}{PARA 0 "> " 0 "" {MPLTEXT 1 0 73 "xhp:=xh*cf+yh*sf: yhp:=-xh *sf+yh*cf: d:=evalf(sqrt((xf-xi)^2+(yf-yi)^2)):" }}{PARA 0 "> " 0 "" {MPLTEXT 1 0 28 "yhp:=evalf(yhp+d/6*sqrt(3)):" }}{PARA 0 "> " 0 "" {MPLTEXT 1 0 51 "x3:=evalf(xhp*cf-yhp*sf): y3:=evalf(xhp*sf+yhp*cf):" }}{PARA 0 "> " 0 "" {MPLTEXT 1 0 43 "[LL[1],[x1,y1],[x3,y3],[x2,y2],LL [2]]: end:" }}}{EXCHG {PARA 0 "> " 0 "" {MPLTEXT 1 0 20 "Koch([[0,0],[ 0,1]]);" }}}{EXCHG {PARA 0 "> " 0 "" {MPLTEXT 1 0 20 "Koch([[0,0],[1,0 ]]);" }}}{EXCHG {PARA 0 "> " 0 "" {MPLTEXT 1 0 20 "Koch([[0,0],[1,1]]) ;" }}}{EXCHG {PARA 0 "> " 0 "" {MPLTEXT 1 0 39 "plot(%,scaling=constra ined,axes=boxed);" }}}{EXCHG {PARA 0 "" 0 "" {TEXT -1 131 "Now we need a driver: we want to take a list of points, and for each pair of subs equent points we want to apply the Koch procedure." }}}{EXCHG {PARA 0 "> " 0 "" {MPLTEXT 1 0 14 "DriveK:=(LL)->" }}{PARA 0 "> " 0 "" {MPLTEXT 1 0 51 "[seq( op(Koch([LL[i],LL[i+1]])),i=1..nops(LL)-1 )]:" }}}{EXCHG {PARA 0 "> " 0 "" {MPLTEXT 1 0 22 "DriveK([[0,0],[1,0]]);" } }}{EXCHG {PARA 0 "> " 0 "" {MPLTEXT 1 0 10 "DriveK(%);" }}}{EXCHG {PARA 0 "> " 0 "" {MPLTEXT 1 0 39 "plot(%,scaling=constrained,axes=box ed);" }}}{EXCHG {PARA 0 "" 0 "" {TEXT -1 230 "The Koch curve is an imp ortant exercise in order to understand fractal dimension, because the \+ increase in length of the figure after each step is obvious: a factor \+ of 4/3 applies for each iteration. Let us verify this numerically:" }} }{EXCHG {PARA 0 "> " 0 "" {MPLTEXT 1 0 94 "Length:=LL->evalf(add(sqrt( (LL[i+1][1]-LL[i][1])^2+(LL[i+1][2]-LL[i][2])^2),i=1..nops(LL)-1));" } }}{EXCHG {PARA 0 "> " 0 "" {MPLTEXT 1 0 28 "Length(Koch([[0,0],[0,1]]) );" }}}{EXCHG {PARA 0 "> " 0 "" {MPLTEXT 1 0 35 "Length((DriveK@@2)([[ 0,0],[0,1]]));" }}}{EXCHG {PARA 0 "> " 0 "" {MPLTEXT 1 0 15 "evalf((4/ 3)^2);" }}}{EXCHG {PARA 0 "> " 0 "" {MPLTEXT 1 0 64 "plot((DriveK@@5)( [[0,0],[1,0]]),scaling=constrained,axes=boxed);" }}}{EXCHG {PARA 0 "" 0 "" {TEXT -1 53 "The capacity dimension of the Koch curve is given as :" }}}{EXCHG {PARA 0 "> " 0 "" {MPLTEXT 1 0 29 "D_Koch:=evalf(log(4)/l og(3));" }}}{EXCHG {PARA 0 "" 0 "" {TEXT -1 53 "Why is this so? Do som e reading on fractal dimension." }}{PARA 0 "" 0 "" {TEXT -1 0 "" }} {PARA 0 "" 0 "" {TEXT -1 59 "How rapidly does the number of points gro w with each layer?" }}}{EXCHG {PARA 0 "> " 0 "" {MPLTEXT 1 0 45 "seq(n ops((DriveK@@i)([[0,0],[1,0]])),i=1..5);" }}}{EXCHG {PARA 0 "" 0 "" {TEXT -1 183 "In order to implement the box counting method we need to pick a pixel size resolution and connect the points. In a plotting pr ogram this is carried out by the graphics engine (in the " }{TEXT 19 4 "plot" }{TEXT -1 109 " command above the default method is to connec t the points by a straight line). Redraw the above graph using " } {TEXT 19 24 "style=point,symbol=point" }{TEXT -1 239 " to observe what the plot engine has done for us. Of course, we can simply iterate the curve to a sufficiently high degree such that at a given smallest box size resolution (e.g., pixel size for a graphics format) all points a re connected." }}{PARA 0 "" 0 "" {TEXT -1 0 "" }}}{EXCHG {PARA 0 "> " 0 "" {MPLTEXT 1 0 31 "K7:=(DriveK@@7)([[0,0],[1,0]]):" }}}{EXCHG {PARA 0 "> " 0 "" {MPLTEXT 1 0 15 "K8:=DriveK(K7):" }}}{EXCHG {PARA 0 "" 0 "" {TEXT -1 88 "The two data sets K7 and K8 allow us to experimen t below. Let us choose the smaller one:" }}}{EXCHG {PARA 0 "> " 0 "" {MPLTEXT 1 0 8 "K_u:=K7:" }}}{EXCHG {PARA 0 "> " 0 "" {MPLTEXT 1 0 17 "k_max:=nops(K_u);" }}}{EXCHG {PARA 0 "" 0 "" {TEXT -1 613 "Choose a m aximum value for the resolution layer (smallest box size) at which the plane with the curve enclosed will be explored. The boxes are squares in our case, because the points that make up the Koch curve are two-d imensional. The total area in which the fractal resides is a rectangle : we can account for that. We are setting up an algorithm analogous to a histogram calculation except that we don't care how many times a gi ven bin has been occupied, we just determine whether it has been hit, \+ and then add up how many hits we get at a certain level of resolution. The baseline for a bin (square) is called " }{TEXT 19 11 "epsilon=1/n " }{TEXT -1 105 ", where n is the resolution layer, and we need to loo k at the log of the total number of box hits versus " }{TEXT 19 14 "lo g(1/epsilon)" }{TEXT -1 71 "; the slope of this graph gives the fracta l dimension (in the limit of " }{TEXT 19 7 "epsilon" }{TEXT -1 81 " go ing to zero). At the lowest layer of resolution (n=1) we use 10 times \+ 5 boxes." }}}{EXCHG {PARA 0 "> " 0 "" {MPLTEXT 1 0 10 "layer:=15;" }}} {EXCHG {PARA 0 "> " 0 "" {MPLTEXT 1 0 8 "LXY:=[]:" }}{PARA 0 "> " 0 " " {MPLTEXT 1 0 25 "for n from 1 to layer do:" }}{PARA 0 "> " 0 "" {MPLTEXT 1 0 26 "for i from 1 to (10*n) do:" }}{PARA 0 "> " 0 "" {MPLTEXT 1 0 25 "for j from 1 to (5*n) do:" }}{PARA 0 "> " 0 "" {MPLTEXT 1 0 10 "P[i,j]:=0;" }}{PARA 0 "> " 0 "" {MPLTEXT 1 0 3 "od:" }}{PARA 0 "> " 0 "" {MPLTEXT 1 0 4 "od: " }}{PARA 0 "> " 0 "" {MPLTEXT 1 0 14 "Pepsilon:=1/n:" }}{PARA 0 "> " 0 "" {MPLTEXT 1 0 13 " Pboxcount:=0:" }}{PARA 0 "> " 0 "" {MPLTEXT 1 0 92 "for k from 1 to k_ max do: ix:=trunc(10*n*(K_u[k][1]+0.5)); iy:=trunc(5*n*(K_u[k][2]+0.2 5));" }}{PARA 0 "> " 0 "" {MPLTEXT 1 0 20 "if (P[ix,iy]=0) then" }} {PARA 0 "> " 0 "" {MPLTEXT 1 0 12 "P[ix,iy]:=1:" }}{PARA 0 "> " 0 "" {MPLTEXT 1 0 25 "Pboxcount:=Pboxcount + 1:" }}{PARA 0 "> " 0 "" {MPLTEXT 1 0 3 "fi:" }}{PARA 0 "> " 0 "" {MPLTEXT 1 0 77 "od: print(\" For res. \",n,\"find \",Pboxcount,\" hits out of \",5*n*10*n,\" boxes \");" }}{PARA 0 "> " 0 "" {MPLTEXT 1 0 67 "LXY:=[op(LXY), [evalf(log10 (1/Pepsilon)),evalf(log10(Pboxcount))]]:" }}{PARA 0 "> " 0 "" {MPLTEXT 1 0 3 "od:" }}}{EXCHG {PARA 0 "> " 0 "" {MPLTEXT 1 0 48 "PL1: =plot(LXY,style=point): plots[display](PL1);" }}}{EXCHG {PARA 0 "" 0 " " {TEXT -1 110 "When the number of boxes is small we have no useful re sults, so let us ignore them in the regression analysis:" }}}{EXCHG {PARA 0 "> " 0 "" {MPLTEXT 1 0 21 "nL:=4; nU:=nops(LXY);" }}}{EXCHG {PARA 0 "> " 0 "" {MPLTEXT 1 0 84 "sol:=fit[leastsquare[[x,y]]]([[seq( LXY[i][1],i=nL..nU)],[seq(LXY[i][2],i=nL..nU)]]);" }}}{EXCHG {PARA 0 " > " 0 "" {MPLTEXT 1 0 28 "FractDim:=coeff(rhs(sol),x);" }}}{EXCHG {PARA 0 "> " 0 "" {MPLTEXT 1 0 67 "PL2:=plot(rhs(sol),x=0.5..1.2,color =blue): plots[display](PL1,PL2);" }}}{EXCHG {PARA 0 "" 0 "" {TEXT 258 11 "Exercise 1:" }}{PARA 0 "" 0 "" {TEXT -1 104 "Push the calculation \+ to higher levels of resolution. What do you find? Then use the larger \+ data sample (" }{TEXT 19 2 "K8" }{TEXT -1 180 ") and repeat the analys is. For a better analysis you want to use a regression fit with uncert ainty estimate of the coefficients. An algorithm for this is provided \+ in the worksheet " }{TEXT 19 16 "DataAnalysis.mws" }{TEXT -1 1 "." }} {PARA 0 "" 0 "" {TEXT -1 0 "" }}}{EXCHG {PARA 0 "" 0 "" {TEXT 259 17 " Duffing Attractor" }}{PARA 0 "" 0 "" {TEXT -1 0 "" }}{PARA 0 "" 0 "" {TEXT -1 0 "" }}{PARA 0 "" 0 "" {TEXT -1 46 "The Duffing oscillator ha s been introduced in " }{TEXT 19 11 "Duffing.mws" }{TEXT -1 625 " and \+ it has been shown that for certain parameter values the limit cycle be comes an interesting structure. Here we carry out a measurement of a c apacity dimension to demonstrate the fractal nature of the attractor, \+ which is therefore called a strange attractor. In maple8 the integrati on in dsolve[numeric] has become fast enough that we can extract usefu l information without having to resort to C-coded routines (such as th e method=dna option, which requires a special library to be loaded). F or maple7 we still recommend the usage of the DNA package; for maple8 \+ it was unavailable yet at the beginning of September, 2002)." }}{PARA 0 "" 0 "" {TEXT -1 0 "" }}{PARA 0 "" 0 "" {TEXT -1 119 "For maple8 (wi thout DNA) we restrict the precision to 13 digits. Going beyond 14 dig its results in a massive slow-down." }}}{EXCHG {PARA 0 "> " 0 "" {MPLTEXT 1 0 20 "restart; Digits:=13;" }}}{EXCHG {PARA 0 "> " 0 "" {MPLTEXT 1 0 43 "unprotect(gamma): with(plots): with(stats):" }}} {EXCHG {PARA 0 "> " 0 "" {MPLTEXT 1 0 32 "gamma:=1/10: eta:=1: delta:= 1/4:" }}}{EXCHG {PARA 0 "> " 0 "" {MPLTEXT 1 0 64 "plot(-1/2*eta*x^2+d elta*x^4/4,x=-3..3,title=\"Potential Energy\");" }}}{EXCHG {PARA 0 "" 0 "" {TEXT -1 47 "We include a periodic driving force (amplitude " } {TEXT 19 2 "F0" }{TEXT -1 24 " and circular frequency " }{TEXT 19 2 "n u" }{TEXT -1 63 ") and start the oscillator from the unstable equilibr ium point:" }}}{EXCHG {PARA 0 "> " 0 "" {MPLTEXT 1 0 33 "nu:=2: F0:=2. 5: DF:=F0*sin(nu*t);" }}}{EXCHG {PARA 0 "> " 0 "" {MPLTEXT 1 0 72 "DE: =diff(x(t),t)=y(t),diff(y(t),t)=-gamma*y(t)+eta*x(t)-delta*x(t)^3+DF; " }}}{EXCHG {PARA 0 "> " 0 "" {MPLTEXT 1 0 59 "IC:=x(0)=0,y(0)=0; #sta rt at the unstable equilibrium point" }}}{EXCHG {PARA 0 "" 0 "" {TEXT -1 71 "In case you haven't looked at Duffing.mws: here is a sample tra jectory:" }}}{EXCHG {PARA 0 "> " 0 "" {MPLTEXT 1 0 63 "#sol:=dsolve(\{ DE,IC\},\{x(t),y(t)\},numeric,output=listprocedure):" }}}{EXCHG {PARA 0 "> " 0 "" {MPLTEXT 1 0 38 "#X:=subs(sol,x(t)): Y:=subs(sol,y(t)):" } }}{EXCHG {PARA 0 "> " 0 "" {MPLTEXT 1 0 57 "#plot('X(t)',t=0..100,nump oints=2000,title=\"Trajectory\");" }}}{EXCHG {PARA 0 "> " 0 "" {MPLTEXT 1 0 74 "#plot(['X(t)','Y(t)',t=0..50],numpoints=1000,title=\" Phase-space Diagram\");" }}}{EXCHG {PARA 0 "" 0 "" {TEXT -1 553 "The P oincare section is obtained by a stroboscopic view at intervals spaced apart by the period of the driving force. An interesting question ari ses: what happens to the Poincare section when we change the phase at \+ which we record the phase-space point? It turns out that one gets a si milar structure. It is possible to record these different sequences, a nd run an animation. We do this before we proceed to compute the fract al dimension, just to convince ourselves that there is no loss of gene rality by fixing the phase at which we record the section." }}}{EXCHG {PARA 0 "> " 0 "" {MPLTEXT 1 0 30 "tau:=evalf(2*Pi/nu); LL:='LL':" }}} {EXCHG {PARA 0 "> " 0 "" {MPLTEXT 1 0 70 "t:='t': sol:=dsolve(\{DE,IC \},\{x(t),y(t)\},numeric,output=listprocedure):" }}}{EXCHG {PARA 0 "> \+ " 0 "" {MPLTEXT 1 0 37 "X:=subs(sol,x(t)): Y:=subs(sol,y(t)):" }}} {EXCHG {PARA 0 "> " 0 "" {MPLTEXT 1 0 13 "dtau:=tau/10;" }}}{EXCHG {PARA 0 "> " 0 "" {MPLTEXT 1 0 112 "for i from 1 to 2000 do: for j fro m 1 to 10 do: t_i:=evalf(i*tau+(j-1)*dtau); LL||j[i]:=[X(t_i),Y(t_i)]: od: od:" }}}{EXCHG {PARA 0 "> " 0 "" {MPLTEXT 1 0 22 "for j from 1 to 10 do:" }}{PARA 0 "> " 0 "" {MPLTEXT 1 0 134 "PL[j]:=plot([seq(LL||j[ i],i=1..2000)],style=point,symbol=point,axes=boxed,view=[-5..5,-8..8], title=\"Poincare Section\",color=blue): od:" }}}{EXCHG {PARA 0 "> " 0 "" {MPLTEXT 1 0 46 "display([seq(PL[j],j=1..10)],insequence=true);" }} }{EXCHG {PARA 0 "" 0 "" {TEXT 260 11 "Exercise 2:" }}{PARA 0 "" 0 "" {TEXT -1 147 "Explore the behaviour of this Duffing dust for small var iations in the drive amplitude and frequency. What is your guess for t he fractal dimension?" }}{PARA 0 "" 0 "" {TEXT -1 0 "" }}}{EXCHG {PARA 0 "" 0 "" {TEXT -1 162 "Now we collect our data sample for fract al dimensions analysis. We record in excess of 100,000 points. This ta kes some time (under 10 minutes on a > 1GHz machine)" }}}{EXCHG {PARA 0 "> " 0 "" {MPLTEXT 1 0 70 "t:='t': sol:=dsolve(\{DE,IC\},\{x(t),y(t) \},numeric,output=listprocedure):" }}{PARA 0 "> " 0 "" {MPLTEXT 1 0 37 "X:=subs(sol,x(t)): Y:=subs(sol,y(t)):" }}}{EXCHG {PARA 0 "> " 0 " " {MPLTEXT 1 0 66 "for i from 1 to 200000 do: t_i:=i*tau; LL[i]:=[X(t_ i),Y(t_i)]: od:" }}}{EXCHG {PARA 0 "" 0 "" {TEXT -1 0 "" }}{PARA 0 "" 0 "" {TEXT -1 145 "To analyze the section we adapt the box counting me thod. We should perform the analysis with varying level of layers and \+ increasing point number " }{TEXT 19 5 "k_max" }{TEXT -1 212 " (the max imum possible value being limited to the pre-computed maximum). We sta rt with 200,000 points, and can repeat the exploration with larger sam ples (requires a re-run of the previous loop for longer times)." }} {PARA 0 "" 0 "" {TEXT -1 0 "" }}{PARA 0 "" 0 "" {TEXT -1 157 "The adju stment to the box counting method concerns a shift in the x- and y-coo rdinates by 5 and 8 units respectively, and a choice of resolution (re ctangle)." }}}{EXCHG {PARA 0 "> " 0 "" {MPLTEXT 1 0 14 "k_max:=200000; " }}}{EXCHG {PARA 0 "> " 0 "" {MPLTEXT 1 0 10 "layer:=20;" }}}{EXCHG {PARA 0 "> " 0 "" {MPLTEXT 1 0 8 "LXY:=[]:" }}{PARA 0 "> " 0 "" {MPLTEXT 1 0 25 "for n from 1 to layer do:" }}{PARA 0 "> " 0 "" {MPLTEXT 1 0 26 "for i from 1 to (10*n) do:" }}{PARA 0 "> " 0 "" {MPLTEXT 1 0 26 "for j from 1 to (12*n) do:" }}{PARA 0 "> " 0 "" {MPLTEXT 1 0 10 "P[i,j]:=0;" }}{PARA 0 "> " 0 "" {MPLTEXT 1 0 3 "od:" }}{PARA 0 "> " 0 "" {MPLTEXT 1 0 4 "od: " }}{PARA 0 "> " 0 "" {MPLTEXT 1 0 14 "Pepsilon:=1/n:" }}{PARA 0 "> " 0 "" {MPLTEXT 1 0 13 " Pboxcount:=0:" }}{PARA 0 "> " 0 "" {MPLTEXT 1 0 80 "for k from 1 to k_ max do: ix:=trunc(n*(LL[k][1]+5)); iy:=trunc(n*(LL[k][2]+8));" }} {PARA 0 "> " 0 "" {MPLTEXT 1 0 20 "if (P[ix,iy]=0) then" }}{PARA 0 "> \+ " 0 "" {MPLTEXT 1 0 12 "P[ix,iy]:=1:" }}{PARA 0 "> " 0 "" {MPLTEXT 1 0 25 "Pboxcount:=Pboxcount + 1:" }}{PARA 0 "> " 0 "" {MPLTEXT 1 0 3 "f i:" }}{PARA 0 "> " 0 "" {MPLTEXT 1 0 78 "od: print(\"For res. \",n,\"f ind \",Pboxcount,\" hits out of \",12*n*10*n,\" boxes\");" }}{PARA 0 " > " 0 "" {MPLTEXT 1 0 67 "LXY:=[op(LXY), [evalf(log10(1/Pepsilon)),eva lf(log10(Pboxcount))]]:" }}{PARA 0 "> " 0 "" {MPLTEXT 1 0 3 "od:" }}} {EXCHG {PARA 0 "> " 0 "" {MPLTEXT 1 0 41 "PL1:=plot(LXY,style=point): \+ display(PL1);" }}}{EXCHG {PARA 0 "" 0 "" {TEXT -1 56 "The numbers and \+ the graph tell us two things (at least):" }}{PARA 0 "" 0 "" {TEXT -1 0 "" }}{PARA 0 "" 0 "" {TEXT -1 98 "1) early on for relatively large b oxes the box count increases with a slope of approximately 1.6; " }} {PARA 0 "" 0 "" {TEXT -1 115 "2) later the slope settles to about 0.6, and epsilon is beginning to be small, i.e., we should trust those res ults." }}{PARA 0 "" 0 "" {TEXT -1 125 "3) clearly, for the last few po ints we are hitting the limit: the number of boxes becomes comparable \+ to the number of points." }}{PARA 0 "" 0 "" {TEXT -1 0 "" }}{PARA 0 " " 0 "" {TEXT -1 26 "So let us extract a slope:" }}}{EXCHG {PARA 0 "> \+ " 0 "" {MPLTEXT 1 0 24 "nL:=10; nU:=nops(LXY)-2;" }}}{EXCHG {PARA 0 "> " 0 "" {MPLTEXT 1 0 84 "sol:=fit[leastsquare[[x,y]]]([[seq(LXY[i][1], i=nL..nU)],[seq(LXY[i][2],i=nL..nU)]]);" }}}{EXCHG {PARA 0 "> " 0 "" {MPLTEXT 1 0 58 "PL2:=plot(rhs(sol),x=0..1.5,color=blue): display(PL1, PL2);" }}}{EXCHG {PARA 0 "> " 0 "" {MPLTEXT 1 0 28 "FractDim:=coeff(rh s(sol),x);" }}}{EXCHG {PARA 0 "" 0 "" {TEXT -1 111 "So it appears as i f the attractor is indeed made of dust, i.e. a set of points of dimens ionality less than one." }}{PARA 0 "" 0 "" {TEXT -1 0 "" }}{PARA 0 "" 0 "" {TEXT -1 429 "If one runs for a million points one can go reasona bly to layer 20, and extract a slope of about 0.5-0.6. The ODE was sol ved using 13 digits, and one has to worry about the inaccuracies that \+ will accumulate in very long time sequences. From this perspective the measurement of the capcity dimension is problematic with the present \+ method. Various methods have been devised to accelerate the convergenc e of the capacity dimension." }}{PARA 0 "" 0 "" {TEXT -1 0 "" }}{PARA 0 "" 0 "" {TEXT -1 59 "Look at research articles in the journal Physic s Letters A:" }}{PARA 0 "" 0 "" {TEXT -1 127 "Phys. Lett. A 141, p. 38 6; A 148, p. 63; A 151, p.43 have O(N*log(N)) accurate methods, and ev en O(N) methods apparently exist." }}{PARA 0 "" 0 "" {TEXT -1 0 "" }} {PARA 257 "" 0 "" {TEXT -1 11 "Exercise 3:" }}{PARA 0 "" 0 "" {TEXT -1 145 "Experiment with different data sets, in particular, confirm th e finding of the small-epsilon regime by using a larger sample of attr actor points." }}{PARA 0 "" 0 "" {TEXT -1 0 "" }}{PARA 0 "" 0 "" {TEXT -1 0 "" }}{PARA 0 "" 0 "" {TEXT -1 51 "It is also illustrative t o look at histograms over " }{TEXT 19 1 "x" }{TEXT -1 4 " or " }{TEXT 19 1 "y" }{TEXT -1 1 ":" }}}{EXCHG {PARA 0 "> " 0 "" {MPLTEXT 1 0 47 " Nh:=1000: for i from 1 to Nh do: Pxh[i]:=0: od:" }}{PARA 0 "> " 0 "" {MPLTEXT 1 0 186 "Nbin:=0: for k from 1 to 200000 do: ix:=trunc(Nh/10* (LL[k][1]+5)); if ix>0 and ix " 0 "" {MPLTEXT 1 0 90 "plot([seq([i,Pxh[i]],i=1..N h)],style=point,symbol=point,axes=boxed,view=[100..900,0..20]);" }}} {EXCHG {PARA 0 "> " 0 "" {MPLTEXT 1 0 90 "plot([seq([i,Pxh[i]],i=1..Nh )],style=point,symbol=point,axes=boxed,view=[100..900,0..10]);" }}} {EXCHG {PARA 0 "> " 0 "" {MPLTEXT 1 0 56 "Nbin:=0: Nh:=1000: for i fro m 1 to Nh do: Pyh[i]:=0: od:" }}{PARA 0 "> " 0 "" {MPLTEXT 1 0 177 "fo r k from 1 to 100000 do: iy:=trunc(Nh/10*(LL[k][2]+7)); if iy>0 and iy " 0 "" {MPLTEXT 1 0 90 "plot([seq([i,Pyh[i]],i=1..Nh)],style=point,symbol=poi nt,axes=boxed,view=[100..900,0..20]);" }}}{EXCHG {PARA 0 "" 0 "" {TEXT -1 584 "The histograms over x or y only, i.e., over position or \+ velocity only are similar in the sense that about 2/3 of the bins are \+ being filled at all, the majority are filled with relatively few point s, and a few bins (we find seven; to reveal them: repeat the plots wit hout a y-range restriction!) get filled very strongly. For an object o f dimensionality greater than one we should not be missing any bins in the 1d histogram. Thus, it is an interesting question to ask whether \+ there are gaps in the histogram. Will these gaps be filled if more poi nts are added to the Poincare section?" }}{PARA 0 "" 0 "" {TEXT -1 0 " " }}{PARA 0 "" 0 "" {TEXT 261 10 "Exercise4:" }}{PARA 0 "" 0 "" {TEXT -1 136 "Scan the histogram to find the gaps. What are the regions of x (or y) that are left out? Find out whether they are filled in eventua lly." }}}{EXCHG {PARA 0 "> " 0 "" {MPLTEXT 1 0 0 "" }}}}{MARK "82" 0 } {VIEWOPTS 1 1 0 3 2 1804 1 1 1 1 }{PAGENUMBERS 0 1 2 33 1 1 }