Resetting Cash-Karp Integrator

2025-07-26 20:34 UTC gpt-4o Open in ChatGPT ↗

You

what would be the best way to reset this integrator, without re instantiating the class?

using System; using UnityEngine;

namespace UnchartedLimbo.ChaosOfSun.ThreeBodyOrbits { ///

/// Cash–Karp adaptive Runge–Kutta integrator (RK45) for systems dy/dt = f(t, y). ///

public class CashKarpIntegrator { ///

/// Signature for your derivative function: given (t, y) fill dydt. ///

public delegate void DerivativeFunc(double t, double[] y, double[] dydt);

    /// <summary>Current independent variable (time).</summary>
    public double T { get; private set; }

    /// <summary>Current state vector (length N).</summary>
    public double[] Y { get; private set; }

    /// <summary>Current step size.</summary>
    public double H { get; private set; }

    readonly double tol;
    readonly double hMin, hMax;
    readonly DerivativeFunc derivs;
    readonly int n;

    // Cash–Karp coefficients
    private static readonly double[] C = {
        0.0,
        1.0 / 5.0,
        3.0 / 10.0,
        3.0 / 5.0,
        1.0,
        7.0 / 8.0 };

    // Cash–Karp coefficients (flat A: 6 rows × 5 cols)
    const int A_ROWS = 6;
    const int A_COLS = 5;
    private static readonly double[] A = new double[A_ROWS * A_COLS]
    {
        // row 0 (unused)
        0.0,       0.0,       0.0,       0.0,       0.0,
        // row 1
        1.0/5.0,    0.0,       0.0,       0.0,       0.0,
        // row 2
        3.0/40.0,   9.0/40.0,   0.0,       0.0,       0.0,
        // row 3
        3.0/10.0,  -9.0/10.0,   6.0/5.0,    0.0,       0.0,
        // row 4
        -11.0/54.0,  5.0/2.0,  -70.0/27.0,  35.0/27.0,  0.0,
        // row 5
        1631/55296.0,
        175/512.0,
        575/13824.0,
        44275/110592.0,
        253/4096.0
    };

    private static readonly double[] B5 = { 37.0 / 378, 0, 250.0 / 621, 125.0 / 594, 0, 512.0 / 1771 };
    private static readonly double[] B4 = { 2825.0 / 27648, 0, 18575.0 / 48384, 13525.0 / 55296, 277.0 / 14336, 1.0 / 4 };

    // workspace
    readonly double[] k1, k2, k3, k4, k5, k6, yTemp, yOut;

    /// <summary>
    /// Constructor.
    /// </summary>
    /// <param name="y0">Initial state (will be copied).</param>
    /// <param name="derivs">Your dy/dt function.</param>
    /// <param name="dtStart">Starting step size.</param>
    /// <param name="tol">Error tolerance.</param>
    /// <param name="dtMin">Minimum allowed step.</param>
    /// <param name="dtMax">Maximum allowed step.</param>
    public CashKarpIntegrator(
        double[] y0,
        DerivativeFunc derivs,
        double dtStart = 0.01,
        double tol = 1e-10,
        double dtMin = 1e-7,
        double dtMax = 0.1
    )
    {
        n = y0.Length;
        Y = new double[n];
        Array.Copy(y0, Y, n);
        this.derivs = derivs;

        T = 0f;
        H = dtStart;
        this.tol = tol;
        hMin = dtMin;
        hMax = dtMax;

        k1 = new double[n];
        k2 = new double[n];
        k3 = new double[n];
        k4 = new double[n];
        k5 = new double[n];
        k6 = new double[n];
        yTemp = new double[n];
        yOut = new double[n];
    }

    /// <summary>
    /// Advance the solution until T >= tEnd (calling Step internally).
    /// </summary>
    public void StepTo(double tEnd)
    {
        while (T < tEnd)
        {
            // clamp last bit so we end exactly at tEnd
            H = Math.Min(H, tEnd - T);
            Step();
        }
    }

