Engineering Math - Chaos Theory

 

 

 

 

Barnsley Fern

 

Barnsley Fern is a fractal set that looks like a shape of a fern. If you plot this and magnify any one part (e.g, the tip of a leaf) you will see the similar shape of the whole leafs wherever you magnify.

 

 

The procedure to plot this as follows.

    i) define four sets of 2D transformation equation as follows. (p stands for the coordinate for the current point, m is a 2x2 matrix that is specially created for this transformation, v is 2x1 vector)

      tf1 = m1 . p + v1

      tf2 = m2 . p + v2

      tf3 = m2 . p + v3

      tf4 = m2 . p + v4

    ii) pick any arbitray points (let's say this point is named as p)

    iii) pick one transformation equation out of the 4 equation defined in step i). The probability of the random selection for each transformation equation should be predefined).

    iv) transform p with the transformation equation selected at step iii).

    v) move the current point to the transformed location and set the new position to the current point (p).

    vi) repeat the process iii) ~ v) and plot the point p on the coordinate at each iterations.

Following is the Octave/Matlab code that I wrote as per the procedure described above.

    m1 = [0.85 0.04;-0.04 0.85]; v1 = [0.0;1.60];

    m2 = [0.2 -0.26;0.23 0.22]; v2 = [0.0;1.60];

    m3 = [-0.15 0.28;0.26 0.24]; v3 = [0.0;0.44];

    m4 = [0.0 0.0;0.0 0.16]; v4 = [0.0;0.0];

     

    prob = [0.85 0.07 0.07 0.01];

     

    Np = 50000;

     

    p_old = [0.5;0.5];

     

    pList = [p_old'];

     

    for i = 1:Np

        r = rand();

        if (r < prob(1))

           p_new = m1 * p_old + v1;

        elseif (r < (prob(1)+prob(2)))

           p_new = m2 * p_old + v2;

        elseif (r < (prob(1)+prob(2)+prob(3)))

           p_new = m3 * p_old + v3;

        else

           p_new = m4 * p_old + v4;

        endif;

        pList = [pList;p_new'];

        p_old = p_new;

    end;

Note : the syntax of if statement of Octave may be a little different from Matlab. I tried this only on Octave, so not sure if this code would work in matlab.