Sk/update basic examples and docstrings

docs/Project.toml

@@ -1,9 +1,14 @@
 [deps]
+BenchmarkTools = "6e4b80f9-dd63-53aa-95a3-0cdb28fa8baf"
 Documenter = "e30172f5-a6a5-5a46-863b-614d45cd2de4"
 OrdinaryDiffEq = "1dea7af3-3e70-54e6-95c3-0bf5283fa5ed"
 Plots = "91a5bcdd-55d7-5caf-9e0b-520d859cae80"
+PositiveIntegrators = "d1b20bf0-b083-4985-a874-dc5121669aa5"
+SparseArrays = "2f01184e-e22b-5df5-ae63-d93ebab69eaf"
 
 [compat]
+BenchmarkTools = "1"
 Documenter = "1"
-OrdinaryDiffEq = "6"
+OrdinaryDiffEq = "6.59"
 Plots = "1"
+SparseArrays = "1.7"

docs/make.jl

@@ -73,6 +73,9 @@ makedocs(modules = [PositiveIntegrators],
          # Explicitly specify documentation structure
          pages = [
              "Home" => "index.md",
+             "Tutorials" => [
+                 "Linear Advection" => "linear_advection.md",
+             ],
              "API reference" => "api_reference.md",
              "Contributing" => "contributing.md",
              "Code of conduct" => "code_of_conduct.md",

docs/src/index.md

@@ -78,8 +78,8 @@ To solve this PDS together with initial values ``u_1(0)=u_2(0)=2`` on the time d
 ```@example LotkaVolterra
 using PositiveIntegrators # load PDSProblem
 
-P(u, p, t) = [2*u[1]  0.0; u[1]*u[2]  0.0] # Production matrix
-d(u, p, t) = [0.0; u[2]] # Destruction vector
+P(u, p, t) = [2*u[1]  0; u[1]*u[2]  0] # Production matrix
+d(u, p, t) = [0; u[2]] # Destruction vector
 
 u0 = [2.0; 2.0] # initial values
 tspan = (0.0, 10.0) # time span
@@ -88,12 +88,10 @@ tspan = (0.0, 10.0) # time span
 prob = PDSProblem(P, d, u0, tspan)
 nothing #hide
 ```
-Now that the problem has been created, we can solve it with any of the methods of [OrdinaryDiffEq.jl](https://docs.sciml.ai/OrdinaryDiffEq/stable/). Here we use the method `Tsit5()`. Please note that [PositiveIntegrators.jl](https://github.com/SKopecz/PositiveIntegrators.jl) currently only provides methods for positive and conservative PDS, see below.
+Now that the problem has been created, we can solve it with any method of [PositiveIntegrators.jl](https://github.com/SKopecz/PositiveIntegrators.jl). In the following, we use the method `MPRK22(1.0)`. In addition, we could also use any method provided by [OrdinaryDiffEq.jl](https://docs.sciml.ai/OrdinaryDiffEq/stable/), but these might possibly generate negative approximations.
 
 ```@example LotkaVolterra
-using OrdinaryDiffEq  #load Tsit5
-
-sol = solve(prob, Tsit5())
+sol = solve(prob, MPRK22(1.0))
 nothing # hide
 ```
 Finally, we can use [Plots.jl](https://docs.juliaplots.org/stable/) to visualize the solution.
@@ -119,7 +117,6 @@ This shows that the sum of the state variables of a conservative PDS remains con
 \sum_{i=1}^N y_i(t) = \sum_{i=1}^N y_i(0)
 ```
 for all times ``t>0``.
-Moreover, a conservative PDS is completely defined by the square matrix ``\mathbf P=(p_{ij})_{i,j=1,\dots,N}``. There is no need to store an additional vector of destruction terms since ``d_{ij} = p_{ji}`` for all ``i,j=1,\dots,N``.
 
 One specific example of a conservative PDS is the SIR model
 ```math
@@ -165,7 +162,7 @@ tspan = (0.0, 100.0); # time span
 prob = ConservativePDSProblem(P, u0, tspan)
 nothing # hide
 ```
-Since the SIR model is not only conservative but also positive, we can use any MPRK scheme from [PositiveIntegrators.jl](https://github.com/SKopecz/PositiveIntegrators.jl) to solve it. Here we use `MPRK22(1.0)`.
+Since the SIR model is not only conservative but also positive, we can use any scheme from [PositiveIntegrators.jl](https://github.com/SKopecz/PositiveIntegrators.jl) to solve it. Here we use `MPRK22(1.0)`.
 Please note that any method from [OrdinaryDiffEq.jl](https://docs.sciml.ai/OrdinaryDiffEq/stable/) can be used as well, but might possibly generate negative approximations.
 
 ```@example SIR
@@ -176,9 +173,10 @@ Finally, we can use [Plots.jl](https://docs.juliaplots.org/stable/) to visualize
 ```@example SIR
 using Plots
 
-plot(sol, legend=:right)
+plot(sol, label = ["S" "I" "R"], legend=:right)
+plot!(sol, idxs = ((t, S, I, R) -> (t, S + I + R), 0, 1, 2, 3), label = "S+I+R") #Plot S+I+R over time.
 ```
-
+We see that there is always a nonnegative number of people in each compartment, while the population ``S+I+R`` remains constant over time.
 
 ## Referencing
 

docs/src/linear_advection.md

@@ -0,0 +1,121 @@
+# Tutorial: Solution of the linear advection equation
+
+This tutorial is about the efficient solution of production-destruction systems (PDS) with a large number of differential equations. 
+We will explore several ways to represent such large systems and assess their efficiency. 
+
+## Definition of the production-destruction system
+
+One example of the occurrence of a PDS with a large number of equations is the space discretization of a partial differential equation. In this tutorial we want to solve the linear advection equation
+
+```math
+\partial_t u(t,x)=-a\partial_x u(t,x),\quad u(0,x)=u_0(x)
+```
+
+with ``a>0``, ``t≥ 0``, ``x\in[0,1]`` and periodic boundary conditions. To keep things as simple as possible, we 
+discretize the space domain as ``0=x_0<x_1\dots <x_{N-1}<x_N=1`` with ``x_i = i Δ x`` for ``i=0,\dots,N`` and ``Δx=1/N``. An upwind discretization of the spatial derivative yields the ODE system
+
+```math
+\begin{aligned}
+&\partial_t u_1(t) =-\frac{a}{Δx}\bigl(u_1(t)-u_{N}(t)\bigr),\\
+&\partial_t u_i(t) =-\frac{a}{Δx}\bigl(u_i(t)-u_{i-1}(t)\bigr),\quad i=2,\dots,N,
+\end{aligned}
+```
+
+where ``u_i(t)`` is an approximation of ``u(t,x_i)`` for ``i=1,\dots, N``.
+This system can also be written as ``\partial_t \mathbf u(t)=\mathbf A\mathbf u(t)`` with ``\mathbf u(t)=(u_1(t),\dots,u_N(t))`` and 
+
+```math
+\mathbf A= \frac{a}{Δ x}\begin{bmatrix}-1&0&\dots&0&1\\1&-1&\ddots&&0\\0&\ddots&\ddots&\ddots&\vdots\\ \vdots&\ddots&\ddots&\ddots&0\\0&\dots&0&1&-1\end{bmatrix}.
+```
+
+In particular the matrix ``\mathbf A`` shows, that there is a single production term and a single destruction term per equation. 
+Furthermore, the system is conservative as ``\mathbf A`` has column sum zero.
+To be precise, the production matrix ``\mathbf P = (p_{i,j})`` of this conservative PDS is given by
+
+```math
+\begin{aligned}
+&p_{1,N}(t,\mathbf u(t)) = \frac{a}{Δ x}u_N(t),\\
+&p_{i,i-1}(t,\mathbf u(t)) = \frac{a}{Δ x}u_{i-1}(t),\quad i=2,\dots,N.
+\end{aligned}
+```
+
+Since the PDS is conservative, we have ``d_{i,j}=p_{j,i}`` and the system is fully determined by the production matrix ``\mathbf P``.
+
+## Solution of the production-destruction system
+
+Now we are ready to define a `ConservativePDSProblem` and to solve this problem with a method of [PositiveIntegrators.jl](https://github.com/SKopecz/PositiveIntegrators.jl) or [OrdinaryDiffEq.jl](https://docs.sciml.ai/OrdinaryDiffEq/stable/). In the following we use ``a=1``, ``N=1000`` and the time domain ``t\in[0,1]``. Moreover, we choose the step function
+
+```math
+u_0(x)=\begin{cases}1, & 0.4 ≤ x ≤ 0.6,\\ 0,& \text{elsewhere}\end{cases}
+```
+
+as initial condition. Due to the periodic boundary conditions and the transport velocity ``a=1``, the solution at time ``t=1`` is identical to the initial distribution, i.e. ``u(1,x) = u_0(x)``.
+
+```@example LinearAdvection
+N = 1000 # number of subintervals
+dx = 1/N # mesh width
+x = LinRange(dx, 1.0, N) # discretization points x_1,...,x_N = x_0
+u0 = @. 0.0 + (0.4 ≤ x ≤ 0.6) * 1.0 # initial solution
+tspan = (0.0, 1.0) # time domain
+
+nothing #hide
+```
+
+As mentioned above, we will try different approaches to solve this PDS and compare their efficiency. These are
+1. an in-place implementation with a dense matrix,
+2. an in-place implementation with a sparse matrix.
+
+### Standard in-place implementation
+
+```@example LinearAdvection
+using PositiveIntegrators # load ConservativePDSProblem
+
+function lin_adv_P!(P, u, p, t)
+    P .= 0.0
+    N = length(u)
+    dx = 1 / N
+    P[1, N] = u[N] / dx
+    for i in 2:N
+        P[i, i - 1] = u[i - 1] / dx
+    end
+    return nothing
+end
+
+prob = ConservativePDSProblem(lin_adv_P!, u0, tspan) # create the PDS
+
+sol = solve(prob, MPRK43I(1.0, 0.5); save_everystep = false)
+
+nothing #hide
+```
+
+```@example LinearAdvection
+using Plots
+
+plot(x,u0)
+plot!(x, last(sol.u))
+```
+
+```@example LinearAdvection
+using SparseArrays
+p_prototype = spdiagm(-1 => ones(eltype(u0), N - 1),
+                      N - 1 => ones(eltype(u0), 1))
+prob_sparse = ConservativePDSProblem(lin_adv_P!, u0, tspan; p_prototype=p_prototype)
+
+sol_sparse = solve(prob_sparse, MPRK43I(1.0, 0.5); save_everystep = false)
+
+nothing #hide
+```
+
+```@example LinearAdvection
+plot(x,u0)
+plot!(x, last(sol_sparse.u))
+```
+
+```@example LinearAdvection
+using BenchmarkTools
+@benchmark solve(prob, MPRK43I(1.0, 0.5); save_everystep = false)
+```
+
+```@example LinearAdvection
+@benchmark solve(prob_sparse, MPRK43I(1.0, 0.5); save_everystep = false)
+```

src/mprk.jl

@@ -122,13 +122,14 @@ The first-order modified Patankar-Euler algorithm for production-destruction sys
 first-order accurate, unconditionally positivity-preserving, and
 linearly implicit.
 
-The scheme was introduced by Burchard et al for conservative production-destruction systems.
+The scheme was introduced by Burchard et al. for conservative production-destruction systems.
 For nonconservative production–destruction systems we use the straight forward extension
 
-``u_i^{n+1} = u_i^n + Δt \\sum_{j, j≠i} \\biggl(p_{ij}^n \\frac{u_j^{n+1}}{u_j^n}-d_{ij}^n \\frac{u_i^{n+1}}{u_i^n}\\biggr) + {\\Delta}t p_{ii}^n - Δt d_{ii}^n\\frac{u_i^{n+1}}{u_i^n}``.
+``u_i^{n+1} = u_i^n + Δt \\sum_{j, j≠i} \\biggl(p_{ij}^n \\frac{u_j^{n+1}}{u_j^n}-d_{ij}^n \\frac{u_i^{n+1}}{u_i^n}\\biggr) + {\\Delta}t p_{ii}^n - Δt d_{ii}^n\\frac{u_i^{n+1}}{u_i^n}``,
 
-The modified Patankar-Euler method requires the special structure of a
-[`PDSProblem`](@ref) or a [`ConservativePDSProblem`](@ref).
+where ``p_{ij}^n = p_{ij}(t^n,\\mathbf u^n)`` and ``d_{ij}^n = d_{ij}(t^n,\\mathbf u^n)``.
+
+The modified Patankar-Euler method requires the special structure of a [`PDSProblem`](@ref) or a [`ConservativePDSProblem`](@ref).
 
 You can optionally choose the linear solver to be used by passing an
 algorithm from [LinearSolve.jl](https://github.com/SciML/LinearSolve.jl)
@@ -137,6 +138,8 @@ You can also choose the parameter `small_constant` which is added to all Patanka
 to avoid divisions by zero. You can pass a value explicitly, otherwise `small_constant` is set to
 `floatmin` of the floating point type used.
 
+The current implementation only supports fixed time steps. 
+
 ## References
 
 - Hans Burchard, Eric Deleersnijder, and Andreas Meister.
@@ -329,14 +332,18 @@ end
     MPRK22(α; [linsolve = ..., small_constant = ...])
 
 A family of second-order modified Patankar-Runge-Kutta algorithms for
-production-destruction systems. Each member of this family is an one-step, two-stage method which is
+production-destruction systems. Each member of this family is an adaptive, one-step, two-stage method which is
 second-order accurate, unconditionally positivity-preserving, and linearly
-implicit. The parameter `α` is described by Kopecz and Meister (2018) and
+implicit. In this implementation the stage-values are conservative as well. 
+The parameter `α` is described by Kopecz and Meister (2018) and
 studied by Izgin, Kopecz and Meister (2022) as well as
-Torlo, Öffner and Ranocha (2022).
+Torlo, Öffner and Ranocha (2022).  
+
+This method supports adaptive time stepping, using the Patankar-weight denominators 
+``σ_i``, see Kopecz and Meister (2018), as first order approximations to estimate the error.
 
 The scheme was introduced by Kopecz and Meister for conservative production-destruction systems.
-For nonconservative production–destruction systems we use the straight forward extension
+For nonconservative production–destruction systems we use a straight forward extension
 analogous to [`MPE`](@ref).
 
 This modified Patankar-Runge-Kutta method requires the special structure of a
@@ -712,12 +719,15 @@ end
 
 A family of third-order modified Patankar-Runge-Kutta schemes for
 production-destruction systems, which is based on the two-parameter family of third order explicit Runge--Kutta schemes.
-Each member of this family is a one-step method with four-stages which is
+Each member of this family is an adaptive, one-step method with four-stages which is
 third-order accurate, unconditionally positivity-preserving, conservative and linearly
 implicit. In this implementation the stage-values are conservative as well.
 The parameters `α` and `β` must be chosen such that the Runge--Kutta coefficients are nonnegative,
 see Kopecz and Meister (2018) for details.
 
+These methods support adaptive time stepping, using the Patankar-weight denominators 
+``σ_i``, see Kopecz and Meister (2018), as second order approximations to estimate the error.
+
 The scheme was introduced by Kopecz and Meister for conservative production-destruction systems.
 For nonconservative production–destruction systems we use the straight forward extension
 analogous to [`MPE`](@ref).
@@ -809,15 +819,18 @@ end
 """
     MPRK43II(γ; [linsolve = ..., small_constant = ...])
 
-A family of third-order modified Patankar-Runge-Kutta schemes for (conservative)
+A family of third-order modified Patankar-Runge-Kutta schemes for
 production-destruction systems, which is based on the one-parameter family of third order explicit Runge--Kutta schemes with
 non-negative Runge--Kutta coefficients.
-Each member of this family is a one-step method with four stages which is
+Each member of this family is an adaptive, one-step method with four stages which is
 third-order accurate, unconditionally positivity-preserving, conservative and linearly
 implicit. In this implementation the stage-values are conservative as well. The parameter `γ` must satisfy
 `3/8 ≤ γ ≤ 3/4`.
 Further details are given in Kopecz and Meister (2018).
 
+This method supports adaptive time stepping, using the Patankar-weight denominators 
+``σ_i``, see Kopecz and Meister (2018), as second order approximations to estimate the error.
+
 The scheme was introduced by Kopecz and Meister for conservative production-destruction systems.
 For nonconservative production–destruction systems we use the straight forward extension
 analogous to [`MPE`](@ref).

src/sspmprk.jl

@@ -3,14 +3,18 @@
     SSPMPRK22(α, β; [linsolve = ..., small_constant = ...])
 
 A family of second-order modified Patankar-Runge-Kutta algorithms for
-production-destruction systems. Each member of this family is a one-step, two-stage method which is
+production-destruction systems. Each member of this family is an adaptive, one-step, two-stage method which is
 second-order accurate, unconditionally positivity-preserving, and linearly
 implicit. The parameters `α` and `β` are described by Huang and Shu (2019) and
 studied by Huang, Izgin, Kopecz, Meister and Shu (2023).
 The difference to [`MPRK22`](@ref) is that this method is based on the SSP formulation of
 an explicit second-order Runge-Kutta method. This family of schemes contains the [`MPRK22`](@ref) family,
 where `MPRK22(α) = SSMPRK22(0, α)` applies.
 
+This method supports adaptive time stepping, using the first order approximations 
+``(σ_i - u_i^n) / τ + u_i^n`` with ``τ=1+(α_{21}β_{10}^2)/(β_{20}+β_{21})``, 
+see (2.7) in Huang and Shu (2019), to estimate the error.
+
 The scheme was introduced by Huang and Shu for conservative production-destruction systems.
 For nonconservative production–destruction systems we use the straight forward extension
 analogous to [`MPE`](@ref).
@@ -408,7 +412,7 @@ end
     SSPMPRK43([linsolve = ..., small_constant = ...])
 
 A third-order modified Patankar-Runge-Kutta algorithm for
-production-destruction systems. This scheme is a one-step, two-stage method which is
+production-destruction systems. This scheme is a one-step, four-stage method which is
 third-order accurate, unconditionally positivity-preserving, and linearly
 implicit. The scheme is described by Huang, Zhao and Shu (2019) and
 studied by Huang, Izgin, Kopecz, Meister and Shu (2023).
@@ -429,6 +433,8 @@ You can also choose the parameter `small_constant` which is added to all Patanka
 to avoid divisions by zero. You can pass a value explicitly, otherwise `small_constant` is set to
 `floatmin` of the floating point type used.
 
+The current implementation only supports fixed time steps. 
+
 ## References
 
 - Juntao Huang, Weifeng Zhao and Chi-Wang Shu.