PositiveIntegrators.jl API

ConservativePDSProblem(P, u0, tspan, p = NullParameters();
                       p_prototype = nothing,
                       analytic = nothing,
                       std_rhs = nothing,
                       linear_invariants = nothing)

A structure describing a conservative system of ordinary differential equation in form of a production-destruction system (PDS). P denotes the function defining the production matrix $P$. The diagonal of $P$ contains production terms without destruction counterparts. u0 is the vector of initial conditions and tspan the time span (t_initial, t_final) of the problem. The optional argument p can be used to pass additional parameters to the function P.

The function P can be given either in the out-of-place form with signature production_terms = P(u, p, t) or the in-place form P(production_terms, u, p, t).

Keyword arguments:

• p_prototype: If P is given in in-place form, p_prototype or copies thereof are used to store evaluations of P. If p_prototype is not specified explicitly and P is in-place, then p_prototype will be internally set to zeros(eltype(u0), (length(u0), length(u0))).
• analytic: The analytic solution of a PDS must be given in the form f(u0,p,t). Specifying the analytic solution can be useful for plotting and convergence tests.
• std_rhs: The standard ODE right-hand side evaluation function callable as du = std_rhs(u, p, t) for the out-of-place form and as std_rhs(du, u, p, t) for the in-place form. Solvers that do not rely on the production-destruction representation of the ODE, will use this function instead to compute the solution. If not specified, a default implementation calling P is used
• linear_invariants: The rows of this matrix contain the linear invariants of the ODE. Certain solvers or callbacks require this matrix. Note that this feature is experimental and its API may change in future releases.

References

• Hans Burchard, Eric Deleersnijder, and Andreas Meister. "A high-order conservative Patankar-type discretisation for stiff systems of production-destruction equations." Applied Numerical Mathematics 47.1 (2003): 1-30. DOI: 10.1016/S0168-9274(03)00101-6

PDSProblem(P, D, u0, tspan, p = NullParameters();
           p_prototype = nothing,
           analytic = nothing,
           std_rhs = nothing,
           linear_invariants = nothing)

A structure describing a system of ordinary differential equations in form of a production-destruction system (PDS). P denotes the function defining the production matrix $P$. The diagonal of $P$ contains production terms without destruction counterparts. D is the function defining the vector of destruction terms $D$ without production counterparts. u0 is the vector of initial conditions and tspan the time span (t_initial, t_final) of the problem. The optional argument p can be used to pass additional parameters to the functions P and D.

The functions P and D can be used either in the out-of-place form with signature production_terms = P(u, p, t) or the in-place form P(production_terms, u, p, t).

Keyword arguments:

• p_prototype: If P is given in in-place form, p_prototype or copies thereof are used to store evaluations of P. If p_prototype is not specified explicitly and P is in-place, then p_prototype will be internally set to zeros(eltype(u0), (length(u0), length(u0))).
• analytic: The analytic solution of a PDS must be given in the form f(u0,p,t). Specifying the analytic solution can be useful for plotting and convergence tests.
• std_rhs: The standard ODE right-hand side evaluation function callable as du = std_rhs(u, p, t) for the out-of-place form and as std_rhs(du, u, p, t) for the in-place form. Solvers that do not rely on the production-destruction representation of the ODE, will use this function instead to compute the solution. If not specified, a default implementation calling P and D is used.
• linear_invariants: The rows of this matrix contain the linear invariants of the ODE. Certain solvers or callbacks require this matrix. Note that this feature is experimental and its API may change in future releases.

References

• Hans Burchard, Eric Deleersnijder, and Andreas Meister. "A high-order conservative Patankar-type discretisation for stiff systems of production-destruction equations." Applied Numerical Mathematics 47.1 (2003): 1-30. DOI: 10.1016/S0168-9274(03)00101-6

prob_pds_bertolazzi

Positive and conservative autonomous nonlinear PDS

