The fern, Lorentz, and title fractals were drawn using the excellent freeware program Fractint. The latest version of this is available from here.
The Sierpinski Gasket transforms were created using The Desktop Fractal Design System from Iterated Systems Inc.
The phase portraits and Poincare sections were created in Euler for IBM OS/2, a freeware linear algebra package, using the scripts listed below. Excuse the rough and ready nature. The lack of comments is partially due to programming constraints.
pendphase.e
comment
Phase diagram of a non-linear pendulum.
endcomment
function pendphase
## pendphase draws a phase space diagram for the non-linear pendulum.
## The function requests parameters and initial conditions.
## It then solves the coupled ODE using a 4th order Runge-Kutta algorithm.
“** PHASE SPACE DIAGRAM OF A NON-LINEAR PENDULUM **”,
” “,
q = input(“Damping strength (q): “);
g = input(“Driver strength (g): “);
f = input(“Driver frequency (f): “);
” “,
“Initial Conditions (>-pi, pi)”,***
v = input(“theta: “);
w = input(“omega: “);
” “,
N = input(“Number of time steps (200): “);***
dt = input(“Step size: “);
k = [0,0;0,0;0,0;0,0];
hold off
plot (0,0)
setplot(-4,4,-4,4)
hold on
for t = 0 to N step dt;
k[1,1] = dt * thetadt(w);
k[1,2] = dt * omegadt(t,v,w,q,g,f);
k[2,1] = dt * thetadt(w+0.5*k[1,2]);
k[2,2] = dt * omegadt(t+dt/2,v+0.5*k[1,1],w+0.5*k[1,2],q,g,f);
k[3,1] = dt * thetadt(w+0.5*k[2,2]);
k[3,2] = dt * omegadt(t+dt/2,v+0.5*k[2,1],w+0.5*k[2,2],q,g,f);
k[4,1] = dt * thetadt(w+k[3,2]);
k[4,2] = dt * omegadt(t+dt,v+k[3,1],w+k[3,2],q,g,f);
v = v + (k[1,1] + 2 * k[2,1] + 2 * k[3,1] + k[4,1])/6;
w = w + (k[1,2]+ 2 * k[2,2] + 2 * k[3,2] + k[4,2])/6;
if v > pi;
v = -v – 2 * pi;
endif;
if v < -pi;
v = v + 2 * pi;
endif;
if w > pi;
w = w – 2 * pi;
endif;
if w < -pi;
w = w + 2 * pi;
endif;
x = v;
y = w;
if t > 20;
xplot(x,y);
endif;
end;
return 0;
endfunction;
function thetadt(omega)
## The time derivative of theta.
## Equal to omega.
return omega;
endfunction;
function omegadt(t,theta,omega,q,g,f)
## The time derivative of omega.
return ((-1/q) * omega – sin(theta) + g * cos(f*t));
endfunction;
pendpoinc.e
comment
Phase diagram of a non-linear pendulum.
endcomment
function pendpoinc
## pendpoinc draws a Poincare section of the non-linear pendulum.
## The function requests parameters and initial conditions.
## It then solves the coupled ODE using a Runge-Kutta algorithm.
“** PHASE SPACE DIAGRAM OF A NON-LINEAR PENDULUM **”,
” “,
q = input(“Damping strength (q): “);
g = input(“Driver strength (g): “);
” “,
“Initial Conditions (>pi)”, v=”input(“theta:” “); w=”input(“omega:” “); ” “, N=”input(“Number” of time steps (>-pi,
200): “);
dt = 3*pi/100;
k = [0,0;0,0;0,0;0,0];
f = 2/3;
hold off
plot (0,0)
setplot(-4,4,-4,4)
hold on
for t = 0 to N step dt;
i = i + dt;
k[1,1] = dt * thetadt(w);
k[1,2] = dt * omegadt(t,v,w,q,g,f);
k[2,1] = dt * thetadt(w+0.5*k[1,2]);
k[2,2] = dt * omegadt(t+dt/2,v+0.5*k[1,1],w+0.5*k[1,2],q,g,f);
k[3,1] = dt * thetadt(w+0.5*k[2,2]);
k[3,2] = dt * omegadt(t+dt/2,v+0.5*k[2,1],w+0.5*k[2,2],q,g,f);
k[4,1] = dt * thetadt(w+k[3,2]);
k[4,2] = dt * omegadt(t+dt,v+k[3,1],w+k[3,2],q,g,f);
v = v + (k[1,1] + 2 * k[2,1] + 2 * k[3,1] + k[4,1])/6;
w = w + (k[1,2]+ 2 * k[2,2] + 2 * k[3,2] + k[4,2])/6;
if v > pi;
v = -v – 2 * pi;
endif;
if v < -pi;
v = v + 2 * pi;
endif;
if w > pi;
w = w – 2 * pi;
endif;
if w < -pi;
w = w + 2 * pi;
endif;
x = v;
y = w;
if t > 20 and i = 3 * pi;
xplot(x,y);
i = 0
endif;
end;
return 0;
endfunction;
function thetadt(omega)
## The time derivative of theta.
## Equal to omega.
return omega;
endfunction;
function omegadt(t,theta,omega,q,g,f)
## The time derivative of omega.
return ((-1/q) * omega – sin(theta) + g * cos(f*t));
endfunction;