Unconditionally Positive and Conservative Third Order Modified Patankar-Runge-Kutta Discretizations of Production-Destruction Systems

TL;DR

This paper develops third order Modified Patankar-Runge-Kutta schemes that are unconditionally positive and conservative for production-destruction systems, extending previous second order methods with proven theoretical conditions and numerical validation.

Contribution

It introduces necessary and sufficient conditions for third order MPRK schemes, enabling their construction independent of specific systems.

Findings

Third order MPRK schemes are unconditionally positive and conservative.

Theoretical conditions for third order accuracy are established.

Numerical experiments confirm the theoretical results.

Abstract

Modified Patankar-Runge-Kutta (MPRK) schemes are numerical methods for the solution of positive and conservative production-destruction systems. They adapt explicit Runge-Kutta schemes to ensure positivity and conservation irrespective of the time step size. The first two members of this class were introduced in [Burchard et. al.: A high-order conservative Patankar-type discretisation for stiff systems of production-destruction equations. Appl. Numer. Math., 47(1):1-30, 2003] and have been successfully applied in a large number of applications. Recently, we introduced a general definition of MPRK schemes and presented a thorough investigation of first and second order MPRK schemes in [Kopecz. S, Meister, A.: On order conditions for modified patankar-runge-kutta schemes. arXiv:1702.04589 [math.NA], 2017.]. A potentially third order Patankar-type method was introduced in [Formaggia L., Scotti, A.: Positivity and conservation properties of some integration schemes for mass action kinetics. SIAM J. Numer. Anal., 49(3):1267-1288, 2011.]. This method uses the MPRK22(1) scheme of [Burchard et. al.: A high-order conservative Patankar-type discretisation for stiff systems of production-destruction equations. Appl. Numer. Math., 47(1):1-30, 2003] as a predictor and a modification of the BDF(3) multistep method as a corrector. This method is at most third order and at least second order accurate. In this paper we continue the work of [Kopecz. S, Meister, A.: On order conditions for modified patankar-runge-kutta schemes. arXiv:1702.04589 [math.NA], 2017.] and present necessary and sufficient conditions for third order MPRK schemes. For the first time, we introduce MPRK schemes, which are third order accurate independent of the specific positive and conservative system under consideration. The theoretical results are confirmed by numerical experiments.

Equations

$$\frac{d{y}_i}{dt}(t)=P_i(\mathbf{y}(t))-D_i(\mathbf{y}(t)),i=1,\dots ,N.$$

$$P_i(\mathbf{y})=\sum _{j=1}^N{p}_{ij}(\mathbf{y}),D_i(\mathbf{y})=\sum _{j=1}^N{d}_{ij}(\mathbf{y}),$$

$$P_i(\mathbf{y})=\sum _{\begin{array}{c}j=1\\ j\ne i\end{array}}^N{p}_{ij}(\mathbf{y}),D_i(\mathbf{y})=\sum _{\begin{array}{c}j=1\\ j\ne i\end{array}}^N{d}_{ij}(\mathbf{y}).$$

$$\frac{d}{dt}\sum_{i=1}^{N}y_{i}=\sum_{i=1}^{N}\bigl{(}P_{i}(\mathbf{y})-D_{i}(\mathbf{y})\bigr{)}=\sum_{i,j=1}^{N}\bigl{(}p_{ij}(\mathbf{y})-d_{ij}(\mathbf{y})\bigr{)}=\sum_{i,j=1}^{N}\bigl{(}\underbrace{p_{ij}(\mathbf{y})-d_{ji}(\mathbf{y})}_{=0}\bigr{)}=0.$$

$${\mathbf{y}}^{n+1}={\mathbf{y}}^n+\Delta t\Phi (t^n,{\mathbf{y}}^n,{\mathbf{y}}^{n+1},\Delta t)$$

$$\sum _{i=1}^N\left(y_i^{n+1}-y_i^n\right)=0$$

$$\begin{array}{rl}{y}^{(k)}& =y^n+\Delta t\sum _{\nu =1}^{k-1}{a}_{k\nu }f(t^n+c_{\nu }\Delta t,y^{(\nu )}),k=1,\dots ,s,\\ y^{n+1}& =y^n+\Delta t\sum _{k=1}^s{b}_{k}f(t^n+c_k\Delta t,y^{(k)}).\end{array}$$

$$\begin{array}[]{c|c}\mathbf{c}&\mathbf{A}\\ \hline\cr&\mathbf{b}\end{array},$$

$$y_i^{(k)}$$

$$y_i^{n+1}$$

$$y_{i}^{n+1}=y_{i}^{n}+\Delta t\sum_{k=1}^{s}b_{k}\sum_{j=1}^{N}\biggl{(}p_{ij}(\mathbf{y}^{(k)})-d_{ij}(\mathbf{y}^{(k)})\frac{y_{i}^{n+1}}{\sigma_{i}}\biggr{)},$$

$$y_i^{n+1}=\frac{y_i^n+\Delta t{\sum }_{k=1}^s{b}_k{\sum }_{j=1}^N{p}_{ij}({\mathbf{y}}^{(k)})}{1+\Delta t{\sum }_{k=1}^s{b}_k{\sum }_{j=1}^N{d}_{ij}({\mathbf{y}}^{(k)})/{\sigma }_i}.$$

$$y_i^{(k)}$$

$$y_i^{n+1}$$

$${\mathbf{M}}^{(k)}{\mathbf{y}}^{(k)}$$

$$\mathbf{M}{\mathbf{y}}^{n+1}$$

$$\begin{array}{rl}{m}_{ii}^{(k)}& =1+\Delta t\sum _{\nu =1}^{k-1}{a}_{k\nu }\sum _{j=1}^N{d}_{ij}({\mathbf{y}}^{(\nu )})/{\pi }_i^{(k)}>0,i=1,\dots ,N,\\ m_{ij}^{(k)}& =-\Delta t\delta \sum _{\nu =1}^{k-1}{a}_{k\nu }{p}_{ij}({\mathbf{y}}^{(\nu )})/{\pi }_j^{(k)}\le 0,i,j=1,\dots ,N,i\ne j,\end{array}$$

$$\begin{array}{rl}{m}_{ii}& =1+\Delta t\sum _{k=1}^s{b}_k\sum _{j=1}^N{d}_{ij}({\mathbf{y}}^{(k)})/{\sigma }_i>0,i=1,\dots ,N,\\ m_{ij}& =-\Delta t\sum _{k=1}^s{b}_k{p}_{ij}({\mathbf{y}}^{(k)})/{\sigma }_j\le 0,i,j=1,\dots ,N,i\ne j.\end{array}$$

$$\frac{d{y}_i}{dt}(t)={\hat{P}}_i(\mathbf{y}(t))-{\hat{D}}_i(\mathbf{y}(t)),i=1,\dots ,N,$$

