3rd and 4th order methods

_typos.toml

@@ -3,5 +3,13 @@
 Strang = "Strang"
 # Citation key for Trotter (Tro:1959:psg) used in docstrings
 Tro = "Tro"
+# David Ketcheson, in the citation key AuzHofKetKoc:2017:psm
+Ket = "Ket"
+# Citation keys for the order barrier: Sheng (She:1989:slp), Suzuki (Suz:1991:gtf)
+She = "She"
+Suz = "Suz"
+# Hansen and Ostermann, in the citation key HanOst:2009:hos
+Han = "Han"
+Ost = "Ost"
 # Splitting variable name
 BA = "BA"

docs/src/api-reference/index.md

@@ -18,6 +18,9 @@ GenericSplitFunction
 LieTrotterGodunov
 StrangMarchuk
 PalindromicPairLieTrotterGodunov
+Ruth3
+Yoshida4
+AdjointPair
 ```
 
 ## Per-node configuration

docs/src/assets/references.bib

@@ -52,3 +52,75 @@ @incollection{Mar:1971:tsm
   year={1971},
   publisher={Academic Press}
 }
+
+@article{Rut:1983:cim,
+  title={A canonical integration technique},
+  author={Ruth, Ronald D.},
+  journal={IEEE Transactions on Nuclear Science},
+  volume={NS-30},
+  number={4},
+  pages={2669--2671},
+  year={1983},
+  doi={10.1109/TNS.1983.4332919},
+  publisher={IEEE}
+}
+
+@article{Yos:1990:cho,
+  title={Construction of higher order symplectic integrators},
+  author={Yoshida, Haruo},
+  journal={Physics Letters A},
+  volume={150},
+  number={5--7},
+  pages={262--268},
+  year={1990},
+  doi={10.1016/0375-9601(90)90092-3},
+  publisher={Elsevier}
+}
+
+@article{AuzHofKetKoc:2017:psm,
+  title={Practical splitting methods for the adaptive integration of nonlinear evolution equations. {P}art {I}: {C}onstruction of optimized schemes and pairs of schemes},
+  author={Auzinger, Winfried and Hofst{\"a}tter, Harald and Ketcheson, David and Koch, Othmar},
+  journal={BIT Numerical Mathematics},
+  volume={57},
+  number={1},
+  pages={55--74},
+  year={2017},
+  doi={10.1007/s10543-016-0626-9},
+  publisher={Springer}
+}
+
+@article{She:1989:slp,
+  title={Solving linear partial differential equations by exponential splitting},
+  author={Sheng, Qin},
+  journal={IMA Journal of Numerical Analysis},
+  volume={9},
+  number={2},
+  pages={199--212},
+  year={1989},
+  doi={10.1093/imanum/9.2.199},
+  publisher={Oxford University Press}
+}
+
+@article{Suz:1991:gtf,
+  title={General theory of fractal path integrals with applications to many-body theories and statistical physics},
+  author={Suzuki, Masuo},
+  journal={Journal of Mathematical Physics},
+  volume={32},
+  number={2},
+  pages={400--407},
+  year={1991},
+  doi={10.1063/1.529425},
+  publisher={AIP}
+}
+
+@article{HanOst:2009:hos,
+  title={High order splitting methods for analytic semigroups exist},
+  author={Hansen, Eskil and Ostermann, Alexander},
+  journal={BIT Numerical Mathematics},
+  volume={49},
+  number={3},
+  pages={527--542},
+  year={2009},
+  doi={10.1007/s10543-009-0236-x},
+  publisher={Springer}
+}

docs/src/devdocs/index.md

@@ -146,8 +146,9 @@ end
 
 This example is written for exactly two operators. A scheme that works for any number
 of them loops over `children` with `Unrolled.@unroll`, as the built-in algorithms in
-`src/solver.jl` do; a plain `for` loop over the heterogeneously typed tuple would be
-type unstable.
+`src/solvers/` do; a plain `for` loop over the heterogeneously typed tuple would be
+type unstable. `advance_one_child!` above is `src/solvers/common.jl`'s
+`_advance_child!`, which the built-in schemes share.
 
 ### Adaptive algorithms
 
@@ -177,9 +178,94 @@ end
 OrdinaryDiffEqCore interface for the estimate; go through them rather than touching
 the `EEst` field, whose location on the integrator is an implementation detail.
 
-See [`PalindromicPairLieTrotterGodunov`](@ref) in `src/solver.jl` for a complete
-example, and [Adaptive time stepping](@ref) for how the two layers of adaptivity
-interact.
+See [`PalindromicPairLieTrotterGodunov`](@ref) in `src/solvers/adjoint_pair.jl` for a
+complete example, and [Adaptive time stepping](@ref) for how the two layers of
+adaptivity interact.
+
+## Coefficient tables
+
+Schemes of order three and above are not written out as passes over the children.
+They are described by a table of coefficients and stepped by one shared traversal,
+which lives in `src/solvers/coefficients.jl`; the schemes themselves
+(`src/solvers/tables.jl`) are then little more than the table plus two interface
+methods. See [Higher order splittings](@ref theory_higher-order) for where the
+tables come from.
+
+```@docs
+OrdinaryDiffEqOperatorSplitting.SplittingCoefficients
+OrdinaryDiffEqOperatorSplitting.coefficients
+OrdinaryDiffEqOperatorSplitting.order
+```
+
+A table-driven scheme is added by giving it a struct, a table, `coefficients`,
+`order`, and an `init_cache` returning the shared
+`SplittingCoefficientsCache` — no stepping code of its own:
+
+```julia
+struct MyThirdOrder{AlgTupleType <: Tuple} <:
+       OrdinaryDiffEqOperatorSplitting.AbstractOperatorSplittingAlgorithm
+    inner_algs::AlgTupleType
+end
+
+const MY_COEFFICIENTS = OrdinaryDiffEqOperatorSplitting.SplittingCoefficients(
+    (7 // 24, 2 // 3), (3 // 4, -2 // 3), (-1 // 24, 1 // 1)
+)
+
+OrdinaryDiffEqOperatorSplitting.coefficients(::MyThirdOrder) = MY_COEFFICIENTS
+OrdinaryDiffEqOperatorSplitting.order(::MyThirdOrder) = 3
+
+function OrdinaryDiffEqOperatorSplitting.init_cache(
+        f::GenericSplitFunction, alg::MyThirdOrder;
+        uprev::AbstractArray, u::AbstractVector,
+)
+    return OrdinaryDiffEqOperatorSplitting.SplittingCoefficientsCache(
+        u, uprev, OrdinaryDiffEqOperatorSplitting.coefficients(alg))
+end
+```
+
+Only the consistency condition is checked when the table is built, so a table that
+constructs successfully can still fail to reach the order it claims; a convergence
+test is the only real check.
+
+Two properties of the traversal are worth knowing when reading or extending it:
+
+- **A zero coefficient is skipped entirely**, synchronization included, since the
+  flow is the identity. `Yoshida4`'s last stage relies on this.
+- **The first flow of a step is always synchronized.** Strang-Marchuk skips that
+  sync because its reverse pass ends on operator 1, leaving that child's buffer
+  current; a general table ends on operator `N`, so operator 1's buffer is stale and
+  the same shortcut would silently corrupt every step after the first.
+
+Because a table with negative coefficients steps some children backward, a
+table-driven scheme also depends on the direction reversal described under
+[Backward sub-steps](@ref devdocs_backward-substeps).
+
+If the scheme's order is odd, wrapping it in `AdjointPair` makes it adaptive for
+free: the adjoint is the same table traversed backwards, so nothing further is
+needed from the scheme.
+
+Any published two- or three-operator table with real coefficients — the `AB`/`ABC`
+tables of [AuzHofKetKoc:2017:psm](@cite) and the collections derived from them — can
+be added this way, needing nothing but the coefficients themselves.
+
+The [complex-coefficient schemes](@ref theory_higher-order) are the exception. The
+traversal would not change for them, but everything downstream of `coefficient * dt`
+would: the sub-problems, their states and their inner integrators would all have to
+be complex, as would the error norms and step size controllers. That is a property of
+the problem being split rather than of the table, so it is not something a table
+alone can opt into.
+
+## [Backward sub-steps](@id devdocs_backward-substeps)
+
+An inner integrator fixes its direction of integration when it is constructed, and
+a scheme with negative coefficients has to step it the other way. Rather than build
+two integrators per child, a sub-step against the child's direction reverses the
+child in place, steps, and reverses it back:
+
+```@docs
+OrdinaryDiffEqOperatorSplitting.reverse_direction!
+OrdinaryDiffEqOperatorSplitting.tstops_and_saveat_heaps
+```
 
 ## Dense output
 

docs/src/topics/adaptivity.md

@@ -6,10 +6,14 @@ Two independent layers of a splitting tree can adapt their step sizes:
   they would outside of this package. Their behavior is governed by the tolerances
   and options they receive.
 - **Splitting nodes** adapt the splitting step itself. This requires an algorithm
-  that produces an error estimate of the splitting error;
-  [`PalindromicPairLieTrotterGodunov`](@ref) is currently the only one. It advances
-  with the average of the two mutually reversed Lie-Trotter sequences and uses half
-  their difference as the local error estimate.
+  that produces an error estimate of the splitting error, which here means running a
+  scheme together with its adjoint: the average of the two is the solution and half
+  their difference is the local error estimate (see
+  [Adjoint pairs](@ref theory_higher-order)).
+  [`PalindromicPairLieTrotterGodunov`](@ref) is that construction for the
+  Lie-Trotter sequence, and [`AdjointPair`](@ref) builds it from any odd-order
+  scheme — `AdjointPair(Ruth3((Tsit5(), Tsit5())))` is a fourth-order adaptive
+  splitting. Both cost two sweeps per step.
 
 By default every node adapts exactly if its own algorithm can: for
 

docs/src/topics/time-integration.md

@@ -165,6 +165,120 @@ $[L_1, L_2] = L_1 L_2 - L_2 L_1$ to cancel, leaving a local truncation error
 of $O(t^3)$ and hence second-order global accuracy. The same argument extends to
 the general $N$-operator palindromic scheme.
 
+## [Higher order splittings](@id theory_higher-order)
+
+Lie-Trotter-Godunov and Strang-Marchuk are both instances of a more general
+construction. A splitting scheme applies the sub-flows in a fixed sequence, each for
+a fixed fraction of the step, so for $N$ operators it is completely described by an
+$S \times N$ table of coefficients $a_{ji}$: stage $j$ advances operator $i$ by
+$a_{ji} \Delta t$. In the linear two-operator case,
+
+```math
+\mathcal{S}(\Delta t) = \prod_{j=1}^{S} e^{a_{jN} \Delta t L_N} \cdots e^{a_{j1} \Delta t L_1} \, .
+```
+
+Lie-Trotter-Godunov is the one-stage table $a = (1, 1)$ and Strang-Marchuk is the
+two-stage table $a = \bigl((\tfrac{1}{2}, 1), (\tfrac{1}{2}, 0)\bigr)$.
+
+Each operator's coefficients must sum to one,
+
+```math
+\sum_{j=1}^{S} a_{ji} = 1 \quad \text{for every } i \, ,
+```
+
+since otherwise the scheme does not even advance every sub-problem by $\Delta t$.
+This is the consistency condition, and it is the only one this package checks when a
+table is constructed. Attaining order $p$ imposes further conditions, one per
+independent commutator up to order $p$, obtained by matching the
+Baker-Campbell-Hausdorff expansion of the product above against that of
+$e^{\Delta t (L_1 + L_2)}$ — the same computation as in the two analyses above,
+carried further.
+
+### Composition: the triple jump
+
+Solving the order conditions directly gets unpleasant quickly. A cheaper route is
+*composition*: build a higher-order scheme out of a symmetric one of lower order.
+If $\mathcal{S}_2$ is any symmetric second-order scheme, then
+
+```math
+\mathcal{S}_4(\Delta t) = \mathcal{S}_2(w_1 \Delta t) \, \mathcal{S}_2(w_0 \Delta t) \, \mathcal{S}_2(w_1 \Delta t)
+```
+
+is symmetric for any weights, and hence of even order. It is of order four as soon
+as the weights satisfy
+
+```math
+2 w_1 + w_0 = 1 \, , \qquad 2 w_1^3 + w_0^3 = 0 \, ,
+```
+
+the first being consistency and the second the cancellation of the third-order term.
+The real solution is $w_1 = 1/(2 - 2^{1/3})$ and $w_0 = -2^{1/3} w_1$, giving
+Yoshida's "triple jump" [Yos:1990:cho](@cite), implemented here as
+[`Yoshida4`](@ref). Writing the three Strang steps out as a flat sequence of flows
+and merging the adjacent flows of the same operator that the composition leaves next
+to each other collapses nine flows to eight — which is exactly the four-stage table
+`Yoshida4` carries, its last stage having a zero second coefficient.
+
+### The order barrier and negative coefficients
+
+Note that $w_0 < 0$ above. This is not an artifact of the construction: no splitting
+scheme with *real* coefficients of order greater than two has all of them positive
+[She:1989:slp,Suz:1991:gtf](@cite). Any third- or higher-order real splitting
+therefore integrates some sub-problem *backward in time* during part of every step,
+which has two practical consequences.
+
+First, the sub-problems must admit a backward flow. For a parabolic sub-problem —
+diffusion, say — the backward evolution is ill-posed and the negative sub-steps are
+violently unstable, so on a reaction-diffusion system the higher-order schemes here
+are not usable on the diffusion operator, however attractive their order. This is
+the reason Strang-Marchuk remains the workhorse despite being only second order.
+
+The barrier is a statement about real coefficients only. Allowing *complex*
+coefficients with positive real part, high-order splittings do exist for analytic
+semigroups [HanOst:2009:hos](@cite): the sub-steps then move along rays into the
+complex time plane rather than backward along the real axis, which keeps a parabolic
+sub-flow well posed.
+
+Second, the implementation has to actually run its sub-integrators backwards. An
+inner integrator fixes its direction of integration at construction, so a negative
+sub-step temporarily reverses it; see the developer documentation for the details.
+
+### Adjoint pairs
+
+The *adjoint* of a scheme is
+
+```math
+\mathcal{S}^*(\Delta t) = \mathcal{S}(-\Delta t)^{-1} \, ,
+```
+
+which for a splitting scheme is simply its whole sequence of flows run in reverse
+order, every coefficient keeping its sign and its operator. A scheme is symmetric
+exactly when $\mathcal{S}^* = \mathcal{S}$, which is why Strang-Marchuk — a
+palindrome — gains an order over Lie-Trotter-Godunov.
+
+If $\mathcal{S}$ has order $p$ with leading local error $C \Delta t^{p+1}$, then
+$\mathcal{S}^*$ has the same order with leading error $(-1)^p C \Delta t^{p+1}$. For
+**odd** $p$ the two signs oppose, so running the pair from the same initial value
+gives, at twice the cost of one scheme,
+
+```math
+\frac{\mathcal{S} + \mathcal{S}^*}{2} \quad \text{of order } p+1 \, ,
+\qquad
+\frac{\mathcal{S} - \mathcal{S}^*}{2} \quad \text{an estimate of the local error of } \mathcal{S} \, ,
+```
+
+the latter being asymptotically correct as $\Delta t \to 0$
+[AuzHofKetKoc:2017:psm](@cite). This is the Milne device applied to a scheme and its
+adjoint, and it is what makes the splitting error itself estimable and hence the
+splitting step adaptive — see [Adaptive time stepping](@ref). The construction is
+[`AdjointPair`](@ref); at $p = 1$, with Lie-Trotter-Godunov as the base, it is the
+pair of mutually reversed sequences implemented directly as
+[`PalindromicPairLieTrotterGodunov`](@ref).
+
+For even $p$ the two leading terms are *equal* rather than opposite: averaging
+cancels nothing and the difference is not an error estimate, which is why
+[`AdjointPair`](@ref) rejects an even-order base scheme.
+
 ## References
 
 ```@bibliography

src/OrdinaryDiffEqOperatorSplitting.jl

@@ -150,11 +150,16 @@ include("function.jl")
 include("config_tree.jl")
 include("problem.jl")
 include("integrator.jl")
-include("solver.jl")
+include("solvers/common.jl")
+include("solvers/coefficients.jl")
+include("solvers/lie_trotter_godunov.jl")
+include("solvers/strang_marchuk.jl")
+include("solvers/tables.jl")
+include("solvers/adjoint_pair.jl")
 include("utils.jl")
 
 export GenericSplitFunction, OperatorSplittingProblem, LieTrotterGodunov, StrangMarchuk,
-    PalindromicPairLieTrotterGodunov
+    PalindromicPairLieTrotterGodunov, Ruth3, Yoshida4, AdjointPair
 export SplitNode, TreeOption
 
 include("precompilation.jl")