ITERATED FUNCTION SYSTEMS (1993 by J. C. Sprott)

Code: (Select All)
100 Rem ITERATED FUNCTION SYSTEMS Ver 1.0 (c) 1993 by J. C. Sprott

110 DefDbl A-Z: Dim A(12)

120 Randomize Timer 'Reseed random number generator

130 Screen 12 'Assume VGA graphics

140 GoSub 300 'Set parameters

150 GoSub 400 'Iterate equations

160 GoSub 500 'Display results

170 GoSub 600 'Test results

180 On T% GOTO 130, 140, 150

190 End



300 Rem Set parameters

310 X = .05: Y = .05: XE = X + .00001: YE = Y

320 For I% = 1 To 12: A(I%) = .1 * (Int(25 * Rnd) - 12): Next I%

330 J0 = Abs(A(1) * A(4) - A(2) * A(3))

340 J1 = Abs(A(7) * A(10) - A(8) * A(9))

350 If J0 + J1 = 0 Or J0 > 1 Or J1 > 1 Then GoTo 320 'Not contracting

360 P = J0 / (J0 + J1)

370 T% = 3: LSUM = 0: N = 1: N1 = 0: N2 = 0

380 XMIN = 1000000#: XMAX = -XMIN: YMIN = XMIN: YMAX = XMAX

390 Return



400 Rem Iterate equations

410 If X <> XE Then If Rnd > P Then R% = 6 Else R% = 0

420 XNEW = A(1 + R%) * X + A(2 + R%) * Y + A(5 + R%)

430 YNEW = A(3 + R%) * X + A(4 + R%) * Y + A(6 + R%)

490 Return



500 Rem Display results

510 If N < 100 Or N > 1000 Then GoTo 560

520 If X < XMIN Then XMIN = X

530 If X > XMAX Then XMAX = X

540 If Y < YMIN Then YMIN = Y

550 If Y > YMAX Then YMAX = Y

560 If N = 1000 Then GoSub 800

570 If X > XL And X < XH And Y > YL And Y < YH And N > 1000 Then PSet (X, Y)

590 Return



600 Rem Test results

610 GoSub 700 'Calculate Lyapunov exponent

620 GoSub 900 'Calculate correlation dimension

630 If N > 21000 Then T% = 2 'Candidate IFS found

640 If Abs(XNEW) + Abs(YNEW) > 1000000# Then T% = 2 'Unbounded

650 If N > 998 And (F < 1 Or L > -.2) Then T% = 2 'Uninteresting

660 If Len(InKey$) Then T% = 0 'User key press

670 X = XNEW: Y = YNEW: N = N + 1

690 Return



700 Rem Calculate Lyapunov exponent

710 XSAVE = XNEW: YSAVE = YNEW: X = XE: Y = YE

720 GoSub 400 'Reiterate equations

730 DLX = XNEW - XSAVE: DLY = YNEW - YSAVE: DL2 = DLX * DLX + DLY * DLY

740 DF = 10000000000# * DL2: RS = 1# / Sqr(DF)

750 XE = XSAVE + RS * (XNEW - XSAVE): YE = YSAVE + RS * (YNEW - YSAVE)

760 XNEW = XSAVE: YNEW = YSAVE

770 LSUM = LSUM + Log(DF): L = .721347 * LSUM / N

790 Return



800 Rem Resize the screen (and discard the first thousand iterates)

810 DX = .1 * (XMAX - XMIN): DY = .1 * (YMAX - YMIN)

820 XL = XMIN - DX: XH = XMAX + DX: YL = YMIN - DY: YH = YMAX + DY

830 If XH - XL < .000001 Or YH - YL < .000001 Then GoTo 890

840 Window (XL, YL)-(XH, YH): Cls

850 Line (XL, YL)-(XH, YH), , B

890 Return



900 Rem Calculate fractal dimension

910 If N < 200 Then GoTo 990 'Wait for transient to settle

920 If N = 200 Then D2MAX = (XMAX - XMIN) ^ 2 + (YMAX - YMIN) ^ 2

930 If N = 200 Or Rnd < .02 Then XS = X: YS = Y 'New reference point

940 DX = XNEW - XS: DY = YNEW - YS

950 D2 = DX * DX + DY * DY

960 If D2 < .001 * D2MAX Then N2 = N2 + 1

970 If D2 < .00001 * D2MAX Then N1 = N1 + 1: F = .434294 * Log(N2 / N1)

990 Return