Chapter 13. Fractal Curves and Dimension

Rovenski Vladimir, Haifa

> restart:

13.1 Sierpinski's Curves

Example 1 . A simple program for plotting with recursion.

> tree := proc(L::algebraic,lev::integer,x0::algebraic, y0::algebraic) local i; global s,p; options remember; s:=s+1; p[s]:=plot([[x0-L,y0+L],[x0,y0],[x0+L,y0+L]]); if lev > 1 then tree(L/2,lev-1,x0-L,y0+L);tree(L/2,lev-1,x0+L,y0+L) fi; RETURN(plots[display]([seq(p[i], i=1..2^lev-1)],scaling=constrained)) end:

> s:=0: tree(100, 3, 0, 0);

>

Replace the angle by a more complicated object

> # p[s]:=plot([x0+L*sqrt(2)*sin(2*t)*cos(t), y0+L*sqrt(2)*sin(2*t)*sin(t), t=-Pi/2..Pi/2], color=green, style=point):

> cactus := proc(L::algebraic,lev::integer,x0::algebraic, y0::algebraic) local i; global s,p; options remember; s:=s+1; p[s]:=plot([x0+L*sqrt(2)*sin(2*t)*cos(t), y0+L*sqrt(2)*sin(2*t)*sin(t), t=-Pi/2..Pi/2], color=green, style=point,symbol=diamond): if lev > 1 then cactus(L/2,lev-1,x0-L,y0+L); cactus(L/2,lev-1,x0+L,y0+L) fi; RETURN(plots[display]([seq(p[i], i=1..2^lev-1)],scaling=constrained)) end:

> s:=0: cactus(100, 3, 0, 0);

Problem 1 .

> serp1:= proc(L::algebraic, lev::integer,x0::algebraic,y0::algebraic) global p,s; options remember; if s=0 then p[0]:=plots[polygonplot]([[x0,y0],[x0+2*L ,y0],[x0,y0+2*L]], style=line) fi; s:=s+1; p[s]:=plots[polygonplot]([[x0,y0+L], [x0+L,y0+L],[x0+L,y0]], color=green); if lev>1 then serp1(L/2,lev-1,x0,y0);serp1(L/2,lev-1,x0+L,y0); serp1(L/2,lev-1,x0,y0+L) fi; RETURN(plots[display]([seq(p[i],i=0..(3^lev-1)/2)],scaling=constrained)) end:

> s:=0: serp1(100, 4, 0, 0);

>

Problem 2 .

> serp2 := proc(a::algebraic,lev::integer,x0::algebraic,y0::algebraic) local i,j; global s,p; options remember; if s=0 then p[0]:=plots[polygonplot]([[x0,y0], [x0,y0+a], [x0+a,y0+a], [x0+a,y0]], style=line) fi; s:=s+1; p[s]:=plots[polygonplot]([[x0+a/3, y0+a/3], [x0+a/3, y0+2*a/3], [x0+2*a/3, y0+2*a/3], [x0+2*a/3, y0+a/3]],color=green); if lev > 1 then for i from 0 to 2 do for j from 0 to 2 do if abs(i-1)+abs(j-1) > 0 then serp2(a/3,lev-1, x0+i*a/3, y0+j*a/3) fi od od fi; RETURN(plots[display]([seq(p[i], i=0..(8^lev-1)/7)], axes=none,scaling=constrained)) end:

> s:=0: serp2(100, 3, 0, 0);

>

13.2 Peano Curves

Problem 1 .

(a)

> restart:

> peano1 := proc(L::algebraic, lev::integer) local A,s,p,i; global Q; Q:=plots[polygonplot]([[L,L/2],[L/2,L],[L,3*L/2], [3*L/2,L]]); s:=1; while s <= lev do A:=L/2^(s-1); p[4]:=plot([[A/2,0], [A/2,A/4]]); p[5]:=plot([[0, A/2], [A/4, A/2]]); p[6]:=plot([[A/2,A/4],[A/4,A/2]], color=white,scaling=constrained); p[0]:=plots[display]([p[6],plottools[scale](Q,1/2,1/2), p[4], p[5]]); for i from 1 to 3 do p[i]:=plottools[rotate](p[0], i*Pi/2) od; Q:=plottools[translate](plots[display]([seq(p[i], i=0..3)]),L,L); s:=s+1 od; RETURN(Q) end:

