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. 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. fun man c z = add (square z) c; We can now investigate a point - say (0.5, 0.5) call that point tim: 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. - (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 - (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. 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. 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. 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. - 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: 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 (* 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. *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 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.