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;
    }
}