|
| 1 | +""" |
| 2 | + GlobalDefectCorrection(alg; gtol=nothing, correction_alg=alg, maxiters=6, |
| 3 | + safety=0.8, companion_abstol, companion_reltol) |
| 4 | +
|
| 5 | +Wrap an ODE algorithm with global error estimation and control based on |
| 6 | +solving for the correction (defect correction). After a dense forward solve |
| 7 | +producing the interpolant `P(t)`, the global error `ε(t) ≈ y(t) - P(t)` is |
| 8 | +estimated by integrating the nonlinear companion ODE |
| 9 | +
|
| 10 | +```math |
| 11 | +ε'(t) = f(P(t) + ε(t), p, t) - P'(t), \\qquad ε(t_0) = 0, |
| 12 | +``` |
| 13 | +
|
| 14 | +which requires no Jacobian, only extra evaluations of `f` on perturbed |
| 15 | +arguments. This is the "solving for the correction" technique of Zadunaisky |
| 16 | +(1976) and Dormand, Duckers and Prince (1984), in the continuous per-step |
| 17 | +dense-output form of Dormand, Lockyer, McGorrigan and Prince (1989); its |
| 18 | +validity under standard local-error step control was proven by Calvo, Higham, |
| 19 | +Montijano and Rández (1996). |
| 20 | +
|
| 21 | +Set the requested absolute endpoint error with `gtol`: |
| 22 | +
|
| 23 | +```julia |
| 24 | +solve(prob, GlobalDefectCorrection(Tsit5(); gtol = 1.0e-6)) |
| 25 | +``` |
| 26 | +
|
| 27 | +The solver then tightens local tolerances until the estimated endpoint global |
| 28 | +error 2-norm is at most `gtol`, like [`GlobalAdjoint`](@ref). Omit `gtol` when |
| 29 | +using the algorithm only with [`global_error_estimate`](@ref). |
| 30 | +`correction_alg` selects the solver for the companion equation, and |
| 31 | +`companion_abstol` / `companion_reltol` control its tolerances. |
| 32 | +
|
| 33 | +This implementation supports forward-time, standard-mass-matrix ODEs with real |
| 34 | +vector states and no callbacks. The wrapped solver must provide a dense |
| 35 | +first-derivative interpolation. |
| 36 | +
|
| 37 | +## References |
| 38 | +
|
| 39 | + - P. E. Zadunaisky, On the estimation of errors propagated in the numerical |
| 40 | + integration of ordinary differential equations, Numerische Mathematik 27 |
| 41 | + (1976). |
| 42 | + - J. R. Dormand, R. R. Duckers and P. J. Prince, Global error estimation with |
| 43 | + Runge-Kutta methods, IMA Journal of Numerical Analysis 4 (1984). |
| 44 | + - J. R. Dormand, M. A. Lockyer, N. E. McGorrigan and P. J. Prince, Global |
| 45 | + error estimation with Runge-Kutta triples, Computers & Mathematics with |
| 46 | + Applications 18 (1989). |
| 47 | +""" |
| 48 | +struct GlobalDefectCorrection{A, CA, G, V} <: GlobalDiffEqAlgorithm |
| 49 | + alg::A |
| 50 | + correction_alg::CA |
| 51 | + gtol::G |
| 52 | + options::V |
| 53 | +end |
| 54 | + |
| 55 | +function GlobalDefectCorrection( |
| 56 | + alg; |
| 57 | + correction_alg = alg, |
| 58 | + gtol = nothing, |
| 59 | + maxiters = 6, |
| 60 | + safety = 0.8, |
| 61 | + companion_abstol = gtol === nothing ? 1.0e-10 : min(gtol / 100, 1.0e-10), |
| 62 | + companion_reltol = gtol === nothing ? 1.0e-8 : min(gtol / 100, 1.0e-8) |
| 63 | + ) |
| 64 | + options = _companion_options( |
| 65 | + gtol, maxiters, safety, companion_abstol, companion_reltol |
| 66 | + ) |
| 67 | + return GlobalDefectCorrection(alg, correction_alg, gtol, options) |
| 68 | +end |
| 69 | + |
| 70 | +SciMLBase.allows_arbitrary_number_types(::GlobalDefectCorrection) = false |
| 71 | +SciMLBase.allowscomplex(::GlobalDefectCorrection) = false |
| 72 | +SciMLBase.isautodifferentiable(::GlobalDefectCorrection) = false |
| 73 | + |
| 74 | +function _correction_rhs(sol) |
| 75 | + prob = sol.prob |
| 76 | + foop = _oop_rhs(prob) |
| 77 | + return let sol = sol, foop = foop, p = prob.p |
| 78 | + function (ε, _, t) |
| 79 | + u = sol(t, continuity = :right) |
| 80 | + du = sol(t, Val{1}, continuity = :right) |
| 81 | + return foop(u + ε, p, t) - du |
| 82 | + end |
| 83 | + end |
| 84 | +end |
| 85 | + |
| 86 | +""" |
| 87 | + global_error_estimate(prob, alg::GlobalDefectCorrection; abstol=1e-6, reltol=1e-3, kwargs...) |
| 88 | +
|
| 89 | +Solve an ODE problem and estimate the 2-norm of its global error at the final |
| 90 | +time by solving for the correction of the dense numerical solution. `abstol` |
| 91 | +and `reltol` control the forward solve; `companion_abstol` and |
| 92 | +`companion_reltol` control the companion error solve. |
| 93 | +""" |
| 94 | +function global_error_estimate( |
| 95 | + prob::SciMLBase.AbstractODEProblem, alg::GlobalDefectCorrection, args...; |
| 96 | + abstol = 1.0e-6, |
| 97 | + reltol = 1.0e-3, |
| 98 | + companion_abstol = alg.options.companion_abstol, |
| 99 | + companion_reltol = alg.options.companion_reltol, |
| 100 | + kwargs... |
| 101 | + ) |
| 102 | + return _companion_error_estimate( |
| 103 | + _correction_rhs, "GlobalDefectCorrection", |
| 104 | + prob, alg.alg, alg.correction_alg, args...; |
| 105 | + abstol, reltol, companion_abstol, companion_reltol, kwargs... |
| 106 | + ) |
| 107 | +end |
| 108 | + |
| 109 | +function SciMLBase.__solve( |
| 110 | + prob::SciMLBase.AbstractODEProblem, alg::GlobalDefectCorrection, args...; |
| 111 | + abstol = something(alg.gtol, 1.0e-6), |
| 112 | + reltol = something(alg.gtol, 1.0e-3), |
| 113 | + kwargs... |
| 114 | + ) |
| 115 | + alg.gtol === nothing && throw( |
| 116 | + ArgumentError("GlobalDefectCorrection requires a positive `gtol` constructor keyword") |
| 117 | + ) |
| 118 | + _validate_tolerances(abstol, reltol, "local") |
| 119 | + haskey(kwargs, :callback) && |
| 120 | + throw(ArgumentError("GlobalDefectCorrection does not currently support callbacks")) |
| 121 | + estimator = (local_abstol, local_reltol) -> global_error_estimate( |
| 122 | + prob, alg, args...; |
| 123 | + abstol = local_abstol, reltol = local_reltol, kwargs... |
| 124 | + ) |
| 125 | + return _refine_to_gtol( |
| 126 | + estimator, prob, alg.alg, alg.gtol, alg.options, args...; |
| 127 | + abstol, reltol, kwargs... |
| 128 | + ) |
| 129 | +end |
0 commit comments