"""Stage 15 - the credit rule in a second domain: ODE step-size control. The root-finding work ended with one transferable object, the credit safeguard, and one honest conclusion: as a one-dimensional root finder the method is subsumed by Chandrupatla 1997. Adaptive step-size control is the natural place to test whether the transferable object transfers, because it has the same shape - propose a bold step, look at the outcome, adjust - and because every practical controller in it is made of constants somebody chose. classical scipy's own: safety 0.9, clamps 0.2 and 10 3 constants PI the same plus two gains 5 constants elementary what the error model alone dictates 0 constants measured the safety factor MEASURED instead of chosen 0 constants Everything except the controller comes from SciPy: the RK45 tableau, the stage combination and the embedded error estimate, via scipy.integrate._ivp.rk. The driver is checked against solve_ivp(RK45) with the same control law before any comparison is made. A difference between two rows can only come from the controller. Cost model: evaluations of the right-hand side, as evaluations of f were in the root-finding work. A run that does not reach the end of the interval is a failure and is reported as one. Accuracy is checked against a reference so that a controller cannot win by being fast and wrong. """ from __future__ import annotations import math from odestep.controllers import ALL from odestep.driver import integrate, validate_against_scipy TOLS = (1e-4, 1e-6, 1e-8, 1e-10) # -------------------------------------------------------------------------- # problems: classics with known or high-accuracy reference endpoints # -------------------------------------------------------------------------- def p_decay(): f = lambda t, y: [-y[0]] return "exponential decay", f, 0.0, 10.0, [1.0], [math.exp(-10.0)] def p_forced(): f = lambda t, y: [-y[0] + math.sin(t)] ex = 0.5 * (math.sin(10.0) - math.cos(10.0) + math.exp(-10.0)) return "y' = -y + sin t", f, 0.0, 10.0, [0.0], [ex] def p_lotka(): def f(t, y): return [1.5 * y[0] - y[0] * y[1], -3.0 * y[1] + y[0] * y[1]] return "Lotka-Volterra", f, 0.0, 15.0, [10.0, 5.0], None def p_vanderpol(mu=5.0): def f(t, y): return [y[1], mu * (1 - y[0] * y[0]) * y[1] - y[0]] return f"van der Pol mu={mu:g}", f, 0.0, 25.0, [2.0, 0.0], None def p_brusselator(): def f(t, y): return [1.0 + y[0] * y[0] * y[1] - 4.0 * y[0], 3.0 * y[0] - y[0] * y[0] * y[1]] return "Brusselator", f, 0.0, 20.0, [1.5, 3.0], None def p_kepler(ecc): """The two-body problem at eccentricity e - DETEST's D-class, the classic hard non-stiff test: the solution has a sharp passage at perihelion.""" def f(t, y): x1, x2, v1, v2 = y r3 = (x1 * x1 + x2 * x2) ** 1.5 return [v1, v2, -x1 / r3, -x2 / r3] y0 = [1.0 - ecc, 0.0, 0.0, math.sqrt((1 + ecc) / (1 - ecc))] return f"Kepler e={ecc:g}", f, 0.0, 20.0, y0, None def p_arenstorf(): """The restricted three-body orbit, the standard showcase for adaptive step control: the step must vary by orders of magnitude along the path.""" mu = 0.012277471 mp = 1.0 - mu def f(t, y): y1, y2, y3, y4 = y d1 = ((y1 + mu) ** 2 + y2 * y2) ** 1.5 d2 = ((y1 - mp) ** 2 + y2 * y2) ** 1.5 return [y3, y4, y1 + 2 * y4 - mp * (y1 + mu) / d1 - mu * (y1 - mp) / d2, y2 - 2 * y3 - mp * y2 / d1 - mu * y2 / d2] y0 = [0.994, 0.0, 0.0, -2.00158510637908252240537862224] return "Arenstorf orbit", f, 0.0, 17.0652165601579625588917206249, y0, y0 PROBLEMS = [p_decay(), p_forced(), p_lotka(), p_vanderpol(), p_brusselator(), p_kepler(0.5), p_kepler(0.9), p_arenstorf()] def reference(f, t0, t1, y0): """A reference endpoint from the same integrator at a far tighter tolerance - not from a different code, so that a disagreement cannot be blamed on a different method.""" from odestep.controllers import Classical r = integrate(f, t0, t1, y0, Classical(), atol=1e-13, rtol=1e-13, max_steps=2000000) return r["y"] if r["reached_end"] else None def scaled_err(y, ref, atol, rtol): """The achieved error in the SAME mixed norm the integrator controls. A purely relative criterion is wrong here and was: the decay problem ends at exp(-10) ~ 4.5e-5, so at atol = 1e-4 the tolerance permits a hundred per cent relative error, and a correct run was being marked as a failure. What a run must be judged against is the accuracy it was asked for.""" if ref is None: return float("nan") s = 0.0 for a, b in zip(y, ref): sc = atol + rtol * abs(b) s += ((a - b) / sc) ** 2 return math.sqrt(s / len(ref)) def main(): print("instrument check: the driver against scipy solve_ivp(RK45), " "same control law") worst = 0.0 for name, mine, ref, d in validate_against_scipy(): rel = abs(mine - ref) / ref worst = max(worst, rel) print(f" {name:10s} ours {mine:>6d} solve_ivp {ref:>6d} " f"{100 * (mine - ref) / ref:+6.1f}% endpoint agrees to {d:.1e}") if worst > 0.05: print(" driver does NOT reproduce the library - stopping") return print(" the tableau, the stages and the error estimate are scipy's;") print(" only the controller is ours\n") refs = {} for name, f, t0, t1, y0, exact in PROBLEMS: refs[name] = exact if exact is not None else reference(f, t0, t1, y0) grand = {c.name: [0, 0, 0] for c in ALL} # nfev, rejected, failures for tol in TOLS: print("=" * 96) print(f"tolerance {tol:g} (atol = rtol = tol)") print("=" * 96) print(f"{'problem':22s}" + "".join(f"{c.name.split(' ')[0]:>15s}" for c in ALL)) print("-" * 96) for name, f, t0, t1, y0, _exact in PROBLEMS: cells = [] for c in ALL: r = integrate(f, t0, t1, y0, c, atol=tol, rtol=tol) err = scaled_err(r["y"], refs[name], tol, tol) # local error control does not bound the global error, so the # bar is generous: a run fails only if it misses the accuracy # it asked for by more than three orders of magnitude, or does # not reach the end at all bad = (not r["reached_end"]) or ( math.isfinite(err) and err > 1000.0) grand[c.name][0] += r["nfev"] grand[c.name][1] += r["rejected"] grand[c.name][2] += bad cells.append(f"{r['nfev']:>9d}{'!' if bad else ' '}" f"{r['rejected']:>4d}") print(f"{name:22s}" + "".join(f"{c:>15s}" for c in cells)) print() print("=" * 96) print("TOTALS over every problem and tolerance " "(nfev, rejected steps, failures)") print("=" * 96) print(f"{'controller':34s}{'constants':>11s}{'nfev':>10s}" f"{'rejected':>11s}{'failures':>11s}") print("-" * 96) for c in ALL: n, rej, bad = grand[c.name] print(f"{c.name:34s}{c.constants:>11d}{n:>10d}{rej:>11d}{bad:>11d}") print("\n'!' marks a run that did not reach the end or missed the " "reference by more\nthan a thousand tolerances - a controller " "cannot win by being fast and wrong.") if __name__ == "__main__": main()