RungeKuttaMethod.cs

using System;
using System.Collections.Generic;

namespace Algorithms.Numeric
{
    /// <summary>
    ///     In numerical analysis, the Runge–Kutta methods are a family of implicit and explicit iterative methods,
    ///     used in temporal discretization for the approximate solutions of simultaneous nonlinear equations.
    ///     The most widely known member of the Runge–Kutta family is generally referred to as
    ///     "RK4", the "classic Runge–Kutta method" or simply as "the Runge–Kutta method".
    ///     </summary>
    public static class RungeKuttaMethod
    {
        /// <summary>
        ///     Loops through all the steps until xEnd is reached, adds a point for each step and then
        ///     returns all the points.
        /// </summary>
        /// <param name="xStart">Initial conditions x-value.</param>
        /// <param name="xEnd">Last x-value.</param>
        /// <param name="stepSize">Step-size on the x-axis.</param>
        /// <param name="yStart">Initial conditions y-value.</param>
        /// <param name="function">The right hand side of the differential equation.</param>
        /// <returns>The solution of the Cauchy problem.</returns>
        public static List<double[]> ClassicRungeKuttaMethod(
            double xStart,
            double xEnd,
            double stepSize,
            double yStart,
            Func<double, double, double> function)
        {
            if (xStart >= xEnd)
            {
                throw new ArgumentOutOfRangeException(
                    nameof(xEnd),
                    $"{nameof(xEnd)} should be greater than {nameof(xStart)}");
            }

            if (stepSize <= 0)
            {
                throw new ArgumentOutOfRangeException(
                    nameof(stepSize),
                    $"{nameof(stepSize)} should be greater than zero");
            }

            List<double[]> points = new();
            double[] firstPoint = { xStart, yStart };
            points.Add(firstPoint);

            var yCurrent = yStart;
            var xCurrent = xStart;

            while (xCurrent < xEnd)
            {
                var k1 = function(xCurrent, yCurrent);
                var k2 = function(xCurrent + 0.5 * stepSize, yCurrent + 0.5 * stepSize * k1);
                var k3 = function(xCurrent + 0.5 * stepSize, yCurrent + 0.5 * stepSize * k2);
                var k4 = function(xCurrent + stepSize, yCurrent + stepSize * k3);

                yCurrent += (1.0 / 6.0) * stepSize * (k1 + 2 * k2 + 2 * k3 + k4);
                xCurrent += stepSize;

                double[] newPoint = { xCurrent, yCurrent };
                points.Add(newPoint);
            }

            return points;
        }
    }
}