Super-Convergent Implicit-Explicit Peer Methods with Variable Step Sizes

TL;DR

This paper develops super-convergent IMEX-Peer methods that maintain high order accuracy with variable step sizes and good stability, enhancing the efficient numerical solution of multi-scale dynamical systems.

Contribution

The paper introduces new super-convergent IMEX-Peer methods of orders 3, 4, and 5, with conditions ensuring high order accuracy under variable step sizes.

Findings

New IMEX-Peer methods of order s+1 for s=2,3,4 stages.

Methods maintain super-convergence with variable step sizes.

Numerical experiments demonstrate improved performance.

Abstract

Dynamical systems with sub-processes evolving on many different time scales are ubiquitous in applications. Their efficient solution is greatly enhanced by automatic time step variation. This paper is concerned with the theory, construction and application of IMEX-Peer methods that are super-convergent for variable step sizes and A-stable in the implicit part. IMEX schemes combine the necessary stability of implicit and low computational costs of explicit methods to efficiently solve systems of ordinary differential equations with both stiff and non-stiff parts included in the source term. To construct super-convergent IMEX-Peer methods which keep their higher order for variable step sizes and exhibit favourable linear stability properties, we derive necessary and sufficient conditions on the nodes and coefficient matrices and apply an extrapolation approach based on already computed stage values. New super-convergent IMEX-Peer methods of order $s+1$ for $s=2,3,4$ stages are given as result of additional order conditions which maintain the super-convergence property independent of step size changes. Numerical experiments and a comparison to other super-convergent IMEX-Peer methods show the potential of the new methods when applied with local error control.

Figures

Tables

Table 1. Table 1: Size of stability regions 𝕊 90 ∘ subscript 𝕊 superscript 90 \mathbb{S}_{90^{\circ}} and 𝕊 0 ∘ subscript 𝕊 superscript 0 \mathbb{S}_{0^{\circ}} , x m ​ a ​ x ​ ( 𝕊 90 ∘ ) subscript 𝑥 𝑚 𝑎 𝑥 subscript 𝕊 superscript 90 x_{max}(\mathbb{S}_{90^{\circ}}) at the negative real axis, y m ​ a ​ x ​ ( 𝕊 0 ∘ ) subscript 𝑦 𝑚 𝑎 𝑥 subscript 𝕊 superscript 0 y_{max}(\mathbb{S}_{0^{\circ}}) at the positive imaginary axis, spectral radius of R − 1 ​ Q superscript 𝑅 1 𝑄 R^{-1}Q , and error constants c i ​ m = | d s + 1 | subscript 𝑐 𝑖 𝑚 subscript 𝑑 𝑠 1 c_{im}=|d_{s+1}| and c e ​ x = | R ​ l s | subscript 𝑐 𝑒 𝑥 𝑅 subscript 𝑙 𝑠 c_{ex}=|R\,l_{s}| for super-convergent IMEX-Peer methods, including those from [ 16 ] .

IMEX-	$\left|{𝕊}_{{90}^{\circ }}\right|$	$x_{max}$	$\left|{𝕊}_{0^{\circ }}\right|$	$y_{max}$	$\rho \left(R^{-1}Q\right)$	$c_{im}$	$c_{ex}$
Peer2s	$2.15$	$-1.41$	$4.47$	$1.21$	$0.128$	${\mathrm{2.37\hspace{0.17em}10}}^{-1}$	${\mathrm{3.23\hspace{0.17em}10}}^{-1}$
Peer2sve	${\mathrm{6.68\hspace{0.17em}10}}^{-5}$	$-{\mathrm{5.68\hspace{0.17em}10}}^{-3}$	$0.14$	$0.36$	$0.863$	${\mathrm{1.94\hspace{0.17em}10}}^{-1}$	${\mathrm{2.83\hspace{0.17em}10}}^{-1}$
Peer3s	$2.67$	$-1.58$	$6.11$	$1.69$	$0.552$	${\mathrm{1.24\hspace{0.17em}10}}^{-1}$	${\mathrm{1.68\hspace{0.17em}10}}^{-1}$
Peer3sv	$0.11$	$-0.25$	$0.55$	$0.43$	$0.254$	${\mathrm{2.29\hspace{0.17em}10}}^{-1}$	${\mathrm{1.43\hspace{0.17em}10}}^{-1}$
Peer4s	$1.07$	$-1.45$	$4.39$	$1.00$	$0.542$	${\mathrm{6.42\hspace{0.17em}10}}^{-2}$	${\mathrm{1.17\hspace{0.17em}10}}^{-1}$
Peer4sve	$1.66$	$-1.68$	$3.11$	$0.92$	$0.118$	${\mathrm{2.02\hspace{0.17em}10}}^{-2}$	${\mathrm{3.37\hspace{0.17em}10}}^{-2}$
Peer4sv	${\mathrm{1.34\hspace{0.17em}10}}^{-3}$	$-{\mathrm{4.05\hspace{0.17em}10}}^{-2}$	$0.63$	$0.67$	$0.632$	${\mathrm{7.47\hspace{0.17em}10}}^{-2}$	${\mathrm{6.75\hspace{0.17em}10}}^{-2}$

Table 2. Table 2: Coefficients of IMEX-Peer3sv and IMEX-Peer4sv which are super-convergent for variable step sizes. Here, E 2 = ( e i ​ j ) subscript 𝐸 2 subscript 𝑒 𝑖 𝑗 E_{2}=(e_{ij}) .

IMEX-Peer3sv, $s=3$
$c_1$	$0.000000000000000$			$p_{11}$	$1.000000000000000$
$c_2$	$0.500000000000000$			$p_{12}$	$0.000000000000000$
$c_3$	$1.000000000000000$			$p_{13}$	$0.000000000000000$
$\gamma$	$0.690969692535085$			$p_{21}$	$1.009534846612963$
$r_{21}$	$0.351562922857064$			$p_{22}$	$-0.000125189884283$
$r_{31}$	$0.346024253990984$			$p_{23}$	$-0.009409656728680$
$r_{32}$	$0.328884660689640$			$p_{31}$	$0.927244072163109$
$e_{21}$	$1.454929231059714$			$p_{32}$	$-0.000247968521087$
$e_{31}$	$-6.099201725139450$			$p_{33}$	$0.073003896357977$
$e_{32}$	$3.157746208382228$
IMEX-Peer4sv, $s=4$
$c_1$	$0.000000000000000$			$p_{11}$	$1.000000000000000$
$c_2$	$-1.598239239549169$			$p_{12}$	$0.000000000000000$
$c_3$	$0.523829503832339$			$p_{13}$	$0.000000000000000$
$c_4$	$1.000000000000000$			$p_{14}$	$0.000000000000000$
$\gamma$	$0.681884472048995$			$p_{21}$	$1.000204745561481$
$r_{21}$	$1.292744499701930$			$p_{22}$	$-0.000195233457439$
$r_{31}$	$1.074957286644128$			$p_{23}$	$-0.000009518220959$
$r_{32}$	$-0.054028162784565$			$p_{24}$	$0.000000006116916$
$r_{41}$	$4.064480810437903$			$p_{31}$	$1.169763235411655$
$r_{42}$	$1.031994574173631$			$p_{32}$	$-0.169740581681421$
$r_{43}$	$-0.534558192336057$			$p_{33}$	$-0.000025123517333$
$e_{21}$	$-0.153830152235951$			$p_{34}$	$0.000002469787099$
$e_{31}$	$0.065444441626366$			$p_{41}$	$1.915153835547942$
$e_{32}$	$-0.976514386415223$			$p_{42}$	$-0.244331567248295$
$e_{41}$	$-0.234155732816782$			$p_{43}$	$-0.671042624270695$
$e_{42}$	$-2.535629358626096$			$p_{44}$	$0.000220355971049$
$e_{43}$	$1.477107513945526$

