Lorenz attractor on a fx-730P

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:

lorenz

The Lorenz attractor: data points were calculated on Casio fx-730P and imported (manually!) to Wolfram Mathematica for 3d plotting

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:

lorentz1

The FX-730P can only display coordinates of one point at a time. To make a 3D-diagram of the attractor, coordinates have to be manually copied into graph plotting software on a more powerful machine…

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
Advertisements
This entry was posted in FX-730P and tagged , , , , , , , , , . Bookmark the permalink.

Leave a Reply

Fill in your details below or click an icon to log in:

WordPress.com Logo

You are commenting using your WordPress.com account. Log Out / Change )

Twitter picture

You are commenting using your Twitter account. Log Out / Change )

Facebook photo

You are commenting using your Facebook account. Log Out / Change )

Google+ photo

You are commenting using your Google+ account. Log Out / Change )

Connecting to %s