A third order Runge Kutta can be used to integrate a system of simultaneous differential equations.
The update formula for the third order Runge Kutta algorithm is shown below.
This algorithm is implemented in the third order Runge Kutta 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
x1 = -4etsin(2t)
x2 = 4etcos(2t)
To solve these equations using the third order Runge Kutta 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 third order Runge Kutta class gives much more accurate results than the second order Runge Kutta classes.
The code for the third order Runge Kutta 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