Table 3. Table 3: Coefficients of IMEX-Peer4sve which is optimally zero-stable, super-convergent for variable step sizes in the explicit part and for constant step sizes in the implicit part. Here, E 2 = ( e i ​ j ) subscript 𝐸 2 subscript 𝑒 𝑖 𝑗 E_{2}=(e_{ij}) .

IMEX-Peer4sve, $s=4$, optimally zero-stable
$c_1$	$-0.868838855210029$			$p_{11}$	$0.000000000000000$
$c_2$	$-0.253884413463736$			$p_{12}$	$0.316402904545681$
$c_3$	$0.754504864110948$			$p_{13}$	$1.127642509582261$
$c_4$	$1.000000000000000$			$p_{14}$	$-0.444045414127942$
$\gamma$	$0.473861788489939$			$p_{21}$	$0.000000000000000$
$r_{21}$	$0.732961380396538$			$p_{22}$	$0.000000000000000$
$r_{31}$	$-2.472299983846101$			$p_{23}$	$-0.017465269321373$
$r_{32}$	$0.077358285702625$			$p_{24}$	$1.017465269321373$
$r_{41}$	$-1.603925020256191$			$p_{31}$	$0.000000000000000$
$r_{42}$	$-2.797576519478004$			$p_{32}$	$0.000000000000000$
$r_{43}$	$-0.278164642408456$			$p_{33}$	$0.000000000000000$
$e_{21}$	$-0.183287385063759$			$p_{34}$	$1.000000000000000$
$e_{31}$	$5.974911797174020$			$p_{41}$	$0.000000000000000$
$e_{32}$	$-2.556627399170977$			$p_{42}$	$0.000000000000000$
$e_{41}$	$2.456065798975378$			$p_{43}$	$0.000000000000000$
$e_{42}$	$-2.032396276261657$			$p_{44}$	$1.000000000000000$
$e_{43}$	$1.255044479285407$

Equations

$$u^{\prime }(t)=F_0(u(t))+F_1(u(t)),$$

$$u^{\prime }(t)=F(u(t)),u(0)=u_0\in \mathbb{V}.$$

$$w_n=(P_n\otimes I)w_{n-1}+{△t}_n(Q_n\otimes I)F(w_{n-1})+{△t}_n(R_n\otimes I)F(w_n)$$

$$w_n=[w_{n,1}^T,\dots ,w_{n,s}^T{]}^T\in {\mathbb{V}}^s,w_{n,i}\approx u(t_n+c_i{△t}_n).$$

$$w_n=P_n{w}_{n-1}+{△t}_n{Q}_{n}F(w_{n-1})+{△t}_n{R}_{n}F(w_n).$$

$$r_n=w(t_n)-P_{n}w(t_{n-1})-{△t}_n{Q}_n{w}^{\prime }(t_{n-1})-{△t}_n{R}_n{w}^{\prime }(t_n).$$

$$w_i(t_{n-1})=u\left(t_n+\frac{c_i-1}{{\sigma }_n}{△t}_n\right),i=1,\dots ,s,$$

$$w(t_n)=$$

$$w(t_{n-1})=$$

$$r_n=\sum _{j\ge 1}{△t}_n^j{d}_{n,j}\otimes u^{(j)}(t_n)$$

$$d_{n,j}=\frac{1}{j!}\left(c^j-\frac{1}{{\sigma }_n^j}{P}_n(c-e{)}^j-\frac{j}{{\sigma }_n^{j-1}}{Q}_n(c-e{)}^{j-1}-j{R}_n{c}^{j-1}\right).$$

$$V_{0}=\big{(}c_{i}^{j-1}\big{)},\qquad V_{1}=\big{(}(c_{i}-1)^{j-1}\big{)},\qquad i,j=1,\ldots,s,$$

$$C{V}_0-\frac{1}{{\sigma }_n}{P}_n(C-I)V_1{S}_n^{-1}-Q_n{V}_{1}D{S}_n^{-1}-R_n{V}_{0}D=0.$$

$$Q_n=\left((C{V}_0-R_n{V}_{0}D)S_n-\frac{1}{{\sigma }_n}{P}_n(C-I)V_1\right)(V_{1}D{)}^{-1},$$

$$w_n=P{w}_{n-1}+{△t}_n{Q}_{n}F(w_{n-1})+{△t}_{n}RF(w_n)$$

$$Q_n=\left((C{V}_0-R{V}_{0}D)S_n-\frac{1}{{\sigma }_n}P(C-I)V_1\right)(V_{1}D{)}^{-1}.$$

$$w_n=(I-z_{n}R{)}^{-1}(P+z_n{Q}_n)w_{n-1}=:M_{im}(z_n,{\sigma }_n)w_{n-1}$$

$$w_n=M_{im}(z_n,{\sigma }_n)M_{im}(z_{n-1},{\sigma }_{n-1})\cdots M_{im}(z_1,{\sigma }_1)w_0.$$

$$\Vert {\epsilon }_n\Vert \le$$

$$\times \underset{0\le t\le t_n}{\max }\Vert u^{(s+1)}(t)\Vert +\mathcal{O}\left({△t}_{max}^{s+1}\right).$$

$${△t}_1^{s+1}\Vert d_{1,s+1}{\Vert }_{\infty }+\dots +{△t}_n^{s+1}\Vert d_{n,s+1}{\Vert }_{\infty }\le (t_{n+1}-t_1){△t}_{max}^s\underset{i=1,\dots ,n}{\max }\Vert d_{i,s+1}{\Vert }_{\infty }.$$

$$\bar{w}(t_i):=w(t_i)-{△t}_i^{s+1}{v}_i\otimes u^{(s+1)}(t_i).$$

$$\begin{array}[]{rll}\bar{r}_{i}&=&\bar{w}(t_{i})-P\bar{w}(t_{i-1})-{\triangle t}_{i}Q_{i}F(\bar{w}(t_{i-1}))-{\triangle t}_{i}RF(\bar{w}(t_{i}))\\[5.69054pt] &=&r_{i}-{\triangle t}_{i}^{s+1}d_{i,s+1}\otimes u^{(s+1)}(t_{i})-T(v_{i-1},v_{i})\otimes u^{(s+1)}(t_{i})+{\cal O}({\triangle t}_{i}^{s+2})\,,\end{array}$$

$$T(v_{i-1},v_i)={△t}_i^{s+1}P{v}_i-{△t}_{i-1}^{s+1}P{v}_{i-1}.$$

$${\sigma }_i^{-j}-{\sigma }_{i-1}^{-j}=\mathcal{O}({△t}_{max})\text{for all}j\ge 1.$$

