|
|
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.ecomment 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; |