# File: fraktale.txt # All procedures were written for Maple V Release 3. # # Original file location: # http://www.math.utsa.edu/mirrors/maple/maplev/fractals/fraktale.ms # # ------------------------------------------------------------------------------------------------ # # Note: Make sure that Digits <= 10 in order to speed up computations. # # Credits: # John Oprea, oprea@math.csuohio.edu, implemetented the J procedure to pass initial # points and the anonymous functions to the fractal procedures. His procedures also # caused me to add further fractals. # # The escape-time fractal procedures are based on a Maple V algorithm taken from the book # 'Maple V - Programming Guide' by M.B. Monagan, K.O. Geddes, G. Labahn and S. Vorkoetter, # Springer Verlag. # # Remark on the first procedure 'mandelbrot' by John Oprea: # The modification on mandelbrot uses the fact that Maple can color things according to a function. # Also, plotting the zero function keeps the graph in a plane. Increasing the values of 'grid' in the # plot3d command results in better resoultion of the fractal. However, this means longer # computation time. 'style=patchnogrid' draws the surface with shaded rectangular patches with # no wireframe grid. This is the best 'style' option for this routine. The procedure mandelbrot returns # the value m. The higher m the more likely the probability that the point (x, y) is part of the # Mandelbrot set. Certainly m depends of the number of iterations (here: 30). It is possible that a # point satisfies the condition abs(z)<2 after n iterations but is outside the range after n+1 iterations. # So to make sure that (x, y) is part of the set you would have to do infinite iterations. m is used in the # 'color=mandelbrot' option to visualize the set: the red area - the Mandelbrot lake - indicates that # it is part of the set. > restart: > mandelbrot:=proc(x, y) > local z, c, m; > c:=evalf(x+y*I); > z:=evalf(x+y*I); > for m from 0 to 30 while abs(z)<2 do > z:=z^2+c > od; > m > end: > plot3d(0, -2 .. 0.7, -1.2 .. 1.2, orientation=[-90,0], grid=[250, 250], style=patchnogrid, > scaling=constrained, color=mandelbrot); > plot3d(0, -1.3 .. -1.25, 0.35 .. 0.4, orientation=[-90,0], grid=[100, 100], style=patchnogrid, scaling=con > strained, color=mandelbrot); > restart: > mandelbrot:=proc(x, y) > local z, m; > z:=evalf(x+y*I); > for m from 0 to 25 while abs(z)<2 do > z:=z^2+(x+y*I) > od; > m > end: > plot3d(mandelbrot, -2 .. 0.5, -1.2 .. 1.2, grid=[50, 50]); > restart: # Julia sets modified by John Oprea, oprea@math.csuohio.edu > julia := proc(c, x, y) > local z, m; > z := evalf(x+y*I); > for m from 0 to 30 while abs(z) < 3 do > z := z^2+c > od; > m > end: > J := proc(d) > global phonyvar; > phonyvar := d; > (x, y) -> julia(phonyvar, x, y) > end: > plot3d(0, -2 .. 2, -1.3 ..1.3, style=patchnogrid, orientation=[-90,0], grid=[150, 150], > scaling=constrained, color=J(-1.25)); > plot3d(0, -1.9 .. 1.9, -1.9 ..1.9, style=patchnogrid, orientation=[-90,0], grid=[100, 100], > scaling=constrained, color=J(-.74543+.11301*I)); > plot3d(0, -1.9 .. 1.9, -1.2 ..1.2, style=patchnogrid, orientation=[-90,0], grid=[150, 150], > scaling=constrained, color=J(-.9+.12*I)); > restart: > # LAMBDAFN > lambdafn_sin:=proc(x, y) > local z, m, p1; > z:=evalf(x+y*I); > p1:=evalf(1+I*0.4); > for m from 0 to 100 while abs(z)<4 do > z:=sin(z)*p1 > od; > m > end: > plot3d(0, -4 .. 4, -3 .. 3, orientation=[-90,0], grid=[50, 50], style=patchnogrid, > scaling=constrained, color=lambdafn_sin); > plot3d(0, -1.5 .. 1.5, -1.1 .. 1.1, orientation=[-90,0], grid=[250, 250], style=patchnogrid, > scaling=constrained, color=lambdafn_sin); > restart: > # BIOMORPH2 > # May 22, 1996 > biomorph2:=proc(x, y) > local z, m, c; > z:=evalf(x+y*I); > c:=z; > for m from 0 to 100 while abs(z)<4 do > z:=cos(z)*z^2+c > od; > m > end: > plot3d(0, -2 .. 2, -1.5 .. 1.5, orientation=[-90,0], grid=[20, 20], style=patchnogrid, > scaling=constrained, color=biomorph2); > restart: > # NEWTON > # Visualization of Newton's method to calculate the complex solution of z^n-1=0 by > # the iteration z(j+1) = z(j) - fn(z)/diff(fn(z)) with fn(z) := z^n-1, and n > 2 the polynomial degree > # of fn (also real values for n are allowed). > solve(z^3-1=0, z); > # The graph shows the starting points from which the method does or does not > # converge to one of the above roots. The documentation of FRACTINT explains the result: > # "The coloring of the plot shows the "basins of attraction" for each root of the polynomial > # -- i.e., an initial guess within any area of a given color would lead you to one of the roots. > # As you can see, things get [fractal] along certain radial lines or "spokes," those being the lines > # between actual roots." > # > # The fractal properties of Netwon's method to calculate complex roots were first > # discovered by John Hubbard of Orsay, France. > newton:=proc(x, y) > local z, m; > z:=evalf(x+y*I); > for m from 0 to 50 while abs(z^3-1) >= 0.001 do > z:=z-(z^3-1)/(3*z^2) > od; > m > end: > plot3d(0, -2 .. 2, -1.5 .. 1.5, orientation=[-90,0], grid=[250, 250], style=patchnogrid, > scaling=constrained, color=newton); > restart: > cmplxmarkmand:=proc(x, y) > # May 25, 1996 > # The original formula 'z:=z^2*c^(p-1)+c' was developed by Mark Peterson. > local c, m, z; > c:=evalf(x+y*I); > z:=evalf(x+y*I); > for m from 0 to 30 while abs(z)<2 do > z:=z^2*c^0.1+c > od; > ln(m) > end: > plot3d(0, -2 .. 0.7, -1.2 .. 1.2, orientation=[-90,0], grid=[250, 250], style=patchnogrid, > scaling=constrained, color=cmplxmarkmand); > restart: > lambda:=proc(x, y) > # Lambda > # May 26, 1996 > # > # Complex version of the population growth equation x(j+1)=r*x(j)*(1-x(j)) > # developed by P. F. Verhulst in 1845 > # Lambda is part of the Julia family > local c, z, m; > c:=0.85+I*0.6; > z:=evalf(x+y*I); > for m from 0 to 30 while abs(z)<3 do > z:=c*z*(1-z) > od; > m > end: > plot3d(0, -1.5 .. 2.5, -1.5 ..1.5, style=patchnogrid, orientation=[-90,0], grid=[25, 25], > scaling=constrained, color=lambda); > restart: > ifs:=proc(imax) > # Iterated Functions Systems > # original program 'Chaos-Spiel fuer ein Farnblatt' written by Heinz-Otto Peitgen, > # Hartmut Juergens and Dietmar Saupe in BASIC, > # published in: 'Bausteine des Chaos - Fraktale', Springer Verlag/Klett-Cotta, > # p. 415, written by the authors mentioned above > # June 06, 1997 > local e1, e2, e3, e4, f1, f2, f3, f4, x, y, xn, yn, z, pts; > pts:=NULL; > e1:=0.5; e2:=0.57; e3:=0.408; e4:=0.1075; > f1:=0; f2:=-0.036; f3:=0.0893; f4:=0.27; > x:=e1; > y:=0; > to imax do > z:=rand()/1e12; > if z<=0.02 then > xn:=e1; > yn:=0.27*y+f1 > elif z<=0.17 then > xn:=-0.139*x+0.263*y+e2; > yn:=0.246*x+0.224*y+f2 > elif z<=0.3 then > xn:=0.17*x-0.215*y+e3; > yn:=0.222*x+0.176*y+f3 > else > xn:=0.781*x+0.034*y+e4; > yn:=-0.032*x+0.739*y+f4 > fi; > pts:=pts, [xn, yn]; > x:=xn; > y:=yn > od; > [pts] > end: > plot(ifs(12500), axes=NONE, style=point, symbol=POINT, scaling=constrained); > restart: > julfn_zsqrd:=proc(x, y) > # Based on a Maple V algorithm taken from the book > # 'Maple V - Programming Guide' by M.B. Monagan, > # K.O. Geddes, G. Labahn and S. Vorkoetter, > # Springer Verlag, modified by John Oprea > # modified by Alexander F. Walz > # May 25, 1996 > local c, z, m; > c:=evalf(-0.5+0.5*I); > z:=evalf(x+y*I); > for m from 0 to 30 while abs(z)<3 do > z:=sin(z)+z^2+c > od; > m > end: > plot3d(0, -1.91 .. 1.37, -1.24 ..1.21, style=patchnogrid, orientation=[-90,0], grid=[250, 250], scaling=co > nstrained, color=julfn_zsqrd); > restart: > mandelbrot_fast:=proc(x, y) > # Based on a Maple V algorithm taken from the book > # 'Maple V - Programming Guide' by M.B. Monagan, > # K.O. Geddes, G. Labahn and S. Vorkoetter, > # Springer Verlag, modified by John Oprea > # modified by Alexander F. Walz > # May 25, 1996 > local xn, xnold, yn, m; > xn:=x; > yn:=y; > for m from 0 to 100 while sqrt(evalhf(xn^2+yn^2)) < 2 do > xnold:=xn; > xn:=evalhf(xn^2-yn^2+x); > yn:=evalhf(2*xnold*yn+y) > od; > m > end: > plot3d(0, -0.845 .. -0.794, 0.172 .. 0.21, orientation=[-90,0], grid=[250, 250], style=patchnogrid, > scaling=constrained, color=mandelbrot_fast); > restart: > chaos:=proc(x1, y1, x2, y2, x3, y3, maxiter, seed1, seed2) > # original procedure written by Tom Williams and modified by John Oprea > local randi, x, y, sx, sy, ir, i, j; > randi:=rand(1 .. 3): > x:=[x1, x2, x3]: > y:=[y1, y2, y3]: > sx:=seed1: > sy:=seed2: > ir:=[sx, sy]: > for i from 1 to maxiter do > j:=randi(): > sx:=(sx+x[j])/2.: > sy:=(sy+y[j])/2.: > ir:=[op(ir), sx, sy] > od: > plot(ir, style=point, symbol=POINT, scaling=constrained) > end: > chaos(0, 0, .5, 1, 1, 0, 9000, .5, .5); > chaos(0, 0, -4, 7, 6, 2, 3000, 1, .5); > restart: > alexerror:=proc(x, iter) > # written by Alexander F. Walz > # May 29, 1996 > # March 02, 1997 > local n, pts, xn, xnold, yn, y; > pts:=[0, 0]; # in case pts is not assigned any value after completion of the for n loop > for y from 0 to 1.1 by 0.01 do > xn := x; > yn := y; > for n to iter while evalhf(xn^2+yn^2) < 4 do > xnold:=xn; > xn := evalhf(xn^2-yn^2+x); > yn := evalhf(2*xnold*yn+y) > od; > if n >= iter then pts := pts, [xn, yn], [xn, -yn] fi > # as a variant try: "if n >= iter then pts := pts, [x, y], [-x, -y] fi;" to compute the Mandelbrot set > od; > pts > end: > i:='i': PLOT(seq(POINTS(alexerror(i/100, 500)), i=-200 .. 70), SYMBOL(POINT), axes=none); > restart: > orbits:=proc(x, y, iter) > # procedure to plot the trajectory a point traverses by being iterated > # through the formula 'z(j+1)=z(j)^2+c' (also called 'orbits') for the > # Mandelbrot set > # written by Alexander F. Walz > # May 30, 1996 > local a, b, xn, xnold, yn, pts; > pts:=NULL; > a:=x; > b:=y; > xn:=x; > yn:=y; > to iter do > xnold:=xn; > xn:=evalf(xn^2-yn^2+a); > yn:=evalf(2*xnold*yn+b); > pts:=pts, [xn, yn] > od; > [pts] > end: > # orbits1.gif 'Spiral' > plot(orbits(-0.61143695190549, -0.36571035906672, 500), style=point, symbol=POINT, > axes=framed, scaling=constrained); > # orbits2.gif > plot(orbits(-0.68181818351150, -0.31095156446099, 300), style=point, symbol=POINT, > axes=framed, scaling=constrained); > # orbits3.gif 'Vortex' > plot(orbits(-0.34555229917169, -0.60430224984884, 300), style=point, symbol=POINT, > axes=framed, scaling=constrained); > # orbits4.gif > plot(orbits(+0.24877810105681, +0.58083451911807, 115), style=line, > axes=framed, scaling=constrained); > # orbits5.gif 'Radiation' > plot(orbits(-0.33382209390402, -0.60821359232068, 500), style=point, symbol=POINT, > axes=framed, scaling=constrained); > # orbits6.gif > plot(orbits(-0.25953079387546, -0.62385896220803, 500), style=point, symbol=POINT, > axes=framed, scaling=constrained); > # orbits7a.gif 'Curved' > plot(orbits(-0.14222874119878, -0.64732701703906, 500), style=point, symbol=POINT, > axes=framed, scaling=constrained); > # orbits7b.gif 'Coil' > plot(orbits(-0.14222874119878, -0.64732701703906, 300), style=line, > axes=framed, scaling=constrained); > # orbits8a.gif 'Spiral Galaxies' > plot(orbits(-0.50977517291904, -0.60039090737700, 500), style=point, symbol=POINT, > axes=framed, scaling=constrained); > # orbits8b.gif 'US Star' > plot(orbits(-0.50977517291904, -0.60039090737700, 300), style=line, > axes=framed, scaling=constrained); > restart: > # LAMBDAFN_FAST > # Using hardware coprocessor > # Based on a Maple V algorithm taken from the book > # 'Maple V - Programming Guide' by M.B. Monagan, > # K.O. Geddes, G. Labahn and S. Vorkoetter, > # Springer Verlag, modified by John Oprea > # modified by Alexander F. Walz > # May 31, 1996 > # > lambdafn_sinfast:=proc(x, y) > local xn, xnold, yn, m; > xn:=x; > yn:=y; > for m from 0 to 100 while evalhf(sqrt(xn^2+yn^2)) < 4 do > xnold:=xn; > xn:=evalhf(sin(xn)*cosh(yn) - cos(xn)*sinh(yn)*0.4); > yn:=evalhf(cos(xnold)*sinh(yn) + sin(xnold)*cosh(yn)*0.4) > od; > m > end: > plot3d(0, -4 .. 4, -3 .. 3, orientation=[-90,0], grid=[250, 250], style=patchnogrid, scaling=constrained, > color=lambdafn_sinfast); > plot3d(0, -1.5 .. 1.5, -1.1 .. 1.1, orientation=[-90,0], grid=[25, 25], style=patchnogrid, scaling=constrai > ned, color=lambdafn_sinfast); > restart: > kamtorus:=proc(angle, step, ende, points) > # written by Alexander F. Walz > # May 31, 1996 > local x, xold, y, orbit, pts; > pts:=NULL; > for orbit from 0 by step to ende do > x:=orbit/3; > y:=x; > to points > while abs(x) < 10^38 or abs(y) < 10^38 do > xold:=x; > x := evalhf(x*cos(angle) + (x^2-y)*sin(angle)); > y := evalhf(xold*sin(angle) - (xold^2-y)*cos(angle)); > pts:=pts, [x, y] > od > od; > [pts] > end: > plot(kamtorus(1.3, 0.05, 1.5, 200), style=point, symbol=POINT, axes=none, scaling=constrained); > restart: > koch := proc(p: numeric) > # SNOWFLAKE > # developed by Swedish mathematician Helge von Koch in 1904 > # algorithm taken from 'Algorithmen fuer Chaos und Fraktale' > # by Dietmar Herrmann, Addison-Wesley > # original program written in BASIC > # p is the number of recursions, should be no more than value 5 > local m, n, k, l, s, h, x, y, pts, t; > h := 3^(-p); > pts := NULL; > x := 0; y := 0; > for n from 0 to (4^p-1) do > m := n; > s := 0; > for l from 0 to p-1 do > t := irem(m, 4); > m := iquo(m, 4); > s := s+irem((t+1), 3) -1 > od; > x := x+cos(Pi*s/3)*h; > y := y+sin(Pi*s/3)*h; > pts := pts, [x, y] > od; > [pts] > end: > plot(koch(5), scaling=constrained, style=LINE, axes=NONE);