$$\begin{array}[]{rll}T(v_{i-1},v_{i})&\!\!\!=\!\!\!&{\triangle t}_{i}^{s+1}P(v_{i}-v_{i-1})+{\cal O}({\triangle t}_{max}){\triangle t}_{i}^{s+1}\\[5.69054pt] &\!\!\!=\!\!\!&{\triangle t}_{i}^{s+1}P(I-P)^{+}(d_{i,s+1}-d_{i-1,s+1})+{\cal O}({\triangle t}_{max}){\triangle t}_{i}^{s+1}={\cal O}({\triangle t}_{max}){\triangle t}_{i}^{s+1}\,,\end{array}$$

$$d_{i,s+1}\in \text{range}(I-P)\text{iff}{v}^T{d}_{i,s+1}=0\text{for all }i=1,\dots ,n.$$

$$d_{s+1}(\sigma )=\frac{1}{(s+1)!}\left(c^{s+1}-\frac{1}{{\sigma }^{s+1}}P(c-e{)}^{s+1}-\frac{s+1}{{\sigma }^s}Q(\sigma )(c-e{)}^s-(s+1)R{c}^s\right),$$

$$v^T{d}_{s+1}(\sigma )=h_0+\sum _{j=1}^s{\tilde{v}}_{s+1-j}^T{\tilde{c}}_{s+1-j}{\sigma }^{-j}+h_{s+1}{\sigma }^{-(s+1)}$$

$$h_0=$$

$${\tilde{v}}^T=$$

$$\tilde{c}=$$

Full text

*Dedicated to the Memory of Willem Hundsdorfer (1954 - 2017)

Keywords: implicit-explicit (IMEX) Peer methods; super-convergence; extrapolation; A-stability; variable step size; local error control

1 Introduction

Many dynamical processes in engineering, physics, chemistry and other areas are modelled by large systems of ordinary differential equations (ODEs) of the form

[TABLE]

where $F_{0}:\mathbb{R}^{m}\rightarrow\mathbb{R}^{m}$ represents the non-stiff or mildly stiff part and $F_{1}:\mathbb{R}^{m}\rightarrow\mathbb{R}^{m}$ gives the stiff part of the equation. Such problems often result from semi-discretized systems of partial differential equations with diffusion, advection and reaction terms. Instead of applying a single explicit or implicit method, an often more appropriate and efficient approach is to use the decomposition of the right-hand side by treating only the $F_{1}$ contribution in an implicit fashion. Thus, favourable stability properties of implicit schemes and the advantage of lower costs for explicit schemes are combined to enhance the overall computational efficiency. Since dynamical systems typically have sub-processes evolving on many different time scales, a good ODE integrator should come with some adaptive error control, making frequent step size changes over its own progress. In smooth regions, a few large steps should speed up the integration, whereas many small steps should be applied in non-smooth terrains. The resulting gains in efficiency can be up to factors of hundreds or more.

IMEX-Peer methods with variable step sizes have been successfully applied by Soleimani, Knoth, and Weiner [17] to fast-wave-slow-wave problems arising in weather prediction. The super-convergent IMEX-Peer methods recently developed by Soleimani and Weiner [18, 19] and Schneider, Lang, and Hundsdorfer [16] can in principle be applied with variable step sizes, but then they might lose their super-convergence property, especially for serious step size changes. Super-convergent explicit Peer methods for variable step sizes have first been constructed by Weiner, Schmitt, Podhaisky, and Jebens [21], exploiting special matrix structures. Another approach to construct such methods is the use of extrapolation as proposed by Schneider, Lang, and Hundsdorfer [16]. This idea goes back to Crouzeix [7] and was also used by Cardone, Jackiewicz, Sandu, and Zhang [5, 6] and later on by Braś, Izzo, and Jackiewicz [2] to construct implicit-explicit general linear and Runge-Kutta methods. The procedure can be easily extended to variable step sizes for IMEX-Peer methods.

In this paper, we use the extrapolation approach to construct new super-convergent IMEX-Peer methods that keep their higher order for variable step sizes and exhibit favourable linear stability properties, including A-stability of the implicit part. Additional order conditions on the nodes and coefficient matrices which maintain the super-convergence property independent of step size changes are derived for implicit, explicit and IMEX-Peer methods. We give formulas for new super-convergent IMEX-Peer methods of order $s+1$ for $s=2,3,4$ stages. Stability regions are computed and compared to those of super-convergent IMEX-Peer methods for constant step sizes from Schneider, Lang, and Hundsdorfer [16]. Eventually, numerical results are presented for a Prothero-Robinson problem, the van der Pol oscillator, a one-dimensional Burgers equation with stiff diffusion and a one-dimensional advection-reaction problem with stiff reactions.

2 Implicit-Explicit Peer Methods with Variable Step Sizes

2.1 Super-convergent implicit Peer methods with variable

step sizes

We apply the so-called Peer methods introduced by Schmitt, Weiner and co-workers [14, 15, 18] to solve initial value problems in the vector space $\mathbb{V}=\mathbb{R}^{m},m\geq 1$,

[TABLE]

The general form of an $s$-stage implicit Peer method with variable step sizes ${\triangle t}_{n}$ is

[TABLE]

with the $m\times m$ identity matrix $I$, the $s\times s$ coefficient matrices $P_{n}=(p_{ij}(\sigma_{n}))$, $Q_{n}=(q_{ij}(\sigma_{n}))$, $R_{n}=(r_{ij}(\sigma_{n}))$, which depend on the step size ratio $\sigma_{n}:={\triangle t}_{n}/{\triangle t}_{n-1}$, and approximations

[TABLE]

Here, $\mathbb{V}^{s}=\mathbb{R}^{ms}$, $t_{n}={\triangle t}_{0}+\ldots+{\triangle t}_{n-1}$, $n\geq 0$, and the nodes $c_{i}\in\mathbb{R}$ are such that $c_{i}\neq c_{j}$ if $i\neq j$, and $c_{s}=1$. Further, $F(w)=[F(w_{i})]\in\mathbb{V}^{s}$ is the application of $F$ to all components of $w\in\mathbb{V}^{s}$. The starting vector $w_{0}=[w_{0,i}]\in\mathbb{V}^{s}$ is supposed to be given, or computed by a Runge-Kutta method, for example.

Peer methods belong to the class of general linear methods introduced by Butcher [3]. All approximations have the same order, which gives the name of the methods. Here, we are interested in A-stable and super-convergent Peer methods with order of convergence $p\!=\!s+1$ even for variable step sizes. For constant step sizes, such methods have been recently constructed by Soleimani and Weiner [18] and Schneider, Lang and Hundsdorfer [16]. In the following, for an $s\times s$ matrix we will use the same symbol for its Kronecker product with the identity matrix as a mapping from the space $\mathbb{V}^{s}$ to itself. Then, (3) simply reads

[TABLE]

In what follows, we discuss requirements and desirable properties for the implicit Peer method (5).

Accuracy. Let $e=(1,\ldots,1)^{T}\in\mathbb{R}^{s}$. We assume pre-consistency, i.e., $P_{n}e\!=\!e$, which means that for the trivial equation $u^{\prime}(t)=0$, we get solutions $w_{n,i}=1$ provided that $w_{0,j}=1,\;j=1,\ldots,s$. The residual-type local errors result from inserting exact solution values $w(t_{n})=[u(t_{n}+c_{i}{\triangle t}_{n})]\in\mathbb{V}^{s}$ in the implicit scheme (5):

[TABLE]

Let $c=(c_{1},\ldots,c_{s})^{T}$ with point-wise powers $c^{j}=(c_{1}^{j},\ldots,c_{s}^{j})^{T}$. Then Taylor expansion with the expressions

