Fifth order Runge Kutta Butcher

by Trent Guidry8. October 2009 04:54
In this article I develop a fifth order Runge Kutta Butcher.
 
The update formula for the fifth order Runge Kutta Butcher algorithm is shown below.
Fifth Order Runge Kutta Butcher
This algorithm is implemented in the Runge Kutta 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 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.
Runge Kutta 5 100
Step size of 0.25.
Runge Kutta 5 250
Step size of 0.5.
Runge Kutta 5 500
Step size of 1.0.
Runge Kutta 5 1000
As can be seen, the Runge Kutta Butcher class gives more accurate results than the fourth order Runge Kutta class.
The code for the Runge Kutta 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];}//k0sx.Value = x0;for (int i = 0; i < rungeKuttaParameters.Length; i++){ks[0, i] = rungeKuttaFunctions[i]();}//k1sx.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]();}//k2sx.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]();}//k3sx.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]();}//k4sx.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]();}//k5sx.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]();}//FinalcurrentX += 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();}}