> peano1(10, 4);

>

(b) We plot the triangular Peano curve by the program peano2(L,lev) .

> peano2 := proc(L, k::algebraic, lev::integer) local f,s,p; global Q; p[0]:=plot([[-L*k/2,0], [-L*k/2,L/2], [L*k/2,L/2], [L*k/2,0]], thickness=1); p[1]:=plots[polygonplot]([[-L,0],[0,L],[L,0]], linestyle=2); Q:=plots[display]([p[0],p[1]]); f:=plottools[transform]((x,y) -> [-x,y]); s:=2; while s < lev+2 do if (s mod 2=0) then p[0]:=plot([[L*(1-(2-k)/2^(s/2)),L/2^(s/2)],[L*(1-(3-k)/2^(s/2+1)),L*(3-k)/2^(s/2+1)]],thickness=1); p[1]:=plot([[L*(1-(2-k)/2^(s/2)),L/2^(s/2)],[L*(1-(2-k)/2^(s/2)),0]],thickness=1,color=white); Q:=plots[display]([p[1],Q,p[0]]) fi; p[2]:=plottools[scale](plottools[rotate](plottools[translate](Q,0,-L),3*Pi/4),1/sqrt(2),1/sqrt(2)); Q:=plots[display]([p[2],f(p[2])],scaling=constrained); s:=s+1 od; RETURN(Q) end: # 0<k<1 defines smoothness of curve.

> peano2(10, 0.7, 6);

(c) Recursion can be used for plotting the image known as the Hilbert's curve.

> peano3 := proc(a,b,c,dab::array,dac::array, lev::integer) local d,e,f,g,h,i,j,k,m,p; global Q; options remember; d[1]:=(a[1]+c[1]-dac[1])/2; d[2]:=(a[2]+c[2]-dac[2])/2; f[1]:=d[1]+dac[1]; f[2]:=d[2]+dac[2]; i[1]:=f[1]+b[1]-a[1]; i[2]:=f[2]+b[2]-a[2]; e[1]:=(a[1]+b[1]-dab[1])/2; e[2]:=(a[2]+b[2]-dab[2])/2; g[1]:=f[1]+e[1]-a[1]; g[2]:=f[2]+e[2]-a[2]; h[1]:=g[1]+dab[1]; h[2]:=g[2]+dab[2]; j[1]:=c[1]+h[1]-f[1]; j[2]:=c[2]+h[2]-f[2]; k[1]:=i[1]-dac[1]; k[2]:=i[2]-dac[2]; m[1]:=e[1]+dab[1]; m[2]:=e[2]+dab[2]; p[1]:=plot([[a[1],a[2]],[f[1],f[2]],[i[1],i[2]],[b[1],b[2]]]); p[2]:=plot([[a[1],a[2]],[b[1],b[2]]],color=white); Q:=plots[display]([p[2],p[1],Q],scaling=constrained); if lev > 0 then peano3(a,d,e,dac,dab,lev-1); peano3(f,g,c,dab,dac,lev-1); peano3(h,i,j,dab,dac,lev-1); dab[1]:=-dab[1]; dab[2]:=-dab[2]; peano3(b,k,m,dac,dab,lev-1); dab[1]:=-dab[1]; dab[2]:=-dab[2] fi; RETURN(Q) end:

> a:=[1,1]: b:=[4,1]: c:=[1,4]: m:=3: N:=2^(m+1)-1: dac:=evalm((c-a)/N): dab:=evalm((b-a)/N): Q:=plot([[a[1], a[2]], [b[1], b[2]]], color=white, thickness=2):

> peano3(a, b, c, dab, dac, m);

Triangle:

> with(plottools): a:=[0,0]: b:=[0,3]: c:=[-3,0]: Q:=plot([[a[1],a[2]],[b[1],b[2]]], color=white):

> m:=3: N:=2^(m+1)-1: dac:=evalm((c-a)/N): dab:=evalm((b-a)/N):

> q[1]:=peano3(a,b,c,dab,dac,m): q[2]:=scale(rotate(translate(q[1],0,-3),Pi/2),4/3,4/3): q[3]:=translate(scale(rotate(q[1], Pi+arcsin(4/5)),5/3,5/3),0,3): plots[display]([q[1], q[2], q[3]],axes=none);

>

13.3 Koch Curves

> restart:

> triad1 := proc(L::algebraic, lev::integer, x0::algebraic,y0::algebraic) global p,s; options remember; s:=s+1; p[s]:=plots[polygonplot]([[x0,y0], [x0,y0+L], [x0+L,y0+L], [x0+L,y0]], color=green); if lev>1 then triad1(L/3, lev-1, x0+L/3, y0-L/3); triad1(L/3,lev-1,x0-L/3,y0+L/3);triad1(L/3, lev-1, x0+L/3, y0+L); triad1(L/3, lev-1, x0+L, y0+L/3) fi; RETURN(plots[display]([seq(p[i], i=1..(4^lev-1)/3)])) end:

> s:=0: triad1(100, 3, 0, 0);

>

> triad2 := proc(a,b, lev::integer) local s,c,d,e,f; global p,Q; options remember; s:=.33; c[1]:=a[1]+(b[1]-a[1])*s; c[2]:=a[2]+(b[2]-a[2])*s; f[1]:=a[1]+(b[1]-a[1])*(1-s); f[2]:=a[2]+(b[2]-a[2])*(1-s); d[1]:=c[1]-(f[2]-c[2]); d[2]:=c[2]+(f[1]-c[1]);e[1]:=d[1]+f[1]-c[1]; e[2]:=d[2]+f[2]-c[2];p[1]:=plot([[c[1],c[2]],[d[1],d[2]],[e[1],e[2]],[f[1],f[2]]]); p[2]:=plot([[c[1],c[2]],[f[1],f[2]]], color=white); Q:=plots[display]([p[2],Q,p[1]], scaling=constrained);if lev>1 then triad2(c,d,lev-1); triad2(d,e,lev-1); triad2(e,f,lev-1) fi; RETURN(Q) end: # 0<s<0.5

> a:=[1,1]: b:=[4,1]: Q:=plot([[a[1],a[2]], [b[1],b[2]]]):

> triad2(a, b, 4);

> with(plottools): setoptions(scaling=constrained, axes=none):

> f:= transform((x,y)->[-x,y]): q||1:=translate(Q,-a[1],-a[2]): q||2:=rotate(f(q||1), Pi/2): q||3:=translate(rotate(q||1, -Pi/2), b[1]-a[1], 0): q||4:=translate(rotate(q||1, Pi), b[1]-a[1], a[1]-b[1]): plots[display]([seq(q||i, i=1..4)],axes=framed);

>

> triad3 := proc(a,b,lev::integer) local s,k,c,d,f,p; global Q; options remember; s:=0.4; k:=2;c[1]:=a[1]+(b[1]-a[1])*s; c[2]:=a[2]+(b[2]-a[2])*s; f[1]:=a[1]+(b[1]-a[1])*(1-s); f[2]:=a[2]+(b[2]-a[2])*(1-s); d[1]:=(a[1]+b[1])/2-k*(f[2]-c[2]); d[2]:=(a[2]+b[2])/2+k*(f[1]-c[1]); p[1]:=plot([[c[1],c[2]],[d[1],d[2]],[f[1],f[2]]]); p[2]:=plot([[c[1],c[2]],[f[1],f[2]]], color=white); Q:=plots[display]([p[2],Q,p[1]], scaling=constrained); if lev>1 then triad3(a,c,lev-1); triad3(c,d,lev-1); triad3(d,f,lev-1); triad3(f,b,lev-1) fi; RETURN(Q) end:

> a:=[1,1]: b:=[4,1]: Q:=plot([[a[1],a[2]],[b[1],b[2]]]):

> triad3(a,b,4);

> setoptions(scaling=constrained, axes=none):

> f:=transform((x,y)->[-x,y]): q||1:=translate(Q,-a[1],-a[2]): q||2:=rotate(f(q||1), 2*Pi/3): q||3:=translate(rotate(q||1,-2*Pi/3), b[1]-a[1],0): plots[display]([seq(q||i, i=1..3)], axes=none);

>

Exercise . Modify programs 13.3.1 and 13.3.2 for getting Koch curves which are inside the triangle and the square.

For obtaining Koch similar curves replace the equilateral triangle in the original

Koch curve by an isosceles triangle or try to find your own "generators".

13.4 Dragon Curve (or Polygon)

> restart: with(plottools): setoptions(scaling=constrained,axes=none):

>