[TABLE]

gives

[TABLE]

from which we obtain

[TABLE]

with

[TABLE]

A pre-consistent method is said to have stage order $q$ if $d_{n,j}=0$ for all $\sigma_{n}$ and $j=1,2,\ldots,q$. With the Vandermonde matrices

[TABLE]

and $C=\text{diag}(c_{1},c_{2},\ldots,c_{s})$, $D=\text{diag}(1,2,\ldots,s)$, and $S_{n}=\text{diag}(1,\sigma_{n},\ldots,\sigma_{n}^{s-1})$, the conditions for having stage order $s$ for the implicit Peer method (5) for variable step sizes are

[TABLE]

Since $V_{1}$ and $D$ are regular, we have the relation

[TABLE]

showing that $Q_{n}$ is uniquely defined by the choice of $P_{n}$, $R_{n}$, the node vector $c$, and the step size ratio $\sigma_{n}$. Moreover, there is an easy way to achieve consistency for any choice of the step sizes ${\triangle t}_{n}$ by setting $P_{n}\equiv P$ and $R_{n}\equiv R$ with constant matrices $P$ and $R$, and recomputing $Q_{n}$ from (14) in each time step. In what follows, we will make use of this simplification and consider implicit Peer methods with variable step sizes ${\triangle t}_{n}$ of the form

[TABLE]

with constant matrices $P$ and $R$, and $Q_{n}$ updated in each time step by

[TABLE]

The matrix $R$ is taken to be lower triangular with constant diagonal $r_{ii}\!=\!\gamma\!>\!0$, $i=1,\ldots,s$, giving singly diagonally implicit methods.

Remark 2.1**.**

Implicit Peer methods of the form (15) that are consistent of order $s$ for constant time steps, i.e., ${\triangle t}_{n}={\triangle t}$ and $Q_{n}=Q$, can be applied in a variable time-step environment without loss of their order of consistency by updating (the original) $Q$ by $Q_{n}$ from (16) in each time step. We will use this modification in the numerical comparisons for our recently developed methods in [16, 19].

Stability. Applying the implicit method (15) to Dahlquist’s test equation $y^{\prime}\!=\!\lambda y$ with $\lambda\in\mathbb{C}$, gives the following recursion for the approximations $w_{n}$:

[TABLE]

with $z_{n}\!:=\!\sigma_{n}\sigma_{n-1}\cdots\sigma_{1}z_{0},\;z_{0}={\triangle t}_{0}\lambda$. Hence,

[TABLE]

The asymptotic behaviour of the matrix product is very difficult to analyse, see e.g. the discussion in Jackiewicz, Podhaisky, and Weiner [11, Sect. 1]. Here, we consider methods that are zero-stable for arbitrary step sizes and A-stable for constant step sizes. Zero stability requires the constant matrix $P\!=\!M_{im}(0,\sigma)$ to be power bounded to have stability for the trivial equation $u^{\prime}(t)\!=\!0$. We will derive methods for which the spectral radius of $M_{im}(z,\sigma)$ satisfies $\rho(M_{im}(z,\sigma))\leq 1$ for all $\sigma\in[\sigma_{min},\sigma_{max}]$ with $0\leq\sigma_{min}<1\leq\sigma_{max}$ and all $z\in\mathbb{C}$ with $\mathrm{Re}(z)\leq 0$. Since for constant step sizes, $M_{im}(\infty,1)\!=\!-R^{-1}Q(1)$ with $Q(1)\neq 0$, A-stability does not imply L-stability. To guarantee good damping properties for very stiff problems, we will aim at having a small spectral radius of $R^{-1}Q(\sigma)$ for $\sigma\in[\sigma_{min},\sigma_{max}]$. Although we cannot prove boundedness of the matrix product in (18) for $n\rightarrow\infty$ and variable step sizes, the methods derived along the design principles described above performed always stable in our numerical applications for various step size patterns.

Super-convergence. Applying the convergence theory from multistep methods, stage order $q\!=\!s$ and zero stability yield convergence of order $p\!=\!s$ for variable step sizes with $0\leq\sigma_{min}<\sigma_{n}<\sigma_{max}$ and ${\triangle t}_{n}\leq{\triangle t}_{max}:=\max_{i=0,\ldots,n}{\triangle t}_{i}$ demonstrated, e.g., in [1, 13]. Here, we are interested in using the degrees of freedom provided by the free parameters in $P$, $R$, and $c$ to have convergence of order $p\!=\!s+1$ without raising the stage order further. This is discussed under the heading super-convergence in the book of Strehmel, Weiner and Podhaisky [20, Sect. 5.3] for non-stiff problems. Similar results for stiff systems were obtained by Hundsdorfer [9]. We follow the approach recently developed in Schneider, Lang, and Hundsdorfer [16] for having an extra order of convergence for Peer methods with constant step sizes to later discuss the property of super-convergence for IMEX-Peer methods based on extrapolation for variable step sizes.

Let $\varepsilon_{n}\!=\!w(t_{n})-w_{n}$ be the global error. Under the standard stability assumption, where products of the transfer matrices are bounded in norm by a fixed constant $K$ (see, e.g., Theorem 2 in [18]), we get the estimate $\|\varepsilon_{n}\|\leq K(\|\varepsilon_{0}\|+\|r_{1}\|+\ldots+\|r_{n}\|)$. Together with stage order $s$, this gives the standard convergence result

[TABLE]

Then we have the following

Theorem 2.1**.**

Assume the implicit Peer method (15) has stage order $s$ and estimate (2.1) holds true for the global error with $\|\varepsilon_{0}\|={\cal O}\left({\triangle t}_{0}^{s}\right)$. Then the method is convergent of order $p\!=\!s$, i.e., the global error satisfies $\varepsilon_{n}={\cal O}\left({\triangle t}_{max}^{s}\right)$. Furthermore, if the initial values are of order $s+1$, $d_{i,s+1}\in\text{range}\,(I-P)$ and ${\triangle t}_{i-1}=\left(1+{\cal O}({\triangle t}_{max})\right){\triangle t}_{i}$ for all $i=1,\ldots,n$, then the order of convergence is $p\!=\!s+1$.

Proof: The first statement follows directly from (2.1) with the estimate

[TABLE]

Suppose that $d_{i,s+1}=(I-P)v_{i}$ with $v_{i}\in\mathbb{R}^{s}$. Since $I-P$ has an eigenvalue zero, $v_{i}$ is not uniquely determined. To fix $v_{i}$, we choose the one with minimum Euclidean norm, i.e., $v_{i}=(I-P)^{+}d_{i,s+1}$ with $(I-P)^{+}$ being the Moore-Penrose inverse. Let now

[TABLE]

Insertion of these modified solution values in the scheme (15) will give modified local errors

[TABLE]

where

[TABLE]

Next, we will show that $T(v_{i-1},v_{i})={\cal O}({\triangle t}_{max}){\triangle t}_{i}^{s+1}$. From the assumption on the step sizes, ${\triangle t}_{i-1}=(1+{\cal O}({\triangle t}_{max})){\triangle t}_{i}$, we deduce $\sigma_{i}^{-1}=1+{\cal O}({\triangle t}_{max})$, which yields

[TABLE]