$${\hat{P}}_i(\mathbf{y})=\{\begin{array}{ll}\mu y_I^{\kappa },& i=J,\\ 0,& \text{otherwise},\end{array}\text{and}{\hat{D}}_i(\mathbf{y})=\{\begin{array}{ll}\mu y_I^{\kappa },& i=I,\\ 0,& \text{otherwise},\end{array}$$

$$\frac{d{y}_I}{dt}=-\mu y_I^{\kappa },\frac{d{y}_J}{dt}=\mu y_I^{\kappa },\frac{d{y}_i}{dt}=0,i\in \{1,\dots ,N\}\setminus \{I,J\}.$$

$$y_I(t)$$

$$y_J(t)$$

$$y_i(t)$$

$$y_I(t)$$

$$y_J(t)$$

$$y_i(t)$$

$${\hat{P}}_i(\mathbf{y})=\sum _{j=1}^N{\hat{p}}_{ij}(\mathbf{y}),{\hat{D}}_i(\mathbf{y})=\sum _{j=1}^N{\hat{d}}_{ij}(\mathbf{y}),$$

$${\hat{p}}_{ij}(\mathbf{y})=\{\begin{array}{ll}\mu y_I^{\kappa },& i=J\text{ and }j=I,\\ 0,& \text{otherwise},\end{array}{\hat{d}}_{ij}(\mathbf{y})=\{\begin{array}{ll}\mu y_I^{\kappa },& i=I\text{ and }j=J,\\ 0,& \text{otherwise},\end{array}$$

$$0\le {\tilde{m}}_{ij},{\tilde{m}}_{ij}^{(k)}\le 1,i,j=1,\dots ,N,$$

$${\xi }_0(\Delta t)+\mu {\xi }_1(\Delta t)+{\mu }^2{\xi }_2(\Delta t)+\cdots +{\mu }^n{\xi }_n(\Delta t)=\mathcal{O}(\Delta t^s)\text{ for all }\mu >0$$

$${\xi }_i(\Delta t)=\mathcal{O}(\Delta t^s),i=0,\dots ,N.$$

$${\xi }_0(\Delta t)+{\mu }_i{\xi }_1(\Delta t)+{\mu }_i^2{\xi }_2(\Delta t)+\cdots +{\mu }_i^n{\xi }_n(\Delta t)=\mathcal{O}(\Delta t^s)$$

$$\mathbf{V}\mathbf{\xi }(\Delta t)=\mathcal{O}(\Delta t^s),$$

Full text

The first two members of this class, the first order MPE scheme and the second order MPRK22(1) scheme, were introduced in [BDM03] and have been successfully applied in a large number of applications. Recently, we introduced a general definition of MPRK schemes and presented a thorough investigation of first and second order MPRK schemes in [KM17].

A potentially third order Patankar-type method was introduced in [FS11]. This method uses the MPRK22(1) scheme of [BDM03] as a predictor and a modification of the BDF(3) multistep method as a corrector. It restricts to the MPRK22(1) approximation, whenever the positivity of the corrector cannot be guaranteed. Hence, this method is at most third order accurate and at least second order accurate.

In this paper we continue the work of [KM17] and present necessary and sufficient conditions for third order MPRK schemes. For the first time, we introduce MPRK schemes, which are third order accurate independent of the specific positive and conservative system under consideration. The theoretical results derived within the first part are subsequently confirmed by numerical experiments for the entire domain of linear and nonlinear as well as nonstiff and stiff systems of differential equations.

1 Introduction

A wide variety of mathematical models for real life problems are given in the form of a system of partial differential equations including stiff production-destruction terms. The development of numerical methods is therefore often based on splitting approaches, where the discretization of the convection and diffusion terms is conducted within a first step and the approximation of the source terms is realized subsequently. Thereby, the time step size of the first part should also be applicable within the second step and even particular properties like conservativity of the source terms and positivity of the constituents have to be maintained independent of the time step size in order to obtain a reliable, appropriate and efficient simulation. Whereas finite volume schemes are well established for the discretization of convection-diffusion equations, the development of unconditionally positivity preserving and conservative methods of higher order for stiff systems of ordinary differential equations is still a challenge. To overcome this gap, the paper is devoted to the derivation and investigation of a class of third order modified Patankar schemes. Therefore, we consider production-destruction systems (PDS) of the form

[TABLE]

By $\mathbf{y}=(y_{1},\dotsc,y_{N})^{T}$ we denote the vector of constituents, which depends on time $t$. Both, the production terms $P_{i}$ and the destruction terms $D_{i}$ are assumed to be non-negative, that is $P_{i},D_{i}\geq 0$ for $i=1,\dots,N$. Furthermore, the production and destruction terms can be written as

[TABLE]

where $d_{ij}(\mathbf{y})\geq 0$ is the rate at which the $i$th constituent transforms into the $j$th component, while $p_{ij}(\mathbf{y})\geq 0$ is the rate at which the $j$th constituent transforms into the $i$th component.

We are interested in PDS which are positive as well as fully conservative.

Definition 1.1**.**

The PDS (1) is called positive, if positive initial values, $y_{i}(0)>0$ for $i=1,\dots,N$, imply positive solutions, $y_{i}(t)>0$ for $i=1,\dots,N$, for all times $t>0$.

Definition 1.2**.**

The PDS (1), (2) is called conservative, if for all $i,j=1,\dots,N$ and $\mathbf{y}\geq 0$, we have $p_{ij}(\mathbf{y})=d_{ji}(\mathbf{y}).$ The system is called fully conservative, if in addition $p_{ii}(\mathbf{y})=d_{ii}(\mathbf{y})=0$ holds for all $\mathbf{y}\geq 0$ and $i=1,\dots,N$.

In the following, we will assume that the PDS (1) is fully conservative, since every conservative PDS can be rewritten as an equivalent fully conservative PDS. For a fully conservative PDS, (2) can be written as

[TABLE]

But for the sake of a simple notation, we will always use the form (2). Examples of positive and conservative PDS, which model academic as well as realistic applications, can be found in Section 3.

If a PDS is conservative the sum of its constituents $\sum_{i=1}^{N}y_{i}(t)$ remains constant in time, since we have

[TABLE]

This motivates the definition of a conservative numerical scheme.

Definition 1.3**.**

Let $\mathbf{y}^{n}$ denote an approximation of $\mathbf{y}(t^{n})$ at time level $t^{n}$. The one-step method

[TABLE]

is called

• •

unconditionally conservative, if

[TABLE]

is satisfied for all $n\in\mathbb{N}$ and $\Delta t>0$.

• •

unconditionally positive, if it guarantees $\mathbf{y}^{n+1}>0$ for all $\Delta t>0$ and $\mathbf{y}^{n}>0$.

An explicit $s$-stage Runge-Kutta method for the solution of an ordinary differential equation $y^{\prime}(t)=f(t,y(t))$ is given by

[TABLE]

The method is characterized by its coefficients $a_{k\nu}$, $b_{k}$, $c_{k}$ for $k=1,\dots,s$, $\nu=1,\dots,k-1$ and can be represented by the Butcher tableau

[TABLE]

with $\mathbf{A}=(a_{k\nu})_{k,\nu=1,\dots,s}$, $\mathbf{c}=(c_{1},\dots,c_{s})^{T}$ and $\mathbf{b}=(b_{1},\dots,b_{s})$. Applied to (1) the method reads

[TABLE]

The idea of the modified Patankar-Runge-Kutta (MPRK) schemes is to adapt explicit Runge-Kutta schemes in such a way that they become positive irrespective of the chosen time step size $\Delta t$, while still maintaining their inherent property to be conservative. One approach to achieve unconditional positivity is the so-called Patankar-trick introduced in [Pat80] as source term linearization in the context of turbulent flow. If we modify (3b) and add a weighting of the destruction terms like

[TABLE]

we obtain

[TABLE]

Thus, if $y_{i}^{n}$, the weights $b_{k}$ for $k=1,\dots,s$ and $\sigma_{i}$ are positive, so is $y_{i}^{n+1}$. The crucial idea of the Patankar-trick is to multiply the destruction terms with weights that comprise $y_{i}^{n+1}$ as a factor themselves.

Weighting only the destruction terms will result in a non-conservative scheme. So the production terms have to be weighted accordingly as well. Since we have $d_{ij}(\mathbf{y})=p_{ji}(\mathbf{y})$, the proper weight for $p_{ij}(\mathbf{y}^{(k)})$ is $y_{j}^{n+1}/\sigma_{j}$.

Definition 1.4**.**

Given a non-negative Runge-Kutta matrix $\mathbf{A}=(a_{ij})_{i,j=1,\dots,s}$, non-negative weights $b_{1},\dotsc,b_{s}$ and $\delta\in\{0,1\}$, the scheme

[TABLE]

for $i=1\dots,N$, is called modified Patankar-Runge-Kutta scheme (MPRK) if

$\pi_{i}^{(k)}$ and $\sigma_{i}$ are unconditionally positive for $k=1,\dots,s$ and $i=1,\dots,N$, 2. 2.

$\pi_{i}^{(k)}$ is independent of $y_{i}^{(k)}$ and $\sigma_{i}$ is independent of $y_{i}^{n+1}$ for $k=1,\dots,s$ and $i=1,\dots,N$.

The weights $1/\sigma_{i}$ and $1/\pi_{i}^{(k)}$ are called Patankar-weights and the denominators $\sigma_{i}$ and $\pi_{i}^{(k)}$ are called Patankar-weight denominators (PWD).

Due to the introduction of the Patankar-weights, $s$ linear systems of size $N\times N$ need to be solved to obtain the stage values and the approximation at the next time level. In consideration of $p_{ii}=d_{ii}=0$ for $i=1,\dots,N$, the scheme (4) can be written in matrix-vector notation as

[TABLE]

with $\mathbf{P}(\mathbf{y}^{n})=(P_{1}(\mathbf{y}^{n}),\dots,P_{N}(\mathbf{y}^{n}))^{T}$ and

[TABLE]

for $k=1,\dots,s$ and

[TABLE]

If $\delta=0$, the matrices $\mathbf{M}^{(k)}$ become diagonal and the production terms appear on the right hand side of (5a).

The following two lemmas of [KM17] show that MPRK schemes, as defined in Definition 1.4, are indeed unconditionally positive and conservative.

Lemma 1.5**.**

A MPRK scheme (4) applied to a conservative PDS is unconditionally conservative. If $\delta=1$, the same holds for all stage values, this is $\sum_{i=1}^{N}(y_{i}^{(k)}-y_{i}^{n})=0$ for $k=1,\dots,s$.

Lemma 1.6**.**

A MPRK scheme (4) is unconditionally positive. The same holds for all the stages of the scheme, this is for all $\Delta t>0$ and $\mathbf{y}^{n}>0$ we have $\mathbf{y}^{(k)}>0$ for $k=1,\dots,s$.

Lemmas 1.5 and 1.6 show that the MPRK schemes as defined in Definition 1.4 possess the desired properties of unconditional positivity and conservation. The only quantities left to choose are the PWDs $\sigma_{i}$ and $\pi_{i}^{(k)}$ for $i=1,\dots,N$ and $k=1,\dots,s$. In [KM17] we introduced the second order MPRK22($a_{21}$) schemes, which use $\pi_{i}=y_{i}^{n}$ and $\sigma_{i}=y_{i}^{n}(y_{i}^{(2)}/y_{i}^{n})^{1/a_{21}}$ for $i=1,\dots,N$.

The MPRK22(1) scheme is equivalent to the original MPRK scheme introduced in [BDM03]. This scheme and the first order modified Patankar-Euler scheme of [BDM03] have been successfully applied to solve physical, biogeochemical and ecosystem models ([BDM05, BBK*+*06, BMZ09, HB10b, HB10a, MB10, WHK13, SD17]), and have also proven beneficial in astrophysics [KM10, Gre16].

In [SB11] it was demonstrated that the MPRK22(1) scheme outperforms standard Runge-Kutta and Rosenbrock methods when solving biogeochemical models without multiple source compounds per system reactions. The same was shown with respect to workload in [BR16], where the Brusselator PDS was solved with different time integration schemes.

In [BBKS07, BRBM08] second order schemes, which ensure conservation in a biochemical sense, were introduced. These schemes require the solution of a non-linear equation in each time step. Other schemes for the same purpose were recently presented in [RB15]. These explicit schemes incorporate the MPRK schemes of [BDM03] to achieve multi-element conservation for stiff problems. A potentially third order Patankar-type scheme was introduced in [FS11]. This method uses the MPRK22(1) scheme a as predictor and a modification of the BDF(3) multistep method as a corrector. It yields the MPRK22(1) approximation, whenever the positivity of the corrector approximation cannot be guaranteed.

Modified Patankar-Runge-Kutta type schemes are also used in the context of partial differential equations. An implicit first order Patankar-type scheme based on a third order SDIRK method was presented in [MO14] and applied to the shallow water equations. In [OH16] Patankar-type Runge-Kutta schemes for linear PDEs were investigated.

In the present paper, we extend the work of [KM17] to third order. We present necessary and sufficient conditions for third order three-stage MPRK schemes, and introduce two families of third order MPRK methods. To our knowledge, this is the first time that third order Patankar-type schemes are presented.

The paper is organized as follows. Section 2 deals with the derivation of conditions for third order three-stage MPRK schemes. In this section also novel third order MPRK schemes are introduced. The test problems of Section 3 are used in Section 4 to show numerical experiments with these novel schemes.

2 Third order MPRK schemes

In this section we assume that all occurring PDS are positive. To prove convergence of the MPRK schemes we investigate the local truncation errors. In doing so we make frequent use of the Landau symbol $\mathcal{O}$ and omit to specify the limit process $\Delta t\to 0$ each time. As customary, we identify $y_{i}^{n}$ and $y_{i}(t^{n})$ for $i=1,\dots,N$ when studying the truncation errors. Furthermore, since we are dealing with positive PDS we assume $y_{i}^{n}>0$ for $i=1,\dots,N$.

To derive necessary conditions that guarantee a certain order of an MPRK scheme, it suffices to consider specific PDS. In this regard, the following family of PDS will be very helpful. Given parameters $I,J\in\{1,\dots,N\}$, $I\neq J$, $\mu>0$ and $\kappa\in\{1,2\}$, we consider

[TABLE]

and initial values $y_{i}(0)=1$ for $i=1,\dots,N$. The PDS can be written in the form

[TABLE]

For $\kappa=1$ the exact solution is given by

[TABLE]

and for $\kappa=2$ the exact solution is given by

[TABLE]

This shows that the PDS is positive. Writing

[TABLE]

with

[TABLE]

we see that the PDS is also fully conservative.

The following lemmas are helpful to find necessary and sufficient conditions for the PWDs of third order MPRK schemes. The first one ensures the boundedness of the MPRK approximations and the stage values.

Lemma 2.1**.**

Let $\mathbf{M}$, $\mathbf{M}^{(k)}$ be given by (7), (6) with $\delta=1$ and $\mathbf{M}^{-1}=(\widetilde{m}_{ij})$, $(\mathbf{M}^{(k)})^{-1}=\widetilde{m}_{ij}^{(k)}$. Then we have

[TABLE]

for $k=1,\dots,s$.

Proof.

See [KM17]. ∎

The next lemma is useful to separate complicated conditions for the PWDs into simpler ones.

Lemma 2.2**.**

The identity

[TABLE]

is equivalent to

[TABLE]

Proof.

Let $0<\mu_{0}<\mu_{1}<\dots<\mu_{n}$. Since (9) is valid for all $\mu>0$, we have

[TABLE]

for $i=0,\dots,n$. This can be rewritten as

[TABLE]

with

[TABLE]

Since all $\mu_{i}$ are distinct for $i=0,\dots,n$, $\mathbf{V}$ is a Vandermonde matrix, and hence, regular. Multiplication of (11) with $\mathbf{V}^{-1}$ yields

[TABLE]

Hence, (10) is satisfied. On the other hand, if (10) is satisfied, so is (9). ∎

As a three-stage MPRK scheme is build on an explicit three-stage Runge-Kutta scheme with non-negative parameters, we must characterize these schemes in some way. It is well known, that an explicit three-stage Runge-Kutta scheme

[TABLE]

is third order accurate, if the conditions

[TABLE]

are satisfied. The last condition particularly implies $a_{21},a_{32},b_{3}\neq 0$.

The following lemma shows, that all explicit three-stage Runge-Kutta schemes of order three can be parameterized by families with at most two free parameters. Later we will use this lemma to characterize all explicit three-stage Runge-Kutta schemes of order three with non-negative parameters.

Lemma 2.3**.**

All explicit third order Runge-Kutta schemes can be parameterized with at most two parameters. The following three cases can occur:

Case I:

[TABLE]

with $\alpha,\beta\neq 0$, $\alpha\neq\beta$, $\alpha\neq\frac{2}{3}$.

Case II:

[TABLE]

with $\gamma\neq 0$.

Case III:

[TABLE]

with $\gamma\neq 0$.

Proof.

See [SWP12, RR01]. ∎

To build an MPRK scheme based on an explicit third order three-stage Runge-Kutta scheme, we must ensure the non-negativity of the occurring parameters. The following lemma characterizes all such Runge-Kutta schemes.

Lemma 2.4**.**

All explicit three-stage third order Runge-Kutta schemes with non-negative parameters can be represented by the following Butcher tableaus:

Case I:

[TABLE]

with

[TABLE]

and

[TABLE]

Case II:

[TABLE]

with

[TABLE]

Proof.

According to Lemma 2.3 we have to distinguish between three different cases.

In case III we have $b_{3}=\gamma$ and $a_{31}=-\frac{1}{4\gamma}$. Thus, $b_{3}>0$ implies $a_{31}<0$ and hence, negative Runge-Kutta parameters are inevitable in this case.

In case II we must restrict $b_{3}=\gamma>0$, such that $a_{31}=\frac{2}{3}-\frac{1}{4\gamma}\geq 0$ and $b_{2}=\frac{3}{4}-\gamma\geq 0$. This is the case for $\frac{3}{8}\leq\gamma\leq\frac{3}{4}$.

In case I things become more technical. First, we need to ensure $a_{21}=\alpha>0$. Since $\beta\neq 0$ and $\beta\neq\alpha$, we must have $a_{32}={\beta(\beta-\alpha)}/{(\alpha(2-3\alpha))}>0$, and due to $\beta=a_{31}+a_{32}$, $\beta>0$ must hold as well. Next, we present conditions for the non-negativity of the remaining Runge-Kutta parameters subject to $\alpha,\beta>0$. From $a_{31}=(3\alpha\beta(1-\alpha)-\beta^{2})/(\alpha(2-3\alpha))\geq 0$ we can conclude

[TABLE]

To ensure $a_{32}={\beta(\beta-\alpha)}/{(\alpha(2-3\alpha))}>0$, we must have

[TABLE]

The requirement $b_{1}=1+(2-3(\alpha+\beta))/(6\alpha\beta)\geq 0$ demands

[TABLE]

To guarantee $b_{2}=(3\beta-2)/(6\alpha(\beta-\alpha))\geq 0$, the conditions

[TABLE]

are necessary. Finally, $b_{3}=({2-3\alpha})/({6\beta(\beta-\alpha)})>0$ implies

[TABLE]

Merging the above conditions, we obtain (13), in which $\alpha_{0}$ denotes the unique solution of $3\alpha(1-\alpha)=(3\alpha-2)/(6\alpha-3)$. The region of feasibility, which contains all pairs $(\alpha,\beta)$ that ensure non-negativity of the Runge-Kutta parameters, is shown in Figure 1.

∎

To see that $\pi_{i}=y_{i}^{n}+\mathcal{O}(\Delta t)$ and $\rho_{i}=y_{i}^{n}+\mathcal{O}(\Delta t)$ for $i=1,\dots,N$ are necessary conditions for the PWDs of a third order MPRK scheme, the following lemma will be helpful.

Lemma 2.5**.**

Given an explicit three-stage Runge-Kutta scheme of order three with non-negative parameters, the nonlinear system

[TABLE]

has the unique positive solution

[TABLE]

Proof.

First, we note that (15) is a solution of (14), owing to (12).

Next, we show that no other solutions exist. If $b_{2}=0$, (12d) becomes $b_{3}(a_{31}+a_{32})=1/2$ and hence (14a) reads $y/2=1/2$, which implies $y=1$. Similar, owing to (12e), $y=1$ is the only positive solution of (14b), hence, we can conclude $x=1$ from (14c). Thus, (15) is the only solution of (14), if $b_{2}=0$.

From now on, we assume $b_{2}\neq 0$. As $b_{3}(a_{31}+a_{32})\neq 0$ as well, since $a_{32}>0$, $b_{3}>0$ and $a_{31}\geq 0$, (14a) represents a line and (14b) represents an ellipse in the $x$-$y$-plane. There are at most two intersections of the line and the ellipse, and thus, the system (14) has at most two solutions. We already know that one of them is (15). To find the hypothetical other one, we assume $y\neq 1$ and compute the intersection of the line (14a) and the hyperbola (14c). Subtraction of (12d) from (14a) yields

[TABLE]

which becomes

[TABLE]

owing to (14c). Division by $1-y\neq 0$ results in

[TABLE]

and thus, we have

[TABLE]

Next, we compute the intersection of the ellipse (14b) and the hyperbola (14c). We obtain

[TABLE]

by subtracting (12e) from (14b), and utilization of (14c) yields

[TABLE]

Owing to $1-y^{2}\neq 0$, as $0<y\neq 1$, we can divide by $1-y^{2}$ and find

[TABLE]

Altogether, owing to (16) and (17), only $b_{2}=b_{3}$ yields a potential second solution of (14). This solution reads

[TABLE]

The remaining question is, if there are any explicit third order Runge-Kutta schemes with non-negative parameters that satisfy $b_{2}=b_{3}$. According to Lemma 2.4, we have to consider two cases to answer this question. In case I, $b_{2}=b_{3}$ can be written as

[TABLE]

This is satisfied if

[TABLE]

which can be reformulated as

[TABLE]

holds true. Thus, $(\alpha,\beta)$ must be a point on the boundary of the circle with center $(1/3,1/3)$ and radius $\sqrt{2}/3$.

Figure 2 shows the feasible region from Lemma 2.4, together with the circle (19). Computing the intersection of the circle (19) and the parabola $3\alpha(1-\alpha)$, yields $\alpha=2/3$. As this value of $\alpha$ is excluded in case I, there is no solution of the system (14) in the situation of case I.

In case II of Lemma 2.4, $b_{2}=b_{3}$ is equivalent to $3/4-\gamma=\gamma$, which is satisfied for $\gamma=3/8$. Due to $a_{21}=a_{31}+a_{32}=2/3$, (18) becomes (15). All things considered, we have shown that (15) is the unique positive solution of (14). ∎

An MPRK scheme (4) with three stages is given by

[TABLE]

for $i=1,\dots,N$. The next theorem gives necessary and sufficient conditions for the Patankar-weights of a third order three stage MPRK scheme.

Theorem 2.6**.**

Given an explicit three-stage third order Runge-Kutta scheme with non-negative weights, the MPRK scheme (20) is of third order, if and only if the conditions

[TABLE]

are satisfied.

Proof.

We use the notation $\phi^{*}$ to represent $\phi(\mathbf{y}^{*})$ for a given function $\phi$. As (20) is an MPRK scheme, all Patankar-weights are positive, i. e. $\pi_{i}>0$, $\rho_{i}>0$ and $\sigma_{i}>0$ for $i=1,\dots,N$.

The Runge-Kutta scheme is of third order and substitution of (12a) and (12b) into (12d) and (12e) shows

[TABLE]

hold true. Furthermore, equation (12f) ensures $a_{21},a_{31},b_{3}>0$.

For a sufficiently smooth function $\phi$ and some $y_{i}^{*}=y_{i}^{n}+\mathcal{O}(\Delta t)$, we can expand $\phi(\mathbf{y}^{*})$ in the form

[TABLE]

in which $\mathbf{H}_{\phi}^{n}$ denotes the Hessian matrix of $\phi$ evaluated at $\mathbf{y}^{n}$. The Taylor expansion of the exact solution of (1) reads

[TABLE]

for $i=1,\dots,N$.

To derive necessary conditions, which allow for third order accuracy, we assume that the MPRK scheme (20) is third order accurate, this is

[TABLE]

for $i=1,\dots,N$. Subtracting (24) from (25) shows $y_{i}^{n+1}-y(t_{i}^{n+1})=\mathcal{O}(\Delta t^{4})$. Utilizing (20i) and (24) this can be written as

[TABLE]

for $i=1,\dots,N$.

From now on, we focus on the solution of the PDS (8) with $I,J\in\{1,\dots,N\}$, $I\neq J$, $\kappa\in\{1,2\}$ and $\mu>0$. Since $\widehat{P}_{I}=0$ and $\widehat{D}_{I}=\widehat{d}_{IJ}=\mu y_{I}^{\kappa}$, it follows that $\partial D_{I}/\partial\mathbf{y}=(\partial D_{I}/\partial y_{I})\mathbf{e}_{I}^{T}=\mu\kappa y_{I}^{\kappa-1}\mathbf{e}_{I}^{T}$, with $\mathbf{e}_{I}$ denoting the $I$th unit column vector and $\mathbf{H}_{D_{I}}=(\partial^{2}\widehat{D}_{I}^{n}/{\partial y_{I}^{2}})\mathbf{e}_{I}\mathbf{e}_{I}^{T}=\mu\kappa(\kappa-1)y_{I}^{\kappa-2}\mathbf{e}_{I}\mathbf{e}_{I}^{T}$. Hence, (26) with $i=I$ becomes

[TABLE]

For $k=1,2$ the destruction terms can be expanded as

[TABLE]

since derivatives of order higher than two vanish. Substituting this into (27), results in

[TABLE]

Owing to (20b), we have

[TABLE]

and from (20f), (2) and (30) we see

[TABLE]

Before we introduce (30) and (2) into (29), we set $\kappa=1$, which implies $(\partial^{2}\widehat{D}_{I}^{n}/{\partial y_{I}^{2}})=\mu\kappa(\kappa-1)y_{I}^{\kappa-2}=0$. Hence, we can drop the terms containing second derivatives in (29) and (2). We can exploit these conditions, when we consider (29) with $\kappa=2$, as some terms can be neglected. Setting $\kappa=1$, we have

[TABLE]

according to (30) and (2), since $\partial^{2}\widehat{D}_{I}^{n}/{\partial y_{I}^{2}}=0$. Substituting this into (29) yields

[TABLE]

A subsequent division by $-\widehat{D}_{I}^{n}=-\mu y_{I}^{n}\neq 0$ and utilization of $\partial\widehat{D}_{I}/\partial y_{I}=\mu\kappa y_{I}^{\kappa-1}=\mu$ results in

[TABLE]

Since $\mu>0$ was chosen arbitrary, we find that this holds true for all $\mu>0$. From Lemma 2.2 we can conclude that

[TABLE]

hold true.

The above equations contain products of Patankar-weights. To find conditions for the PWDs, we determine the limits of the Patankar-weights. In this regard, equation (32) shows $y_{I}^{n+1}/\sigma_{I}\to 1$. Substitution of this into (33) and (34) yields

[TABLE]

and

[TABLE]

Next, we show that none of the Patankar-weights ${y_{I}^{(2)}}/{\pi_{I}}$ and ${y_{I}^{(3)}}/{\rho_{I}}$ can tend to infinity. To do so, we must consider two cases. If $b_{2}=0$, (22a) becomes $b_{3}(a_{31}+a_{32})=1/2$, so we can conclude $y_{I}^{(3)}/\rho_{I}\to 1$ from (35) and thus, $y_{I}^{(2)}/\pi_{I}\to 1$ from (36). If $b_{2}>0$, both terms on the left hand side of (35) are positive, since $a_{21},a_{32},b_{3}>0$ and $a_{31}\geq 0$, hence, $y_{I}^{(2)}/\pi_{I}\not\to\infty$ and ${y_{I}^{(3)}}/{\rho_{I}}\not\to\infty$. Consequently, owing to (36), we find that none of the Patankar-weights ${y_{I}^{(2)}}/{\pi_{I}}$ or ${y_{I}^{(3)}}/{\rho_{I}}$ can tend to zero, as this would require the other weight to tend to infinity. Denoting by $\Gamma_{I}^{(2)}$ and $\Gamma_{I}^{(3)}$ the limits of $y_{I}^{(2)}/\pi_{I}$ and $y_{I}^{(3)}/\rho_{I}$, we have

[TABLE]

with $\Gamma_{I}^{(2)},\Gamma_{I}^{(3)}>0$ and

[TABLE]

Now we consider the case $\kappa=2$ in (29). From (2) and (37) we see

[TABLE]

which implies

[TABLE]

Together with (30) we find

[TABLE]

and

[TABLE]

Substituting this and (32) into (29) yields

[TABLE]

Owing to (33) and (34), we find

[TABLE]

which shows

[TABLE]

due to (32) and (37). Utilization of (37) results in an additional equation

[TABLE]

containing $\Gamma_{I}^{(2)}$ and $\Gamma_{I}^{(3)}$. To determine the values $\Gamma_{I}^{(2)}$ and $\Gamma_{I}^{(3)}$ we can use Lemma 2.5, which shows that the system (38), (39) and (40) has the unique solution

[TABLE]

Together with (32) and (34), this leads to

[TABLE]

Substituting this into (30) and (2), we find $y_{I}^{(2)}=y_{I}^{n}+\mathcal{O}(\Delta t)$ and $y_{I}^{(3)}=y_{I}^{n}+\mathcal{O}(\Delta t)$. Hence, we obtain

[TABLE]

from which we conclude

[TABLE]

In an analogous manner, we can conclude

[TABLE]

from (41). Since $I$ was chosen arbitrary, we can let it run from 1 to $N$ and find that (21a) and (21b) are necessary conditions.

The same is true for the other equations derived above. In particular,

[TABLE]

hold true. These equations are helpful to find a concise representation of (33). On account of (20b) and (45), we have

[TABLE]

for $i=1,\dots,N$. Similar, (20f) and (46) show

[TABLE]

for $i=1,\dots,N$. According to (47), we have $\mathbf{y}^{(2)}-\mathbf{y}^{n}=\mathcal{O}(\Delta t)$, and from (23) we see

[TABLE]

Hence,

[TABLE]

follows from (48). Substituting this into (33), and taking account of (43) and (44), results in

[TABLE]

Again, $I$ was chosen arbitrary, and letting it run from $1$ to $N$, shows that condition (21c) is necessary. Finally, analogous to (42), (32) and (25) show

[TABLE]

By letting $I$ run from $1$ to $N$, we see that also condition (21d) is necessary.

Now we show that the conditions (21) are sufficient, to make (20) a third order MPRK scheme. We start our investigation with the choice $\delta=1$. The MPRK scheme (20) can be written in the form of three linear systems

[TABLE]

Since $\delta=1$, utilizing Lemma 2.1 yields $(\mathbf{M}^{(2)})^{-1}=\mathcal{O}(1)$, $(\mathbf{M}^{(3)})^{-1}=\mathcal{O}(1)$ and $\mathbf{M}^{-1}=\mathcal{O}(1)$. Thus, we can conclude $\mathbf{y}^{(2)}=\mathcal{O}(1)$, $\mathbf{y}^{(3)}=\mathcal{O}(1)$ and $\mathbf{y}^{n+1}=\mathcal{O}(1)$. Together with conditions (21a), (21b), and (21d), this results in

[TABLE]

for $i=1,\dots,N$, since $y_{i}^{n}>0$. The boundedness of the Patankar-weights (49) shows that (20b) yields

[TABLE]

for $i=1,\dots,N$. This allows us to use (23) to expand $p_{ij}(\mathbf{y}^{(2)})$ and $d_{ij}(\mathbf{y}^{(2)})$ in the form

[TABLE]

for $i,j=1,\dots,N$. Substituting this into (20f), and taking account of (50), we find

[TABLE]

as well. Now, we show that (52) and (54) are valid for $\delta=0$ as well. In this case, owing to (20b), we have

[TABLE]

for $i=1,\dots,N$. Since $1/\pi_{i}=\mathcal{O}(1)$ according to (21a), we can conclude $y_{i}^{(2)}=\mathcal{O}(1)$ and thus $y_{i}^{(2)}/\pi_{i}=\mathcal{O}(1)$ for $i=1,\dots,N$. Utilizing this in (20b) we find

[TABLE]

as before in (52). Hence, (53) is holds true for $\delta=0$ as well. This, together with (21b), shows

[TABLE]

and in addition $y_{i}^{(3)}/\rho_{i}=\mathcal{O}(1)$ for $i=1,\dots,N$. Consequently, we have

[TABLE]

for $i=1,\dots,N$, as in (54). The remaining part of the proof is independent of the value of $\delta$.

Owing to (52) and (54), (23) shows

[TABLE]

for $i,j=1,\dots,N$ and $k=2,3$. Utilizing this and (51), the solution on the next time level (20i) satisfies

[TABLE]

for $i=1,\dots,N$. Hence, we can conclude

[TABLE]

from (21d). Inserting this and (55) into (20i) shows

[TABLE]

for $i=1,\dots,N$, since $b_{1}+b_{2}+b_{3}=1$ according to (12c). Now, we can conclude

[TABLE]

from (21d). Introducing this relation and (55) into (20i) yields

[TABLE]

for $i=1,\dots,N$. From (52) and (21a) we can conclude

[TABLE]

thus, (20b) shows

[TABLE]

for $i=1,\dots,N$. Similar, (50) and (20f) imply

[TABLE]

for $i=1,\dots,N$. Thus, inserting this and (55) into (20f) shows

[TABLE]

for $i=1,\dots,N$. Finally, substitution of (57) and (59) into (56) results in

[TABLE]

since $b_{2}a_{21}+b_{3}(a_{31}+a_{32})=1/2$ due to (12). Hence, we even have

[TABLE]

for $i=1,\dots,N$.

This enables the proof of the third order accuracy of the MPRK scheme. Substitution of (60) and (55) into (20i) yields

[TABLE]

Taking account of (57) and (59), and using $b_{2}a_{21}^{2}+b_{3}(a_{31}+a_{32})^{2}=1/3$, this can be written in the form

[TABLE]

It remains to expand $b_{2}(\mathbf{y}^{(2)}-\mathbf{y}^{n})+b_{3}(\mathbf{y}^{(3)}-\mathbf{y}^{n})$ up to $\mathcal{O}(\Delta t^{3})$. Therefore, we use (20f) and (55) to see

[TABLE]

for $i=1,\dots,N$. Insertion of (57) and (58) shows

[TABLE]

for $i=1,\dots,N$. Utilization of (20b) results in

[TABLE]

for $i=1,\dots,N$, since $b_{3}a_{21}a_{32}=1/6$ owing to (12f). Substitution of (57) and (59) into (62) in combination with (21c) yields

[TABLE]

for $i=1,\dots,N$. Inserting this into (61) results in

[TABLE]

for $i=1,\dots,N$. A comparison with (24) completes the proof. ∎

The following theorem defines a family of third order MPRK schemes. It is based on the idea to use a second order MPRK22($\alpha$) scheme of [KM17] to compute the PWDs $\sigma_{i}$, as condition (21d) of Theorem 2.6 shows that $\sigma_{i}$ must be a second order approximation of $y_{i}(t^{n+1})$ for $i=1,\dots,N$.

Theorem 2.7**.**

Given an explicit three-stage third order Runge-Kutta scheme with non-negative weights, the MPRK scheme (20) is of third order, if we choose

[TABLE]

Proof.

We need to verify that the choice of PWDs (63) satisfies conditions (21) of Theorem 2.6. Therefore, we make repeatedly use of the statements of the proof of Theorem 2.6. We also use Newton’s generalized binomial theorem111The theorem can be deduced from the binomial series $\sum_{k=0}^{\infty}\binom{s}{k}x^{k}=(1+x)^{s}$, which is convergent for $\lvert x\rvert<1$. See for instance [How01]., which states that

[TABLE]

with

[TABLE]

holds true for $s\in\mathbb{R}$, if $x>0$ and $\lvert y/x\rvert<1$. The theorem implies

[TABLE]

and

[TABLE]

for $s,\eta\in\mathbb{R}$, since $y_{i}^{n}>0$.

First, we note that condition (21a) is clearly satisfied by (63a). This allows us to conclude

[TABLE]

for $i=1,\dots,N$, along the same lines as in (57) of Theorem 2.6. Introducing this into (63b) shows

[TABLE]

for $i=1,\dots,N$. Thus, condition (21b) holds true as well.

Next, we verify condition (21c). From (66) and (67) we find

[TABLE]

for $i=1,\dots,N$. Defining $f(\Delta t)=1/(\xi+\Delta t\eta)$ for some constants $\xi$ and $\eta$, we can conclude $f(\Delta t)=f(0)+f^{\prime}(0)\Delta t+\mathcal{O}(\Delta t^{2})=1/\xi-\eta/\xi^{2}+\mathcal{O}(\Delta t^{2})$, and hence,

[TABLE]

for $i=1,\dots,N$. Consequently,

[TABLE]

for $i=1,\dots,N$. Substituting this and (63a) into condition (21c) shows

[TABLE]

for $i=1,\dots,N$. Hence, condition (21c) holds true. Finally, (20b) and (63f) with PWDs (63c) form the MPRK22($a_{21}$) scheme of [KM17]. As this is a second order scheme, condition (21d) is satisfied as well. ∎

The family of schemes introduced in Theorem 2.7 can be written in the form

[TABLE]

with $p=3a_{21}(a_{31}+a_{32})b_{3}$, $q=a_{21}$, $\beta_{2}=1/(2a_{21})$ and $\beta_{1}=1-\beta_{2}$ for $i=1,\dots,N$. We denote the members of this family, which derive from case I in Lemma 2.4, by MPRK43I($\alpha$,$\beta$) if $\delta=1$ and by MPRK43Incs($\alpha$,$\beta$) if $\delta=0$. If the method comes from case II in Lemma 2.4, we denote it by MPRK43II($\gamma$) if $\delta=1$ and by MPRK43IIncs($\gamma$) if $\delta=0$.

To our knowledge, this is the first time that third order MPRK schemes are presented. A third order Patankar type scheme based on a BDF method was presented in [FS11].

As the schemes (68) incorporate the MPRK22($a_{21}$) scheme, we must restrict $a_{21}$ to $a_{21}\geq 1/2$. Hence, the permissible Runge-Kutta parameters are given by the Butcher tableaus of Lemma 2.4 with the additional restriction $\alpha\geq 1/2$ in case I.

The MPRK scheme (68) can be understood as a four stage MPRK scheme with corresponding Butcher tableau

[TABLE]

The extra stage to compute the PWDs $\sigma_{i}$ requires no additional function evaluations, but nevertheless an additional linear system needs to be solved. It is a future concern to prove or disprove, whether the construction of a third order three stage MPRK scheme is possible.

In the numerical experiments of Section 4 we consider six specific MPRK43 schemes. Choosing $\alpha=1$ and $\beta=1/2$ in case I of Lemma 2.4 yields the Butcher tableau

[TABLE]

The corresponding scheme MPRK43I($1$,$1/2$) is the only MPRK43 scheme with $p=q=1$ and it uses the original MPRK22(1) scheme from [BDM03], which is based on Heun’s method, to compute the PWDs $\sigma_{i}$. Setting $\alpha=1/2$ and $\beta=3/4$ in case I of Lemma 2.4, we obtain the Butcher tableau

[TABLE]

and the method MPRK43I($1/2$,$3/4$) with $p=q=1/2$. This method uses the MPRK22($1/2$) scheme, which is adapted from the midpoint method, to compute the PWDs $\sigma_{i}$. Case II of Lemma 2.4 with $\gamma=1/2$ provides the Butcher tableau

[TABLE]

The associated MPRK scheme MPRK43II($2/3$) with $p=q=2/3$. It employs Ralston’s method MPRK22(2/3) to calculate the PWDs $\sigma_{i}$. Besides the above three schemes, we also use their counterparts with $\delta=0$ in the numerical experiments of Section 4.

Of course, many other schemes are members of the family (68). Here we made the restriction $p=q$ to allow for the same PWDs in (68c) and (68d). We also selected schemes, which use MPRK22($\alpha$) schemes to compute the PWDs $\sigma_{i}$, that were investigated in [KM17].

The numerical experiments, presented in Section 4, will confirm the third order accuracy of the MPRK43 schemes. Additionally, numerical solutions of the Robertson problems will show that these schemes have the ability to integrate stiff PDS.

3 Test problems

For our numerical experiments, we consider the same test cases as in [KM17]. A simple linear test problem for which the analytical solution is known, two non-stiff nonlinear test problems and the stiff Robertson problem.

Linear test problem

The simple linear test case is given by

[TABLE]

with a constant parameter $a$ and initial values $y_{1}(0)=y_{1}^{0}$ and $y_{2}(0)=y_{2}^{0}$. We can write the right hand side in the form (2) with

[TABLE]

and $p_{ii}(\mathbf{y})=d_{ii}(\mathbf{y})=0$ for $i=1,2$. The system describes exchange of mass between to constituents. The analytical solution is

[TABLE]

with the asymptotic solution

[TABLE]

The system is conservative and we get

[TABLE]

In the numerical simulations of Section 4 we use $a=5$ and initial values $y_{1}^{0}=0.9$ and $y_{2}^{0}=0.1$. The solution is approximated on the time interval $[0,1.75]$.

Nonlinear test problem

The non-stiff nonlinear test problem reads

[TABLE]

with initial conditions $y_{i}(0)=y_{i}^{0}$ for $i=1,2,3$. To express the right hand side in the form (2) we can use

[TABLE]

and $p_{ij}(\mathbf{y})=d_{ij}(\mathbf{y})=0$ for all other combinations of $i$ and $j$.

The system represents a biogeochemical model for the description of an algal bloom, that transforms nutrients ($y_{1}$) via phytoplankton ($y_{2}$) into detritus ($y_{3}$). In the numerical simulations of Section 4 we use the initial conditions $y_{1}^{0}=9.98$, $y_{2}^{0}=0.01$ and $y_{3}^{0}=0.01$. The solution is approximated on the time interval $[0,30]$.

Original Brusselator test problem

As another non-stiff nonlinear test case we consider the original Brusselator problem [LN71, HNW93]

[TABLE]

with constant parameters $k_{i}$ and initial values $y_{i}(0)=y_{i}^{0}$ for $i=1,\dots,6$. The system can be written in the form (2), setting

[TABLE]

and $p_{ij}(\mathbf{y})=d_{ji}(\mathbf{y})=0$ for all other combinations of $i$ and $j$.

In the numerical simulations of Section 4 we set $k_{i}=1$ for $i=1,\dots,6$ and the initial values $y_{1}(0)=y_{2}(0)=10$, $y_{3}(0)=y_{4}(0)=\mathtt{eps}\approx 2.2204\cdot 10^{-16}$, and $y_{5}(0)=y_{6}(0)=0.1$. The time interval of interest is $[0,6]$.

Robertson test problem

To demonstrate the practicability of MPRK schemes in the case of stiff systems, we apply the schemes to the Robertson test case, which is given by

[TABLE]

with initial values $y_{i}(0)=y_{i}^{0}$ for $i=1,2,3$. For this problem the production and destruction rates (2) are given by

[TABLE]

and $p_{ij}(\mathbf{y})=d_{ij}(\mathbf{y})=0$ for all other combinations of $i$ and $j$.

We use the initial values $y_{1}(0)=y_{1}^{0}=1-2\mathtt{eps}$ and $y_{2}^{0}=y_{3}^{0}=\mathtt{eps}\approx 2.2204\cdot 10^{-16}$ in the numerical simulations of Section 4.

In this problem the reactions take place on very different time scales, the time interval of interest is $[10^{-6},10^{10}]$. Therefore, a constant time step size is not appropriate. In the numerical simulations we use $\Delta t_{i}=4^{i-1}\Delta t_{0}$ with $\Delta t_{0}=10^{-6}$ in the $i$th time step. The small initial time step size $\Delta t_{0}$ is chosen to obtain an adequate resolution of $y_{2}$.

4 Numerical results

In this section, we confirm the theoretical convergence order of the novel MPRK43 schemes. We also show numerical approximations of MPRK43 schemes applied to the stiff Robertson problem (72), the nonlinear test problem (70) and the Brusselator problem (71).

To visualize the order of the MPRK schemes we use a relative error $E$ taken over all time steps and all constituents:

[TABLE]

where $M$ denotes the number of executed time steps. To compute the error $E$ we need to know the analytic solution, which is known for the linear test case, but not for the other test problems. Hence, we computed a reference solution, using the Matlab functions ode45 for the non-stiff nonlinear problems and ode23s for the Robertson problem. In both cases we utilized the tolerances $\mathtt{AbsTol}=\mathtt{RelTol}=10^{-10}$.

Convergence order

Figure 3 shows error plots of six MPRK43 schemes applied to the linear test problem (69), the nonlinear test problem (70) and the Brusselator (71). In all cases the third order accuracy is confirmed. Moreover, Figure 3(a) shows that MPRK43I(1,1/2) is less accurate than MPRK43IIncs(1,1/2) and MPRK43I(1/2,3/4) is more accurate than MPRK43Incs(1/2), when applied to the linear test problem. Hence, we cannot make a general statement, whether to choose $\delta=0$ or $\delta=1$ in the MPRK43 schemes.

Stiff problems

Figure 4 shows numerical approximations of six MPRK43 schemes applied to the stiff Robertson problem (72). As mentioned, the time step size in the $k$th time step was chosen as $\Delta t_{k}=4^{k-1}\Delta t_{0}$ with initial time step size $\Delta t_{0}=10^{-6}$. Hence, only 29 time steps are necessary to traverse the time interval $[10^{-6},10^{10}]$. The small initial time step was chosen to obtain an adequate resolution of the component $y_{2}$ in the starting phase. To visualize the evolution of $y_{2}$, it is multiplied by $10^{4}$.

All six schemes generate adequate solutions. For this test problem, the variants with conservative stage values (left column) can be seen to be more accurate than those with non-conservative stage values (right column). But the overall accuracy is excellent with regard to the fairly large time steps in use.

In [KM17] we reported that some MPRK22ncs schemes generate oscillations, when applied to the Robertson problem. In case of the MPRK43 schemes, we did not encounter any issues.

Like in [KM17], we additionally show numerical solutions of the six MPRK43 schemes applied to the nonlinear test problem (70) and the Brusselator (71) in Figures 5 and 6.

5 Summary and Outlook

In this paper we have extended the work of [KM17] to third order by deriving necessary and sufficient conditions for three-stage third order schemes. We also introduced the MPRK43I and MPRK43II schemes, which, to our knowledge, are the first third order Patankar-type schemes presented in literature. These schemes can be regarded as four-stage third order MPRK schemes and it is a future research topic, to investigate the construction of three-stage third order MPRK schemes. In addition, the search for other possible PWDs is of interest as well.

The numerical experiments have shown that the MPRK43 schemes are capable of integrating stiff ODEs, such as the Robertson problem. However, in absence of a thorough analysis of truncation errors and stability, we cannot make statements which schemes of the MPRK43 family are preferable. This is a future research topic as well.