<< Runge Kutta Butcher | Latest version of the source code >>

Runge Kutta Fehlberg

by Trent Guidry 9. October 2009 07:33

In this post I develop the RungeKutta54Fehlberg class, which implements a fifth/fourth order Runge Kutta Fehlberg with adaptive step sizing.

The update formula for the fifth/fourth order Runge Kutta Fehlberg algorithm is shown below.

k₀ = f(x_i, y_i)

k₁ = f(x_i + 1/4 h, y_i + 1/4 k₀h)

k₂ = f(x_i + 3/8 h, y_i + (3/32 k₀ + 9/32 k₁)h)

k₃ = f(x_i + 12/13 h, y_i + (1932/2197 k₀ – 7200/2197 k₁ + 7296/2197 k₂)h)

k₄ = f(x_i + h, y_i + (439/216 k₀ - 8 k₁ +3680/513 k₂ -845/4104 k₃)h)

k₅ = f(x_i + 1/2 h, y_i + (-8/27 k₀ + 2 k₁ - 3544/2565 k₂ + 1859/4104 k₃ - 11/40 k₄)h)

y_i+1 = y_i + (25/216 k₀ + 1408/2565 k₂ + 2197/4104 k₃ – 1/5 k₄)h

z_i+1 = z_i + (16/135 k₀ + 6656/12825 k₂ + 28561/56430 k₃ – 9/50 k₄ + 2/55 k₄)h

Where y is a fourth order Runge Kutta and z is a fifth order Runge Kutta. An estimate of the error can be obtained by subtracting the two values obtained. If the error exceeds a specified threshold, the results can be recalculated using a smaller step size.

The approach to estimating the new step size is shown below and is the same as the one at the link Runge Kutta Fehlberg Method.

h_new = h_old(εh_old/2|z_i+1– y_i+1|)¼

This algorithm is implemented in the Runge Kutta 54 Fehlberg class, which is given at the bottom of this post.

Unlike the other Runge Kutta classes posted previously, this class has a number of different constructors. These constructors allow for the specification of the maximum error and also whether to use the fifth order Runge Kutta or the fourth order Runge Kutta to compute the new values.

As an example of using this class, I will use the following equations.

x₁’ = x₁ – 2x₂

x₂’ = 2x₁ + x₂

Where x₁(0) = 0, x₂(0) = 4

These equations have an analytical solution given below.

x₁ = -4e^tsin(2t)

x₂ = 4e^tcos(2t)

For this example, I will set the step size to the full range and let the estimated error control the actual steps used in the algorithm. To solve these equations using the Runge Kutta 54 Fehlberg class, with a maximum error of 0.001 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 RungeKutta54Fehlberg(0.001);
double[][] dRes = ode.Integrate(parameters, T, functions, 3.3, 0.1);

Maximum error of 0.001.

Runge Kutta Fehlberg 0001

Maximum error of 0.01.

Runge Kutta Fehlberg 0010

Maximum error of 0.1.

Runge Kutta Fehlberg 0100

Maximum error of 1.0.

Runge Kutta Fehlberg 1000