The definition of $d_{i,s+1}$ in (11) gives the polynomial representation $d_{i,s+1}=\sum_{j=0,\ldots,s+1}a_{j}\sigma^{-j}$ with $\sigma$-independent $a_{j}\in\mathbb{R}^{s}$ (see also (26) and (27) for more details). Hence, we have $d_{i,s+1}-d_{i,s}={\cal O}({\triangle t}_{max})$. Using ${\triangle t}_{i-1}^{s+1}=(1+{\cal O}({\triangle t}_{max})){\triangle t}_{i}^{s+1}$, we conclude that

[TABLE]

which, due to (10), reveals $\bar{r}_{i}\!=\!{\cal O}({\triangle t}_{max}{\triangle t}_{i}^{s+1})$ in (21). This yields, in the same way as above, $\|\bar{\varepsilon}_{n}\|\!=\!\|\bar{w}(t_{n})-w_{n}\|\!\leq\!K\|\varepsilon_{0}\|+{\cal O}({\triangle t}_{max}^{s+1})$. Since $\|\bar{\varepsilon}_{n}-\varepsilon_{n}\|\!\leq\!{\triangle t}_{n}^{s+1}\|v_{n}\|_{\infty}\|u^{(s+1)}(t_{n})\|$ and $\|\varepsilon_{0}\|\!=\!{\cal O}({\triangle t}_{0}^{s+1})$, this shows convergence of order $s+1$ for the global errors $\varepsilon_{n}$. $\Box$

Recall that the range of $I-P$ consists of the vectors that are orthogonal to the null space of $I-P^{T}$. If the method is zero-stable, then this null space has dimension one. Therefore, up to a constant there is a unique vector $v\in\mathbb{R}^{s}$ such that $(I-P^{T})v=0$. Then we have

[TABLE]

Since $d_{i,s+1}$ depends on $\sigma_{i}$, these equations have to be satisfied for all $\sigma_{i}$. In the following, we will drop the index $i$ and examine $v^{T}d_{s+1}(\sigma)$ as a function of $\sigma$. From (11), we find

[TABLE]

where $Q(\sigma)$ is taken from (16) with $\sigma_{n}=\sigma$. Replacing $Q$, using the definition of $S_{n}$, and separating all powers of $\sigma$, we eventually get the polynomial representation

[TABLE]

with the $\sigma$-independent coefficients

[TABLE]

Here, $\tilde{D}:=(s+1)D^{-1}-I$. Note that we have used the relation $v^{T}P=v^{T}$ to eliminate $P$ in (31). With $h_{j}:=\tilde{v}_{s+1-j}^{T}\,\tilde{c}_{s+1-j}$ for $j=1,\ldots,s$, condition (25) can be fulfilled by adding the $s+2$ additional equations $h_{j}\equiv 0$ to the consistency conditions in order to achieve super-convergence for variable step sizes. The special structure of the coefficients $\tilde{c}_{s+1-j}$ allows the following statement.

Lemma 2.1**.**

Assume $c_{1},\ldots,c_{s-1}<1$ with $s\geq 2$, $c_{i}\neq c_{j}$ for $i\neq j$, and $c_{s}=1$. Then, $\tilde{c}_{1}=0$ and $\tilde{c}_{2},\ldots,\tilde{c}_{s}\neq 0$.

Proof: The conditions on $c_{i}$ guarantee the regularity of the Vandermonde matrix $V_{1}$. Let $x_{i}:=c_{i}-1$, $i=1\ldots,s$. Then, we have $x^{s}=(c-e)^{s}$ and $V_{1}=(x_{i}^{j-1})$, $i,j=1,\ldots,s$. From (30), we deduce $V_{1}\tilde{c}=x^{s}$. The choice $c_{s}=1$ yields $x_{s}=0$ and hence $\tilde{c}_{1}=0$ from the last equation. Further, assumption $c_{i}<1$ gives $x_{1},\ldots,x_{s-1}<0$. This allows division by $x_{i}$, resulting in the linear equations

[TABLE]

Now, let us consider the polynomial of order $s-1$,

[TABLE]

Then, $p(x_{i})=0$ is the $i$-th row of system (32) and hence $x_{1},\ldots,x_{s-1}$ are the $s-1$ roots of $p$, i.e., $p(x)=(x-x_{1})\cdots(x-x_{s-1})$. The theorem of Vieta shows

[TABLE]

Since $x_{i}<0$ for all $i=1,\ldots,s-1$, we observe that all products in the sum have the same number of factors which have one and the same sign, i.e., the sums cannot vanish. More precisely, $\text{sgn}(\kappa_{j})=(-1)^{j}$ and hence $\text{sgn}(\tilde{c}_{s+1-j})=-1$, which proves the statement. $\Box$

We are now ready to formulate additional simplified conditions for the super-convergence of implicit Peer methods when they are applied with variable step sizes.

Theorem 2.2**.**

Assume the implicit Peer method (15) has stage order $s$ and estimate (2.1) holds true for the global error with $\|\varepsilon_{0}\|={\cal O}\left({\triangle t}_{0}^{s+1}\right)$. Let ${\triangle t}_{i-1}=(1+{\cal O}({\triangle t}_{max})){\triangle t}_{i}$ for all $i=1,\ldots,n$. Then the method is convergent of order $p\!=\!s+1$, i.e., the global error satisfies $\varepsilon_{n}={\cal O}\left({\triangle t}_{max}^{s+1}\right)$, if for all $v\in\mathbb{R}^{s}$ with $(I-P^{T})v=0$, the following additional conditions are satisfied:

[TABLE]

Proof: Condition $v^{T}d_{s+1}(\sigma)=0$ requires $h_{j}=0$ for $j=0,\ldots,s+1$ in (27). Observe that $h_{s}=0$ is always satisfied since $c_{s}=1$ and hence $\tilde{c}_{1}=0$. The property $\tilde{c}_{j}\neq 0$ for $j=2,\ldots,s$, leads to $\tilde{v}^{T}_{j}=0$, which is equivalent to $v^{T}(c^{j}-jRc^{j-1})=0$. The remaining conditions follow directly from $h_{0}=h_{s+1}=0$. $\Box$

2.2 Super-convergent explicit Peer methods for variable steps sizes

Super-convergent explicit Peer methods for variable step sizes with a special structure of the matrix $P$ have first been constructed by Weiner, Schmitt, Podhaisky, and Jebens [21]. A convenient way to construct such methods for more general $P$ is the use of extrapolation as proposed by Schneider, Lang, and Hundsdorfer [16]. This idea goes back to Crouzeix [7] and was also used by Cardone, Jackiewicz, Sandu, and Zhang [6] to construct implicit-explicit diagonally implicit multistage integration methods. The procedure can be easily extended to variable step sizes.

Assume that all approximations $w_{n,j}$ obtained from method (15) have stage order $s$. Then, we can use $w_{n-1}$ and most recent values $w_{n,j}$, $j=1,\ldots,i-1$, already available for the computation in the $i$-th stage, to extrapolate $F(w_{n})$ by

[TABLE]

where $\tau_{n}=\max{({\triangle t}_{n-1},{\triangle t}_{n})}$ and the $s\times s$-matrices $E_{1,n}$ and $E_{2,n}$ of extrapolation coefficients depend on the step size ratio $\sigma_{n}$. Here, $E_{2,n}$ is a strictly lower triangular matrix. Replacing $F(w_{n})$ in (15) gives the explicit method

[TABLE]

Note that $RE_{2,n}$ is strictly lower triangular since $R$ is lower triangular. We will discuss consistency and super-convergence of this explicit method.

