'Program REGION3D.BAS plots chaotic regions for 3-D system Sprott & Horton '(c) 1999 by J. C. Sprott 'Compile with PowerBASIC for DOS 3.5 '\$error all on SCREEN 12 defext a - z shared a, c randomize timer line(0, 0) - (639, 479), ,b call insolve(1, 1, 638, 478) while inkey\$ <> chr\$(27): wend end SUB calculate(xp%, yp%, c%) h = 0.1 'Iteration step size dr = 1e-8 'Perturbation for LE calculation a = (479 - yp%) / 479 c = 10^(4*xp% / 639) x=0: y=0: z=c1*x xe = x + dr: ye = y: ze = z chaos% = 0 iter& = 0: lsum = 0: ub% = 0 for j% = 1 to 32000 call rk4(x, y, z, h) call rk4(xe, ye, ze, h) dx = xe - x: dy = ye - y: dz = ze - z dsq = dx * dx + dy * dy + dz * dz rs = sqr(dsq) / dr if rs > 0 then xe = x + dx / rs: ye = y + dy / rs: ze = z + dz / rs if j% > 2000 then lsum = lsum + log(rs): incr iter& end if if abs(x) > 1e5 then ub% = -1: exit for next j% if ub% = 0 then le! = lsum / (iter& * h) c% = 7 'limit cycle if le! > .01 then c% = 15 'chaotic if le! < -.01 then c% = 0 'fixed point else c% = 8 'unbounded end if END SUB def fnx(x,y,z) = y def fny(x,y,z) = z def fnz(x,y,z) = -a*z-Y+1-c*(1+tanh(x)) SUB rk4 (x, y, z, h) 'Advances (x, y, z) with a fourth order Runge-Kutta method 'User must provide the functions fnx, fny, fnz, and a step size h k1x = h * fnx(x, y, z) k1y = h * fny(x, y, z) k1z = h * fnz(x, y, z) k2x = h * fnx(x + .5 * k1x, y + .5 * k1y, z + .5 * k1z) k2y = h * fny(x + .5 * k1x, y + .5 * k1y, z + .5 * k1z) k2z = h * fnz(x + .5 * k1x, y + .5 * k1y, z + .5 * k1z) k3x = h * fnx(x + .5 * k2x, y + .5 * k2y, z + .5 * k2z) k3y = h * fny(x + .5 * k2x, y + .5 * k2y, z + .5 * k2z) k3z = h * fnz(x + .5 * k2x, y + .5 * k2y, z + .5 * k2z) k4x = h * fnx(x + k3x, y + k3y, z + k3z) k4y = h * fny(x + k3x, y + k3y, z + k3z) k4z = h * fnz(x + k3x, y + k3y, z + k3z) x = x + (k1x + 2 * (k2x + k3x) + k4x) / 6 y = y + (k1y + 2 * (k2y + k3y) + k4y) / 6 z = z + (k1z + 2 * (k2z + k3z) + k4z) / 6 END SUB SUB insolve (x1%, y1%, x2%, y2%) 'Paints the screen at successively higher resolution d%=min(x2%-x1%,y2%-y1%) d%=2^int(log2(d%)) while d%>0 d2%=2*d% for x%=x1% to x2%-d%+1 step d% old%=x% 30 then tanh = sgn(x) else tanh = 1 - 2 / (EXP(2 * x) + 1) end if END FUNCTION