> dragon1 := proc(m::integer, k::algebraic) local b,i; global c,p; p[1]:=plot([[0,0],[1-k,0],[1,k],[1,1]]); c[1]:=circle([1,0],0.04); for i from 1 to m-1 do b:=[2^(i/2)*cos(i*Pi/4),2^(i/2)*sin(i*Pi/4)]; c[i+1]:=circle(b,0.05); p[i+1]:=plots[display]([plot([[b[1],b[2]-k],[b[1],b[2]],[b[1]-k,b[2]]], color=white), p[i],plot([[b[1],b[2]-k],[b[1]-k,b[2]]]),translate(rotate(translate(p[i],-b[1],-b[2]),-Pi/2),b[1],b[2])]) od; RETURN(plots[display]([p[m], seq(c[i], i=1..m)],axes=framed,scaling=constrained)) end:

> dragon1(10, 1/4);

> plots[display]([p[7], rotate(p[7],Pi/2), rotate(p[7], -Pi/2), rotate(p[7],Pi)],scaling=constrained,axes=framed);

>

> dragon2 := proc(k::algebraic, N::integer) local t,i,q1,q2,q3,q4,delta; global p; q2:=[k,0]; q3:=[1-k,0]; delta:=evalm(q3-q2); p[0]:=plot([[q2[1],q2[2]], [q3[1],q3[2]]]); for i from 1 to N do if (i mod 2=0) then t[i]:=t[i/2] else t[i]:=(i mod 4)-2 fi; q4:=evalm(q3+k/(1-2*k)*delta); delta:=evalm(t[i]*[delta[2],-delta[1]]); q1:=evalm(q3); q2:=evalm(q4+k/(1-2*k)*delta); q3:=evalm(q2+delta); p[i]:=plot([[q1[1],q1[2]], [q2[1],q2[2]], [q3[1],q3[2]]]) od; RETURN(plots[display]([seq(p[i], i=0..N)])) end:

> dragon2(0.2, 128);

>

13.5 Menger Curve

Plot the figure based on the cube of zero rank hexahedron([x0,y0,z0], a0) .

> restart: with(plots): with(plottools): x0:=1: y0:=2: z0:=3: a0:=.8:

Warning, existing definition for changecoords has been overwritten

> display(cutout(hexahedron([x0, y0, z0], a0), 1/3));

Elements of and can be obtained as follows:

> N2:={$3..421} minus {seq(21*i+2, i=1..20)}: nops(N2); N3:={$3..8842} minus ({seq(421*j+2, j=0..20)} union {seq(seq(421*j+21*i+3, i=1..20), j=0..20)}): nops(N3);

>

> menger1 := proc(L::algebraic, lev::integer, x0,y0,z0) local i,j,k,f; global p,s; options remember; s:=s+1;f:=0.4; p[s]:=cutout(hexahedron([x0,y0,z0],L),f); if s=2 then p[1]:=p[2] fi; if lev > 1 then for i from -1 to 1 do for j from -1 to 1 do for k from -1 to 1 do if (k=0 and abs(i)+abs(j)>1) or (k<>0 and abs(i)+abs(j)>0) then menger1(L/3,lev-1, x0+2*i*L/3, y0+2*j*L/3, z0+2*k*L/3) fi od od od fi; RETURN(plots[display]([seq(p[i], i=1..(20^lev-1)/19)],orientation=[43,60],scaling=constrained)) end:

> s:=0: menger1(100, 2, 0, 0, 0);

>

> menger2 := proc(L::algebraic, lev::integer, x0, y0, z0) local i,j,k; global p,s; options remember; if s=0 then p[0]:=hexahedron([x0,y0,z0],L,style=line) fi; for i from -1 to 1 do for j from -1 to 1 do for k from -1 to 1 do if abs(i)+abs(j)+abs(k)=1 then s:=s+1; p[s]:=hexahedron([x0+2*i*L/3, y0+2*j*L/3, z0+2*k*L/3], L/3) fi od od od; if lev > 1 then for i from -1 to 1 do for j from -1 to 1 do for k from -1 to 1 do if abs(i)+abs(j)+abs(k)>1 then menger2(L/3,lev-1,x0+2*i*L/3, y0+2*j*L/3, z0+2*k*L/3) fi od od od fi; RETURN(plots[display]([seq(p[i], i=0..6*(20^lev-1)/19)],scaling=constrained)) end:

> s:=0: menger2(100, 3, 0, 0, 0);

>