13,147,540 members (32,145 online)
alternative version

#### Stats

29.2K views
36 bookmarked
Posted 26 Sep 2009

# The Butterfly Effect

, 26 Sep 2009
 Rate this:
An investigation of Chaos Theory Attractors, using a 4th order Runge-Kutta solver

## Introduction

I've always been fascinated with the idea of using computers for solving big problems, so rather than write yet another tool or technique for building crud applications, I hope this article helps to show some of the more fun things you can do with C# and .NET.

## Background

When writing a Brief History of Time, Stephen Hawking was advised that the readability of his book was inversely proportional to the number of equations that he presented. Unfortunately, the background to this code is quite math heavy and a full treatment of it can be viewed here at my website. I suspect the audience on this site is savvy enough to get to grips with at least the basics.

Many dynamic systems in classical physics can be described fully using Ordinary Differential Equations of 2 types. Initial Value problems, where the start conditions tend to have profound effects on the solutions. As an example, imagine a bowling ball in 10-pin bowling. Where the pins end up is highly dependent on the throwers initial input. Boundary Values usually deal with the same situations but have additional constraints imposed upon them. For instance, waves crashing against walls.

During the 70s, Ed Lorenz was studying cloud formation using computational techniques to solve these types of problems. His differential equations dealt with fluid flow in the atmosphere that he was able to split into 3 dimensions:

• dx/dt = s(y - x) where s = Prandtl number
• dy/dt = x(p - z) - y where p = Rayleigh number
• dz/dt = xy - Bz

By using a method known as the 4th Order Runge-Kutta routine, he was able to plot the solutions and model them as they evolved with time. Whilst evaluating the effect of changing the initial conditions, Lorenz noticed that the shapes plotted would vary quite radically with very small changes. Even more, he found a range of 'stable' initial values that would plot quite beautiful and symmetric shapes. Known as the Lorenz attractor, it is often mistakenly known as the butterfly of the 'Butterfly Effect'.

A Runge Kutta routine, essentially boils down to sequentially generating a co-ordinate plot in each dimension. If you are familiar with concepts such as the Taylor, Binomial and Fibonacci series, then the Runge-Kutta routine should be easy to understand. Solving for a time-step instance n of Y:

Yn+1 = Yn + k1/6 + k2/3 + k3/3 + k4/6 + 0 (h^5)

where:

k1 = hf'(Xn, Yn)
k2 = hf'(Xn + h/2, Yn + k1/2)
k3 = hf'(Xn + h/2, Yn + k2/2)
k4 = hf'(Xn + h, Yn + k3)

h, is a constant representing the time step length.
f', represents the differential component in terms of X and Y that Lorenz derived at the beginning of the article :)

OK, enough maths and time for some code!

We can rewrite Lorenz's equations to solve for each co-ordinate in space-time using C# thus:

private float CalculateLorenzIComponent(float x, float y, float t)
{
float I = t * ((-_sigma * x) + (_sigma * y));
return I;
}
private float CalculateLorenzJComponent(float x, float y, float z, float t)
{
float K = t * ((_rho * x) - y - (x * z));
return K;
}
private float CalculateLorenzKComponent(float x, float y, float z, float t)
{
float L = t * ((x * y) - (_beta * z));
return L;
}

But this is a time-step calculation, so we need a looping function that increments our t variable and calculates the k1, k2, k3 and k4 terms for larger equations. What's more, we would like to be able to visualize the solution in a graph so we should create some kind of list of drawable objects for WPF to then render. The driver for this is the LorenzRoutine() function below, which with not a little imagination is a good candidate for some refactoring with delegates or generics.

for (int n = 0; n < _iterations; n++)
{
//1st order approximation
I1 = CalculateLorenzIComponent(X1, Y1, T1);
J1 = CalculateLorenzJComponent(X1, Y1, Z1, T1);
K1 = CalculateLorenzKComponent(X1, Y1, Z1, T1);

//2nd order approximation
I2 = CalculateLorenzIComponent(X1 + (h / 2) * I1, Y1 + (h / 2) * J1, T1 + (h / 2));
J2 = CalculateLorenzJComponent(X1 + (h / 2) *
I1, Y1 + (h / 2) * J1, Z1 + (h / 2) * K1, T1 + h / 2);
K2 = CalculateLorenzKComponent(X1 + (h / 2) *
I1, Y1 + (h / 2) * J1, Z1 + (h / 2) * K1, T1 + (h / 2));

//3rd order approximation
I3 = CalculateLorenzIComponent(X1 + (h / 2) * I2, Y1 + (h / 2) * J2, T1 + h / 2);
J3 = CalculateLorenzJComponent(X1 + (h / 2) *
I2, Y1 + (h / 2) * J2, Z1 + (h / 2) * K1, T1 + (h / 2));
K3 = CalculateLorenzKComponent(X1 + (h / 2) *
I2, X1 + (h / 2) * J2, Z1 + (h / 2) * K1, T1 + (h / 2));

//4th order approximation
I4 = CalculateLorenzIComponent(X1 + (h / 2) * I3, Y1 + (h / 2) * J3, T1 + (h / 2));
J4 = CalculateLorenzJComponent(X1 + (h / 2) *
I3, Y1 + (h / 2) * J3, Z1 + (h / 2) * K1, T1 + (h / 2));
K4 = CalculateLorenzKComponent(X1 + (h / 2) *
I3, X1 + (h / 2) * J3, Z1 + (h / 2) * K1, T1 + (h / 2));

//Taylor Series Expansion in 3 dimensions
float X2 = X1 + (h / 6) * (I1 + 2 * I2 + 2 * I3 + I4);
float Y2 = Y1 + (h / 6) * (J1 + 2 * J2 + 2 * J3 + J4);
float Z2 = Z1 + (h / 6) * (K1 + 2 * K2 + 2 * K3 + J4);

if (n > 0)
{
BuildVector3D(new Point3D(X1, Y1, Z1),
new Point3D(X2, Y2, Z2), 1, 0x000000, 100, 0x000000, 100);
}

X1 = X2;
Y1 = Y2;
Z1 = Z2;
}

