In this post I develop the Runge Kutta 2 Improved Polygon class, which implements a second order Runge Kutta Improved Polygon.
The update formula for the Runge Kutta Improved Polygon algorithm is shown below.
k0 = f(xi, yi)
k1 = f(xi + 1/2 h, yi + 1/2 hk0)
yi+1 = yi + k1h
This algorithm is implemented in the Runge Kutta 2 Improved Polygon 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 2 Improved Polygon 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 RungeKutta2ImprovedPolygon();
double[][] dRes = ode.Integrate(parameters, T, functions, 3.3, 0.1);
A graph of these results are given below.
As can be seen, the Runge Kutta 2 Improved Polygon class gives much more accurate results than the Euler class with a step size of 0.1.
The code for the Runge Kutta 2 Improved Polygon class is given below. The code for the Rung Kutta Base class was posted previously.
public class RungeKutta2ImprovedPolygon : RungeKuttaBase
{
public RungeKutta2ImprovedPolygon()
:
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[2, 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]();
}
//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[1, i] * currentStep;
row[i + 1] = rungeKuttaParameters[i];
}
returnCollection.Add(row);
}
return returnCollection.ToArray();
}
}
Subscribe