190 lines
7.5 KiB
Python
190 lines
7.5 KiB
Python
"""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()
|