Accuracy. Taylor expansion with exact values $F(w(t_{n}))$ gives for the residual-type error vector

[TABLE]

Then, the residual-type local error of the explicit Peer method (38) reads

[TABLE]

with

[TABLE]

We can achieve stage order $s$, if the underlying implicit Peer method has stage order $s$, i.e., $d_{n,j}=0$ for all $\sigma_{n}$ and $j=1,\ldots,s$, and if we choose

[TABLE]

with a constant $s\times s$-matrix $E_{2}$ and $S_{n}=\text{diag}(1,\sigma_{n},\ldots,\sigma_{n}^{s-1})$ as defined above. This gives $l_{n,j}=0$ for all $\sigma_{n}$ and $j=0,\ldots,s-1$ and eventually $r_{n}={\cal O}({\triangle t}_{n}^{s+1})$.

Super-convergence. Under standard stability assumptions as for the implicit method, we derive the global error estimate for the explicit Peer method defined in (38),

[TABLE]

Analogously, we have

Theorem 2.3**.**

Assume the implicit Peer method (15) has stage order $s$ and estimate (2.2) holds true for the global error with $\|\varepsilon_{0}\|={\cal O}\left({\triangle t}_{0}^{s}\right)$. Then the explicit method (38) is convergent of order $p\!=\!s$, i.e., the global error satisfies $\varepsilon_{n}={\cal O}\left({\triangle t}_{max}^{s}\right)$. Furthermore, if the initial values are of order $s+1$, $(d_{i,s+1}+Rl_{i,s})\in\text{range}\,(I-P)$ and ${\triangle t}_{i-1}=(1+{\cal O}({\triangle t}_{max})){\triangle t}_{i}$ for all $i=1,\ldots,n$, then the order of convergence is $p\!=\!s+1$.

Proof: Replacing $d_{i,s+1}$ by $d_{i,s+1}+Rl_{i,s}$ in the proof of Theorem 2.1 gives the desired result. $\Box$

Thus, super-convergence for variable step sizes is achieved if for all $i=1,\ldots,n,$ it holds

[TABLE]

If the underlying implicit method is already super-convergent, the conditions simplify to $v^{T}Rl_{i,s}=0$. Next, we will study the $l_{i,s}$ as functions of $\sigma$ and derive sufficient conditions for order $s+1$.

From (41) and (42), we get

[TABLE]

The investigation of the product $v^{T}(d_{s+1}(\sigma)+Rl_{s}(\sigma))$ yields the following

Theorem 2.4**.**

Assume the explicit Peer method (38) has stage order $s$ and estimate (2.2) holds true for the global error with $\|\varepsilon_{0}\|={\cal O}\left({\triangle t}_{0}^{s+1}\right)$. Let ${\triangle t}_{i-1}=(1+{\cal O}({\triangle t}_{max})){\triangle t}_{i}$ for all $i=1,\ldots,n$. Then the method is convergent of order $p\!=\!s+1$, i.e., the global error satisfies $\varepsilon_{n}={\cal O}\left({\triangle t}_{max}^{s+1}\right)$, if for all $v\in\mathbb{R}^{s}$ with $(I-P^{T})v=0$, the following additional conditions are satisfied:

[TABLE]

Proof: The proof follows the same way as demonstrated in the proof of Theorem 2.1. The coefficients of $\sigma^{-s},\ldots,\sigma^{-1}$ are again expressed as products $\tilde{v}_{1}\tilde{c}_{1},\ldots,\tilde{v}_{s}\tilde{c}_{s}$ with

[TABLE]

Due to $\tilde{c}_{j}\neq 0$ for $j=2,\ldots,s$, we have $\tilde{v}^{T}_{j}=0$ and hence $v^{T}(c^{j}-jRE_{2}c^{j-1})=0$. The other condition remains unchanged. $\Box$

We would like to conclude with the following observation: If we start with a super-convergent implicit Peer method for variable step sizes, i.e., the additional conditions in Theorem 2.2 are already fulfilled, then (46) disappears and (47) changes to $v^{T}R(E_{2}-I)c^{j-1}=0$. This can be rewritten to $v^{T}R(E_{2}-I)CV_{0}=0$. Since $R(E_{2}-I)$ and $V_{0}$ are regular matrices, $C$ must be singular to satisfy (47) for $v\neq 0$. That means, one of the nodes $c_{i}$ must be zero, because we always assume $c_{i}\neq c_{j}$. We will discuss this point later.

2.3 Super-convergent IMEX-Peer methods with variable step sizes

We now apply the implicit and explicit methods (15) and (38) to systems of the form

[TABLE]

where $F_{0}$ will represent the non-stiff or mildly stiff part, and $F_{1}$ gives the stiff part of the equation. The resulting IMEX scheme is

[TABLE]

where ${\hat{Q}}_{n}=Q_{n}+RE_{1,n}$, ${\hat{R}}=RE_{2}$, and extrapolation is used only on $F_{0}$. Combining the local consistency analysis for both the explicit and implicit method, the residual-type local errors for the IMEX-Peer methods have the form

[TABLE]

Super-convergence. In order to construct super-convergent IMEX-Peer methods of order $s+1$ for variable step sizes, we have to impose consistency of order $s$ and ensure that for all $v\in\mathbb{R}^{s}$ with $(I-P^{T})v=0$ it holds

[TABLE]

for all $\sigma$. We have the following

Theorem 2.5**.**

Let the $s$-stage implicit Peer method (15) defined by the coefficients $(c,P,Q_{n},R)$, with $Q_{n}$ from (16), be zero-stable and suppose its stage order is equal to $s$. Let the initial values satisfy $w_{0,i}-u(t_{0}+c_{i}{\triangle t}_{0})={\cal O}({\triangle t}_{0}^{s+1})$, $i=1,\ldots,s,$ and ${\triangle t}_{i-1}=\left(1+{\cal O}({\triangle t}_{max})\right){\triangle t}_{i}$, $i=1,\ldots,n$. Then the IMEX-Peer method (49) is convergent of order $s+1$, i.e., the global error satisfies $\varepsilon_{n}={\cal O}({\triangle t}_{max}^{s+1})$, if for all $v\in\mathbb{R}^{s}$ with $(I-P^{T})v=0$, the following additional conditions are satisfied:

[TABLE]

Proof: Suppose $d_{i,s+1}=(I-P)v_{d,i}$ and $R\,l_{i,s}=(I-P)v_{l,i}$ with $v_{d,i},v_{l,i}\in\mathbb{R}^{s}$. Again, we fix these vectors by setting $v_{d,i}=(I-P)^{+}d_{i,s+1}$ and $v_{l,i}=(I-P)^{+}Rl_{i,s}$ with $(I-P)^{+}$ being the Moore-Penrose inverse. Let now

[TABLE]

Inserting these modified values in (49) gives the modified residual-type local errors

[TABLE]

which can be rearranged to

[TABLE]

Then, Taylor expansions yields

[TABLE]

with $T(\cdot,\cdot)$ and $r_{i}$ as defined in (22) and (50), respectively. The same arguments as in the proof of Theorem 2.1 show $\bar{r}_{i}={\cal O}({\triangle t}_{max}{\triangle t}_{i}^{s+1})$ and eventually the convergence of order $s+1$ for the global errors $\varepsilon_{n}=w(t_{n})-w_{n}$. $\Box$

