Scheme code for the Lorenz attractor.

lorenz.rkt

;;;; Scheme code to plot the chaotic behaviour of the Lorenz attractor.
;;;; Date: 09/01/2016

#lang racket
(require plot)

;;; The parameters for the attractor (these values were used by Lorenz to
;;; demonstrate the chaotic behaviour.
(define rho 28)
(define sigma 10)
(define beta 8/3)

;;; The Lorenz attractor is described by a coupled system of
;;; three ordinary differential equations:
(define (dxdt x)
  (match-let (((list x1 x2 x3) (vector->list x)))
    (vector (* sigma (- x2 x1))
            (- (* x1 (- rho x3)) x2)
            (- (* x1 x2) (* beta x3)))))

;;; We use the 4th-order Runge-Kutta method to solve the
;;; ODE system with a step size of h=0.01.
(define h 0.01)
(define (runge-kutta t-final x-init)
  (define (vect-mult c v) (vector-map (lambda (vi) (* c vi)) v))
  (define (vect-add v1 v2 . vs)
    (apply vector-map (append (list + v1 v2) vs)))
  (let ((xs (list x-init)))
    (for ((tk (in-range 0 t-final h)))
      (let* ((xk (car xs))
             (f1 (vect-mult h (dxdt xk)))
             (f2 (vect-mult h (dxdt (vect-add xk (vect-mult 1/2 f1)))))
             (f3 (vect-mult h (dxdt (vect-add xk (vect-mult 1/2 f2)))))
             (f4 (vect-mult h (dxdt (vect-add xk f3)))))
        (set! xs 
              (cons 
               (vect-add xk 
                         (vect-mult 1/6 (vect-add f1 
                                                  (vect-mult 2 f2) 
                                                  (vect-mult 2 f3) 
                                                  f4))) 
               xs))))
    (reverse xs)))

;;; We then plot the solution values in phase-space.
(parameterize ((plot-width 400)
               (plot-height 400)
               (plot-font-size 11)
               (plot-x-label "x")
               (plot-y-label "y")
               (plot-z-label "z")
               (plot3d-angle 132)
               (plot3d-altitude 16))
  (plot3d
   (lines3d (runge-kutta 60 (vector 0 1 0)) #:color 3)))