## 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 readibility 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 atmostphere 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 4^{th} 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 symetric 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 familar 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 visualise 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++)
{
I1 = CalculateLorenzIComponent(X1, Y1, T1);
J1 = CalculateLorenzJComponent(X1, Y1, Z1, T1);
K1 = CalculateLorenzKComponent(X1, Y1, Z1, T1);
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));
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));
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));
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:

## History

- 26/09/2009 - Initial creation