The $2s+1$ additional conditions (52)-(54) are quite demanding. We have already mentioned the fact that (54) requests that one of the nodes $c_{i}$, $i\neq s$, must be zero. In this case, the method delivers two vectors, $w_{n-1,s}$ and $w_{n,i}$ with a certain $i$, that approximate $u(t_{n})$. We note that the difference of these approximations is used in the extrapolation process as an additional degree of freedom. The matrix $E_{1,n}$ in (42) is still well defined. However, it is not always possible to construct such methods at all or with good stability properties in particular. In many practical applications, it might be sufficient that the explicit method has the property of super-convergence for variable step sizes and the implicit method is only super-convergent for constant step sizes. We have constructed such methods as well. They have to fulfill the following additional conditions for all $v\in\mathbb{R}^{s}$ with $(I-P^{T})v=0$ and for all $\sigma$:

[TABLE]

Due to the second condition for $\sigma=1$, the first one can be replaced by the often simpler requirement $v^{T}Rl_{s}(1)=0$. Using Theorem 2.4 and the definition of $Rl_{s}$, we find the explicit relations

[TABLE]

Compared to (52)-(54), the number of conditions has been significantly reduced. Moreover, since condition (54) disappeared, the restriction $c_{i}=0$ for a certain $i$ is no longer necessary.

2.4 Stability of super-convergent IMEX-Peer methods

We consider the usual split scalar test equation

[TABLE]

with complex parameters $\lambda_{0}$ and $\lambda_{1}$. Applying an IMEX-Peer method (49) to (63) gives the recursion

[TABLE]

with $z_{i}^{(n)}={\triangle t}_{n}\lambda_{i}$, $i=0,1$. As for the implicit method itself, an analysis of matrix products formed by $M_{1}M_{2}\cdots M_{n}$ would be far too complicated. Therefore, we restrict ourselves to constant step sizes and require

[TABLE]

with $z_{i}={\triangle t}\lambda_{i}$, $i=0,1$. Then, the stability regions of the IMEX-Peer method applied with constant step sizes are defined by the sets

[TABLE]

in the left-half complex plane for $\alpha\in[0^{\circ},90^{\circ}]$. Further, we define the stability region of the corresponding explicit method (with constant step sizes) as

[TABLE]

with the stability matrix $M(z_{0},0)=(I-z_{0}{\hat{R}})^{-1}(P+z_{0}{\hat{Q}})$. Efficient numerical algorithms to compute $\mathbb{S}_{\alpha}$ and $\mathbb{S}_{E}$ are extensively described in [6, 12].

Our goal is to construct IMEX-Peer methods for which $\mathbb{S}_{E}$ is large and $\mathbb{S}_{E}\backslash\mathbb{S}_{\alpha}$ is as small as possible for angles $\alpha$ that are close to $90^{\circ}$. We will construct super-convergent IMEX-Peer methods with A-stable implicit part for constant step sizes, i.e., the stability region $\mathbb{S}_{90^{\circ}}$ is non-empty. Concerning variable step sizes, we follow the design principles already stated in the stability discussion in Section 2.1.

2.5 Practical Issues

Starting procedure. In order to execute the first step of the IMEX-Peer method (49), we have to choose $t_{1}$, ${\triangle t}_{0}$, ${\triangle t}_{1}$, and need to approximate the $s$ initial values $w_{0,i}\approx u(t_{1}-(1-c_{i}){\triangle t}_{0})$. For this, we apply a suitable integration method with continuous output, e.g. a Runge-Kutta or BDF scheme, on the interval $[t_{0},t_{0}+\tau]$ with $\tau>0$. The accuracy of the continuous numerical solution $\tilde{w}(t)$ can be controlled by standard step size control or by choosing $\tau$ sufficiently small. Denoting the minimum and maximum component of the node vector $c$ by $c_{min}$ and $c_{max}$, respectively, we require

[TABLE]

This linear system for $t_{1}$ and ${\triangle t}_{0}$ has the unique solution

[TABLE]

The initial values are now taken from

[TABLE]

Note that $w_{0,i}=u_{0}$ for index $i$ with $c_{i}=c_{min}$. Eventually, we set ${\triangle t}_{1}={\triangle t}_{0}$.

Step size selection. We extend the approach proposed by Soleimani, Knoth, and Weiner in [17] to locally approximate ${\triangle t}_{n}^{s}u^{(s)}(t_{n})$, which mimics the leading error term of an embedded solution of order $s-1$. Let $F=F_{0}+F_{1}$ and define

[TABLE]

with $\alpha$ and $\beta$ determined through

[TABLE]

where $e_{s}^{T}=(0,\ldots,0,1)$ and $\delta\in[0,1]$ is chosen as a weighting factor. Then Taylor expansion of the exact solution shows the desired property:

[TABLE]

In our numerical experiments, we have discovered that the use of old function values, i.e., $\delta=0$ in (72), works quite reliable for stiff and very stiff problems. For mildly stiff problems, the choice $\delta=1$ often leads to a slightly better performance. For our examples in Section 4, we will present results for $\delta=0$.

The new step size is computed by

[TABLE]

with the weighted relative maximum error

[TABLE]

In order to reach the time end point $T$ with a step of averaged normal length, we adjust after each step size ${\triangle t}_{new}$ to ${\triangle t}_{new}=(T-t_{n})/\lfloor(1+(T-t_{n})/{\triangle t}_{new})\rfloor$.

Given an overall tolerance $TOL$, the step is accepted and the computation is continued with ${\triangle t}_{n+1}={\triangle t}_{new}$, if $err\leq TOL$. Otherwise, the step is rejected and repeated with ${\triangle t}_{n}={\triangle t}_{new}$.

3 Construction of super-convergent IMEX-Peer methods with variable

step sizes

3.1 The case $s=2$

First, we have a negative result. With $c_{1}=0$, $c_{2}=1$, and pre-consistency $Pe=e$, the coefficient matrices are

[TABLE]

The first condition (52) for super-convergence reads $(-1/2,0)\,v=0$, which gives, up to scaling, $v=(0,1)^{T}$. Then, (53) reduces to $1-2\gamma=1-3\gamma=0$, which is not possible for any $\gamma$.

Next we try to find methods that satisfy (60)-(62) with $c_{1}\neq 0$. There are indeed candidates with $c_{1}=2/3$, $p_{2}=0$, $e_{12}=3/(4\gamma)$, and $r_{21}=3/4-2\gamma$. The remaining parameters $p_{1}$ and $\gamma$ are chosen such that the implicit part is A-stable and the stability regions of the IMEX-method are optimized. Good results are obtained for the following method:

[TABLE]

We will refer to this method as IMEX-Peer2sve.

3.2 The cases $s=3$ and $s=4$

In order to construct super-convergent methods for variable step sizes, we have to satisfy conditions (52)-(54) for all $v\in\mathbb{R}^{s}$ with $(I-P^{T})v=0$ and one of the nodes $c_{i}$ being zero. A surprisingly simple choice is $c_{1}=0$ and $v=e_{1}$, which yields the validity of (53) and (54). Then, equation (52) yields one condition for the remaining nodes. We find $c_{2}=0.5$ for $s=3$ and $c_{3}=(5c_{2}-1)/(10c_{2}-5)$ for $s=4$. Furthermore, the first row of $P$ is $e_{1}$, which goes along with pre-consistency. The value of $c_{3}$ and the remaining coefficients of $P$, $R$ and $E_{2}$ are chosen in such a way that the implicit Peer methods are A-stable and the IMEX-Peer methods exhibit good stability properties and moderate error constants. This has been done using the Matlab-routine fminsearch, where we included the desired properties in the objective function and used random start values for the remaining degrees of freedom. Different combinations of weights in the objective function have been employed to select promising candidates which were then tested in various problems. We will refer to the methods finally selected as IMEX-Peer3sv and IMEX-Peer4sv.

