Scale error estimate by dt²/2 for second-order ODEs

The probabilistic error estimate is the standard deviation of the
predicted constraint residual. For second-order ODEs the constraint
operates on u'', so converting to a local error in u requires scaling
by dt²/2 (Taylor series), not just dt as for first-order ODEs.

This aligns with the approach in pnkraemer/probdiffeq#852 and
significantly improves step size selection for second-order problems
like Pleiades.

Author: nathanaelbosch

src/perform_step.jl

@@ -188,7 +188,11 @@ function compute_scaled_error_estimate!(integ, cache)
     @unpack err_tmp = cache
     t = integ.t + integ.dt
     err_est_unscaled = estimate_errors!(cache)
-    err_est_unscaled .*= integ.dt
+    if integ.f isa DynamicalODEFunction # second-order ODE
+        err_est_unscaled .*= integ.dt^2 / 2
+    else
+        err_est_unscaled .*= integ.dt
+    end
     if integ.f isa DynamicalODEFunction # second-order ODE
         DiffEqBase.calculate_residuals!(
             err_tmp,