# lorenz-circuit.R
#
# Simulate the Lorenz circuit described in [1].
#
# K. Myneni, 2010-02-03, Creative Consulting for Research & Education
#
# Notes:
#
# 1) This script is written for the free, multi-platform, R computing environment.
# See [2] to obtain R, and [3] for notes on using R. Simply load this script
# using,
#
# source("lorenz-circuit.R")
#
# from within R, and two graphical windows displaying figures similar to
# fig. 2 and fig. 3, from ref. [1] will be computed and displayed.
#
# 2) The add-on package for R, odesolve, must be installed. See ref. [3] for
# instructions. You may also refer to the documentation provided with the
# odesolve package to understand how to setup and solve ordinary differential
# equations.
#
# References:
#
# [1] N. J. Corron, 2010, A Simple Circuit Implementation of a Chaotic Lorenz System,
# http://ccreweb.org/documents/physics/chaos/LorenzCircuit3.html
#
# [2] http://www.r-project.org
#
# [3] K. Myneni, Notes on Using R, http://ccreweb.org/documents/programming/R-notes.html
#
## Load the add-on R package, odesolve, for solving ordinary differential
## equations.
require("odesolve")
## Define the differential equations to be solved (eqns. 4 from ref. [1])
derivs <- function( t, x, p ) {
xdot = p[1]*(-x[1] + x[2])
ydot = (p[2]/3)*(3 - x[3])*x[1] - x[2]
zdot = -p[3]*x[3] + x[1]*x[2]
return( list( c( xdot, ydot, zdot ) ) )
}
## Set parameters used in [1]
sigma = 10
R = 30
b = 8/3
params = c( sigma, R, b )
## Pick some initial values of the variables, near the attractor
x0 = c( -1, 1, 0 )
## Set up an array of time values for the solution
times <- seq(0, 1000, by = 0.01)
## Solve the differential equations, and create arrays, t, x, y, z
## from the solution. Skip past the transient region, where the system
## is moving from its initial point onto the attractor.
out <- lsoda( x0, times, derivs, params )
n = length(out[,1])
range = 5000:7000
t = out[range, 1]
x = out[range, 2]
y = out[range, 3]
z = out[range, 4]
## Make a plot to show the time series of x, y, and z, similar to
## figure 2 in [1].
par(mfrow=c(3,1))
plot(t, x, type="l", col="blue")
plot(t, y, type="l", col="red")
plot(t, z, type="l", col="green")
## Open a new graphics output window and plot two views of the
## attractor, as shown in figure 3 of [1], using more points from
## the output to show the structure of the attractor more clearly.
x11()
range2 = 5000:10000
t = out[range2, 1]
x = out[range2, 2]
y = out[range2, 3]
z = out[range2, 4]
par(mfrow=c(1,2))
plot(x, y, type="l", col="blue")
plot(x, z, type="l", col="red")