**Who said chaos can’t be modelled on an ancient machine with RAM limited to just several kilobytes? Execution time may be hours or weeks, but with a bit of careful programming, fractals and strange attractors will be calculated to the same resolution as on more capable graphic workstations.**

The code at the end of this post is an attempt to solve the system of differential equations of the Lorenz attractor on a Casio fx-730P programmable calculator. The machine is low on computational power but has a fair amount of memory for programs and data (8 kb!), so I decided on a quite sophisticated numerical algorithm that offers good accuracy at a relatively low computational cost. The method is the Cash-Karp variant of the fifth-order Runge-Kutta method with adaptive stepsize control. Cash-Karp is an all-round method widely used in scientific computing. The downside is that the code is quite long; with all control and output routines it takes almost 2 kilobytes on the little Casio!

Before going into details of the code, let’s take a look at the results of the program:

The image shows the first 5000 data points of the Lorenz attractor calculated on the fx-730P. Points were joined with straight line segments and 3D-plotted in Wolfram Mathematica. The exact Mathematica command was this:

lines = Table[{Hue[.45 + i/20000], Line[{attr[[i]], attr[[i+1]]}]}, {i, 1, Length[attr] - 1}]; Graphics3D[{Thickness[.002], lines}, Background -> Black, AspectRatio -> .65, Axes -> True]

The “attr” variable is a list of 5000 points in 3D-space. This list was manually created from the output of the Casio program!

Now for some details about the Casio program that was used to calculate this. The program lets the user specify the initial value for the calculation. The initial point doesn’t really matter because the same attractor is reached from whatever point you begin, but to produce the image shown above, I started at (0,1,1). The initial stepsize has to be provided – it again doesn’t matter too much because the algorithm will adjust the stepsize according to the prescribed accuracy. I took 0.01. The absolute error tolerance (accuracy) of the method is hard-coded in a variable called L (line number 4170). I found that the accuracy of 10^-6 is achieved at a stepsize so large that there would be surprisingly long gaps between calculated points in a 3D-plot. It shows how accurate the Cash-Karp method is, but to generate enough points for a smooth 3D-plot, an upper limit is set on the stepsize (line number 4180, variable V).

When executed, the program displays one calculated point at a time and waits until EXE is pressed. It gives you an opportunity to write down the coordinates on paper – the fx-730P doesn’t have enough memory to store all the results and the coordinates are overwritten after each step! One point takes around 12 seconds to calculate, unless the prescribed accuracy is not reached – in which case the calculation is repeated with a shorter stepsize.

The program asks in the beginning for the number of points to be calculated. The next time the program is executed, it lets you continue the calculation from the last point found previously, or start over with a new initial condition. This is how an output may look like:

The three numbers on the display are, well, x- y- and z- coordinates of a newly calculated point along the attractor. Square brackets at the beginning of the line provide some diagnostics about the stepsize:

- “[ ]” is shown when the stepsize was unchanged or increased from the previous step;
- “[M]” indicates that the stepsize was forced to the hard-coded maximum;
- “[*]” and a beep means that the stepsize was reduced by the adaptive algorithm.

Below is the complete code of the Lorenz attractor solver. When stored in calculator’s memory, the code takes up 1898 bytes (excluding REM-comments at the end). The program uses several one- and two-dimensional arrays, which when allocated, occupy another 504 bytes of RAM.