The BuildVector3D() function very simply adds a new 3-Dimensional vector onto a list of objects that is later used in RenderScreen() to draw each graph element to the AttractorCanvas WPF Control.

public void RenderScreen(System.Windows.Controls.Canvas c, float originX, float originY)
{
System.Windows.Media.SolidColorBrush scb
= new SolidColorBrush(System.Windows.Media.Colors.Black);

PathGeometry pathGeometry = new PathGeometry();
PathFigureCollection pathFigureCollection = new PathFigureCollection();

int n = 0;
foreach (Vector3D v in this.drawableObjects)
{
PathFigure pathFigure = new PathFigure();

if (pathFigureCollection.Count > 0)
{
var oldPathFigure = pathFigureCollection[n - 1];
if (oldPathFigure != null)
pathFigure.StartPoint =
}
float[] point1 = MatrixVectorMultiply(TransformMatrix, v.Origin.ToArray());
float[] point2 = MatrixVectorMultiply(TransformMatrix, v.Destination.ToArray());

point1[0] / (1 - (point1[2] / FocalLength)) + originX,
point1[1] / (1 - (point1[2] / FocalLength)) + originY);
(point2[0] / (1 - (point2[2] / FocalLength))) + originX,
point2[1] / (1 - (point2[2] / FocalLength)) + originY);

n++;
}
pathGeometry.Figures = pathFigureCollection;
Path path = new Path();
path.Stroke = scb;
path.StrokeThickness = 1;
path.Data = pathGeometry;
}

RenderScreen()'s main responsibility is to take our graph of vectors and perform all necessary steps to paint the graph on the canvas. We are planning to render a scene in WPF, for which we will choose a Quadric Bezier segment to try and attain some smoothness. Next, because we wish to rotate the image we perform a matrix multiplication on the vector in 3-dimensions via the MatrixVectorMultiply method, a discussion of which is probably worth an article in itself. Finally to cater for zoom, both points in the Bezier segment are the subject of a focal length calculation to expand or contract the segment's length.

## Using the Code

The rest of the code in the download is really just plumbing to get the graphs painted to the screen and allow the user to vary the initial conditions, rotate and regenerate the image. Of course, due to the real time nature of the screen rendering, be careful not to specify too many iterations or ask the app to rotate the graph too quickly else your PC will begin to choke.

## Points of Interest

I have provided a second example, which I originally wrote in the MathCAD functional language back in 1999 to show the Rossler Attractor. Hopefully this should demonstrate how the solver should be refactored and made more generic for these types of problems. Rossler's equations model chemical kinetics and give rise to an aesthetically pleasing spiral that seems to be hanging half off an unseen table ledge.

Further readings can be found at:

If you're into Silverlight and want see how to pass event notifications back and forth between it and JavaScript, then see here for the web version.

## History

• 26/09/2009 - Initial creation

## Share

 Technical Lead Alpha Integralis Limited United Kingdom
CatchExAs aka Nick Wilton runs Alpha Integralis, a software development consultancy to companies in the City of London.

Main interests include exploring Algo-Trading, Quant Development and Physics using C# and Java.

www.nickwilton.info/Blog

## You may also be interested in...

 First Prev Next
 thanks kartalkartal10-Jan-13 7:56 kartalkartal 10-Jan-13 7:56
 Re: thanks CatchExAs10-Apr-14 6:47 CatchExAs 10-Apr-14 6:47
 My vote of 5 Filip D'haene26-May-11 9:30 Filip D'haene 26-May-11 9:30
 solving differential equations is a standard problem headmyshoulder27-Sep-09 7:40 headmyshoulder 27-Sep-09 7:40
 Re: solving differential equations is a standard problem CatchExAs27-Sep-09 10:49 CatchExAs 27-Sep-09 10:49
 Last Visit: 31-Dec-99 18:00     Last Update: 24-Sep-17 4:38 Refresh 1