Robert Dodier pisze:
> Try changing the definition of f to f(z, c) := expand (z*z + c);
> I think that could speed up your computations.
> It might also fix the stack overflow problem.
Yes it works great, thx.
>
> Some other random comments --- I recommend constructing a function
> instead of a program, so you can run it repeatedly with different
> parameters. I also recommend separating the plotting stuff from the
> computational stuff (i.e. compute list of points in one function and the
> plot in another).
>
> Hope this helps
Thx. You are right, of course I'm programming in the simplest way,
almost without functions.
Step from batch file to functions is good but needs time and pefect
understanding of notation ( which I do not have ).
Adam
Here is quick version :
========== code =====================
/* compute c for given radius, angle and l ( lambda=multiplier of fixed
point) */
epsilon:-0.001;
radius:1+epsilon;
angle:(3-sqrt(5))/2; /* in radians */
l:radius * exp(%i*angle*2*%pi);
c:l/2 - (l*l)/4;
trigsimp(%);
rectform(%);
c:float(%), numer;
/* define function ( map) for dynamical system z(n+1)=f(zn,c) */
f(z,c):=expand (z*z + c); /* expand speed up computations and fix the
stack overflow problem. Robert Dodier */
/* maximal number of iterations */
iMax:1000; /* to big couses bind stack overflow */
EscapeRadius:10;
/* define z-plane ( dynamical ) */
zxMin:-0.8;
zxMax:0.2;
zyMin:-0.2;
zyMax:0.8;
/* resolution is proportional to number of details and time of drawing */
iXmax:1000;
iYmax:1000;
/* compute critical point */
zcr:rhs(solve(diff(f(z,c),z,1)));
/* save critical point to 2 lists */
xcr:makelist (realpart(zcr), i, 1, 1); /* list of re(z) */
ycr:makelist (imagpart(zcr), i, 1, 1); /* list of im(z) */
/* ------------------- compute forward orbit of critical point ----------*/
z:zcr; /* first point */
orbit:[z];
for i:1 thru iMax step 1 do
block
(
z:f(z,c),
if abs(z)>EscapeRadius then return(i) else orbit:endcons(z,orbit)
);
/*-------------- save orbit to draw it later on the screen
----------------------------- */
/* save the z values to 2 lists */
xx:makelist (realpart(zcr), i, 1, 1); /* list of re(z) */
yy:makelist (imagpart(zcr), 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)
);
/* draw reversed orbit of beta using draw package */
load(draw);
draw2d(
title= concat(" dynamical plane for f(z,c):=z*z+",string(c)),
user_preamble = "set key outside top",
terminal = 'screen,
pic_width = iXmax,
pic_height = iYmax,
yrange = [zyMin,zyMax],
xrange = [zxMin,zxMax],
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),
key = "critical point",
color =blue,
points_joined = false,
points(xcr,ycr)
);