tutorial: construction of B-series and average vector field example (#68

)

* construct B-series from a given function

* add average vector field method to tutorials

* format

* set version to v0.1.25

* add test

* format

* fix typo

* fix typo

* fix typo

* fix typo

* fix formatting of jldoctest

* fix doctest

* tutorial on creating B-series

* fix formatting of jldoctest

* fix tutorial

* modified equation of AVF

Project.toml

@@ -1,7 +1,7 @@
 name = "BSeries"
 uuid = "ebb8d67c-85b4-416c-b05f-5f409e808f32"
 authors = ["Hendrik Ranocha <[email protected]> and contributors"]
-version = "0.1.24"
+version = "0.1.25"
 
 [deps]
 Latexify = "23fbe1c1-3f47-55db-b15f-69d7ec21a316"
@@ -18,5 +18,5 @@ OrderedCollections = "1"
 Polynomials = "2.0.23, 3"
 Reexport = "1"
 Requires = "1"
-RootedTrees = "2.10.4"
+RootedTrees = "2.12"
 julia = "1.6"

docs/make.jl

@@ -52,6 +52,7 @@ makedocs(modules = [BSeries],
              # "Introduction" => "introduction.md",
              "Tutorials" => [
                  "tutorials/bseries_basics.md",
+                 "tutorials/bseries_creation.md",
                  "tutorials/modified_equations.md",
                  "tutorials/modifying_integrators.md",
                  "tutorials/symbolic_computations.md",

docs/src/tutorials/bseries_basics.md

@@ -5,14 +5,14 @@ methods, generically or when applied to a specific ordinary differential equatio
 
 
 ```@example bseries-basics
-# Load the packages we will use.  
+# Load the packages we will use.
 # These must first be installed using: import Pkg; Pkg.add("package_name")
 using BSeries
 using Latexify  # Only needed for some pretty-printing cells below using `latexify`
 import SymPy; sp=SymPy;
 ```
 
-# B-Series for a generic ODE
+## B-Series for a generic ODE
 
 First we specify the Butcher coefficients of the RK method.
 This can include symbolic expressions and parameterized families of methods.
@@ -53,9 +53,11 @@ We can also print out the B-series coefficients this way:
 coeffs4
 ```
 
-In this form, the rooted trees are printed as level sequences.  The corresponding coefficients are on the right.
+In this form, the rooted trees are printed as level sequences.
+The corresponding coefficients are on the right.
 
-# Exact series and local error
+
+## Exact series and local error
 
 We can also get the B-series of the exact solution:
 
@@ -92,7 +94,8 @@ This confirms again the accuracy of the method, and shows us that we
 can eliminate one of the leading error terms completely if we take
 $\alpha=3/4$ (this is known as Ralston's method, or sometimes as Heun's method).
 
-# B-Series for a specific ODE
+
+## B-Series for a specific ODE
 
 Next, let us define an ODE.  We'll consider the Prothero-Robinson problem:
 
@@ -144,7 +147,8 @@ And their difference, which is the local error:
 expr = sp.simplify(evaluate(ff, u, h, coeffs4) - evaluate(ff, u, h,coeffs_ex))[1]
 ```
 
-# B-series for a generic RK method
+
+## B-series for a generic RK method
 
 We can also examine just the elementary differentials, without specifying a RK method:
 

docs/src/tutorials/bseries_creation.md

@@ -0,0 +1,155 @@
+# [Creating B-series](@id tutorial-bseries-creation)
+
+We have already seen some ways of creating B-series in the
+[basic tutorial](@ref tutorial-bseries-basics). In this tutorial, we look at
+ways to obtain the B-series of different time integration methods.
+
+
+## B-series for Runge-Kutta methods
+
+[BSeries.jl](https://github.com/ranocha/BSeries.jl) and
+[RootedTrees.jl](https://github.com/SciML/RootedTrees.jl) provide
+the type `RungeKuttaMethod` as wrapper of Butcher coefficients `A, b, c` of
+Runge-Kutta methods.
+
+```@example ex:RK4
+using BSeries
+
+# classical RK4
+A = [0 0 0 0
+      1//2 0 0 0
+      0 1//2 0 0
+      0 0 1 0]
+b = [1 // 6, 1 // 3, 1 // 3, 1 // 6]
+rk = RungeKuttaMethod(A, b)
+```
+
+Instead of passing the Butcher coefficients explicitly
+to [`bseries`](@ref), you can also pass the wrapper struct. In fact, this is
+the preferred method.
+
+```@example ex:RK4
+series = bseries(rk, 5)
+```
+
+We can check that the classical Runge-Kutta method is indeed fourth-order accurate.
+
+```@example ex:RK4
+series - ExactSolution(series)
+```
+
+
+## B-series for additive Runge-Kutta methods
+
+[BSeries.jl](https://github.com/ranocha/BSeries.jl) and
+[RootedTrees.jl](https://github.com/SciML/RootedTrees.jl) also support additive
+Runge-Kutta methods via the wrapper `AdditiveRungeKuttaMethod`. For example,
+we can write the Störmer-Verlet method as additive Runge-Kutta method following
+Table II.2.1 of Hairer, Lubich, and Wanner (2002).
+
+```@example ex:SV
+using BSeries
+
+As = [
+    [0 0; 1//2 1//2],
+    [1//2 0; 1//2 0],
+]
+bs = [
+    [1 // 2, 1 // 2],
+    [1 // 2, 1 // 2],
+]
+ark = AdditiveRungeKuttaMethod(As, bs)
+```
+
+We can create the B-series as usual, truncated to order 3.
+
+```@example ex:SV
+series = bseries(ark, 3)
+```
+
+For additive Runge-Kutta methods like this, we use colored rooted trees.
+Again, we can check the order of accuracy by comparing the coefficients to
+the exact solution.
+
+```@example ex:SV
+series - ExactSolution(series)
+```
+
+
+### References
+
+- Ernst Hairer, Gerhard Wanner, Christian Lubich.
+  Geometric numerical integration.
+  Springer, 2002.
+
+
+## [B-series for the average vector field method](@id tutorial-bseries-creation-AVF)
+
+Consider the autonomous ODE
+
+```math
+u'(t) = f\bigl( u(t) \bigr).
+```
+
+The average vector field method
+
+```math
+u^{n+1} = u^{n} + \Delta t \int_0^1 f\bigl(\xi u^{n+1} + (1 - \xi) u^{n}\bigr) \mathrm{d} \xi
+```
+
+introduced by McLachlan, Quispel, and Robidoux (1999) is a second-order accurate
+method. Quispel and McLaren (2008) discovered that it is indeed a B-series method.
+Its coefficients are given explicitly as
+
+```math
+\begin{align*}
+b(.) &= 1, \\
+b([t_1, ..., t_n]) &= b(t_1)...b(t_n) / (n + 1).
+\end{align*}
+```
+
+by Celledoni, McLachlan, McLaren, Owren, Quispel, and Wright (2009). We can
+implement this up to order 5 in BSeries.jl as follows.
+
+```@example ex:AVF
+using BSeries
+
+series = bseries(5) do t, series
+    if order(t) in (0, 1)
+        return 1 // 1
+    else
+        v = 1 // 1
+        n = 0
+        for subtree in SubtreeIterator(t)
+            v *= series[subtree]
+            n += 1
+        end
+        return v / (n + 1)
+    end
+end
+```
+
+We can check that this method is second-order accurate by comparing it to
+the B-series of the exact solution, truncated at the same order.
+
+```@example ex:AVF
+series - ExactSolution(series)
+```
+
+### References
+
+- Robert I. McLachlan, G. Reinout W. Quispel, and Nicolas Robidoux.
+  "Geometric integration using discrete gradients."
+  Philosophical Transactions of the Royal Society of London.
+  Series A: Mathematical, Physical and Engineering Sciences 357,
+  no. 1754 (1999): 1021-1045.
+  [DOI: 10.1098/rsta.1999.0363](https://doi.org/10.1098/rsta.1999.0363)
+- G. Reinout W. Quispel, and David Ian McLaren.
+  "A new class of energy-preserving numerical integration methods."
+  Journal of Physics A: Mathematical and Theoretical 41, no. 4 (2008): 045206.
+  [DOI: 10.1088/1751-8113/41/4/045206](https://doi.org/10.1088/1751-8113/41/4/045206)
+- Elena Celledoni, Robert I. McLachlan, David I. McLaren, Brynjulf Owren,
+  G. Reinout W. Quispel, and William M. Wright.
+  "Energy-preserving Runge-Kutta methods."
+  ESAIM: Mathematical Modelling and Numerical Analysis 43, no. 4 (2009): 645-649.
+  [DOI: 10.1051/m2an/2009020](https://doi.org/10.1051/m2an/2009020)

docs/src/tutorials/modified_equations.md

@@ -72,7 +72,7 @@ using LaTeXStrings, Plots
 
 fig = plot(xguide=L"$q$", yguide=L"$p$")
 default(linewidth=2)
-plot!(fig, sol_ref, vars=(2, 1), label="Reference solution")
+plot!(fig, sol_ref, idxs=(2, 1), label="Reference solution")
 scatter!(fig, last.(sol_euler.u), first.(sol_euler.u),
          label="Explicit Euler, dt = $dt")
 plot!(fig, xlims=(0.0, 9.0), ylims=(0.0, 5.0))
@@ -124,7 +124,7 @@ function solve_modified_equation(ode, truncation_orders, dt)
     modified_f, _ = build_function(series, u_sym, expression=Val(false))
     modified_ode = ODEProblem((u, params, t) -> modified_f(u), ode.u0, ode.tspan)
     modified_sol = solve(modified_ode, Tsit5())
-    plot!(fig, modified_sol, vars=(2, 1),
+    plot!(fig, modified_sol, idxs=(2, 1),
           label="Modified ODE, order $(truncation_order-1)")
   end
   fig
@@ -199,7 +199,7 @@ sol_ref = solve(ode, Tsit5(), abstol=1.0e-9, reltol=1.0e-9)
 
 fig = plot(xguide=L"$q$", yguide=L"$p$")
 default(linewidth=2)
-plot!(fig, sol_ref, vars=(1, 2), label="Reference solution")
+plot!(fig, sol_ref, idxs=(1, 2), label="Reference solution")
 
 savefig(fig, "lotka_volterra_reference.svg"); nothing # hide
 ```
@@ -318,7 +318,7 @@ for truncation_order in 2:3
   modified_f, _ = build_function(series, u_sym, expression=Val(false))
   modified_ode = ODEProblem((u, params, t) -> modified_f(u), ode.u0, ode.tspan)
   modified_sol = solve(modified_ode, Tsit5(), abstol=1.0e-9, reltol=1.0e-9)
-  plot!(fig, modified_sol, vars=(1, 2),
+  plot!(fig, modified_sol, idxs=(1, 2),
         label="Modified equation, order $(truncation_order-1)")
 end
 
@@ -366,7 +366,7 @@ sol_ref = solve(ode, Tsit5(), abstol=1.0e-9, reltol=1.0e-9)
 
 fig = plot(xguide=L"$q$", yguide=L"$v$")
 default(linewidth=2)
-plot!(fig, sol_ref, vars=(2, 1), label="Reference solution")
+plot!(fig, sol_ref, idxs=(2, 1), label="Reference solution")
 scatter!(fig, last.(sol_baseline.u), first.(sol_baseline.u),
          label="Störmer-Verlet, dt = $dt")
 
@@ -436,7 +436,7 @@ for truncation_order in 3:2:5
   modified_f, _ = build_function(series, u_sym, expression=Val(false))
   modified_ode = ODEProblem((u, params, t) -> modified_f(u), [q0, v0], ode.tspan)
   modified_sol = solve(modified_ode, Tsit5(), abstol=1.0e-9, reltol=1.0e-9)
-  plot!(fig, modified_sol, vars=(1, 2),
+  plot!(fig, modified_sol, idxs=(1, 2),
         label="Modified equation, order $(truncation_order-1)")
 end
 fig
@@ -452,6 +452,39 @@ the fourth-order modified equation is even a bit more accurate than the
 second-order one.
 
 
+## Modified equation of the average vector field method
+
+Here, we reproduce Example 1 of [^CelledoniMcLachlanOwrenQuispel2010].
+First, we create the B-series of the average vector field method as described
+in the [tutorial on creating B-series](@ref tutorial-bseries-creation-AVF).
+
+```@example ex:AVF-mod-eq
+using BSeries
+
+series = bseries(5) do t, series
+    if order(t) in (0, 1)
+        return 1 // 1
+    else
+        v = 1 // 1
+        n = 0
+        for subtree in SubtreeIterator(t)
+            v *= series[subtree]
+            n += 1
+        end
+        return v / (n + 1)
+    end
+end
+```
+
+Next, we compute the coefficients of its modified equation.
+
+```@example ex:AVF-mod-eq
+coefficients = modified_equation(series)
+```
+
+Remember that the coefficients of the B-series need to be divided by the
+symmetry of the rooted trees to get the final expressions.
+
 
 ## References
 
@@ -470,3 +503,9 @@ second-order one.
   Algebraic Structures of B-series.
   Foundations of Computational Mathematics
   [DOI: 10.1007/s10208-010-9065-1](https://doi.org/10.1007/s10208-010-9065-1)
+
+[^CelledoniMcLachlanOwrenQuispel2010]:
+  Elena Celledoni, Robert I. McLachlan, Brynjulf Owren, and G. R. W. Quispel (2010)
+  Energy-preserving integrators and the structure of B-series.
+  Foundations of Computational Mathematics
+  [DOI: 10.1007/s10208-010-9073-1](https://doi.org/10.1007/s10208-010-9073-1)

docs/src/tutorials/modifying_integrators.md

@@ -64,7 +64,7 @@ using LaTeXStrings, Plots
 
 fig = plot(xguide=L"$q$", yguide=L"$p$")
 default(linewidth=2)
-plot!(fig, sol_ref, vars=(2, 1), label="Reference solution")
+plot!(fig, sol_ref, idxs=(2, 1), label="Reference solution")
 scatter!(fig, last.(sol_euler.u), first.(sol_euler.u),
          label="Explicit Euler, dt = $dt")
 plot!(fig, xlims=(0.0, 4.0), ylims=(0.0, 2.5))
@@ -109,7 +109,7 @@ for truncation_order in 2:4
   modified_f, _ = build_function(series, u_sym, expression=Val(false))
   modified_ode = ODEProblem((u, params, t) -> modified_f(u), ode.u0, tspan)
   modified_sol_euler = solve(modified_ode, Euler(), dt=dt)
-  plot!(fig, modified_sol_euler, vars=(2, 1),
+  plot!(fig, modified_sol_euler, idxs=(2, 1),
         label="Euler, modified ODE order $(truncation_order-1)")
 end
 plot!(fig, xlims=(0.0, 4.0), ylims=(0.0, 2.5))