1 REM Lorenz attractor 2 REM solved with Runge-Kutta 5th order 3 REM For Casio FX-730P 10 INPUT "[N]ew [C]ontinue",M$ 20 IF M$="N" THEN GOSUB 4001 30 INPUT "How many steps",P 40 FOR I=0 TO 2 50 X(I)=Y(I) 60 NEXT I 70 FOR K=1 TO P 80 U$="[ ]": IF S=V THEN U$="[M]" 90 GOSUB 9001 100 FOR I=0 TO 2 110 G(I)=F(I) 120 NEXT I 130 GOSUB 2001 140 M=0 150 FOR I=0 TO 2 160 IF ABS(E(I))>M THEN M= ABS(E(I)) 170 NEXT I 180 M=M/L 190 IM M<=1 THEN 250 200 U$="[*]": BEEP 210 IF M>6561 THEN S=.1*S 220 IF M<=6561 THEN S=.9*S*M^-.25 230 IF T+S=T THEN PRINT "stepsize underflow" 240 GOTO 130 250 T=T+S 260 FOR I=0 TO 2 270 Y(I)=X(I) 280 NEXT I 290 IF M<=1.89E-4 THEN S=5*S 300 IF M>1.89E-4 THEN S=.9*S*M^-.2 310 IF ABS(S)> ABS(V) THEN S=V 320 PRINT U$; RND(X(0),-4); RND(X(1),-4); RND(X(2),-4) 330 NEXT K 999 END 2001 REM RK4 algorithm routine 2002 REM Cash-Karp Runge-Kutta step 2010 FOR I=0 TO 2 2020 K(0,I)=G(I) 2030 X(I)=Y(I)+.2*S*K(0,I) 2040 NEXT I 2050 GOSUB 9001 2060 FOR I=0 TO 2 2070 K(1,I)=F(I) 2080 X(I)=Y(I)+S*(.075*K(0,I)+.225*K(1,I)) 2090 NEXT I 2100 GOSUB 9001 2110 FOR I=0 TO 2 2120 K(2,I)=F(I) 2130 X(I)=Y(I)+S*(.3*K(0,I)-.9*K(1,I)+1.2*K(2,I)) 2140 NEXT I 2150 GOSUB 9001 2160 FOR I=0 TO 2 2170 K(3,I)=F(I) 2180 X(I)=Y(I)+S*(C(0)*K(0,I)+C(1)*K(1,I)) 2190 X(I)=X(I)+S*(C(2)*K(2,I)+C(3)*K(3,I)) 2200 NEXT I 2210 GOSUB 9001 2220 FOR I=0 TO 2 2230 K(4,I)=F(I) 2240 X(I)=Y(I)+S*(C(4)*K(0,I)+C(5)*K(1,I)) 2250 X(I)=X(I)+S*(C(6)*K(2,I)+C(7)*K(3,I)+C(8)*K(4,I)) 2260 NEXT I 2270 GOSUB 9001 2280 FOR I=0 TO 2 2290 X(I)=Y(I)+S*(C(9)*K(0,I)+C(10)*K(2,I)+C(11)*K(3,I)) 2300 X(I)=X(I)+S*C(12)*F(I) 2310 E(I)=S*(C(13)*K(0,I)+C(14)*K(2,I)+C(15)*K(3,I)) 2320 E(I)=E(I)+S*(C(16)*K(4,I)+C(17)*F(I)) 2330 NEXT I 2340 RETURN 4001 REM Initialize memory 4010 CLEAR 4020 DIM X(3): DIM F(3): DIM G(3) 4030 DIM Y(3): DIM E(3) 4040 DIM K(5,3): DIM C(18) 4050 C(0)=-11/54:C(1)=2.5:C(2)=-70/27:C(3)=35/27 4060 C(4)=1631/55296:C(5)=175/512:C(6)=575/13824 4070 C(7)=44275/110592:C(8)=253/4096:C(9)=37/378 4080 C(10)=250/621:C(11)=125/594:C(12)=512/1771 4090 C(13)=C(9)-2825/27648:C(14)=C(10)-18575/48384 4100 C(15)=C(11)-13525/55296:C(16)=-277/14336 4110 C(17)=C(12)-.25 4120 INPUT "Initial stepsize",S 4130 INPUT "X=",Y(0) 4140 INPUT "Y=",Y(1) 4150 INPUT "Z=",Y(2) 4160 T=0 4170 L=1E-6 4180 V=.01 4190 RETURN 9001 REM Derivatives of the Lorenz attractor ODEs 9010 F(0)=-10*(X(0)-X(1)) 9020 F(1)=28*X(0)-X(1)-X(0)*X(2) 9030 F(2)=X(0)*X(1)-8*X(2)/3 9040 RETURN 10001 REM Variables 10002 REM S stepsize (changed after each step) 10003 REM P number of steps requested 10004 REM L error tolerance 10005 REM T independent variable (incremented after each step) 10006 REM M max absolute error of a step 10007 REM V largest step allowed 10008 REM Y(3) dependent variables x, y, z 10009 REM X(3) temporary storage for new values of Y() 10010 REM F(3) rhs of the attractor ODEs 10011 REM C(18) Cash-Karp coefficients 10012 REM E(3) Error estimate (calculated after each step) 10013 REM K(5,3) Runge-Kutta increments 10014 REM G(3) Temp storage for F() at the beginning of a step 10015 REM U$ Set to "[*]" and beeps if stepsize reduced 10016 REM U$ Set to "[M]" if working at max allowed stepsize