|
Example - Lorenz Equations
* The "\" symbol has been used at the end of lines to indicate that the line was too long for the display and has been
wrapped to the line below in order to fit within the width of the web page. These lines would need to be joined in the normal
program editor.
I. Introduction
This program demonstrates the solutions of the Lorenz equations, which model convective currents in fluids and demonstrate chaotic behaviour.
A finite element array is used to perform the simulation.
II. Mathematics Involved
The Lorenz system of differential equations is  , where a , b and c are constants. If the values a=28 , b=8/3 , c=10 are used (the default values in the program), the famous Lorenz butterfly is generated when the solution is viewed in the x-z plane. The solution is chaotic because a small variation in the initial conditions produces a large variation in the subsequent behaviour.
Usually, b and c are held constant and a is varied. If this is the case then the following scenarios apply:
a<1: solution decays rapidly to the origin
a>1: solution converges to one of two fixed points
For larger values of a, some solutions approach a strange attractor (like the “butterfly”).
As a gets even bigger, the nonlinearity is increased and the solutions are more likely to approach a strange attractor than converge to a fixed point.
III. How the Program Works
1) Initialisation of the Finite Element Array
The differential equations are easily translated into a one-dimensional finite element array:
'set up the finite element array
IF EXISTS ("lorenz") THEN Lorenz.Delete()
GLOBAL f AS FiniteElmtArray
f.Name = "lorenz"
f.dimensions(1) = 1
DIM fa AS FiniteElmtAttribute
f.add("x")
f.add("y")
f.add("z")
GLOBAL err%:err% = 0.05
SET fa = f.x
fa.SetInitialValue(xinit%)
fa.Order = 1
fa.SetFormula("-c%*(x()-y())",1)
fa.AllowableErrorFactor = err%
fa.AllowableErrorMin = err%
SET fa = f.y
fa.SetInitialValue(yinit%)
fa.Order = 1
fa.SetFormula("a%*x()-y()-x()*z()",1)
fa.AllowableErrorFactor = err%
fa.AllowableErrorMin = err%
SET fa = f.z
fa.SetInitialValue(zinit%)
fa.Order = 1
fa.SetFormula("b%*(x()*y()-z())",1)
fa.AllowableErrorFactor = err%
fa.AllowableErrorMin = err%
f.AdvanceTo(0)
The range of the graph is then calculated (roughly) by running the finite element array for 15 seconds and using the maximum and minimum values of x and z over the period as the limits for the plot. This is also done whenever the view is changed to a different combination of axes, e.g. x and y.
2) The advance() Procedure
If either of the Go buttons is clicked, the program advances the finite element array, then displays the next part of the solution curve as a line (by calling the DrawLine() procedure). To save on memory and display time, only one Line object is actually used: it is moved each time but its old location is not refreshed, only its new one. Please note that this is not a generally recommended approach to displaying things on the screen. It was chosen because there is no way to easily limit the number of line objects that may eventually be created. It suffers from a significant problem in that if any object on the screen covers up part of the display, the result will be a blank area (colored white) since the objects will no longer exist for the purpose of redrawing the parts of the display that have been covered up.
SUB advance()
'advance the finite element array
'change enabled/disabled buttons
DIM tt%,oldx%,oldy%,oldz%,newx%,newy%,newz%
oldx% = f.x()
oldy% = f.y()
oldz% = f.z()
tt% = f.t + 0.002
IF Lor.tlab.Caption <> "t=" + STR$ (tt%,"999.9") THEN
'only change time label if the time has changed
'when expressed to one decimal place
Lor.tlab.Caption = "t=" + STR$ (tt%,"999.9")
Lor.Page.SetDirtyRect(220,500,280,520)
Lor.Refresh()
END IF
f.AdvanceTo(tt%)
newx% = f.x()
newy% = f.y()
newz% = f.z()
CALL DrawLine(oldx%,oldy%,oldz%,newx%,newy%,newz%,tt%)
END SUB
The results of the program after some 165 seconds are as shown below:
The program has two Go buttons to enable the user to view the effects of chaos. For example, if the simulation is run for 10 seconds with the default initial conditions, then stopped, then the x0 value changed by 0.0001 and the simulation restarted in red, it can be seen that the trajectories start off very close together.
However, after only a fairly short period of time they are following quite different paths.
|
