In this post I develop the RungeKutta3 class, which implements a third order Runge Kutta.
The update formula for the third order Runge Kutta algorithm is shown below.
k0 = f(xi, yi)
k1 = f(xi + 1/2 h, yi + 1/2 hk0)
k2 = f(xi + h, yi + (-k0 + 2k1)h)
yi+1 = yi + (1/6 (k0 + 4 k1 + k2))h
This algorithm is implemented in the Runge Kutta 3 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 Runge Kutta 3 class, with a step size of 0.1 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 RungeKutta3();
double[][] dRes = ode.Integrate(parameters, T, functions, 3.3, 0.1);
Step size of 0.1.
Step size of 0.25.
Step size of 0.5.
As can be seen, the Runge Kutta 3 class gives much more accurate results than the second order Runge Kutta classes.
The code for the Runge Kutta 3 class is given below. The code for the Runge Kutta Base class was posted previously.
public class RungeKutta3 : RungeKuttaBase
{
public RungeKutta3()
:
base()
{
}
public override double[][] Integrate(Parameter[] rungeKuttaParameters, Parameter x, Func<double>[] rungeKuttaFunctions, double xEnd, double step)
{
Debug.Assert(rungeKuttaParameters.Length == rungeKuttaFunctions.Length);
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] = x;
for (int i = 0; i < rungeKuttaParameters.Length; i++)
{
row[i + 1] = rungeKuttaParameters[i];
}
returnCollection.Add(row);
double currentX = xStart;
double[,] ks = new double[3, rungeKuttaParameters.Length];
double x0 = xStart;
double[] y0 = new double[rungeKuttaParameters.Length];
double currentStep = step;
while (currentX < xEnd)
{
if (currentX + currentStep > xEnd) { currentStep = xEnd - currentX; }
x0 = currentX;
for (int i = 0; i < rungeKuttaParameters.Length; i++)
{
y0[i] = rungeKuttaParameters[i];
}
//k0s
x.Value = x0;
for (int i = 0; i < rungeKuttaParameters.Length; i++)
{
ks[0, i] = rungeKuttaFunctions[i]();
}
//k1s
x.Value = x0 + currentStep / 2.0;
for (int i = 0; i < rungeKuttaParameters.Length; i++)
{
rungeKuttaParameters[i].Value = y0[i] + ks[0, i] * currentStep / 2.0;
}
for (int i = 0; i < rungeKuttaParameters.Length; i++)
{
ks[1, i] = rungeKuttaFunctions[i]();
}
//k2s
x.Value = x0 + currentStep;
for (int i = 0; i < rungeKuttaParameters.Length; i++)
{
rungeKuttaParameters[i].Value = y0[i] + (-ks[0, i] + 2.0 * ks[1, i]) * currentStep;
}
for (int i = 0; i < rungeKuttaParameters.Length; i++)
{
ks[2, i] = rungeKuttaFunctions[i]();
}
//Final
currentX += currentStep;
row = new double[rungeKuttaParameters.Length + 1];
row[0] = currentX;
for (int i = 0; i < rungeKuttaParameters.Length; i++)
{
rungeKuttaParameters[i].Value = y0[i] + ((ks[0, i] + 4.0 * ks[1, i] + ks[2, i]) / 6.0) * currentStep;
row[i + 1] = rungeKuttaParameters[i];
}
returnCollection.Add(row);
}
return returnCollection.ToArray();
}
}
Subscribe