Click here to Skip to main content
Click here to Skip to main content

The Butterfly Effect

, 26 Sep 2009
Rate this:
Please Sign up or sign in to vote.
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 Smile | :)

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();
        QuadraticBezierSegment quadraticBezierSegment = new QuadraticBezierSegment();
    
        if (pathFigureCollection.Count > 0)
        {
            var oldPathFigure = pathFigureCollection[n - 1];
            if (oldPathFigure != null)
                pathFigure.StartPoint = 
        ((QuadraticBezierSegment)oldPathFigure.Segments[0]).Point2;
        }
        float[] point1 = MatrixVectorMultiply(TransformMatrix, v.Origin.ToArray());
        float[] point2 = MatrixVectorMultiply(TransformMatrix, v.Destination.ToArray());
   
        quadraticBezierSegment.Point1 = new Point(
                point1[0] / (1 - (point1[2] / FocalLength)) + originX,
                point1[1] / (1 - (point1[2] / FocalLength)) + originY);
        quadraticBezierSegment.Point2 = new Point(
                (point2[0] / (1 - (point2[2] / FocalLength))) + originX,
                point2[1] / (1 - (point2[2] / FocalLength)) + originY);
     
        pathFigure.Segments.Add(quadraticBezierSegment);
    
        pathFigureCollection.Add(pathFigure);
        n++;
    }
    pathGeometry.Figures = pathFigureCollection;
    Path path = new Path();
    path.Stroke = scb;
    path.StrokeThickness = 1;
    path.Data = pathGeometry;
    c.Children.Add(path);
}

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

License

This article, along with any associated source code and files, is licensed under The Code Project Open License (CPOL)

About the Author

CatchExAs
Technical Lead Alpha Integralis Limited
United Kingdom 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
Follow on   LinkedIn

Comments and Discussions

 
Questionthanks Pinmemberkartalkartal10-Jan-13 7:56 
AnswerRe: thanks PinprofessionalCatchExAs10-Apr-14 6:47 
GeneralMy vote of 5 PinmemberFilip D'haene26-May-11 9:30 
Generalsolving differential equations is a standard problem Pinmemberheadmyshoulder27-Sep-09 7:40 
GeneralRe: solving differential equations is a standard problem PinmemberCatchExAs27-Sep-09 10:49 

General General    News News    Suggestion Suggestion    Question Question    Bug Bug    Answer Answer    Joke Joke    Rant Rant    Admin Admin   

Use Ctrl+Left/Right to switch messages, Ctrl+Up/Down to switch threads, Ctrl+Shift+Left/Right to switch pages.

| Advertise | Privacy | Mobile
Web04 | 2.8.140721.1 | Last Updated 27 Sep 2009
Article Copyright 2009 by CatchExAs
Everything else Copyright © CodeProject, 1999-2014
Terms of Service
Layout: fixed | fluid