    /// <summary>
    /// Perform a single adaptive step: compute new Y, update T and H.
    /// </summary>
    void Step()
    {
        derivs(T, Y, k1);
        ComputeStage(1, k1, k2);
        ComputeStage(2, k1, k3);
        ComputeStage(3, k1, k4);
        ComputeStage(4, k1, k5);
        ComputeStage(5, k1, k6);

        double errMax = 0.0;
        for (int i = 0; i < n; i++)
        {
            double yi5 = Y[i] + H * (
                B5[0]*k1[i] +
                B5[1]*k2[i] +
                B5[2]*k3[i] +
                B5[3]*k4[i] +
                B5[4]*k5[i] +
                B5[5]*k6[i]
            );
            double yi4 = Y[i] + H * (
                B4[0]*k1[i] +
                B4[1]*k2[i] +
                B4[2]*k3[i] +
                B4[3]*k4[i] +
                B4[4]*k5[i] +
                B4[5]*k6[i]
            );
            yOut[i] = yi5;
            double err = Math.Abs(yi5 - yi4);

            if (err > errMax)
                errMax = err;
        }

        double scale = errMax > 0.0
            ? 0.9 * Math.Pow(tol / errMax, 0.2)
            : 5.0;

        double hNew = Math.Clamp(scale * H, hMin, hMax);

        if (errMax <= tol || H <= hMin)
        {
            // either “good enough” or “we’re as small as we dare go” ⇒ accept
            Array.Copy(yOut, Y, n);
            T += H;
            H = hNew;
        }
        else
        {
            // too big error but we can still shrink: retry
            H = hNew;
        }
    }

    /// <summary>
    /// Perform one of the 5 intermediate RK stages.
    /// </summary>
    private void ComputeStage(int stageIndex, double[] kPrev, double[] kOut)
    {
        double tStage = T + (A[stageIndex * A_COLS + 0] /*unused*/ + /*but*/0f) * H;
        // actually C array handles times, so you'd still use C[stageIndex], but omitted here
        for (int j = 0; j < n; j++)
        {
            double sum = 0.0;
            for (int m = 0; m < stageIndex; m++)
            {
                double coeff = A[stageIndex * A_COLS + m];
                double kj = (m == 0) ? kPrev[j] : GetStageKValue(m, j);
                sum += coeff * kj;
            }
            yTemp[j] = Y[j] + H * sum;
        }
        derivs(tStage, yTemp, kOut);
    }

    /// <summary>
    /// Returns the j-th component of the m-th stage array:
    /// m=1→k2, 2→k3, 3→k4, 4→k5, 5→k6.
    /// </summary>
    private double GetStageKValue(int m, int j)
    {
        return m switch
        {
            1 => k2[j],
            2 => k3[j],
            3 => k4[j],
            4 => k5[j],
            5 => k6[j],
            _ => throw new ArgumentOutOfRangeException($"Invalid RK stage index: {m}. expected 1..5")
        };
    }

}

}

ChatGPT

To reset the CashKarpIntegrator without re-instantiating the class, you’d ideally want a method that:

• Resets the time T to 0
• Resets the current state Y to a new initial state y0
• Resets the step size H (optional, but useful)

✅ Recommended Approach

Add a Reset method to the class like this:

/// <summary>
/// Resets the integrator to a new initial state and time.
/// </summary>
/// <param name="y0">New initial state vector (will be copied).</param>
/// <param name="t0">New initial time (default 0).</param>
/// <param name="dtStart">Optional new step size (if you want to change it).</param>
public void Reset(double[] y0, double t0 = 0.0, double? dtStart = null)
{
    if (y0.Length != n)
        throw new ArgumentException($"Expected length {n}, got {y0.Length}", nameof(y0));

    Array.Copy(y0, Y, n);
    T = t0;

    if (dtStart.HasValue)
        H = Math.Clamp(dtStart.Value, hMin, hMax);
}

🔧 Usage Example

double[] newInitialState = { 1.0, 0.0, 0.0 }; // Example
integrator.Reset(newInitialState, t0: 0.0, dtStart: 0.01);

---

🔒 Why this is safe and efficient:

• No need to reallocate or garbage collect internal buffers (k1, k2, …, yOut).
• Keeps the object’s memory footprint stable and performant.
• Respects your encapsulation and avoids class-level restructuring.

Let me know if you want it to reset more internal state or reset entirely back to constructor values.