'Program ABSCHAOS.BAS makes bifurcation diagram for chaotic flow with abs term '(c) 1998 by J. C. Sprott 'Compile with PowerBASIC for DOS 3.5 defext a-z screen 12 amin=-.8 amax=-.5 xmin=-1 xmax=2 x=0: y=0: z=0 h=.002 tmax&=1000000 line(10,0)-(630,460),,b locate 1,1: print trim\$(str\$(csng(xmax))); locate 15,1: print"X"; locate 29,1: print trim\$(str\$(csng(xmin))); locate 30,1: print trim\$(str\$(csng(amin))); locate 30,41: print"A"; locate 30,81-len(trim\$(str\$(csng(amax)))): print trim\$(str\$(csng(amax))); xticks%=6: yticks%=6 for i%=1 to xticks%-1 xp%=10+620*i%/xticks% line(xp%,1)-(xp%,11) line(xp%,449)-(xp%,459) next i% for i%=1 to yticks%-1 yp%=460-460*i%/yticks% line(11,yp%)-(21,yp%) line(619,yp%)-(629,yp%) next i% for i%=11 to 628 a=amin+(amax-amin)*(i%-10)/620 for t&=0 to tmax& call rk4(x,y,z,h) if abs(x)>1e3 then x=0: y=0: z=0: exit for if t&>tmax&/2 and xxolder then call extremum (xolder, xold, x, dum, xe) j%=460-460*(xe-xmin)/(xmax-xmin) pset(i%,j%) end if xolder=xold xold=x next t& if inkey\$=chr\$(27) then end next i% while inkey\$="": wend end def fnx(x,y,z)=y def fny(x,y,z)=z def fnz(x,y,z)=a*z-y+abs(x)-1 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 extremum (yl, yo, yh, xe, ye) 'Given three values yl @ x=-1, yo @ x=0 and yh @ x=1 'Calculate the extremum ye and the value xe at which ye occurs 'Fits the data to a parabola y = yo + bx + cx^2 b = (yh - yl) / 2 c = (yh + yl) / 2 - yo IF c = 0 THEN 'Second derivative is zero (no extremum) xe = 0 ye = yo ELSE xe = -b / (2 * c) ye = yo - b * b / (4 * c) END IF END SUB