Fix Hermite fallback for ExplicitRK without B_interp

When B_interp is nothing (e.g. default ExplicitRK with Dormand-Prince),
fall back to Hermite interpolation instead of throwing. Also skip the
custom _ode_addsteps! in this case so k contains only the 2 entries
Hermite expects. For the out-of-place interpolant, compute bi weights
inline to support ForwardDiff.Dual Θ values.

Co-Authored-By: Chris Rackauckas <[email protected]>
Co-Authored-By: Claude Opus 4.5 <[email protected]>

Author: claude

lib/OrdinaryDiffEqExplicitRK/src/explicit_rk_perform_step.jl

@@ -330,25 +330,25 @@ function generic_rk_interpolant(Θ, dt, y₀, k, B_interp, bi; idxs = nothing, o
 
     inv_dt_factor = order <= 1 ? one(dt) : inv(dt)^(order - 1)
 
-    # Fill pre-allocated bi buffer with polynomial weights
-    for i in 1:nstages
-        bi[i] = eval_poly_derivative(Θ, @view(B_interp[i, :]), order)
-    end
-
+    # Compute weights inline (Θ may be a ForwardDiff.Dual, so we cannot
+    # store into the pre-allocated Float64 buffer here).
+    b1 = eval_poly_derivative(Θ, @view(B_interp[1, :]), order)
     return if isnothing(idxs)
-        interp_sum = k[1] * bi[1]
+        interp_sum = k[1] * b1
         for i in 2:nstages
-            interp_sum = interp_sum + k[i] * bi[i]
+            bval = eval_poly_derivative(Θ, @view(B_interp[i, :]), order)
+            interp_sum = interp_sum + k[i] * bval
         end
         if order == 0
             y₀ + dt * interp_sum
         else
             interp_sum * inv_dt_factor
         end
     else
-        interp_sum = k[1][idxs] * bi[1]
+        interp_sum = k[1][idxs] * b1
         for i in 2:nstages
-            interp_sum = interp_sum + k[i][idxs] * bi[i]
+            bval = eval_poly_derivative(Θ, @view(B_interp[i, :]), order)
+            interp_sum = interp_sum + k[i][idxs] * bval
         end
         if order == 0
             y₀[idxs] + dt * interp_sum

lib/OrdinaryDiffEqExplicitRK/src/interpolants.jl

@@ -16,6 +16,10 @@ const ExplicitRKCacheTypes = Union{ExplicitRKCache, ExplicitRKConstantCache}
         always_calc_begin = false, allow_calc_end = true,
         force_calc_end = false
     )
+    # Only compute all stage vectors when B_interp is available for generic RK
+    # interpolation. Without B_interp, the default Hermite addsteps suffices.
+    isnothing(cache.B_interp) && return nothing
+
     (; A, c, stages) = cache
 
     if length(k) < stages || always_calc_begin
@@ -37,6 +41,10 @@ end
         always_calc_begin = false, allow_calc_end = true,
         force_calc_end = false
     )
+    # Only compute all stage vectors when B_interp is available for generic RK
+    # interpolation. Without B_interp, the default Hermite addsteps suffices.
+    isnothing(cache.tab.B_interp) && return nothing
+
     (; kk, tmp, tab) = cache
     (; A, c, stages) = tab
 
@@ -94,6 +102,7 @@ end
 end
 
 # Generate interpolant methods for derivative orders 0-3
+# Only dispatch when B_interp is available; otherwise fall through to Hermite default.
 for order in 0:3
     @eval begin
         @muladd function _ode_interpolant(
@@ -102,6 +111,17 @@ for order in 0:3
                 idxs, T::Type{Val{$order}}, differential_vars
             )
             B_interp = get_B_interp(cache)
+            if isnothing(B_interp)
+                # Fall back to Hermite interpolation (the generic default)
+                dv = OrdinaryDiffEqCore.interpolation_differential_vars(
+                    differential_vars, y₀, idxs
+                )
+                return OrdinaryDiffEqCore.hermite_interpolant(
+                    Θ, dt, y₀, y₁, k,
+                    Val{cache isa OrdinaryDiffEqMutableCache},
+                    idxs, T, dv
+                )
+            end
             bi = get_bi(cache)
             return generic_rk_interpolant(Θ, dt, y₀, k, B_interp, bi; idxs = idxs, order = $order)
         end
@@ -112,6 +132,16 @@ for order in 0:3
                 idxs, T::Type{Val{$order}}, differential_vars
             )
             B_interp = get_B_interp(cache)
+            if isnothing(B_interp)
+                # Fall back to Hermite interpolation (the generic default)
+                dv = OrdinaryDiffEqCore.interpolation_differential_vars(
+                    differential_vars, y₀, idxs
+                )
+                return OrdinaryDiffEqCore.hermite_interpolant!(
+                    out, Θ, dt, y₀, y₁, k,
+                    idxs, T, dv
+                )
+            end
             bi = get_bi(cache)
             return generic_rk_interpolant!(out, Θ, dt, y₀, k, B_interp, bi; idxs = idxs, order = $order)
         end