In this post I develop the Runge Kutta 5 Butcher class, which implements a fifth order Runge Kutta Butcher.
The update formula for the fifth order Runge Kutta Butcher algorithm is shown below.
k0 = f(xi, yi)
k1 = f(xi + 1/4 h, yi + 1/4 k0h)
k2 = f(xi + 1/4 h, yi + (1/8 k0 + 1/8 k1)h)
k3 = f(xi + 1/2 h, yi + (-1/2 k1 + k2)h)
k4 = f(xi + 3/4 h, yi + (3/16 k0 + 9/16 k3)h)
k5 = f(xi + h, yi + (-3/7 k0 + 2/7 k1 + 12/7 k2 - 12/7 k3 + 8/7 k4)h)
yi+1 = yi + ((7k0 + 32k2 + 12k3 + 32k4 + 7k5) /90)h
This algorithm is implemented in the Runge Kutta 5 Butcher 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 5 Butcher 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 RungeKutta5Butcher();
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.
Step size of 1.0.
As can be seen, the Runge Kutta 5 Butcher class gives more accurate results than the fourth order Runge Kutta class.
The code for the Runge Kutta 5 Butcher class is given below. The code for the Runge Kutta Base class was posted previously.
public class RungeKutta5Butcher : RungeKuttaBase
{
public RungeKutta5Butcher()
:
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[6, 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 / 4.0;
for (int i = 0; i < rungeKuttaParameters.Length; i++)
{
rungeKuttaParameters[i].Value = y0[i] + (ks[0, i] / 4.0) * currentStep;
}
for (int i = 0; i < rungeKuttaParameters.Length; i++)
{
ks[1, i] = rungeKuttaFunctions[i]();
}
//k2s
x.Value = x0 + currentStep / 4.0;
for (int i = 0; i < rungeKuttaParameters.Length; i++)
{
rungeKuttaParameters[i].Value = y0[i] + (ks[0, i] / 8.0 + ks[1, i] / 8.0) * currentStep;
}
for (int i = 0; i < rungeKuttaParameters.Length; i++)
{
ks[2, i] = rungeKuttaFunctions[i]();
}
//k3s
x.Value = x0 + currentStep / 2.0;
for (int i = 0; i < rungeKuttaParameters.Length; i++)
{
rungeKuttaParameters[i].Value = y0[i] + (-ks[1, i] / 2.0 + ks[2, i]) * currentStep;
}
for (int i = 0; i < rungeKuttaParameters.Length; i++)
{
ks[3, i] = rungeKuttaFunctions[i]();
}
//k4s
x.Value = x0 + currentStep * 3.0 / 4.0;
for (int i = 0; i < rungeKuttaParameters.Length; i++)
{
rungeKuttaParameters[i].Value = y0[i] + (ks[0, i] * 3.0 / 16.0 + ks[3, i] * 9.0 / 16.0) * currentStep;
}
for (int i = 0; i < rungeKuttaParameters.Length; i++)
{
ks[4, i] = rungeKuttaFunctions[i]();
}
//k5s
x.Value = x0 + currentStep;
for (int i = 0; i < rungeKuttaParameters.Length; i++)
{
rungeKuttaParameters[i].Value = y0[i] + (-ks[0, i] * 3.0 / 7.0 + ks[1, i] * 2.0 / 7.0 + ks[2, i] * 12.0 / 7.0 - ks[3, i] * 12.0 / 7.0 + ks[4, i] * 8.0 / 7.0) * currentStep;
}
for (int i = 0; i < rungeKuttaParameters.Length; i++)
{
ks[5, 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] + ((7.0 * ks[0, i] + 32.0 * ks[2, i] + 12.0 * ks[3, i] + 32.0 * ks[4, i] + 7.0 * ks[5, i]) / 90.0) * currentStep;
row[i + 1] = rungeKuttaParameters[i];
}
returnCollection.Add(row);
}
return returnCollection.ToArray();
}
}
Subscribe