Diversion : Mandelbrot sets
The Mandelbrot set is a set of complex numbers. To find if a point is
in the Mandelbrot set we work out the Mandelbrot function for that point
then repeatedly apply function over and over starting at zero. If the
resulting sequence does not "explode" then the original point is in the
Mandelbrot set.
We can represent a complex number as a pair of reals and define the
functions square and add.
<![CDATA[fun add (r1,i1) (r2,i2) = (r1+r2,i1+i2):real*real;
fun square(r,i) = (r*r-i*i, 2.0*r*i);
]]>The Mandelbrot function is simply z ²+c.
<![CDATA[fun man c z = add (square z) c;
]]>We can now investigate a point - say (0.5, 0.5) call that point tim:
<![CDATA[val tim = (0.5,0.5);
]]>The Mandelbrot function for that point is (man tim), we apply this
function to zero and notice that the first application brings us to tim again.
<![CDATA[- (man tim) (0.0,0.0);
val it = (0.5, 0.5) : real * real
]]>The label it hold the most recently computed value, we can feed that
back into man tim
<![CDATA[- (man tim) it;
val it = (0.5,1.0): real * real
- man tim it;
val it = (~0.25,1.5):real*real
- man tim it;...
]]>If you continue you will find the values get bigger and bigger until
an overflow occurs. This means that tim is not in the Mandelbrot set. We
shall try point (0.3, 0.1) called jane.
<![CDATA[val jane = (0.3,0.1);
man jane (0.0,0.0);
man jane it;
man jane it;
...
]]>You should find that jane's sequence settles down to around
(0.32..., 0.28...) after a few goes.
The point (~0.5,0.5) is bob. You will notice that bob takes a long
time to settle down.
<![CDATA[val bob = (~0.5,0.5);
man bob (0.0,0.0);
man bob it;
]]>man bob it; ...
Copy and paste - if you have it - will speed things up, but not enough...
Rather than apply (man bob) just once each time we can use the twice
function to speed things up. We also define zero to save typing.
<![CDATA[val zero = (0.0,0.0);
fun twice f = f o f;
fun t256 f = twice(twice(twice twice)) f; (* apply 256 times *)
]]>We can now launch the point bob through 256 iterations then look
at it to see if it has settled down into a pattern.
<![CDATA[- t256 (man bob) zero;
val it = (~0.40855...,0.268...)
- man bob it;
val it = (~0.40511..,0.2806...)
...
]]>As you can, the point is still bobbing about in the second/third
decimal place.
Further experiments
Observe the behaviour of the following points, some of them settle
down almost immediately some must be put through a few hundred iterations
first, some never seem to exhibit any kind of predictable behaviour. In
each case you are advised to go through a few hundred iterations then
try a few function applications one at a time.
Points:
<![CDATA[val fast = (0.2,0.2);
val vacillate = (~1.0,0.0);
val bound = (0.15625,0.5625);
val tri = (~0.1,0.7);
]]>All very well Andrew, but I can't put this stuff on a T-shirt
True. To get pictures of the fabulous Mandelbrot set try copying
in the following
<![CDATA[
(* Mandelbrot set investigator New Jersey ML *)
(* A Cumming Napier University *)
(* Start off with some useful complex number functions *)
fun square(x,y) = (x*x-y*y,2.0*x*y);
fun add (x,y) (u,v) = (x+u,y+v):real*real;
fun scalar (s:real) (x,y) = (s*x,s*y);
fun dot (x1,y1) (x2,y2) = (x1*x2,y1*y2):real*real;
fun sub p q = add p (scalar ~1.0 q);
fun dist p q = let
fun r(x,y)=sqrt(x*x+y*y)
in r(sub p q) end;
val zero = (0.0,0.0);
(* Mandelbrot functions *)
fun man c z = add (square z) c;
(* Multiple applications of a function *)
fun twice f = f o f;
fun thrice f = f o f o f;
fun quice f = twice twice f;
fun t256 f = twice(twice(twice twice)) f; (* apply 256 times *)
(* Catagorise the point type 1,2,3,4,* or " " *)
(* Complete 200 iterations then test for period of 1,2,3 or 4 *)
fun cat p = let
val small = 0.001;
val a= t256 (man p) zero;
in if dist a (man p a)<0.001 then "1"
else if dist a (twice (man p) a)<small then "2"
else if dist a (thrice (man p) a)<small then "3"
else if dist a (quice (man p) a)<small then "4"
else "*"
end handle Overflow => " ";
fun for (r1:real) r2 d f = if (0.0 > (r2-r1)*d) then ""
else (f r1)^(for (r1+d) r2 d f);
fun line x1 x2 y = (for x1 x2 ((x2-x1)/78.0) (fn i=> cat(i,y)))^"\n";
fun box((x1,y1),(x2,y2))= for y2 y1 ((y1-y2)/24.0) (fn i=> line x1 x2 i);
fun K a b = b;
val ibox=( (~2.0,~1.0),(1.0,1.0));
(* Box shifting stuff *)
fun zoom(p, q) = (add (scalar 0.75 p) (scalar 0.25 q),
add (scalar 0.25 p) (scalar 0.75 q));
fun shift v (p,q) = let val w = dot v (sub q p) in
(add w p, add w q) end;
val right = shift (0.5,0.0);
val left = shift (~0.5,0.0);
val up = shift (0.0,0.5);
val down = shift (0.0,~0.5);
fun doit x = K (print(box x)) x;
doit ibox;
]]>
The function doit takes the bottom left and top right corners of a box and
prints the set to the screen - the box itself is the output, you can shift
the box left, right, up and down or zoom in.
<![CDATA[
*3
*33333
33*
* *211111111111*
*111111111111111111444
**1111111111111111111111*
*1111111111111111111111111**
** ** * *111111111111111111111111111
222222222* *1111111111111111111111111111
*22222222222**111111111111111111111111111
* **444222222222222221111111111111111111111111
*22222222222**111111111111111111111111111
222222222* *1111111111111111111111111111
** ** * *111111111111111111111111111
*1111111111111111111111111**
**1111111111111111111111*
*111111111111111111444
* *211111111111*
33*
*33333
*3
]]>
Try typing in the following
<![CDATA[
doit ibox;
doit (zoom ibox);
doit (zoom it);
doit (up(left it));
]]>If you have a very fast machine or plenty of time you can produce
colour BMP pictures of the Mandelbrot set.
Prize puzzle - find a point which has order 131. Submit the point
together with ML code demonstrating it's order by mail to andrew. The
prize is either a large cash sum (£5) or hearty congratulations.