/* Unit ( package ) for Maxima http://maxima.sourceforge.net for discrete dynamical system based on complex quadratic polynomial fc(z):=z*z+c version : 2008.12.21 licence : GNU Public License usage : ( save as qdynamics.mac in directory (for example) C:\Program Files\Maxima-5.16.3\share\maxima\5.16.3\share\dynamics ) load("qdynamics"); examples of use are in the description of procedures ------------------------------ tested in : wxMaxima 0.7.6 http://wxmaxima.sourceforge.net Maxima 5.16.3 http://maxima.sourceforge.net Using Lisp GNU Common Lisp (GCL) GCL 2.6.8 (aka GCL) Distributed under the GNU Public License. Dedicated to the memory of William Schelter. ---------------------------------------------- bugs : Maxima 5.11.0 http://maxima.sourceforge.net Using Lisp SBCL 1.0.2 Distributed under the GNU Public License. See the file COPYING. Dedicated to the memory of William Schelter. This is a development version of Maxima. The function bug_report() provides bug reporting information. (%i1) load("qdynamics20081103.mac"); ... (%o1) qdynamics20081103.mac (%i2) centers:GiveCenters(6); Maxima encountered a Lisp error: Error during processing of --eval option "(cl-user::run)": invalid number of arguments: 2 Automatically continuing. To reenable the Lisp debugger set *debugger-hook* to nil. -------------------------------------------------------- */ disp("qdynamics.mac 2008.12.21 ; unit for discrete dynamical system based on quadratic polynomial")$ disp(" for drawing it uses draw package = interface to gnuplot by Mario Rodriguez Riotorto http://www.telefonica.net/web2/biomates"); disp("uses jtroot3 Maxima package (Jenkins-Traub polynomial root finder by Raymond Toy for GiveCentersP"); disp("uses cpoly.lisp needs Revision 1.16 - Wed Dec 17 22:02:30 2008 UTC from cvs http://maxima.cvs.sourceforge.net/viewvc/maxima/maxima/src/cpoly.lisp"); /* ======================= definitions ============================================================== */ /* basic funtion = monic and centered complex quadratic polynomial http://en.wikipedia.org/wiki/Complex_quadratic_polynomial */ f(z,c):=z*z+c $ /* iterated function */ F(i, z, c) := if i=1 then f(z,c) else f(F(i-1, z, c),c) $ /* */ G(i,z,c):=F(i, z, c)-z $ /* inverse functions to f(z,c) */ f_inverse_plus(z,c):=sqrt(z-c)$ f_inverse_minus(z,c):=-sqrt(z-c)$ /* ======================== functions ===================================================================== */ /* gives angle of complex number in turns modulo 1 360 deg = 1 turn = 2*%pi radians */ cturn(z):=mod(float(carg(z)/(2*%pi)),1)$ /* gives distancve between 2 complex points */ GiveDistance(c1,c2):=abs(abs(c1)-abs(c2))$ /* ---------------------- IsEqual ------------------ function that compares 2 complex numbers if they are similar with precision=eps then returns true else returns false */ IsEqual(c1,c2,eps):= block( if GiveDistance(c1,c2)<=eps then return( true) else return (false) )$ /* ------------- log2 ------------- */ log2(x):=log(x)/log(2)$ /* ------------ critical_orbit --------------------------------- this function computes forward orbit of critical point z=0 of complex quadratic polynomial f(z,c) for given parameter c gives a list of iMax points z as a result example of use : orbit:critical_orbit(-1,100); period:GivePeriodOfOrbit(orbit,0.001); DrawOrbit(orbit); DrawOrbitDistance(orbit); */ critical_orbit(c,iMax):= block( /* iMax : Number of iterations , for most cases iMax:100 */ ER:1000000, z:0+0*%i, /* first point = critical point */ orbit:[z], /* compute forward orbit */ i:0, if (abs(z)>ER) then return(orbit), loop, z:rectform(f(z,c)), orbit:endcons(z,orbit), i:i+1, if ((abs(z)1 then block ( loop, i:i-1, if i>0 and not IsEqual(orbit[iLast],orbit[i],eps) then go(loop), /* ? */ period:iLast-i ), return(period) )$ /* shorter version of GivePeriodOfC usefull for testing centers : centers:GiveCenters(4); a:map(gpc,centers); */ gpc(c):=GivePeriodOfC(c,0.01)$ /* ------------------- GivePeriodOfOrbit --------------------------------------- Gives period of periodic orbit {z1,..,zn} for given parameter c example of use : orbit:critical_orbit(-1,100); period:GivePeriodOfOrbit(orbit,0.01); testing centers : gpc(c):=GivePeriodOfC(c,0.01); a:map(gpc,centers); */ GivePeriodOfOrbit(orbit,eps):= block ( period:0, iLast:length(orbit), i:iLast, if iLast>1 then block ( loop, i:i-1, if i>1 and not IsEqual(orbit[iLast],orbit[i],eps) then go(loop), /* ? */ period:iLast-i ), return(period) )$ /* ---------- GiveFixedPoint --------------------------- For f(z,c) there are 2 fixed points z1 and z2, typically one is attracting=alfa, one repelling=beta see : http://en.wikipedia.org/wiki/Periodic_points_of_complex_quadratic_mappings */ GiveFixedPoint(c):= float(rectform(solve([z*z+c=z],[z])))$ /*----------------Give repellor -------------------------*/ GiveRepellingFixedPoint(c):= block( [fixed:GiveFixedPoint(c)], /* Find which is repelling */ if (abs(2*rhs(fixed[1]))<1) then return(rhs(fixed[2])) else return(rhs(fixed[1])) )$ /* ------------- GiveAttractor -----------------------*/ GiveAttractingFixedPoint(c):= block( fixed:GiveFixedPoint(c), /* Find which is repelling */ if (abs(2*rhs(fixed[1]))<1) then return(rhs(fixed[2])) else return(rhs(fixed[1])) )$ /* ------------- GiveCenters ----------------------- Computes list of complex points c that are centers of hyperbolic components of period p of Mandelbrot set Uses allroots, which works up to period 6 ( for 7 not all centers have period 7 ) This function uses an algorithm by Robert Munafo example of use : centers:GiveCenters(4); testing centers : a:map(gpc,centers); */ GiveCenters(p):= block( /* if cc could not be computed then gives empty list */ [cc:[],f:ifactors(p)], /* ------------ prime numbers ------------------------------------- */ if (p=1) then g:G(p,0,c) elseif primep(p) then g:G(p,0,c)/G(1,0,c) /* prime > 1 */ /* ------------ composite numbers ------------------------------------- */ /* PowerOfPrime, for example 9=3*3 */ elseif length(f)=1 then g:G(p,0,c)/G(p/2,0,c) /* Product Of 2 different Primes */ elseif length(f)=2 and f[1][2]+f[2][2]=2 then g:G(p,0,c)*G(1,0,c)/(G(f[1][1],0,c)*G(f[2][1],0,c)) /* p=q*r*r, a prime times the square of another prime */ elseif length(f)=2 and f[1][2]+f[2][2]=3 /* (Zn(p)*Zn(r))/(Zn(q*r)*Zn(r*r)) */ then g:G(p,0,c)*G(r(p),0,c)/(G(p/r(p),0,c)*G(r(p)*r(p),0,c)), cc:allroots(expand(%i*g)=0,c), cc:map(rhs,cc),/* remove string "c=" */ return(cc) )$ /* works up to period 8 fpprec: 16 is good up to period 6, 32 for period 7, 64 for period 8 load(cpoly); needs Revision 1.16 - Wed Dec 17 22:02:30 2008 UTC from cvs http://maxima.cvs.sourceforge.net/viewvc/maxima/maxima/src/cpoly.lisp. example of use : load(qdynamics;); load(cpoly); fpprec:16; /* standard value */ centers:GiveCenters_bf(3); DrawPointsT(centers,"centers of components of period 3",screen); length(centers); */ GiveCenters_bf(p):= block( /* if cc could not be computed then gives empty list */ [cc:[],f:ifactors(p)], /* ------------ prime numbers ------------------------------------- */ if (p=1) then g:G(p,0,c) elseif primep(p) then g:G(p,0,c)/G(1,0,c) /* prime > 1 */ /* ------------ composite numbers ------------------------------------- */ /* PowerOfPrime, for example 9=3*3 */ elseif length(f)=1 then g:G(p,0,c)/G(p/2,0,c) /* Product Of 2 different Primes */ elseif length(f)=2 and f[1][2]+f[2][2]=2 then g:G(p,0,c)*G(1,0,c)/(G(f[1][1],0,c)*G(f[2][1],0,c)) /* p=q*r*r, a prime times the square of another prime */ elseif length(f)=2 and f[1][2]+f[2][2]=3 /* (Zn(p)*Zn(r))/(Zn(q*r)*Zn(r*r)) */ then g:G(p,0,c)*G(r(p),0,c)/(G(p/r(p),0,c)*G(r(p)*r(p),0,c)), cc:bfallroots(expand(%i*g)=0), cc:map(rhs,cc),/* remove string "c=" */ return(cc) )$ /* ------------- GiveCentersP ----------------------- Computes list of complex points c that are centers of hyperbolic components of period p of Mandelbrot set Uses polyroots, which works up to period 7 ( for fpprec : 32 in 340 sec ), time for period 6 is 35 sec for period 8 only some roots are period 8 centers This function uses an algorithm by Robert Munafo example of use : fpprec : 32; (Default value: 16, which not work for period 7 ) centers:GiveCentersP(4); testing centers : a:map(gpc,centers); DrawPoints(centers,"centers of components of period 4"); */ GiveCentersP(p):= block( load("jtroot3.mac"), /* jtroot3 Maxima package (Jenkins-Traub polynomial root finder) by Raymond Toy */ maperror:false, /* if cc could not be computed then gives empty list */ [cc:[],f:ifactors(p)], /* ------------ prime numbers ------------------------------------- */ if (p=1) then g:G(p,0,c) elseif primep(p) then g:G(p,0,c)/G(1,0,c) /* prime > 1 */ /* ------------ composite numbers ------------------------------------- */ /* PowerOfPrime, for example 9=3^2 or 4=2^2 */ elseif length(f)=1 then g:radcan(G(p,0,c)/G(p/2,0,c)) /* Product Of 2 different Primes */ elseif length(f)=2 and f[1][2]+f[2][2]=2 then g:radcan(G(p,0,c)*G(1,0,c)/(G(f[1][1],0,c)*G(f[2][1],0,c))) /* p=q*r*r, a prime times the square of another prime (Zn(p)*Zn(r))/(Zn(q*r)*Zn(r*r)) */ elseif length(f)=2 and f[1][2]+f[2][2]=3 then g:G(p,0,c)*G(r(p),0,c)/(G(p/r(p),0,c)*G(r(p)*r(p),0,c)), cc:polyroots(g,c), return(cc) )$ /* ================================================================================================== ================================ draw functions ================================================== ================================================================================================== */ /* draws a list of complex (2D) points example of use : centers:GiveCenters(4); DrawPoints(centers,"centers of components of period 4"); */ DrawPoints(list,sTitle):= block( /* save the orbit to 2 lists */ xx: map(realpart, list), yy: map(imagpart, list), /* draw using gnuplot thru Maxima draw package */ load(draw), /* ( interface to gnuplot ) by Mario Rodriguez Riotorto http://www.telefonica.net/web2/biomates */ draw2d( title= sTitle, user_preamble = "", terminal = 'screen, pic_width = 1000, pic_height = 1000, xlabel = "re ", ylabel = "im", point_type = filled_circle, points_joined = false, color =red, point_size = 0.5, points(xx,yy) ) )$ /* example of use : centers:GiveCenters(6); DrawPointsT(centers,"centers of components of period 6",screen); DrawPointsT(centers,"centers of components of period 6",png); look for a file in directory ( example ) C:\Program Files\Maxima-5.16.3\wxMaxima */ DrawPointsT(list,sTitle,sTerminal):= block( /* save the orbit to 2 lists */ xx: map(realpart, list), yy: map(imagpart, list), /* draw using gnuplot thru Maxima draw package */ load(draw), /* ( interface to gnuplot ) by Mario Rodriguez Riotorto http://www.telefonica.net/web2/biomates */ draw2d( file_name = "p", pic_width=1000, pic_height= 1000, title= sTitle, user_preamble = "", terminal = concat(sTerminal), pic_width = 1000, pic_height = 1000, xlabel = "re ", ylabel = "im", point_type = filled_circle, points_joined = false, color =red, point_size = 0.5, points(xx,yy) ) )$ /* --------------------------- DrawOrbit ------------------------------------------- */ /* draws orbit on parameter plane on screen using draw package example of use : orbit:critical_orbit(-1,100); DrawOrbit(orbit); */ DrawOrbit(orbit):= block( /* save the orbit to 2 lists */ xx:makelist (realpart(orbit[1]), i, 1, 1), /* list of re(z) */ yy:makelist (imagpart(orbit[1]), i, 1, 1), /* list of im(z) */ for i:2 thru length(orbit) step 1 do block ( xx:cons(realpart(orbit[i]),xx), yy:cons(imagpart(orbit[i]),yy) ), load(draw), /* ( interface to gnuplot ) by Mario Rodriguez Riotorto http://www.telefonica.net/web2/biomates */ draw2d( title= " Orbit ", user_preamble = "", terminal = 'screen, pic_width = 1000, pic_height = 1000, xlabel = "Z.re ", ylabel = "Z.im", point_type = filled_circle, points_joined = false, color =red, point_size = 0.5, /*key = "critical orbit",*/ points(xx,yy) ) )$ /* ------------------- DrawOrbitDistance ---------------------------------------*/ /* Draws distnce between points of orbit distance = abs(abs(orbit[n])abs(orbit[n-1])) see image : http://commons.wikimedia.org/wiki/Image:Distance.png example of use : orbit:critical_orbit(-1,100); DrawOrbitDistance(orbit); */ DrawOrbitDistance(orbit):= block( distance:GiveDistance(orbit[1],orbit[2]), /* distance between points */ /* save the orbit to 2 lists */ xx:makelist (1, i, 1, 1), /* list of number of iteration */ yy:makelist (distance, i, 1, 1), /* list of distance */ for n:3 thru length(orbit) step 1 do block ( xx:cons(n,xx), yy:cons(GiveDistance(orbit[n],orbit[n-1]),yy) ), load(draw), /* ( interface to gnuplot ) by Mario Rodriguez Riotorto http://www.telefonica.net/web2/biomates */ draw2d( title= " distance between points zn and z(n+1) where z(n+1)=f(zn)= zn*zn+c and z0=0 is critical point ", user_preamble = "", terminal = 'screen, pic_width = 1000, pic_height = 1000, xlabel = "iteration ", ylabel = "distance", point_type = filled_circle, points_joined = false, color =red, point_size = 0.5, points(xx,yy) ) )$