refactor(ode): bundle Fehlberg and step contexts

src/main/java/org/spaceroots/mantissa/ode/AdaptiveStepsizeIntegrator.java

@@ -4,7 +4,7 @@
  * Base class for integrators that adapt their step size while solving ordinary differential
  * equations.
  *
- * <p>The class centralizes step bounds, tolerance bookkeeping and initialization so concrete
+ * <p>The class centralizes step bounds, tolerance bookkeeping, and initialization so concrete
  * algorithms can focus on method-specific interpolation and local error estimation. For each state
  * component {@code i}, the error threshold is {@code absTol_i + relTol_i * max(|y_m|, |y_{m+1}|)};
  * scalar tolerances are broadcast when vector forms are not provided. A tentative step is accepted
@@ -68,16 +68,16 @@ protected AdaptiveStepsizeIntegrator(
   /**
    * Build an integrator with per-component tolerances and step size bounds.
    *
-   * <p>Vector tolerances let callers adapt admissible error to the magnitude of each state
-   * coordinate. Arrays are used as provided; callers must ensure their length matches the problem
-   * dimension. The default step handler performs no action until replaced.
+   * <p>Vector tolerances let callers adapt admissible error to the scale of each state coordinate.
+   * Arrays are used as provided; callers must ensure their length matches the problem dimension.
+   * The default step handler performs no action until replaced.
    *
    * @param minStep minimal allowed step magnitude; positive even when integrating backward.
    * @param maxStep maximal allowed step magnitude; positive and not smaller than {@code minStep}.
-   * @param vecAbsoluteTolerance absolute error tolerance for each component; must align with state
-   *     dimension and contain non-negative entries.
-   * @param vecRelativeTolerance relative error tolerance for each component; must align with state
-   *     dimension and contain non-negative entries.
+   * @param vecAbsoluteTolerance absolute error tolerance for each component; must align with the
+   *     state dimension and contain non-negative entries.
+   * @param vecRelativeTolerance relative error tolerance for each component; must align with the
+   *     state dimension and contain non-negative entries.
    */
   protected AdaptiveStepsizeIntegrator(
       double minStep,
@@ -106,13 +106,12 @@ protected AdaptiveStepsizeIntegrator(
    * Set an explicit initial step size to be reused by the next integration run.
    *
    * <p>A positive value inside {@code [minStep, maxStep]} is honored directly by {@link
-   * #initializeStep(FirstOrderDifferentialEquations, boolean, int, double[], double, double[],
-   * double[], double[], double[])}. Negative or out-of-range values instruct the integrator to
-   * estimate an initial step from the derivatives instead. The stored value is not persisted across
-   * calls to {@link #resetInternalState()}.
+   * #initializeStep(StepInitializationContext)}. Negative or out-of-range values instruct the
+   * integrator to estimate an initial step from the derivatives instead. The stored value is not
+   * persisted across calls to {@link #resetInternalState()}.
    *
-   * @param initialStepSize desired starting step; must be positive and within configured bounds, or
-   *     it is ignored in favor of automatic estimation.
+   * @param initialStepSize the desired starting step must be positive and within configured bounds,
+   *     or it is ignored in favor of automatic estimation.
    */
   @SuppressWarnings("unused")
   public void setInitialStepSize(double initialStepSize) {
@@ -176,73 +175,48 @@ public void addSwitchingFunction(
    * step is derived from the ratio of scaled state and derivative norms, refined with a trial Euler
    * step, and finally clamped to the configured minimum and maximum step sizes.
    *
-   * @param equations system of first-order differential equations providing derivative evaluation;
-   *     must not be {@code null}.
-   * @param forward {@code true} to integrate toward increasing time, {@code false} to integrate
-   *     backward.
-   * @param order order of the integration method that will consume the step; influences the scaling
-   *     of the heuristic.
-   * @param scale per-component scaling factors applied to the state vector; must match {@code y0}
-   *     length and contain only non-zero values.
-   * @param t0 starting independent variable value associated with {@code y0}.
-   * @param y0 state vector at {@code t0}; not modified by this method.
-   * @param yDot0 first derivative of {@code y0}; computed by {@code equations} and reused in place.
-   * @param y1 workspace receiving a trial state after an Euler prediction; length must match {@code
-   *     y0}.
-   * @param yDot1 workspace receiving the derivative at the trial state; length must match {@code
-   *     y0}.
+   * @param context bundled initialization parameters for the step size heuristic.
    * @return signed initial step size bounded by configured limits and oriented according to {@code
-   *     forward}.
-   * @throws DerivativeException if {@code equations} fails while evaluating derivatives at the
-   *     trial point.
+   *     context.forward()}.
+   * @throws DerivativeException if derivative evaluation fails while probing the trial point.
    */
-  public double initializeStep(
-      FirstOrderDifferentialEquations equations,
-      boolean forward,
-      int order,
-      double[] scale,
-      double t0,
-      double[] y0,
-      double[] yDot0,
-      double[] y1,
-      double[] yDot1)
-      throws DerivativeException {
+  public double initializeStep(StepInitializationContext context) throws DerivativeException {
 
     if (initialStep > 0) {
-      // use the user provided value
-      return forward ? initialStep : -initialStep;
+      // use the user-provided value
+      return context.forward() ? initialStep : -initialStep;
     }
 
-    // very rough first guess : h = 0.01 * ||y/scale|| / ||y'/scale||
+    // very rough first guess: h = 0.01 * ||y/scale|| / ||y'/scale||
     // this guess will be used to perform an Euler step
     double ratio;
     double yOnScale2 = 0;
     double yDotOnScale2 = 0;
-    for (int j = 0; j < y0.length; ++j) {
-      ratio = y0[j] / scale[j];
+    for (int j = 0; j < context.y0().length; ++j) {
+      ratio = context.y0()[j] / context.scale()[j];
       yOnScale2 += ratio * ratio;
-      ratio = yDot0[j] / scale[j];
+      ratio = context.yDot0()[j] / context.scale()[j];
       yDotOnScale2 += ratio * ratio;
     }
 
     double h =
         ((yOnScale2 < 1.0e-10) || (yDotOnScale2 < 1.0e-10))
             ? 1.0e-6
             : (0.01 * Math.sqrt(yOnScale2 / yDotOnScale2));
-    if (!forward) {
+    if (!context.forward()) {
       h = -h;
     }
 
     // perform an Euler step using the preceding rough guess
-    for (int j = 0; j < y0.length; ++j) {
-      y1[j] = y0[j] + h * yDot0[j];
+    for (int j = 0; j < context.y0().length; ++j) {
+      context.y1()[j] = context.y0()[j] + h * context.yDot0()[j];
     }
-    equations.computeDerivatives(t0 + h, y1, yDot1);
+    context.equations().computeDerivatives(context.t0() + h, context.y1(), context.yDot1());
 
     // estimate the second derivative of the solution
     double yDDotOnScale = 0;
-    for (int j = 0; j < y0.length; ++j) {
-      ratio = (yDot1[j] - yDot0[j]) / scale[j];
+    for (int j = 0; j < context.y0().length; ++j) {
+      ratio = (context.yDot1()[j] - context.yDot0()[j]) / context.scale()[j];
       yDDotOnScale += ratio * ratio;
     }
     yDDotOnScale = Math.sqrt(yDDotOnScale) / h;
@@ -253,16 +227,16 @@ public double initializeStep(
     double h1 =
         (maxInv2 < 1.0e-15)
             ? Math.max(1.0e-6, 0.001 * Math.abs(h))
-            : Math.pow(0.01 / maxInv2, 1.0 / order);
+            : Math.pow(0.01 / maxInv2, 1.0 / context.order());
     h = Math.min(100.0 * Math.abs(h), h1);
-    h = Math.max(h, 1.0e-12 * Math.abs(t0)); // avoids cancellation when computing t1 - t0
+    h = Math.max(h, 1.0e-12 * Math.abs(context.t0()));
     if (h < getMinStep()) {
       h = getMinStep();
     }
     if (h > getMaxStep()) {
       h = getMaxStep();
     }
-    if (!forward) {
+    if (!context.forward()) {
       h = -h;
     }
 
@@ -387,8 +361,7 @@ public double getMaxStep() {
   /**
    * Optional user-supplied initial step size; a negative value signals that automatic estimation
    * should be used instead. The value is mutable between runs and is consumed by {@link
-   * #initializeStep(FirstOrderDifferentialEquations, boolean, int, double[], double, double[],
-   * double[], double[], double[])}.
+   * #initializeStep(StepInitializationContext)}.
    */
   private double initialStep;
 
@@ -401,7 +374,7 @@ public double getMaxStep() {
 
   /**
    * Relative error tolerance applied uniformly when vector tolerances are not provided. It scales
-   * the acceptable error proportionally to the magnitude of each state component during adaptive
+   * the acceptable error proportionally to the scale of each state component during adaptive
    * control.
    */
   protected double scalRelativeTolerance;
@@ -429,7 +402,7 @@ public double getMaxStep() {
 
   /**
    * Aggregates registered switching functions and coordinates event detection. It controls sampling
-   * frequency, root finding and state changes triggered by events encountered during the adaptive
+   * frequency, root finding, and state changes triggered by events encountered during the adaptive
    * integration process.
    */
   protected SwitchingFunctionsHandler switchesHandler;

src/main/java/org/spaceroots/mantissa/ode/DormandPrince54Integrator.java

@@ -71,11 +71,7 @@ public class DormandPrince54Integrator extends RungeKuttaFehlbergIntegrator {
   public DormandPrince54Integrator(
       double minStep, double maxStep, double scalAbsoluteTolerance, double scalRelativeTolerance) {
     super(
-        true,
-        C,
-        A,
-        B,
-        new DormandPrince54StepInterpolator(),
+        new RungeKuttaFehlbergMethod(true, C, A, B, new DormandPrince54StepInterpolator()),
         minStep,
         maxStep,
         scalAbsoluteTolerance,
@@ -104,11 +100,7 @@ public DormandPrince54Integrator(
       double[] vecAbsoluteTolerance,
       double[] vecRelativeTolerance) {
     super(
-        true,
-        C,
-        A,
-        B,
-        new DormandPrince54StepInterpolator(),
+        new RungeKuttaFehlbergMethod(true, C, A, B, new DormandPrince54StepInterpolator()),
         minStep,
         maxStep,
         vecAbsoluteTolerance,
@@ -145,11 +137,11 @@ public int getOrder() {
    * Computes the normalized root-mean-square local error estimate for a tentative step.
    *
    * <p>The implementation combines the Dormand–Prince error coefficients with the stage derivatives
-   * to form the embedded fourth-order estimate, scales each component by the greater magnitude of
-   * the start and end state, and divides by either scalar or per-component tolerances. The returned
+   * to form the embedded fourth-order estimate, scales each component by the greater scale of the
+   * start and end state, and divides by either scalar or per-component tolerances. The returned
    * value is compared against {@code 1.0}; values above one trigger a rejected step and step-size
    * reduction, while values below one permit acceptance. The calculation is deterministic and
-   * agnostic to integration direction.
+   * agnostic to the integration direction.
    *
    * @param yDotK derivatives for all Runge–Kutta stages; index {@code [k][j]} stores stage {@code
    *     k} for component {@code j}.

src/main/java/org/spaceroots/mantissa/ode/DormandPrince853Integrator.java

@@ -9,8 +9,8 @@
  * the normalized local error below one. It automatically estimates an initial step, enforces
  * user-provided minimum and maximum magnitudes, and reuses the final derivative of each accepted
  * step (the <i>first same as last</i> property) to save one evaluation on the next step. Continuous
- * output is available through a step interpolator so callers can sample intermediate states without
- * retriggering derivative evaluations.
+ * output is available through a step interpolator, so callers can sample intermediate states
+ * without retriggering derivative evaluations.
  *
  * <p>Use this class when problems are non-stiff, require tight error control, or benefit from
  * higher-order dense output. Instances are mutable and not thread-safe; confine one integrator to a
@@ -213,19 +213,15 @@ public class DormandPrince853Integrator extends RungeKuttaFehlbergIntegrator {
    *     runs; values smaller than this may be used only for the final partial step.
    * @param maxStep positive upper bound for the absolute step size used during adaptation; must not
    *     be smaller than {@code minStep} and is honored for both forward and backward directions.
-   * @param scalAbsoluteTolerance uniform absolute error threshold applied to each component when
+   * @param scalAbsoluteTolerance a uniform absolute error threshold applied to each component when
    *     normalizing local error estimates; must be non-negative.
-   * @param scalRelativeTolerance uniform relative error threshold scaled by the largest magnitude
-   *     of the start and end state values for each component; must be non-negative.
+   * @param scalRelativeTolerance a uniform relative error threshold scaled by the largest scale of
+   *     the start and end state values for each component; must be non-negative.
    */
   public DormandPrince853Integrator(
       double minStep, double maxStep, double scalAbsoluteTolerance, double scalRelativeTolerance) {
     super(
-        true,
-        C,
-        A,
-        B,
-        new DormandPrince853StepInterpolator(),
+        new RungeKuttaFehlbergMethod(true, C, A, B, new DormandPrince853StepInterpolator()),
         minStep,
         maxStep,
         scalAbsoluteTolerance,
@@ -249,20 +245,15 @@ public DormandPrince853Integrator(
    * @param vecAbsoluteTolerance absolute error thresholds per state component used to scale local
    *     truncation error estimates; array length must match the problem dimension.
    * @param vecRelativeTolerance relative error thresholds per state component multiplied by the
-   *     larger magnitude of the start and end values; array length must match the problem
-   *     dimension.
+   *     larger scale of the start and end values; array length must match the problem dimension.
    */
   public DormandPrince853Integrator(
       double minStep,
       double maxStep,
       double[] vecAbsoluteTolerance,
       double[] vecRelativeTolerance) {
     super(
-        true,
-        C,
-        A,
-        B,
-        new DormandPrince853StepInterpolator(),
+        new RungeKuttaFehlbergMethod(true, C, A, B, new DormandPrince853StepInterpolator()),
         minStep,
         maxStep,
         vecAbsoluteTolerance,
@@ -299,7 +290,7 @@ public int getOrder() {
    *
    * <p>The method combines two embedded error estimators (orders five and three) to produce a
    * single scalar error ratio. Each state component is scaled by either scalar or vector
-   * tolerances, using the larger magnitude of the start and end values to normalize derivatives. A
+   * tolerances, using the larger scale of the start and end values to normalize derivatives. A
    * return value above one signals that the step should be rejected and retried with a smaller
    * size; values at or below one allow the step to be accepted and possibly grown for the next
    * attempt. The computation does not modify the provided arrays.
@@ -310,8 +301,8 @@ public int getOrder() {
    *     dimension and is used for scaling tolerances.
    * @param y1 state estimate at the end of the candidate step produced by the high-order formula;
    *     length matches {@code y0} and is used for scaling tolerances.
-   * @param h signed size of the trial step currently under evaluation; absolute value is used when
-   *     computing the returned ratio.
+   * @param h the signed size of the trial step currently under evaluation; absolute value is used
+   *     when computing the returned ratio.
    * @return normalized error ratio; values greater than one trigger step rejection and reduction,
    *     while values at or below one allow acceptance.
    */