The code for the Runge Kutta 54 Fehlberg class is given below. The code for the Runge Kutta Base class was posted previously.

    public class RungeKutta54Fehlberg : RungeKuttaBase
    {
        private double[] _maximumNumberOfErrors;
        private double _maxAllErrors = double.NaN;
        private bool _useFifthOrderEstimate = true;

        public RungeKutta54Fehlberg()
            : base()
        {
        }

        public RungeKutta54Fehlberg(bool useFifthOrderEstimate)
        {
            _useFifthOrderEstimate = useFifthOrderEstimate;
        }

        public RungeKutta54Fehlberg(double maxAllErrors)
        {
            _maxAllErrors = maxAllErrors;
        }

        public RungeKutta54Fehlberg(double maxAllErrors, bool useFifthOrderEstimate)
        {
            _maxAllErrors = maxAllErrors;
            _useFifthOrderEstimate = useFifthOrderEstimate;
        }

        public RungeKutta54Fehlberg(double[] maxErrors)
        {
            _maximumNumberOfErrors = maxErrors;
        }

        public RungeKutta54Fehlberg(double[] maxErrors, bool useFifthOrderEstimate)
        {
            _maximumNumberOfErrors = maxErrors;
            _useFifthOrderEstimate = useFifthOrderEstimate;
        }

        public override double[][] Integrate(Parameter[] rungeKuttaParameters, Parameter x, Func<double>[] rungeKuttaFunctions, double xEnd, double step)
        {
            Debug.Assert(rungeKuttaParameters.Length == rungeKuttaFunctions.Length);
            if (_maximumNumberOfErrors == null && !double.IsNaN(_maxAllErrors))
            {
                _maximumNumberOfErrors = new double[rungeKuttaParameters.Length];
                for (int i = 0; i < _maximumNumberOfErrors.Length; i++) { _maximumNumberOfErrors[i] = _maxAllErrors; }
            }

if (_maximumNumberOfErrors != null) { Debug.Assert(_maximumNumberOfErrors.Length == rungeKuttaParameters.Length); }

            double[] dErrors = new double[rungeKuttaParameters.Length];
            double[] ys = new double[rungeKuttaParameters.Length];
            double[] zs = new double[rungeKuttaParameters.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 j = 0; j < rungeKuttaParameters.Length; j++)
                {
                    y0[j] = rungeKuttaParameters[j];
                }

                //k0s
                x.Value = x0;
                for (int j = 0; j < rungeKuttaParameters.Length; j++)
                {
                    ks[0, j] = rungeKuttaFunctions[j]();
                }

                //k1s
                x.Value = x0 + currentStep / 4.0;
                for (int j = 0; j < rungeKuttaParameters.Length; j++)
                {
                    rungeKuttaParameters[j].Value = y0[j] + ks[0, j] * currentStep / 4.0;
                }
                for (int j = 0; j < rungeKuttaParameters.Length; j++)
                {
                    ks[1, j] = rungeKuttaFunctions[j]();
                }

                //k2s
                x.Value = x0 + currentStep * 3.0 / 8.0;
                for (int j = 0; j < rungeKuttaParameters.Length; j++)
                {
                    rungeKuttaParameters[j].Value = y0[j] + (ks[0, j] * 3.0 / 32.0 + ks[1, j] * 9.0 / 32.0) * currentStep;
                }
                for (int j = 0; j < rungeKuttaParameters.Length; j++)
                {
                    ks[2, j] = rungeKuttaFunctions[j]();
                }

                //k3s
                x.Value = x0 + currentStep * 12.0 / 13.0;
                for (int j = 0; j < rungeKuttaParameters.Length; j++)
                {
                    rungeKuttaParameters[j].Value = y0[j] + (ks[0, j] * 1932.0 / 2197.0 - ks[1, j] * 7200.0 / 2197.0 + ks[2, j] * 7296.0 / 2197.0) * currentStep;
                }
                for (int j = 0; j < rungeKuttaParameters.Length; j++)
                {
                    ks[3, j] = rungeKuttaFunctions[j]();
                }

                //k4s
                x.Value = x0 + currentStep;
                for (int j = 0; j < rungeKuttaParameters.Length; j++)
                {
                    rungeKuttaParameters[j].Value = y0[j] + (ks[0, j] * 439.0 / 216.0 - ks[1, j] * 8.0 + ks[2, j] * 3680.0 / 513.0 - ks[3, j] * 845.0 / 4104.0) * currentStep;
                }
                for (int j = 0; j < rungeKuttaParameters.Length; j++)
                {
                    ks[4, j] = rungeKuttaFunctions[j]();
                }

                //k5s
                x.Value = x0 + currentStep / 2.0;
                for (int j = 0; j < rungeKuttaParameters.Length; j++)
                {
                    rungeKuttaParameters[j].Value = y0[j] + (-ks[0, j] * 8.0 / 27.0 + ks[1, j] * 2.0 - ks[2, j] * 3544.0 / 2565.0 + ks[3, j] * 1859.0 / 4104.0 - ks[4, j] * 11.0 / 40.0) * currentStep;
                }
                for (int j = 0; j < rungeKuttaParameters.Length; j++)
                {
                    ks[5, j] = rungeKuttaFunctions[j]();
                }

                //Calculate estimated errors.
                for (int j = 0; j < rungeKuttaParameters.Length; j++)
                {
                    ys[j] = y0[j] + (ks[0, j] * 25.0 / 216.0 + ks[2, j] * 1408.0 / 2565.0 + ks[3, j] * 2197.0 / 4104.0 - ks[4, j] * 1.0 / 5.0) * currentStep;
                    zs[j] = y0[j] + (ks[0, j] * 16.0 / 135.0 + ks[2, j] * 6656.0 / 12825.0 + ks[3, j] * 28561.0 / 56430.0 - ks[4, j] * 9.0 / 50.0 + ks[5, j] * 2.0 / 55.0) * currentStep;
                    dErrors[j] = Math.Abs(ys[j] - zs[j]);
                }
                bool bErrorExceeded = false;
                double minErrorRatio = double.PositiveInfinity;
                if (_maximumNumberOfErrors != null)
                {
                    for (int j = 0; j < rungeKuttaParameters.Length; j++)
                    {
                        if (dErrors[j] > _maximumNumberOfErrors[j])
                        {
                            bErrorExceeded = true;
                        }
                        if (minErrorRatio > Math.Abs(_maximumNumberOfErrors[j] / (2.0 * (ys[j] - zs[j]))))
                        {
                            minErrorRatio = Math.Abs(_maximumNumberOfErrors[j] / (2.0 * (ys[j] - zs[j])));
                        }
                    }
                }

                if (!bErrorExceeded)
                {
                    //Final
                    currentX += currentStep;
                    row = new double[rungeKuttaParameters.Length + 1];
                    row[0] = currentX;
                    for (int j = 0; j < rungeKuttaParameters.Length; j++)
                    {
                        rungeKuttaParameters[j].Value = (_useFifthOrderEstimate) ? zs[j] : ys[j];
                        row[j + 1] = rungeKuttaParameters[j];
                    }
                    returnCollection.Add(row);
                }
                else
                {
                    x.Value = x0;
                    for (int j = 0; j < rungeKuttaParameters.Length; j++)
                    {
                        rungeKuttaParameters[j].Value = y0[j];
                    }
                }
                if (_maximumNumberOfErrors != null)
                {
                    double newStep = currentStep * Math.Pow(minErrorRatio * currentStep, 0.25);
                    if (!bErrorExceeded || newStep < currentStep)
                    {
                        currentStep *= Math.Pow(minErrorRatio * currentStep, 0.25);
                    }
                    else
                    {
                        currentStep /= 2.0;
                    }
                    if (currentStep > step) { currentStep = step; }
                }
            }
            return returnCollection.ToArray();
        }
    }

Tags: Runge Kutta, Differential Equations, Numerical Methods, C#

Runge Kutta ButcherIn this post I develop the RungeKutta5Butcher class, which implements a fifth order Runge Kutta Butc...Fourth order Runge KuttaIn this post I develop the RungeKutta4 class, which implements a fourth order Runge Kutta.Third order Runge KuttaIn this post I develop the Runge Kutta 3 class, which implements a third order Runge Kutta.

Add comment

Name*
E-mail*
Website
Country Country flag

b i u quote

Comment
Preview

Notify me when new comments are added

Runge Kutta Fehlberg

TextBox

Tag cloud