We have also constructed a $4$-stage IMEX-Peer method, denoted by IMEX-Peer4sve, with the property that the explicit method is super-convergent for variable step sizes and the implicit method is only super-convergent for constant step sizes. In this case, conditions (60)-(62) must be satisfied, where $c_{i}$, $i=1,2,3$, are still free parameters. We set $v=e_{s}$, which gives (61) since then $v^{T}(C-I)=0$. The additional degrees of freedom in the nodes allow us to achieve greater stability regions and smaller error constants compared to IMEX-Peer4sv. The method found is optimally zero-stable, i.e., one eigenvalue of $P$ equals one (due to pre-consistency) and the others are zero.

The coefficients of all new methods for $c$, $P$, $R$, and $E_{2}$ are given in Table 2 and Table 3. Values for the stability regions as well as other constants are collected in Table 1. More details on the stability regions are shown in Figure 1. Obviously, the new property of super-convergence for variable step sizes comes with significantly smaller stability regions, except for IMEX-Peer4sve which even slightly improves $S_{90^{\circ}}$ of IMEX-Peer4s.

4 Numerical examples

We will present results for two ODE and two PDE problems. In order to guarantee that errors of the initial values do not affect the computations, unknown initial values as well as reference solutions $Y$ at the final time are computed by ode15s from Matlab with sufficiently high tolerances. In the comparisons, the global errors are computed by $err=\max_{i}|Y_{i}-\hat{Y}_{i}|/(1+|Y_{i}|)$, where $\hat{Y}$ is the numerical approximation.

All calculations have been done with Matlab-Version R2017a on a Latitude 7280 with an i5-7300U Intel processor at 2.7 GHz.

4.1 Prothero-Robinson Problem

In order to study the rate of convergence under stiffness and changing step sizes, we consider the Prothero-Robinson type equation used in [16, 17],

[TABLE]

where $t\in[0,5]$. The first term is treated explicitly and the second implicitly. Initial values are taken from the analytic solution $y(t)\!=\!(\cos(t),\sin(t))$. For constant step sizes ${\triangle t}=0.05/i$, $i=1,\ldots,6$, we consider the $\sigma$-dependent sequences

[TABLE]

with ${\triangle t}_{1}=2{\triangle t}/(1+\sigma)$ and $N=T/{\triangle t}$. Results for $\sigma=1.0,1.1,1.2$ are shown in Figure 2. Since the $4$-stage methods become instable for $\sigma=1.2$, these results are omitted. One can nicely see that all new methods keep their order of convergence observed for constant step sizes and, therefore, perform quite robust with respect to changing the step size. This is, of course, not the case for the methods that are only super-convergent for constant step sizes.

4.2 Van der Pol Oscillator

Next we consider the well known stiff van der Pol oscillator

[TABLE]

with $y_{1}(0)=2$, $y_{2}(0)=0$, and $t\in[0,2]$. The first term is treated explicitly and the second implicitly. This singularly perturbed problem challenges any code and its efficient solution requires a step size adaptation over several orders of magnitude, see e.g. [8] and the discussions therein. The tolerances are $atol=rtol=10^{-3-i}$, $i=0,1,\ldots,4$ and the calculations are started with initial step $\tau=atol$ for all methods. The results are shown and discussed in Figure 3.

4.3 Burgers Problem

The first PDE problem is taken from [4], see also [19] for further numerical results with super-convergent IMEX-Peer methods. We consider

[TABLE]

with initial value $u(0,x)=\sin(\pi(x+1))$ and homogeneous Dirichlet boundary conditions. The source term is defined through

[TABLE]

The spatial discretization is done by finite differences with $\triangle x=1/2500$. We treat the diffusion implicitly and all other terms explicitly.

We have used tolerances $atol=rtol=10^{-2-i}$, $i=0,1,\ldots,5$ and initial step sizes $\tau=\sqrt{atol}$. The results are plotted and discussed in Figure 4.

4.4 Linear Advection-Reaction Problem

A second PDE problem for an accuracy test is the linear advection-reaction system from [10]. The equations are

[TABLE]

for $0<x<1$ and $0<t\leq 1$, with parameters

[TABLE]

and with the following initial and boundary conditions:

[TABLE]

Note that there are no boundary conditions for $v$ since $\alpha_{2}$ is set to be zero.

Fourth-order finite differences on a uniform mesh consisting of $m=400$ nodes are applied in the interior of the domain. At the boundary, we can take third-order upwind biased finite differences, which here does not affect an overall accuracy of four [10] and gives rise to a spatial error of $1.5\,10^{-5}$. In the IMEX setting, the reaction is treated implicitly and all other terms explicitly.

We have used tolerances $atol=rtol=10^{-3-i}$, $i=0,1,\ldots,5$ and an initial step size $\tau=10^{-3}$ for all runs. The results are plotted and discussed in Figure 4.

5 Conclusion

We have developed a new class of $s$-stage super-convergent IMEX-Peer methods with A-stable implicit part, which maintain their super-convergence order of $s+1$ for variable step sizes. A-stability is important to solve problems with function contributions that have large imaginary eigenvalues in the spectrum of their Jacobian. Applying the idea of extrapolation and studying the $\sigma$-dependent coefficients in the local error representations, we first derived additional conditions for implicit and explicit Peer methods, which are then combined to state $2s+1$ corresponding conditions for IMEX-Peer methods. An interesting theoretical result is that one of the nodes must be zero. Such methods exist for $s\!>\!2$. We designed new methods for $s=3,4$. However, the new property of super-convergence for variable step sizes reduces the scope for achieving good stability properties, resulting in significantly smaller stability regions compared to the super-convergent IMEX-Peer methods from [16]. We also constructed methods for $s=2,4$ having an explicit part that is super-convergent for variable step sizes, whereas the implicit part is only super-convergent for constant steps. In all cases, we employed the Matlab-routine fminsearch with varying objective functions and starting values to find suitable methods with stability regions as large as possible, good damping properties for very stiff problems and small error constants.

We have implemented our newly designed methods with local error control based on linear combinations of old function evaluations to approximate the leading error term of an embedded solution of order $s\!-\!1$. From our observations made for four numerical examples, we can draw the following conclusions: (i) The new methods perform quite robust with respect to changing the step size and, as expected, show their theoretical order at the same time. (ii) For problems that demand a fast step size adaptation over several orders of magnitudes, like the van der Pol oscillator, the new methods have the potential to perform better. (iii) For problems that can be integrated with moderate step size changes, like the Burgers problem, super-convergence for constant step sizes is still sufficient to profit from the additional order and possibly from the larger stability regions.

6 Acknowledgement

J. Lang was supported by the German Research Foundation within the collaborative research center TRR154 “Mathematical Modeling, Simulation and Optimisation Using the Example of Gas Networks” (DFG-SFB TRR154/2-2018, TP B01) and the Graduate Schools Computational Engineering (DFG GSC233) and Energy Science and Engineering (DFG GSC1070).