One of my favourite topics in math is fractals. For this post, I drew the Hilbert Curve fractal in octave.
The Hilbert Curve is a space-filling curve: although it is one-dimensional (a series of line segments) itself, it covers a two-dimensional area. I recently learned about this fractal from a post on Andrea Hawksley's blog.
To draw the Hilbert Curve, I wrote a script in octave to implement the Lindenmayer system algorithm from Wikipedia. The results are shown below, after 5 and 7 iterations, respectively. The script I wrote is included at the end of this post.
With its space-filling properties, the Hilbert Curve has some commonalities with tiling, another mathematical topic that really interests me. I really like the aesthetics of this fractal.
%% Hilbert.m - this is a script to draw the space-filling Hilbert curve. % From Wikipedia, the algorithm is as follows: % A = - B F + A F A + F B - % B = + A F - B F B - F A + % Where: F indicates move forward, - is a left turn (which actually should % be positive, using the right hand rule), + is a right turn, and A and B % are place-holders that are expanded as indicated above on each iteration. % (In the script below, numbers 1 - 5 are used for A, B, -, +, and F, % respectively) % This script has 3 parts: 1. generate the instruction list to the depth % (number of iterations) desired. 2. calculate the coordinates of each % vertex on the curve. 3. plot the curve. %% Generate the instruction list A = [ 3 2 5 4 1 5 1 4 5 2 3 ]; B = [ 4 1 5 3 2 5 2 3 5 1 4 ]; Depth = input("How many iterations? \n"); InstList = A; Iter = 1; while Iter < Depth InstListNew = ; %start with a blank list to fill in for Item = 1:length(InstList) switch (InstList(Item)) case 1 InstListNew = [ InstListNew, A ]; %expand the instruction list where appropriate case 2 InstListNew = [ InstListNew, B ]; otherwise InstListNew = [ InstListNew, InstList(Item) ]; endswitch end InstList = InstListNew; Iter = Iter + 1; end %% Calculate the vertex coordinates R3 = [ 0 -1; 1 0 ]; %rotation matrix for left turn ("3" in instruction list) R4 = [ 0 1; -1 0 ]; %rotation matrix for right turn ("4" in instruction list) dir = [ 1; 0 ]; %initial direction vector points along the x-axis posx = 0; %initial position is (0,0) posy = 0; xpts = [ posx ]; %initialize x and y vectors, plus distance tally and a step counter ypts = [ posy ]; distance = [ 0 ]; step = 1; for Index = 1:length(InstList) switch (InstList(Index)) case 3 dir = R3*dir; %rotate left case 4 dir = R4*dir; %rotate right case 5 posx = posx + dir(1); %move forward in the direction currently faced. posy = posy + dir(2); step = step + 1; xpts(step) = posx; %record current position ypts(step) = posy; distance(step) = (step - 1); %record cumulative distance travelled endswitch end %% Draw the curve plot(xpts,ypts,'b-'); %plot the recorded points xlim([-1, (max(xpts)+1)]); ylim([-1, (max(ypts)+1)]); %It might also be of interest to plot x and/or y vs the distance travelled % (commented out for now) %plot(distance,xpts,'g-'); %plot(distance,ypts,'r-'); %Note: save the figure using the command 'print('filename.png','-dpng')'