\[\begin{aligned} \mathbf{u}'=\begin{pmatrix}2 &-1 &-1\\-1 &2 &-1\\-1& -1& 2\end{pmatrix}\begin{pmatrix}5u_2u_3/(10^{-2} + (u_2u_3)^2) + u_2u_3/(10^{-16} + u_2u_3(10^{-8} + u_2u_3))\\ 10u_1u_3^2\\ 0.1(u_3 - u_2 - 2.5)^2u_1u_2\end{pmatrix} \end{aligned}\]

with initial value $\mathbf{u}_0 = (0.0, 1.0, 2.0)^T$ and time domain $(0.0, 1.0)$. There is one independent linear invariant, e.g. $u_1+u_2+u_3 = 3.0$.

References

• Enrico Bertolazzi. "Positive and conservative schemes for mass action kinetics." Computers and Mathematics with Applications 32 (1996): 29-43. DOI: 10.1016/0898-1221(96)00142-3

prob_pds_brusselator

Positive and conservative autonomous nonlinear PDS

\[\begin{aligned} u_1' &= -u_1,\\ u_2' &= -u_2u_5,\\ u_3' &= u_2u_5,\\ u_4' &= u_5,\\ u_5' &= u_1 - u_2u_5 + u_5^2u_6 - u_5,\\ u_6' &= u_2u_5 - u_5^2u_6, \end{aligned}\]

with initial value $\mathbf{u}_0 = (10.0, 10.0, 0.0, 0.0, 0.1, 0.1)^T$ and time domain $(0.0, 20.0)$. There are two independent linear invariants, e.g. $u_1+u_4+u_5+u_6 = 10.2$ and $u_2+u_3 = 10.0$.

References

• Luca Bonaventura, and Alessandro Della Rocca. "Unconditionally Strong Stability Preserving Extensions of the TR-BDF2 Method." Journal of Scientific Computing 70 (2017): 859 - 895. DOI: 10.1007/s10915-016-0267-9

prob_pds_linmod

Positive and conservative autonomous linear PDS

\[\begin{aligned} u_1' &= u_2 - 5u_1,\\ u_2' &= 5u_1 - u_2, \end{aligned}\]

with initial value $\mathbf{u}_0 = (0.9, 0.1)^T$ and time domain $(0.0, 2.0)$. There is one independent linear invariant, e.g. $u_1+u_2 = 1$.

References

• Hans Burchard, Eric Deleersnijder, and Andreas Meister. "A high-order conservative Patankar-type discretisation for stiff systems of production-destruction equations." Applied Numerical Mathematics 47.1 (2003): 1-30. DOI: 10.1016/S0168-9274(03)00101-6

prob_pds_linmod_inplace

Same as prob_pds_linmod but with in-place computation.

prob_pds_minmapk

Positive and nonconservative autonomous nonlinear PDS

\[\begin{aligned} u_1' &= k_6u_6-k_7u_1-k_1u_1u_2 +k_2u_4,\\ u_2' &= k_5u_3-k_1u_2u_2,\\ u_3' &= k_2u_4-k_3u_1u_3+k_4u_5-k_5u_3,\\ u_4' &= k_1u_1u_2-k_2u_4,\\ u_5' &= k_3u_1u_3-k_4u_5,\\ u_6' &= k_7u_1-k_6u_6, \end{aligned}\]

with constants

\[\begin{aligned} k_1 &=\frac{100}{3}, & k_2 &=\frac{1}{3}, & k_3 &=50,\\ k_4 &=0.5, & k_5 &=\frac{10}{3} , & k_6 &= 0.1,\\ k_7 &= 0.1. \end{aligned}\]

The initial value is $\mathbf{u}_0 = (0.1, 0.175, 0.15, 1.15, 0.81, 0.5)^T$ and the time domain $(0, 200)$. There are two independent linear invariants, e.g. $u_1+u_4+u_6=1.75$ and $u_2+u_3+u_4+u_5 =2.285$.

References

• Sergio Blanes, Arieh Iserles, and Shev Macnamara. "Positivity preserving methods for ordinary differential equations." ESAIM: Mathematical Modelling and Numerical Analysis 56 (2022): 1843–1870. DOI: 10.1051/m2an/2022042
• Otto Hadač, František Muzika, Vladislav Nevoral, Michal Přibyl, and Igor Schreiber "Minimal oscillating subnetwork in the Huang-Ferrell model of the MAPK cascade." PLoS ONE 12 (2017): e0178457. DOI: 10.1371/journal.pone.0178457

prob_pds_nonlinmod

Positive and conservative autonomous nonlinear PDS

\[\begin{aligned} u_1' &= -\frac{u_1u_2}{u_1 + 1.0},\\ u_2' &= \frac{u_1u_2}{u_1 + 1.0} - 0.3u_2,\\ u_3' &= 0.3 u_2, \end{aligned}\]

with initial value $\mathbf{u}_0 = (9.98, 0.01, 0.01)^T$ and time domain $(0.0, 30.0)$. There is one independent linear invariant, e.g. $u_1+u_2+u_3 = 10.0$.

References

• Hans Burchard, Eric Deleersnijder, and Andreas Meister. "A high-order conservative Patankar-type discretisation for stiff systems of production-destruction equations." Applied Numerical Mathematics 47.1 (2003): 1-30. DOI: 10.1016/S0168-9274(03)00101-6

prob_pds_npzd

Positive and conservative autonomous nonlinear PDS

\[\begin{aligned} u_1' &= 0.01u_2 + 0.01u_3 + 0.003u_4 - \frac{u_1u_2}{0.01 + u_1},\\ u_2' &= \frac{u_1u_2}{0.01 + u_1}- 0.01u_2 - 0.5( 1 - e^{-1.21u_2^2})u_3 - 0.05u_2,\\ u_3' &= 0.5(1 - e^{-1.21u_2^2})u_3 - 0.01u_3 - 0.02u_3,\\ u_4' &= 0.05u_2 + 0.02u_3 - 0.003u_4 \end{aligned}\]

with initial value $\mathbf{u}_0 = (8.0, 2.0, 1.0, 4.0)^T$ and time domain $(0.0, 10.0)$. There is one independent linear invariant, e.g. $u_1+u_2+u_3+u_4 = 15.0$.

References

• Hans Burchard, Eric Deleersnijder, and Andreas Meister. "Application of modified Patankar schemes to stiff biogeochemical models for the water column." Ocean Dynamics 55 (2005): 326-337. DOI: 10.1007/s10236-005-0001-x

prob_pds_robertson

Positive and conservative autonomous nonlinear PDS

\[\begin{aligned} u_1' &= -0.04u_1+10^4 u_2u_3,\\ u_2' &= 0.04u_1-10^4 u_2u_3-3⋅10^7 u_2^2,\\ u_3' &= 3⋅10^7 u_2^2, \end{aligned}\]

with initial value $\mathbf{u}_0 = (1.0, 0.0, 0.0)^T$ and time domain $(0.0, 10^{11})$. There is one independent linear invariant, e.g. $u_1+u_2+u_3 = 1.0$.

References

• Ernst Hairer, Gerd Wanner. "Solving Ordinary Differential Equations II - Stiff and Differential-Algebraic Problems." 2nd Edition, Springer (2002): Section IV.1.

prob_pds_sir

Positive and conservative autonomous nonlinear PDS

\[\begin{aligned} u_1' &= -2u_1u_2,\\ u_2' &= 2u_1u_2 - u_2,\\ u_3' &= u_2, \end{aligned}\]

with initial value $\mathbf{u}_0 = (0.99, 0.005, 0.005)^T$ and time domain $(0.0, 20.0)$. There is one independent linear invariant, e.g. $u_1+u_2+u_3 = 1.0$.

References

• Ronald E. Mickens, and Talitha M. Washington. "NSFD discretizations of interacting population models satisfying conservation laws." Computers and Mathematics with Applications 66 (2013): 2307-2316. DOI: 10.1016/j.camwa.2013.06.011

prob_pds_stratreac

Positive and nonconservative autonomous nonlinear PDS

\[\begin{aligned} u_1' &= r_5 - r_6 - r_7,\\ u_2' &= 2r_1 - r_2 + r_3 - r_4 + r_6 - r_9 + r_{10} - r_{11},\\ u_3' &= r_2 - r_3 - r_4 - r_5 - r_7 - r_8,\\ u_4' &= -r_1 -r_2 + r_3 + 2r_4+r_5+2r_7+r_8+r_9,\\ u_5' &= -r_8+r_9+r_{10}-r_{11},\\ u_6' &= r_8-r_9-r_{10}+r_{11}, \end{aligned}\]

with reaction rates

\[\begin{aligned} r_1 &=2.643⋅ 10^{-10}σ^3 u_4, & r_2 &=8.018⋅10^{-17}u_2 u_4 , & r_3 &=6.12⋅10^{-4}σ u_3,\\ r_4 &=1.567⋅10^{-15}u_3 u_2 , & r_5 &= 1.07⋅ 10^{-3}σ^2u_3, & r_6 &= 7.11⋅10^{-11}⋅ 8.12⋅10^6 u_1,\\ r_7 &= 1.2⋅10^{-10}u_1 u_3, & r_8 &= 6.062⋅10^{-15}u_3 u_5, & r_9 &= 1.069⋅10^{-11}u_6 u_2,\\ r_{10} &= 1.289⋅10^{-2}σ u_6, & r_{11} &= 10^{-8}u_5 u_2, \end{aligned}\]

where

\[\begin{aligned} T &= t/3600 \mod 24,\quad T_r=4.5,\quad T_s = 19.5,\\ σ(T) &= \begin{cases}1, & T_r≤ T≤ T_s,\\0, & \text{otherwise}.\end{cases} \end{aligned}\]

The initial value is $\mathbf{u}_0 = (9.906⋅10^1, 6.624⋅10^8, 5.326⋅10^{11}, 1.697⋅10^{16}, 4⋅10^6, 1.093⋅10^9)^T$ and the time domain $(4.32⋅ 10^{4}, 3.024⋅10^5)$. There are two independent linear invariants, e.g. $u_1+u_2+3u_3+2u_4+u_5+2u_6=(1,1,3,2,1,2)\cdot\mathbf{u}_0$ and $u_5+u_6 = 1.097⋅10^9$.

References

• Stephan Nüsslein, Hendrik Ranocha, and David I. Ketcheson. "Positivity-preserving adaptive Runge-Kutta methods." Communications in Applied Mathematics and Computer Science 16 (2021): 155-179. DOI: 10.2140/camcos.2021.16.155

MPE([linsolve = ..., small_constant = ...])

The first-order modified Patankar-Euler algorithm for production-destruction systems. This one-step, one-stage method is first-order accurate, unconditionally positivity-preserving, and linearly implicit.

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}$,

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 or a ConservativePDSProblem.

You can optionally choose the linear solver to be used by passing an algorithm from LinearSolve.jl as keyword argument linsolve. You can also choose the parameter small_constant which is added to all Patankar-weight denominators 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. "A high-order conservative Patankar-type discretisation for stiff systems of production-destruction equations." Applied Numerical Mathematics 47.1 (2003): 1-30. DOI: 10.1016/S0168-9274(03)00101-6

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 adaptive, one-step, two-stage method which is second-order accurate, unconditionally positivity-preserving, and linearly 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).

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 a straight forward extension analogous to MPE.

This modified Patankar-Runge-Kutta method requires the special structure of a PDSProblem or a ConservativePDSProblem.

You can optionally choose the linear solver to be used by passing an algorithm from LinearSolve.jl as keyword argument linsolve. You can also choose the parameter small_constant which is added to all Patankar-weight denominators to avoid divisions by zero. You can pass a value explicitly, otherwise small_constant is set to floatmin of the floating point type used.

References

• Hans Burchard, Eric Deleersnijder, and Andreas Meister. "A high-order conservative Patankar-type discretisation for stiff systems of production-destruction equations." Applied Numerical Mathematics 47.1 (2003): 1-30. DOI: 10.1016/S0168-9274(03)00101-6
• Stefan Kopecz and Andreas Meister. "On order conditions for modified Patankar-Runge-Kutta schemes." Applied Numerical Mathematics 123 (2018): 159-179. DOI: 10.1016/j.apnum.2017.09.004
• Thomas Izgin, Stefan Kopecz, and Andreas Meister. "On Lyapunov stability of positive and conservative time integrators and application to second order modified Patankar-Runge-Kutta schemes." ESAIM: Mathematical Modelling and Numerical Analysis 56.3 (2022): 1053-1080. DOI: 10.1051/m2an/2022031
• Davide Torlo, Philipp Öffner, and Hendrik Ranocha. "Issues with positivity-preserving Patankar-type schemes." Applied Numerical Mathematics 182 (2022): 117-147. DOI: 10.1016/j.apnum.2022.07.014

SSPMPRK22(α, β; [linsolve = ..., small_constant = ...])

A family of second-order modified Patankar-Runge-Kutta algorithms for 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 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 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.

This modified Patankar-Runge-Kutta method requires the special structure of a PDSProblem or a ConservativePDSProblem.

You can optionally choose the linear solver to be used by passing an algorithm from LinearSolve.jl as keyword argument linsolve. You can also choose the parameter small_constant which is added to all Patankar-weight denominators to avoid divisions by zero. You can pass a value explicitly, otherwise small_constant is set to floatmin of the floating point type used.

References

• Juntao Huang and Chi-Wang Shu. "Positivity-Preserving Time Discretizations for Production–Destruction Equations with Applications to Non-equilibrium Flows." Journal of Scientific Computing 78 (2019): 1811–1839 DOI: 10.1007/s10915-018-0852-1
• Juntao Huang, Thomas Izgin, Stefan Kopecz, Andreas Meister and Chi-Wang Shu. "On the stability of strong-stability-preserving modified Patankar-Runge-Kutta schemes." ESAIM: Mathematical Modelling and Numerical Analysis 57 (2023):1063–1086 DOI: 10.1051/m2an/2023005

MPRK43I(α, β; [linsolve = ..., small_constant = ...])

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 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.

These modified Patankar-Runge-Kutta methods require the special structure of a PDSProblem or a ConservativePDSProblem.

You can optionally choose the linear solver to be used by passing an algorithm from LinearSolve.jl as keyword argument linsolve. You can also choose the parameter small_constant which is added to all Patankar-weight denominators to avoid divisions by zero. You can pass a value explicitly, otherwise small_constant is set to floatmin of the floating point type used.

References

• Stefan Kopecz and Andreas Meister. "Unconditionally positive and conservative third order modified Patankar–Runge–Kutta discretizations of production–destruction systems." BIT Numerical Mathematics 58 (2018): 691–728. DOI: 10.1007/s10543-018-0705-1

MPRK43II(γ; [linsolve = ..., small_constant = ...])

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 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.

These modified Patankar-Runge-Kutta methods require the special structure of a PDSProblem or a ConservativePDSProblem.

You can optionally choose the linear solver to be used by passing an algorithm from LinearSolve.jl as keyword argument linsolve. You can also choose the parameter small_constant which is added to all Patankar-weight denominators to avoid divisions by zero. To display the default value for data type type evaluate MPRK43II(gamma).small_constant_function(type), where type can be, e.g., Float64.

References

• Stefan Kopecz and Andreas Meister. "Unconditionally positive and conservative third order modified Patankar–Runge–Kutta discretizations of production–destruction systems." BIT Numerical Mathematics 58 (2018): 691–728. DOI: 10.1007/s10543-018-0705-1

SSPMPRK43([linsolve = ..., small_constant = ...])

A third-order modified Patankar-Runge-Kutta algorithm for 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). The difference to MPRK43I or MPRK43II is that this method is based on the SSP formulation of an explicit third-order Runge-Kutta method.

The scheme was introduced by Huang, Zhao and Shu for conservative production-destruction systems. For nonconservative production–destruction systems we use the straight forward extension analogous to MPE.

This modified Patankar-Runge-Kutta method requires the special structure of a PDSProblem or a ConservativePDSProblem.

You can optionally choose the linear solver to be used by passing an algorithm from LinearSolve.jl as keyword argument linsolve. You can also choose the parameter small_constant which is added to all Patankar-weight denominators to avoid divisions by zero. To display the default value for data type type evaluate SSPMPRK43. small_constant_function(type), where type can be, e.g., Float64.

The current implementation only supports fixed time steps.

References

• Juntao Huang, Weifeng Zhao and Chi-Wang Shu. "A Third-Order Unconditionally Positivity-Preserving Scheme for Production–Destruction Equations with Applications to Non-equilibrium Flows." Journal of Scientific Computing 79 (2019): 1015–1056 DOI: 10.1007/s10915-018-0881-9
• Juntao Huang, Thomas Izgin, Stefan Kopecz, Andreas Meister and Chi-Wang Shu. "On the stability of strong-stability-preserving modified Patankar-Runge-Kutta schemes." ESAIM: Mathematical Modelling and Numerical Analysis 57 (2023):1063–1086 DOI: 10.1051/m2an/2023005

MPDeC(K; [nodes = :gausslobatto, linsolve = ..., small_constant = ...])

A family of arbitrary order modified Patankar-Runge-Kutta algorithms for production-destruction systems. Each member of this family is an adaptive, one-step method which is Kth order accurate, unconditionally positivity-preserving, and linearly implicit. The integer K must be chosen to satisfy 2 ≤ K ≤ 10. Available node choices are Lagrange or Gauss-Lobatto nodes, with the latter being the default. These methods support adaptive time stepping, using the numerical solution obtained with one correction step less as a lower-order approximation to estimate the error. The MPDeC schemes were introduced by Torlo and Öffner (2020) for autonomous conservative production-destruction systems and further investigated in Torlo, Öffner and Ranocha (2022).

For nonconservative production–destruction systems we use a straight forward extension analogous to MPE. A general discussion of DeC schemes applied to non-autonomous differential equations and using general integration nodes is given by Ong and Spiteri (2020).

The MPDeC methods require the special structure of a PDSProblem or a ConservativePDSProblem.

You can optionally choose the linear solver to be used by passing an algorithm from LinearSolve.jl as keyword argument linsolve. You can also choose the parameter small_constant which is added to all Patankar-weight denominators to avoid divisions by zero. You can pass a value explicitly, otherwise small_constant is set to 1e-300 in double precision computations or floatmin of the floating point type used.

References

• Davide Torlo and Philipp Öffner. "Arbitrary high-order, conservative and positivity preserving Patankar-type deferred correction schemes." Applied Numerical Mathematics 153 (2020): 15-34. DOI: 10.1016/j.apnum.2020.01.025

• Davide Torlo, Philipp Öffner, and Hendrik Ranocha. "Issues with positivity-preserving Patankar-type schemes." Applied Numerical Mathematics 182 (2022): 117-147. DOI: 10.1016/j.apnum.2022.07.014

• Benjamin W. Ong and Raymond J. Spiteri. "Deferred Correction Methods for Ordinary Differential Equations." Journal of Scientific Computing 83 (2020): Article 60. DOI: 10.1007/s10915-020-01235-8

isnegative(sol::ODESolution)

Returns true if sol.u contains negative elements.

Please note that negative values may occur when plotting the solution, depending on the interpolation used.

See also isnonnegative.

isnonnegative(u)

Negation of isnegative.

rel_max_error_tend(sol, ref_sol)

Returns the relative maximum error between sol and ref_sol at time sol.t[end].

rel_max_error_overall(sol, ref_sol)

Returns the maximum of the relative maximum errors between sol and ref_sol over all time steps.

rel_l1_error_tend(sol, ref_sol)

Returns the relative l1 error between sol and ref_sol at time sol.t[end].

rel_l2_error_tend(sol, ref_sol)

Returns the relative l2 error between sol and ref_sol at time sol.t[end].

work_precision_adaptive(prob, algs, labels, abstols, reltols, alg_ref;
                        adaptive_ref = false,
                        abstol_ref = 1e-14,
                        reltol_ref = 1e-13,
                        compute_error = rel_max_error_tend,
                        seconds = 2,
                        numruns = 20,
                        kwargs...)

Returns a dictionary to create work-precision diagrams. The problem prob is solved by each algorithm in algs for all tolerances defined in abstols and reltols. For the respective tolerances the error and computing time are stored in the dictionary. If the solve is not successful for the given tolerances, then (Inf, Inf) is stored in the dictionary. The strings in the array labels are used as keys of the dictionary. The reference solution used for error computations is computed with the algorithm alg_ref. Additional keyword arguments are passed on to solve.

Keyword arguments:

• adaptive_ref: If true the refenerce solution is computed adaptively with tolerances abstol_ref and reltol_ref. Otherwise $10^5$ steps are used.
• abstol_ref: See adaptive_ref.
• reltol_ref: See adaptive_ref.
• compute_error(sol::ODESolution, ref_sol::ODESolution): A function to compute the error between sol and ref_sol.
• seconds: If the measured computing time of a single solve is larger than seconds, then this computing time is stored in the dictionary.
• numruns: If the measured computing time of a single solve is less than or equal to seconds, then numruns solves are performed and the median of the respective computing times is stored in the dictionary.

work_precision_adaptive(prob, algs, labels, abstols, reltols, alg_ref;
                        adaptive_ref = false,
                        abstol_ref = 1e-14,
                        reltol_ref = 1e-13,
                        compute_error = rel_max_error_tend,
                        seconds = 2,
                        numruns = 20,
                        kwargs...)

Adds work-precision data to the dictionary dict, which was created with work_precion_fixed_adaptive. See work_precision_adaptive for the meaning of the inputs.

work_precision_fixed(prob, algs, labels, dts, alg_ref;
                     compute_error = rel_max_error_tend,
                     seconds = 2,
                     numruns = 20)

Returns a dictionary to create work-precision diagrams. The problem prob is solved by each algorithm in algs for all the step sizes defined in dts. For each step size the error and computing time are stored in the dictionary. If the solve is not successful for a given step size, then (Inf, Inf) is stored in the dictionary. The strings in the array labels are used as keys of the dictionary. The reference solution used for error computations is computed with the algorithm alg_ref.

Keyword arguments:

• compute_error(sol::ODESolution, ref_sol::ODESolution): Function to compute the error between sol and ref_sol.
• seconds: If the measured computing time of a single solve is larger than seconds, then this computing time is stored in the dictionary.
• numruns: If the measured computing time of a single solve is less than or equal to seconds, then numruns solves are performed and the median of the respective computing times is stored in the dictionary.

work_precision_fixed!(dict, prob, algs, labels, dts, alg_ref;
                      compute_error = rel_max_error_tend,
                      seconds = 2,
                      numruns = 20)
)

Adds work-precision data to the dictionary dict, which was created with work_precion_fixed. See work_precision_fixed for the meaning of the inputs.