In this post I develop the Euler class, which implements a first order Runge Kutta.
The Euler algorithm is an extremely simple, although very inaccurate, algorithm with an update formula shown below.
k0 = f(xi, yi)
yi+1 = yi + k0h
This algorithm is implemented in the Euler class, which is given at the bottom of this post.
As an example of using this class, I will use the following equations.
x1’ = x1 – 2x2
x2’ = 2x1 + x2
Where x1(0) = 0, x2(0) = 4
These equations have an analytical solution given below.
x1 = -4etsin(2t)
x2 = 4etcos(2t)
To solve these equations using the Euler class, with a step size of 0.01 from 0.0 to 3.3, the code below is used.
Parameter T =
new Parameter(0.0);
Parameter X1 =
new Parameter(0.0);
Parameter X2 =
new Parameter(4.0);
Parameter[] parameters =
new Parameter[] { X1, X2 };
Func<
double>[] functions =
new Func<
double>[]
{
() => X1 - 2.0 * X2,
() => 2.0 * X1 + X2
};
RungeKuttaBase ode = new Euler();
double[][] dRes = ode.Integrate(parameters, T, functions, 3.3, 0.1);
A graph of these results are given below.
Euler step size 0.01
Euler step size 0.10
Euler step size 0.25
Euler step size 0.50
As can be seen from the above graphs, extremely small step sizes are needed to get reasonably accurate results from an Euler algorithm.
The code for the Euler class is given below. The code for the Runge Kutta Base class was posted previously.
public class Euler : RungeKuttaBase
{
public Euler()
:
base()
{
}
public override double[][] Integrate(Parameter[] rungeKuttaParameters, Parameter x, Func<double>[] rungeKuttaFunctions, double xEnd, double step)
{
double xStart = x;
double stepSize = (xEnd - xStart) / step;
int integerStepSize = (int)stepSize;
if (stepSize > (double)integerStepSize) { integerStepSize++; }
integerStepSize++;
Collection<double[]> returnCollection = new Collection<double[]>();
double[] row = new double[rungeKuttaParameters.Length + 1];
row[0] = xStart;
for (int i = 0; i < rungeKuttaParameters.Length; i++)
{
row[i + 1] = rungeKuttaParameters[i];
}
returnCollection.Add(row);
double currentX = x;
double[] dCurrentFs = new double[rungeKuttaParameters.Length];
double currentStep = step;
while (currentX < xEnd)
{
if (currentX + currentStep > xEnd) { currentStep = xEnd - currentX; }
x.Value = currentX;
for (int j = 0; j < rungeKuttaParameters.Length; j++)
{
dCurrentFs[j] = rungeKuttaFunctions[j]();
}
row = new double[rungeKuttaParameters.Length + 1];
for (int j = 0; j < rungeKuttaParameters.Length; j++)
{
rungeKuttaParameters[j].Value = rungeKuttaParameters[j] + dCurrentFs[j] * currentStep;
row[j + 1] = rungeKuttaParameters[j];
}
currentX += currentStep;
row[0] = currentX;
returnCollection.Add(row);
}
return returnCollection.ToArray();
}
}
Subscribe