TL;DR
This paper introduces high-order, conservative, and positivity-preserving Patankar-type deferred correction schemes for solving production-destruction systems, extending the applicability of MPRK methods beyond third order with proven properties.
Contribution
The authors develop a novel high-order deferred correction scheme incorporating Patankar methods, ensuring conservation and positivity for PDS problems, which was not previously achieved.
Findings
Schemes can be constructed up to arbitrarily high order.
Numerical simulations confirm theoretical properties.
Methods maintain positivity and conservation at discrete level.
Abstract
Production-destruction systems (PDS) of ordinary differential equations (ODEs) are used to describe physical and biological reactions in nature. The considered quantities are subject to natural laws. Therefore, they preserve positivity and conservation of mass at the analytical level. In order to maintain these properties at the discrete level, the so-called modified Patankar-Runge-Kutta (MPRK) schemes are often used in this context. However, up to our knowledge, the family of MPRK has been only developed up to third order of accuracy. In this work, we propose a method to solve PDS problems, but using the Deferred Correction (DeC) process as a time integration method. Applying the modified Patankar approach to the DeC scheme results in provable conservative and positivity preserving methods. Furthermore, we demonstrate that these modified Patankar DeC schemes can be constructed up to…
| Notation | Meaning |
|---|---|
| Number of constituents and dimension of the ODE system | |
| Index for constituents | |
| Value of the th constituents | |
| Vector of all the constituents | |
| Number of time intervals | |
| Index for a timestep | |
| Timestep | |
| Variables at timestep | |
| Number of subtimeintervals in a timeinterval | |
| Index for subtimesteps | |
| Subtimestep | |
| Variable at subtimestep | |
| Vector of variables at all subtimestep | |
| Number of iterations of the DeC procedure | |
| Index of the iteration | |
| Variables for timestep at the subtimestep and iteration | |
| Vector of variables for all subtimesteps at the iteration | |
| First order operator of DeC procedure | |
| High order operator of DeC procedure | |
| High order operator of mPDeC procedure | |
| Solution of the system . |
Peer Reviews
No public reviews on file for this paper yet. If you reviewed it on a platform where reviews are public (OpenReview, ICLR, NeurIPS, ICML), you can paste yours below so the community can read it here.
Code & Models
Videos
No videos yet. Explain this paper in a talk, walkthrough, or lecture? Add one.
Arbitrary high-order, conservative and positive preserving Patankar-type deferred correction schemes
Philipp Öffner and Davide Torlo
Institute of Mathematics, University of Zurich, Switzerland Corresponding author: P. Öffner, [email protected] author: D. Torlo, [email protected]
(August 6th, 2019)
Abstract
Production-destruction systems (PDS) of ordinary differential equations (ODEs) are used to describe physical and biological reactions in nature. The considered quantities are subject to natural laws. Therefore, they preserve positivity and conservation of mass at the analytical level.
In order to maintain these properties at the discrete level, the so-called modified Patankar-Runge-Kutta (MPRK) schemes are often used in this context. However, up to our knowledge, the family of MPRK has been only developed up to third order of accuracy. In this work, we propose a method to solve PDS problems, but using the Deferred Correction (DeC) process as a time integration method. Applying the modified Patankar approach to the DeC scheme results in provable conservative and positivity preserving methods. Furthermore, we demonstrate that these modified Patankar DeC schemes can be constructed up to arbitrarily high order. Finally, we validate our theoretical analysis through numerical simulations.
1 Introduction
The modelling of geobiochemical processes or ecosystems leads often to systems of ordinary differential equations (ODEs) which can be formulated in the so-called production-destruction systems (PDS) as described in [5, 10] for example. To guarantee the physical and chemical laws, the quantities have to fulfil several conditions like positivity and conservation.
The applied numerical method should not violate these conditions and big efforts have been devoted to designing conservative and positivity preserving schemes, since classical approaches like Runge Kutta (RK) schemes do not guarantee these properties.
In [4] the authors suggest modified Patankar-type methods of first and second order which verify the desired properties, i.e., conservation and positivity. Recently, further extensions were done to construct modified Patankar-Runge-Kutta (MPRK) schemes of second and third order [14, 16, 15, 12, 11]. As the name suggests, all these schemes use, as a basic procedure, the Runge-Kutta method, which has been modified by weighting the production and destruction terms as suggested in [20]. Thanks to these weighting coefficients, the schemes are forced to maintain positivity of the variables and to conserve some quantities of interest. However, the described and constructed schemes are, up to our knowledge, at most third order accurate.
In this paper, we present a way to construct arbitrary high-order, positivity preserving, numerically robust and conservative schemes for PDS. Differently from previous schemes, we do not start building our schemes on RK methods. We consider the Deferred Correction (DeC) procedure, a high order time integration technique, and we modify it, in order to obtain a positivity preserving, conservative and arbitrary high-order scheme. Moreover, we provide a proof of the desired properties.
The paper is organised as follows.
In section 2 we introduce the production-destruction systems and we give a short introduction about the so-called Patankar trick and how it was applied in [5] to construct a modified Patankar-type scheme starting from the explicit Euler method. Afterwards, in section 3, we introduce the Deferred Correction (DeC) method and we discuss conservation and positivity for this classical formulation. In section 4, we build the main core of this work, explaining our modification of DeC through the Patankar trick (mPDeC) and we prove that the obtained mPDeC schemes are positive preserving, conservative and arbitrary high-order accurate. In section 5, we validate our theoretical investigations, considering three different benchmark problems, which are also discussed in different literature references, as [4, 14]. Finally, we give a summary and an outlook for possible extensions.
2 Production–Destruction Systems
In this paper we consider production-destruction systems (PDS) of the form
[TABLE]
where represents the vectors of constituents, denotes the time and the initial condition. Moreover, and represent the production and destruction rates of the -th constituent and both terms are assumed to be non-negative, i.e, for . These systems rise naturally to describe geochemical processes as it is described in [4, 5] and we recapitulate their notations and definitions in this section.
The production and destruction terms can also be written in a matrix form as follows
[TABLE]
where each term and are Lipschitz continuous functions and may depend linearly or non–linearly on . Furthermore, the term describes the rate of change from the -th to the -th constituent while is the rate at which the -th constituent is transformed into the -th.
We are interested in (fully) conservative and positive production–destruction systems. To clarify these expressions we repeat the definitions from [14].
Definition 2.1**.**
*The PDS (1) is called positive if positive initial values for imply positive solutions, for for all times .
The PDS (1) is called conservative if at any time , we have that*
[TABLE]
is fulfilled. In the analytic form (1), the conservation property (3) is equivalent to the following relation for the matrix representation (2)
[TABLE]
Moreover, the system is called fully conservative if additionally holds for all and .
As it is described in [14] every conservation PDS can be written in a fully conservative formulation. We can rewrite the two terms of (4) into one matrix of exchanging quantities defined as
[TABLE]
Clearly, from property (4), we have that . With this notation, let us define the total exchange rate for the -th constituent as
[TABLE]
A numerical method suited to solve a conservative and positive PDS (1) should mimic, at the discrete level, the continuous setting properties. For a one-step methods, we can introduce the discrete analogues of definitions (2.1).
Definition 2.2**.**
Let denote the approximation of at the time level . A one-step method
[TABLE]
with process function , is called
- •
unconditionally conservative* if for all and *
[TABLE]
holds;
- •
unconditionally positive* if for all and , we have that .*
Example 2.3**.**
Let us consider as an example the explicit Euler method. The method is defined by
[TABLE]
It is conservative since
[TABLE]
holds. Conversely, the explicit Euler method is not unconditionally positive. Consider a conservative and positive PDS (1) where we assume that the right hand side is not identical zero. Then, there exists a such that . Since the PDS is conservative, we can at least find one constituent , where . Choosing
[TABLE]
we obtain
[TABLE]
This demonstrates the violation of the positivity for the explicit Euler method for unbounded timesteps .
To build an unconditionally positive numerical scheme, Patankar had the idea in [20] of weighting the destruction term in the original explicit Euler methods with the following coefficient
[TABLE]
Hence, the scheme (13) is unconditionally positive, but the conservation relation is violated. In [4] a modification of the Patankar scheme (13) was presented, resulting in an unconditionally positive and conservative method. It is defined as follows.
[TABLE]
The scheme is implicit and can be solved inverting the mass matrix in the system where is
[TABLE]
The construction of the mass matrix must follow substantial prescriptions in order to preserve the positivity of the scheme, as suggested in [15].
Remark 2.4**.**
*Extensions of the modified Patankar scheme (14) to Runge-Kutta schemes were proposed in [14, 15] and further developed in [11, 12]. Special focus lies in the weighting of the production and destruction terms as it is investigated for example in [16] and references therein. Families of second and third order modified Patankar-Runge-Kutta (MPRK) schemes can be found in the mentioned literature. We do not provide the definition of MPRK because the modified Patankar scheme (14) already gives us the basic idea for the new methods we want to propose. We will prove that these methods are positivity preserving, conservative and arbitrary high-order.
3 Deferred Correction Methods
There are various approaches to solve numerically an ODE. A first approach is given by finite differences, where the derivative in time is replaced by differences of states in different timesteps. Backward and forward Euler are examples of this kind of strategy. Another approach would be to reformulate the ODE by integrating it in time. With different quadrature formulas and approximation techniques one can obtain various Runge-Kutta methods (explicit and implicit ones), see [9, 24] for details. However, we follow a different approach in this paper.
We start our investigation with the Deferred Correction (DeC) method introduced in [7]. In its original formulation, it is an explicit, arbitrary high order method for ODEs. Further extensions of DeC can be found in the literature, including semi-implicit approaches as in [19]. However, in this work we will not consider the semi-implicit framework. Instead, we will focus on the explicit DeC approach used by Abgrall in [1]. In our opinion, his notation describes DeC in a more compact way than in previous works [7, 6, 17].111We like to mention that Abgrall focused on DeC as a time integration scheme in the context of finite element methods. Applying a classical RK method, a dense mass matrix has to be inverted and Abgrall wanted to avoid this. By using a DeC scheme, instead, he showed that a mass matrix free approach is possible [1]. Nevertheless, the main idea is always the same and it is based on the Picard-Lindelöf theorem in the continuous setting. The theorem states the existence and uniqueness of solutions for ODEs. The classical proof makes use of the so-called Picard iterations to minimize the error and to prove convergence. The foundation of DeC relies on mimicking the Picard iterations at the discrete level. The approximation error decreases with several iteration steps. For the description of DeC, Abgrall introduces two operators: and .
Here, the operator represents a low-order easy-to-solve numerical scheme, e.g. the explicit Euler method, and is a high-order operator that can present difficulties in its practical solution, e.g. an implicit RK scheme. The DeC method can be written as a combination of these two operators.
Given a timeinterval we subdivide it into subintervals , where and and we mimic for every subinterval the Picard–Lindelöf theorem for both operators and . We drop the dependency on the timestep for subtimesteps and substates as denoted in Figure 1.
Then, the operator is given by
[TABLE]
Here, the term denotes an interpolation polynomial of order evaluated at the points . In particular, we use Lagrange polynomials , where and for any . Using these properties, we can actually compute the integral of the interpolants, thanks to a quadrature rule in the same points with weights . We can rewrite
[TABLE]
The operator represents an order numerical scheme if set equal to zero, i.e., . Unfortunately, the resulting scheme is implicit and, further, the terms may be non-linear. Because of this, the only formulation is not explicit and more efforts have to be made to solve it.
For this purpose, we introduce a simplification of the operator. Instead of using a quadrature formula at the points we evaluate the integral in equation (16) applying the left Riemann sum. The resulting operator is given by the forward Euler discretization for each state in the timeinterval, i.e.,
[TABLE]
with coefficients .
To simplify the notation and to describe DeC, we introduce the vector of states for the variable at all subtimesteps222We provide a table with all definitions and notations in the appendix A.
[TABLE]
Now, the DeC algorithm uses a combination of the and operators to provide an iterative procedure. The aim is to recursively approximate , the numerical solution of the scheme, similarly to the Picard iterations in the continuous setting. The successive states of the iteration process will be denoted by the superscript , where is the iteration index, e.g. . The total number of iterations (also called correction steps in the following) is denoted by . To describe the procedure, we have to refer to both the -th subtimestep and the -th iteration of the DeC algorithm. We will indicate the variable by . Finally, the DeC method can be written as
DeC Algorithm
[TABLE]
where is the number of iterations that we want to compute. Using the procedure (21), we need, in particular, as many iterations as the desired order of accuracy, i.e., .
Notice that, in every step, we solve the equations for the unknown variables which appears only in the formulation, the operator that can be easily inverted. Conversely, is only applied to already computed predictions of the solution . Therefore, the scheme 21 is completely explicit and arbitrary high order as stated in [1] with the following proposition.
Proposition 3.1**.**
Let and be two operators defined on , which depend on the discretization scale , such that
- •
* is coercive with respect to a norm, i.e., independent of , such that for any we have that*
[TABLE]
- •
* is Lipschitz with constant uniformly with respect to , i.e., for any *
[TABLE]
We also assume that there exists a unique such that . Then, if , the DeC is converging to and after iterations the error is smaller than .
Remark 3.2**.**
The DeC procedure is naturally conservative if is conservative, but it is not positivity preserving if is positivity preserving. Indeed, the coefficients of the operator can be negative and spoil the positivity of the scheme. This is one of the points that make us modify the classical DeC into the scheme that we propose in this work.
Remark 3.3**.**
Any DeC scheme can be interpreted as a RK scheme [6]. The main difference between RK and DeC is that the latter gives a general approach to the time discretization and does not require a specification of the coefficients for every order of accuracy. On the other side, by rewriting a DeC method as a RK scheme, it requires a number of stages equal to , which is bigger than classical RK stages. However, one can notice that every subtimestep is independent of another, so one can compute sequentially the corrections and in parallel the subtimesteps, obtaining a computational cost of just corrections.
Example 3.4**.**
For clarity, we provide here an example of a second order DeC scheme. To get this order of accuracy, we need DeC iterations and one subtimestep . Reminding that , the method (21) for the first step reads
[TABLE]
Substituting this term into the first correction steps leads finally to
[TABLE]
The correction step is not modifying the initial subtimestep. Therefore, with , we get
[TABLE]
where . This scheme coincides with the strong stability preserving Runge-Kutta method of second order [8].
Remark 3.5**.**
*Before we modify our DeC framework, we want to give some final remarks.
The presented DeC approach is not the most general version. In our description we always include both endpoints in the point distribution of the subtimesteps, i.e., and . However, this is not necessary, as it is already described in [7], where also Gauss-Legendre nodes are applied. Then, the approximation at the endpoint is done via extrapolation. Nevertheless, we do not consider in this work this class of point distribution.
Secondly, instead of using the explicit Euler method in , explicit high-order RK methods can also be applied. In principle, this yields a faster increase of the order of accuracy in the iterative procedure, but it has been shown, that it leads also to some problems of smoothness of the error behaviour as it is described in [6], which results in a drop down of the expected accuracy order. However, we will consider this approach in future research.*
4 Modified Patankar Deferred Correction Scheme
In this section, we are going to propose a positivity preserving, conservative and arbitrary high-order scheme, that will be denoted as modified Patankar Deferred Correction (mPDeC).
The DeC procedure (21) serves us as a starting point to construct this scheme, and, thanks to its structure, we will be able to prove the hypotheses of Proposition 3.1. This yield us directly the desired order condition for our modified DeC scheme without performing a specific Taylor expansion for every order of accuracy. We will adapt DeC in such a way to obtain all the properties we are interested in.
The conservation can be easily guaranteed by the consistency of the two operators, i.e., by the consistency of the two schemes described by and by .
Conversely, more effort is required to produce a positivity preserving scheme. For this purpose, we follow the ideas of Patankar [20] and Burchard et al. [4] of weighting the destruction and production terms in the scheme. Their aim is to obtain a mass matrix shaped as in the modified Patankar scheme (14) where all of the positive terms are collected on the diagonal, while the negative terms are put in the non-diagonal entries. This will guarantee that the mass matrix is diagonally dominant by columns, with positive diagonal values, and, thus, its inverse will be positive. Therefore, we introduce some coefficients similar to the ones proposed in (14).
Finally, as we have seen in the example (2.3), an explicit scheme is not positivity preserving and the investigation in [4, 16, 11] support our decision to modify the DeC scheme in order to get a fully implicit method.
Because of all the above mentioned considerations, we came to the conclusion of modifying the operator, to make it fully implicit. In particular, it has to depend on both the previous and the current corrections of the DeC procedure. We redefine it as follows.
[TABLE]
where if and if .
Remark 4.1**.**
*The modification of the scheme is done only through the coefficients on both the production and the destruction terms. The fact that these coefficients depend on the new correction means that we are modifying the mass matrix of the whole DeC correction step.
These coefficients allow to choose in which term of the mass matrix we want to put each term and , according to the sign of the coefficient. The pseudo-algorithm 1 provides the construction steps of the mass matrix, see B. There, it is straightforward to see that the diagonal terms are all positive and the off–diagonal are all negative. The index takes care of the sign of the destruction and production terms which are added in the mass matrix. It is inspired by the explanation given in [14, Remark 2.5], that states that, when negative entries in the Butcher Tableau of the RK scheme appear, one has to interchange the destruction terms with the production ones to guarantee the positivity preserving property. With the function we are taking this into account. In our opinion, it is complicated and unclear to investigate higher order RK schemes properties because of these exchanges depending on the Butcher Tableau. While, with this DeC approach, we can in few lines generalize every order scheme.
Moreover, it is helpful to notice that the coefficients that we are using to modify the contributions, namely , are converging to 1 as the iteration index of the DeC increases. In subsection 4.2 we will make this statement more precise and we will study how fast these coefficients converge to 1.*
Most of the terms in the operator will cancel out through the iteration process, therefore we keep the operator as presented in the original DeC (18).
[TABLE]
Now, we propose the modified Patankar DeC scheme as follows.
mPDeC Algorithm
[TABLE]
One can notice that, using the fact that initial states are identical for any correction , the DeC correction step (24) can be rewritten for , and into
[TABLE]
We keep both formulations (24) and (25) to prove different properties. The DeC formulation (24) will help us to demonstrate the accuracy order of the scheme whereas formulation (25) will be used to prove conservation and positivity. Before we start to prove these properties, we give a small example to get used to the formulation (24). Furthermore, we like to mention that although the algorithm (24) seems quite complex, it is actually easy to implement. We put a small pseudo-code in the B and we refer to the repository 333https://git.math.uzh.ch/abgrall_group/deferred-correction-patankar-scheme for a Julia version of the code.
Example 4.2**.**
We give a small example of the constructed method, applying the DeC approach at second order of accuracy, as already considered in example 3.4, i.e., DeC iterations and one subtimestep . In this case, we recall that and that .
The method (24) for the first step reads
[TABLE]
where the last step is obtained considering, again the fact that for the iteration all the states coincide. Collecting the mass matrix terms as in (15), one can solve the previous equation for . Substituting this term into the second iteration step leads finally to
[TABLE]
The correction step has no effect on the initial subtimestep. Therefore, we get with :
[TABLE]
where . This scheme coincides with a modified Runge Kutta Patankar scheme of second order as it is presented in [14] .
4.1 Conservation and positivity of modified Patankar DeC
In this section, we are proving that the proposed scheme is unconditionally conservative and positivity preserving.
Theorem 4.3**.**
The mPDeC scheme in (25) is unconditionally conservative for all substages, i.e.,
[TABLE]
for all and .
Proof.
Using formulation (25), we can easily see that
[TABLE]
To get this result, we have just used the definition of the scheme (25) in (27) and the property (4) of the production and destruction operators in (28). In the last step, we have exchanged the sums over and .
∎
To demonstrate the positivity of the scheme, we introduce some preliminary results.
Lemma 4.4**.**
The mass matrix of every correction step of the mPDeC scheme described in (25) is diagonal dominant by columns.
Proof.
At each step we are solving an implicit linear system where the mass matrix is given by
[TABLE]
Under the assumption that and are always positive, it is straightforward to see that all the terms of the sum of are positive by construction and that all the terms of the sum of the non-diagonal terms for are negative. Moreover, we can demonstrate that
[TABLE]
by showing
[TABLE]
where we have used the property of the and matrices to obtain the previous computation. Finally, this proves that the mass matrix is diagonally dominant by columns. ∎
Using Lemma 4.4 we prove the following theorem.
Theorem 4.5**.**
The mPDeC scheme defined in (25) is positivity preserving, i.e., if then , for all and .
Proof.
Using lemma 4.4, we can prove that the inverse of any mass matrix obtained from the DeC iterations is positive, i.e., . The proof follows the path of what was proposed in [14]. Using the Jacobi method, we can converge to with iterative matrices for , where
[TABLE]
Here, is the identity and is the diagonal of . If we denote the iteration matrix as , we can see that it has spectral radius smaller than one, since is diagonally dominant. This means that the Jacobi method is convergent to . Now, since and from previous lemma 4.4 and, by induction, also is positive, we can say that will be positive.
∎
4.2 Convergence order
To prove that the solution of the mPDeC procedure is high-order accurate, we mimic the proof of the original DeC convergence as in [1]. We denote by the solution of the operator, i.e., . This solution coincides with the solution of the classical operator defined in (16).
We want to prove that for each iteration step the following inequalities are fulfilled:
[TABLE]
which implies that for each iteration step we obtain one order of accuracy more than the previous iteration. After iterations we, finally, get
[TABLE]
To prove that the inequalities (34) and (36) are valid, we have to demonstrate the following
the coercivity of the operator (as in the inequality (34)) 2. 2.
the Lipschitz inequality for operator used in (36) 3. 3.
the high-order accuracy of the operator , i.e., .
Let us start with the coercivity lemma.
Lemma 4.6** (Coercivity of ).**
Given any , there exists a positive , such that, the operator verifies
[TABLE]
Proof.
We remind that the beginning states coincide for all the variables, i.e., . So, the evolution part simplifies in the two operators and we get the following relation
[TABLE]
This proves that with constant the equation (38) holds.
∎
Before proving Lipschitz continuity, we need two lemmas. The first one proves that each stage of the scheme is a first order approximation of the previous timestep.
Lemma 4.7**.**
For every subtimestep and correction , there exists a matrix , such that
[TABLE]
holds. Moreover, , where does not depend on .
Proof.
For any and , the equation (25) tells us that the mass matrix can be written as where does not depend on , but only on and the production–destruction functions. It is defined as
[TABLE]
This leads to an inverse
[TABLE]
Now, we can define by
[TABLE]
So, we can write
[TABLE]
∎
With the following lemma, we prove that the mPDeC process generates a Cauchy sequence similar to the continuous Picard iterations. Moreover, at each iteration, we differ from the previous step by an error of one order of accuracy more. We will drop the dependency on the subtimestep , as all the relations hold for all of them.
Lemma 4.8**.**
Let and verifying Lemma 4.7, then
[TABLE]
holds where are constants independent from .
Proof.
We prove the lemma by induction.
For , equation (43) follows directly from Lemma 4.7, i.e., .
Given , as induction hypothesis, (43) holds for , i.e.,
[TABLE]
where and denotes the th row of the matrix . We can prove that (43) is verified also for . Using Lemma (4.7), we obtain
[TABLE]
Inserting the induction step (44) we get
[TABLE]
Here, denotes the Hadamard product and . The induction step is evaluated for every entry . Using the regularity of , we expand its Taylor series in for every constituent . Thanks again to the result of Lemma (4.7), we can write
[TABLE]
where is the point of the Lagrange form of the remainder of the Taylor expansion. Hence, we can proceed as follows
[TABLE]
which finally proves equation (43) for . ∎
Now, let us prove on the Lipschitz continuity of the operator .
Lemma 4.9** (Lipschitz continuity of ).**
Let fulfil Lemma 4.7. Then, the operator is Lipschitz continuous with constant , i.e.,
[TABLE]
Proof.
Now, we apply Lemma (4.8) to substitute the new operator (22) with the original one of the classical DeC (17) adding an error of order to the operator. We get another order from the time integration, such that
[TABLE]
and, trivially, holds. Together, we obtain
[TABLE]
Now, we have to take care about the different variables in the operators. Let us start studying the operator . We are focusing on each line of the schemes for an arbitrary subtimestep. The difference is given by
[TABLE]
Now, we can compute the difference of the two terms
[TABLE]
In last step, we have used the regularity of the solutions and and the fact that brings an error of order zero times given by the time integration. Then, we have used the Lipschitz continuity of the functions .
Overall, the constant depends on the operators and and the Lemma is proven.
∎
Finally, we need to show that the solution of the operator is an -order accurate solution. This is given directly by the definition of the operator (22), since it is an -order accurate approximation of the original problem (1) when the two input coincide and, thus, the modification coefficients become 1 and the operator becomes the original one (17).
Theorem 4.10** (Convergence of mPDeC).**
Let and be the operators defined in (23) and (22) respectively. The mPDeC procedure (24) gives an approximation solution with order of accuracy equal to .
Proof.
With Lemma 4.6 we proved the coercivity of the operator , which verifies the inequality in (34). The definition of the mPDeC scheme (24) gives us the equality (35) and the Lipschitz continuity lemma 4.9 proves the inequality (36). Moreover, we know that is an -accurate approximation of the exact solution.
So, overall, we have
[TABLE]
∎
All the desired properties (unconditionally positivity, unconditionally conservation and high-order accuracy) are fulfilled by the proposed scheme.
5 Numerics
In this section, we validate our theoretical investigation of section 4 considering some test cases from [14, 4]. We focus here only on systems of ordinary differential equations (ODE) (stiff and non-stiff). However, the mPDeC schemes can be in general used as time-integration methods for a semidiscrete formulation of partial differential equations, where the spatial discretization is already provided by RD, DG, FR, (c.f. [2, 3, 23]) or your favourite space discretization method.
As part of future research, we will consider these schemes in real applications like non-equilibrium flows or shallow water equations as it was already done, for example, for MPRK together with a WENO approach in [11] or a DG one in [18]. In this work we focus on systems of ODEs. In all the numerical tests, we applied the mPDeC approach on equidistant subtimestep points distributions.
5.1 Linear Model
We start by considering a simple linear test case proposed in [4, 18]. The initial value problem for the PDS is given by
[TABLE]
The initial values of (51) are positive and we can rewrite the right hand side of the ODE system in a PDS format as follows
[TABLE]
and for . The system describes the exchange of mass between two constituents. The analytical solution is given by
[TABLE]
The problem is considered on the time interval and, analogously to [4], we use in the simulations. In Figure (2) we plot the analytical solution (dotted, blue line) and the approximated solutions using 2-nd (solid line, green) and 5-th (dash-dotted line, black) order mPDeC methods. The purple lines represent the sum of the constituents and they are constantly equal to 1, since the methods are conservative. Qualitatively, we see that the 5-th order method approximates better the analytical solutions. Furthermore, to verify the order of convergence of the methods, we consider also the error behavior of the different order schemes. Differently from Kopecz and Meister [14, 15] instead of calculating the relative errors, we compute the absolute discrete errors taken over all the timesteps and all the constituents:
[TABLE]
After a comparison between the final time error and the one proposed (53), we do not observe much discrepancy. Therefore, we will provide only results obtained with (53).
In Figure 3, the left picture shows the error decay for mPDeC schemes at different discretization scales . In the right picture, we plot the slope of the error decay for different orders of accuracy. These graphs demonstrate the high-order accuracy of the proposed methods and the expected convergence rates, validating the theoretical results. It is also possible to test the scheme with higher order of accuracy. However, we have notice a reduction of the order as we reach orders higher than , probably due to Runge phenomena. These are well known issues that arise also with the usual DeC methods [7] using equidistant points distribution in the subtimesteps. A possible solution of this problem can be the usage of Gauß-Lobatto nodes as point distributions. This and stability investigations will be part of future research.
5.2 Nonlinear test problem
In this next subsection, we consider the nonlinear test problem
[TABLE]
with initial condition . As before, this problem was proposed in [14]. The PDS system in the matrix formulation can be expressed by
[TABLE]
and for all other combinations of and . This system (54) is used to describe an algal bloom, that transforms nutrients via phytoplankton into detritus . In our test, we consider the time interval and . We calculate the reference solution with the strong stability preserving Runge-Kutta method 10 stages 4th order introduced by Ketcheson [13], further investigated in [22] and implemented in Julia, see [21] for details.
In Figure 4, the 6-th order mPDeC (black, dash-dotted lines) approximates very precisely the reference solution. The 2-nd order method (solid line, green) shows the same structure as the reference solution but it exhibits a severe error. However, the approximated second order solution is comparable with the results obtained in [14]. We see again that the conservation property is fulfilled in the purple lines.
Since we lack of an analytical solution, in the error plots, we compare successive errors between two refinements of the time mesh
[TABLE]
Here, the subscript indicates the number of equispaced timesteps used to subdivide the total time interval. The results are presented in Figure 5. As for the linear case, we can see that the error decay fulfils the expected behavior and that the order of accuracy tends to the correct one. The slight decrease of the slope function in the right picture using sixth order can be explained by the fact that the error values are close to machine precision in that area and this causes the deprecation of the slope.
These plots verify our theoretical investigations from section 4.
5.3 Robertson Test case
In the last test case, we prove the robustness of the mPDeC schemes in presence of stiff problems. The proposed test is the Robertson problem for a chemical reaction system. It consists of
[TABLE]
with initial conditions .444To avoid the division by zero in the mPDeC scheme, we slightly modify the initial condition in the practical implementation, i.e., with . The time interval of interest is . The PDS for (56) reads
[TABLE]
and zero for the other combinations.
In the Robertson test case, the numerical scheme has to deal with several time scales. Therefore, a constant time step size is not suitable for this purpose.
Following again the literature [14], we use increasing time steps with , where indicates the -th timestep. To make the small values visible on the plot, we multiply it by . As a comparison, we calculate the reference solution (dotted, blue line) using the function Rodas4555 A 4-th order A-stable stiffly stable Rosenbrock method with a stiff-aware 3rd order interpolant. from Julia, where we split the time-interval into 55 subdomains and we solve it on every subdomain with relative tolerance and absolute tolerance . We plot again a second order (green, solid lines) and fifth order (black, dashed-dotted lines) approximations generated by the mPDeC methods and, as it can be seen in figure 6, the designed methods produce reliable and robust results for this kind of stiff problems. As always, the conservation and the positivity properties are fulfilled.
Finally, we can say that the simulations run in this section express the quality of the mPDeC schemes. Moreover, they show that all the targeted properties are obtained even for very problematic test cases.
6 Summary and Outlook
In this paper, we presented a way to build positivity preserving, conservative and arbitrary high-order numerical schemes for production-destructions systems of equations. We adapted the idea of [4] to build modified Patankar type schemes to the Deferred Correction method as an underlying scheme. By altering the and operators using the modified Patankar trick, we were able to obtain schemes with the desired properties. We proved that the proposed modified Patankar DeC (mPDeC) schemes are arbitrary high-order, conservative and positivity preserving. In numerical simulations, we confirmed our theoretical considerations with various test cases.
However, further research can be pursued in this direction. As it was investigated in [15, 14] for families of MPRK, it is possible to study the accuracy and the stability of the method varying the weightings of the production-destruction terms of the schemes. In the spirit of the work [14], a change of the weighting of the Patankar modification in the and operators should be easily applicable to the mPDeC schemes and theoretical investigations will be considered in future research, in particular regarding the stability conditions. Also the distribution of the subtimesteps between and plays a big role on stability and accuracy of the scheme. Many choices are valid and the possible influence of the properties of the method must be carefully analysed. This idea is already work in progress for the classical DeC approach and will be extended to the mPDeC version in the future. Finally, we want to apply and analyse this type of schemes in context of partial differential equations. Here, we focus on applications and problems as described in [11, 12, 18]. As one can see, there are still many open questions and tasks for the mPDeC schemes and we are looking forward to continue our work in this field.
Appendix A Notation
We provide a small table 1 with the notation of symbols used along the paper. Even if some of the notations are ambiguous, the used indices should always clarify the referred meaning. We prefer to keep this notation to keep fluid the reading.
Appendix B Algorithm
We present a pseudo-code for the creation of the mass matrix in Algorithm 1 and one for the mPDeC algorithm in box Algorithm 2. Both algorithms are very simple. The first one consists of 3 loops: 2 for the constituents and one for the subtimesteps and an if statement. The second one consists of 3 nested loops: one for timesteps , one for corrections of the DeC algorithm and one for the subtimesteps .
Acknowledgements
P. Öffner has been funded by the the SNF project (Number 175784).
Davide Torlo is supported by ITN ModCompShock project funded by the European Union’s Horizon 2020 research and innovation program under the Marie Sklodowska-Curie grant agreement No 642768.
The reference list from the paper itself. Each links out to its DOI / PubMed record.
- 1[1] R. Abgrall. High order schemes for hyperbolic problems using globally continuous approximation and avoiding mass matrices. Journal of Scientific Computing , 73(2):461–494, Dec 2017.
- 2[2] R. Abgrall, E. l. Meledo, and P. Öffner. On the connection between residual distribution schemes and flux reconstruction. ar Xiv preprint ar Xiv:1807.01261 , 2018.
- 3[3] R. Abgrall and D. Torlo. Asymptotic preserving deferred correction residual distribution schemes. ar Xiv preprint ar Xiv:1811.09284 , 2018.
- 4[4] H. Burchard, E. Deleersnijder, and A. Meister. A high-order conservative patankar-type discretisation for stiff systems of production–destruction equations. Applied Numerical Mathematics , 47(1):1–30, 2003.
- 5[5] H. Burchard, E. Deleersnijder, and A. Meister. Application of modified patankar schemes to stiff biogeochemical models for the water column. Ocean Dynamics , 55(3-4):326–337, 2005.
- 6[6] A. Christlieb, B. Ong, and J.-M. Qiu. Integral deferred correction methods constructed with high order runge-kutta integrators. Mathematics of Computation , 79(270):761–783, 2010.
- 7[7] A. Dutt, L. Greengard, and V. Rokhlin. Spectral Deferred Correction Methods for Ordinary Differential Equations. BIT Numerical Mathematics , 40(2):241–266, 2000.
- 8[8] S. Gottlieb, D. I. Ketcheson, and C.-W. Shu. Strong stability preserving Runge-Kutta and multistep time discretizations . World Scientific, 2011.
