1 Introduction↩︎

In this paper, we consider systems of ordinary differential equations (ODEs) representing space discretisations of partial differential equations (PDEs) of the form \[\label{eq:main} \begin{eqnarray} \dot{y} &=& F(y), \\ F(y) &:=& F_{D}(y)+F_{A}(y)+F_{R}(y). \end{eqnarray}\tag{1}\] where \(F_{D}(y),F_{A}(y),F_{R}(y) \in \mathbb{R}^{d}\) represent diffusion, advection and reaction terms respectively. Typically 1 arises from the spatial discretisation of an advection-diffusion-reaction problem in \(N\) spatial dimensions, and assuming a spatial grid of size \(\Delta x\), the dimension of the ODE grows like \(d =\mathcal{O}((\Delta x)^{-N})\).

The numerical time integration of 1 is very challenging since this is a stiff ODE. Following [1], we call a differential equation stiff when a standard explicit numerical integrator such as a standard explicit Runge-Kutta method faces a severe time-step restriction. In particular, in the case of 1 the eigenvalues of the Jacobian of \(F_{D}\) are typically distributed along the negative real axis in an interval that grows as \([-\mathcal{O}(\Delta x^{-2}),0]\) (for a symmetric diffusion operator), which implies that the time step \(\Delta t\) that can be used by the explicit Euler method suffers from the CFL \(\Delta t \leq C \Delta x^{2}\). Furthermore, the eigenvalues of the Jacobian of \(F_{A}\) are typically distributed along the imaginary axis in an interval that grows as \([- i\mathcal{O}(\Delta x^{-1}), i \mathcal{O}(\Delta x^{-1})]\), while the eigenvalues of the Jacobian of the reaction term \(F_{R}\) even though are usually not related to \(\Delta x\) can have a ratio \(\max_{j}|\mathcal{R} \lambda_{j}|/\min_{j}|\mathcal{R} \lambda_{j}|\) that can be very large or vary over several orders of magnitude when modeling multiscale reaction terms.

There are two different approaches to dealing with the issue of stiffness within the class of Runge-Kutta methods. The first relates to the design of implicit methods, with a prime example being the implicit Euler method. These methods have the advantage of achieving stability unconditionally with respect to the stepsize, which comes with the cost of solving a system of nonlinear equations. When the dimension of the corresponding differential equation is low, the nonlinear system can be solved efficiently using Newton’s method [2]. As the dimension of the differential equation increases, diagonally implicit Runge-Kutta methods are attractive as they lead to a triangular system which is easier to solve using lower-dimensional Jacobian matrices using an appropriate partitioning of the coefficient matrix as proposed in [3]. Nevertheless, when one deals with differential equations like 1 , dealing with the corresponding nonlinear equations (even in the triangular form) can be cumbersome as it requires specialised linear algebra techniques, where the use of preconditioners raises issues with both memory access and storage on parallel supercomputers.

An alternative approach to using implicit methods for solving stiff differential equations is to use explicit stabilised methods [4]. In their basic form, these are first-order explicit methods that can take larger time steps than the standard Euler method (for an easily computable numerical cost). This property together with the simplicity of implementation has made explicit stabilised methods the method of choice for diffusion-dominated problems. In particular, the second-order Runge Kutta Chebyshev (RKC) method was very successful for large dimensional parabolic problems, while methods with nearly optimal extended stability domain sizes were proposed as the second-order ROCK2 method [5] and the fourth-order ROCK4 method [6]. In the spirit of implicit-explicit methods (IMEX) [7], [8], the versatility of explicit stabilised methods allowed many extensions to partitioned Runge-Kutta methods for addressing problems with additional advection and possibly stiff reaction terms. Examples include Implicit Runge Kutta Chebyshev (IRKC) [9] to deal with stiff reaction terms and Partitioned Runge Kutta Chebyshev (PRKC) [10] to include advection terms. In addition, the Partitioned Implicit Orthogonal Runge Kutta Chebyshev method (PIROCK) [11] combined these two approaches to propose a second-order method for solving advection-diffusion reaction equations. In particular, PIROCK has recently emerged as an efficient alternative to classical splitting approaches for solar atmosphere simulations in astrophysics [12] and for large-scale gyrokinetic simulations in plasma physics [13], [14]. Another recent use of explicit stabilised methods is for solving optimisation problems. More precisely, in [15] by exploiting the good stability properties of explicit stabilised methods a new method called the Runge–Kutta Chebyshev descent (RKCD) was proposed and showed to perform similarly to the state of the art for strongly convex optimisation problems.

A major limitation of methods such as IRKC, PRKC, and PIROCK that use an explicit stabilised method to deal with the \(F_{D}\) term in 1 , is that they are at most second-order accurate in \(\Delta t\). Indeed, while the construction of the PIROCK method which is based on a splitting into three terms needed to solve 12 order conditions, similar construction of higher order methods would be cumbersome due to the dramatic increase of order conditions for partitioned Runge-Kutta methods. On the other hand, implicit methods can achieve much higher order of accuracy in \(\Delta t\), but as discussed above the implementation as well as the storage requirements of the Newton method and variants in high dimensions become an issue. In this paper, inspired by the success of RKCD as an optimisation method, we revisit diagonally implicit methods and their implementation. In particular, instead of solving the nonlinear system with variants of Newton’s method, we introduce a new variant of the RKCD method that only outputs first-order information. This method is simple to implement due to its explicit nature which circumvents the need for specialised linear algebra tools. Furthermore, it leads to rigorous estimates of the computational complexity of the proposed method. In particular, we can show that

1. In the pure linear diffusion case (\(F_{A}=F_{R}=0\) in 1 ), we prove that one time step of an SDIRK method¹ (of arbitrary order) by an iterative method requires \(\mathcal{O}(\sqrt{\Delta t}/\Delta x)\) evaluations of \(F_{D}\) avoiding the standard CFL condition of standard explicit Runge Kutta methods. This essentially implies that we are able to implement a higher order diagonally implicit Runge-Kutta method at the cost of an explicit stabilised method.

2. For a class of linear advection-diffusion PDE (\(F_{R}=0\)), we prove that in spite of relatively large Peclet numbers the proposed partitioning allows to implement one step of an SDIRK method with a similar complexity to the one of the pure diffusion case discussed at the previous point.

3. For various nonlinear diffusion-advection-reaction problems, numerical experiments illustrate the robustness of the approach where we show that the explicit stabilised iterations converge in the case of an order 4 SDIRK method with error control, a state-of-the-art Runge-Kutta method proposed in [1] and considered in this paper to illustrate the new approach.

The paper is organized as follows. In Section 2 we revisit the ideas of numerical stability of Runge-Kutta methods, while also discussing the basics of explicit stabilised methods. In Section 3 we revisit implicit Runge-Kutta methods and discuss the main idea behind solving the nonlinear system using an explicit stabilised algorithm. In addition, we discuss a specific application of the new explicit implementation of implicit methods in the context of advection-diffusion-reaction equations and present an algorithm for doing so. A detailed convergence analysis in the case of pure linear diffusion and linear diffusion-advection is presented in Section 4. Finally, in Section 5 we present a number of numerical computations that illustrate the theoretical convergence analysis in the pure linear diffusion case, and the linear diffusion-advection case. In addition, applying the new explicit stabilised implementation exSDIRK4 of an order 4 SDIRK method to several nonlinear advection-diffusion-reaction problems illustrates the wide applicability and versatility of our approach.

2 Preliminaries↩︎

In this paper, we are interested in the numerical integration of systems of ODEs of the form \[\label{eq:ode} \dot{y}(t)=g(y(t)),\qquad y(0)=y_0, \quad t>0\tag{2}\] where \(y(t) \in \mathbb{R}^{d}\) and the vector field \(g: \mathbb{R}^{d} \mapsto \mathbb{R}^{d}\) is assumed smooth. In this section, we recall the main tools for stability and accuracy analysis of Runge-Kutta methods.

2.1 Stability function and stability domain↩︎

Let \(\lambda\in \mathbb{C}\) and consider the Dahlquist test equation \[\label{eq:testode} \dot{y}(t) = \lambda y(t), \quad y(0)=y_0\in \mathbb{R}.\tag{3}\] In the case where \(\lambda \in \mathbb{C}^-=\{z \in \mathbb{C};\,\text{Re}(z) \leq 0\}\), equation 3 has a bounded solution for \(t \rightarrow +\infty\). The simplest Runge-Kutta method for solving 2 is the explicit Euler method which when applied to 3 yields the recurrence \[y_{n+1}=(1+\lambda \Delta t)y_{n}\] More generally, a Runge-Kutta method with step size \(\Delta t\) applied to 3 produces the following recurrence \[\label{eq:rec95dal} y_{n}=R(\Delta t\lambda)y_{n-1}=R(\Delta t \lambda)^ny_0,\tag{4}\] where \(R(z)\) is called the stability function and the stability domain is defined by \[\mathcal{S}=\{z\in\mathbb{C}\,;\,|R(z)|\leq 1\}.\] Note that \(z \in \mathcal{S}\) ensures that \(y_{n}\) in 4 remains bounded as \(n \rightarrow \infty\) similarly to the true solution. In the case of the explicit Euler we have that \(\mathcal{S}_{\text{E}}=\{|z+1|\leq 1 \}\) which corresponds to a disc of radius one, centred at \(-1\) on the complex plane. For stiff problems where \(|\lambda| \gg 1\) this implies a severe time-step restriction namely \(\Delta t < 2/|\lambda|\) when \(\lambda <0\). In particular, for spatial discretization of the heat equation eigenmodes are typically of size \(\mathcal{O}(\Delta x ^{-2})\) which would in turn imply the CFL condition \(\Delta t < c \Delta x^{2}\).

To avoid such severe time-step restrictions one solution is to use implicit methods. The simplest one of those is the implicit Euler method which when applied to 3 yields the recurrence \[y_{n+1}=\frac{1}{1-\lambda \Delta t} y_{n}.\] In this case we have \(\mathcal{S}_{\text{IE}}=\{|z-1| \geq 1 \}\) which contains \(\mathbb{C}^{-}\). This property is known as A-stability and can only hold in the case of implicit methods. A-stability means that there is no restriction on the time-step size for eigenvalues with negative real parts contrary to the explicit Euler method and other explicit methods. Nevertheless, in the case of severely stiff dissipative problems, a desirable additional property to A-stability is \(R(\infty)=0\) which is known as L-stability [1].

2.2 Higher order Implicit Runge-Kutta methods↩︎

The implicit Euler is an example of an \(L\)-stable integrator meaning that is particularly suitable for integrating stiff equations. However, from the point of view of accuracy, it is a first-order method. To address this issue families of \(L\)-stable Runge-Kutta methods have been proposed [1] with an arbitrarily higher order of accuracy. Important examples are the families of the Radau and Lobatto Runge-Kutta methods are that achieve orders \(2m-1\) and \(2m-2\) respectively for any number \(m\) of internal stages [16]. In particular, Radau with \(m=1\) reduces to the implicit Euler method, while in the case where \(m=2\) we obtain the following implicit Runge-Kutta method for integrating 2 , \[\label{eq:radau2} \begin{eqnarray} Y_{1}&=& y_{n}+\frac{5\Delta t}{12 }g\left(Y_{1}\right)-\frac{\Delta t}{12}g(Y_{2}), \\ Y_{2}&=& y_{n}+\frac{3\Delta t}{4 }g\left(Y_{1}\right)+\frac{\Delta t}{4}g(Y_{2}), \\ y_{n+1}&=& Y_{2}. \end{eqnarray}\tag{5}\] For a general \(m\)-stage implicit Runge-Kutta method with coefficients \(a_{ij},b_j, i,j=1,\ldots,m\), given the numerical solution \(y_n\) at time \(t_{n}\) one needs to solve the following nonlinear system in dimension \(m\cdot d\) with unknowns the internal stages \(Y_i\in\mathbb{R}^d,i=1,\ldots,m\), \[\label{eq:nonlinRK} Y_i= y_n + \Delta t \sum_{j=1}^m a_{ij} g(Y_j),\quad i=1,\ldots,m,\tag{6}\] to obtain the solution \(y_{n+1}\) at time \(t_{n+1}\) by \[\label{eq:nonlinRKfinal} y_{n+1}= y_n + \Delta t \sum_{i=1}^m b_i g(Y_i).\tag{7}\] It is common to represent Runge-Kutta methods using a Butcher tableau \[\begin{array}{c|c} c&\mathcal{A}\\ \hline & b \end{array}\] where \(\mathcal{A}\) is the matrix of coefficients \(\{a_{ij}\}_{i,j=1}^{m}\), \(b\) is the weights vector \(\{b_i\}_{i=1}^{m}\), and the vector \(c\) is such that \(c_i=\sum_{j=1}^{m} a_{ij}\). In compact matrix and tensor notations, 6 can be written as \[\label{eq:nonlinearimp} \boldsymbol{Y} = \mathbb{1}_m \otimes y_n + \Delta t (\mathcal{A}\otimes I_d) g[\boldsymbol{Y}]\tag{8}\] where \(\mathcal{A}\in \mathbb{R}^{m\times m}\) is the Runge-Kutta matrix of coefficients \(a_{ij}\) and \[\label{eq:summary} \mathbb{1}_m = \begin{pmatrix} 1\\ 1\\ \vdots\\ 1\\ \end{pmatrix}\in\mathbb{R}^m,\quad \boldsymbol{Y}=\begin{pmatrix} Y_1\\ Y_2\\ \vdots\\ Y_m\\ \end{pmatrix}\in\mathbb{R}^{m \cdot d},\quad g[\boldsymbol{Y}] = \begin{pmatrix} g(Y_1)\\ g(Y_2)\\ \vdots\\ g(Y_m)\\ \end{pmatrix}\in\mathbb{R}^{m \cdot d}\tag{9}\] In order to implement a general \(m\)-stage implicit Runge-Kutta method, for example, an \(m\)-stage Radau method, one employs Newton’s method or variants of it. These methods make use of the Jacobian of the system which is of size \(m d\times md\). Diagonally implicit Runge-Kutta methods (DIRK) are an attractive alternative to reduce this computational cost. This is because in this case, the matrix \(\mathcal{A}\) in 8 is lower triangular, which implies that the Jacobians used by the Newton method or its variants to solve 8 are of dimension \(d\times d\). Furthermore, when all the diagonal coefficients are equal, the method is referred to as singly diagonally implicit (SDIRK), which may offer the potential to further reduce computational costs [1]. A particular diagonal implicit method of interest in this paper is SDIRK4 as presented in [1] and whose Butcher tableau is given by \[\label{eq:sdirk4} \begin{array} {c|ccccc} \\ & \frac{1}{4}\\[1mm] &\frac{1}{2} &\frac{1}{4} \\[1mm] & \frac{17}{50}& \frac{-1}{25}&\frac{1}{4} \\[1mm] & \frac{371}{1360} &\frac{-137}{2720} &\frac{15}{544} &\frac{1}{4}\\[1mm] &\frac{25}{24} &\frac{-49}{48} &\frac{125}{16} &\frac{-85}{12} &\frac{1}{4}\\[1mm] \hline y_1&\frac{25}{24} &\frac{-49}{48} &\frac{125}{16} &\frac{-85}{12} &\frac{1}{4}\\[1mm] \hat{y}_1&\frac{59}{48} &\frac{-17}{96} &\frac{225}{32} &\frac{-85}{12} &0\\[1mm] err&-\frac{3}{16} &\frac{-27}{32} &\frac{25}{32} & 0 &\frac{1}{4} \end{array}.\tag{10}\]

2.3 Explicit stabilised Runge-Kutta methods↩︎

Even with the dimensionality reduction of diagonally implicit Runge-Kutta methods in terms of the Jacobians used for solving the nonlinear equations, implicit methods in high dimensions become prohibitively expensive. On the other hand, explicit schemes avoid the solution of such high-dimensional systems, but have bounded stability domains. For this reason, standard explicit methods face severe step size restrictions when integrating stiff problems. Explicit stabilised Runge-Kutta methods represent a good compromise, they are fully explicit schemes with extended stability domains over the negative real axis.

Explicit stabilised Runge-Kutta methods are constructed to have (nearly) maximal stability domain sizes along the negative real axis constrained by the fact that they remain explicit. In particular, for an \(s\)-stage method the stability polynomial would be of degree \(s\). The goal is then to find, for a given integer \(p\) (the order of the method), a polynomial of degree \(s\) (the number of internal stages of the method) satisfying \[R_s(z)=1+z+\frac{z^2}{2}+\dots+\frac{z^p}{p!}+\alpha_{p+1}z^{p+1}+\dots+\alpha_sz^s,\] such that \(|R_s(z)|\leq1\) for \(z\in [-L_s,0]\) with \(L_s>0\) is as large as possible. For \(p=1\), the optimal value of \(L_s\) is \(2s^2\) [1], [4]. The corresponding polynomial is the shifted Chebyshev polynomial \(R_s(z)=T_s(1+z/s^2)\) where \(T_s(\cos\theta) = \cos(s\theta)\) is the first kind Chebyshev polynomial of degree s. This polynomial oscillates between \(-1\) and \(1\) which introduces singularities in the stability domain and prevents convergence for some eigenvalues. For this reason, we modify the stability function in the following way \[\label{eq:stabcheb} R_s(z)=\frac{T_s(\omega_0+\omega_1z)}{T_s(\omega_0)},\tag{11}\] where \(\omega_0=1+\eta/s^2\), \(\omega_1=T_s(\omega_0)/T_s'(\omega_0)\), and \(\eta\) is called the damping parameter. The polynomial 11 now oscillates between \(-\alpha_s(\eta)\) and \(\alpha_s(\eta)\) for every \(z\in[-L_{s,\eta},-\delta_\eta]\), where \[\label{eq:alphaL} \alpha_s(\eta)=\frac{1}{T_s(\omega_0)}<1, \quad L_{s,\eta}=\frac{1+\omega_0}{\omega_1},\tag{12}\] where \(\delta_\eta\) is a positive real close to zero (\(R_s(0)=1\)). Using the well-known three-term recurrence relation of the first kind Chebyshev polynomials \(T_{j+1}(z)=2zT_j(z)-T_{j-1}(z)\), for \(j\geq2\) (with \(T_0(z)=1\) and \(T_1(z)=z\)), we can derive an efficient and stable implementation (in terms of memory and rounding errors) of the so-called Chebyshev method or the first order Runge-Kutta-Chebyshev (RKC) method with stability function 11 , as analyzed by [17] (also for second order explicit stabilised methods). When applied to non-linear problems where the Jacobian of \(g\) has negative real eigenvalues, this method allows for taking much larger time steps than standard explicit methods would allow. Algorithm 1 summarises one step of the first order \(s\)-stage RKC method with step-size \(\Delta t\) for the solution of 2 .

Note that there are three parameters needed for implementing this Algorithm, the damping \(\eta\), the number of stages \(s\) as well as the step-size \(\Delta t\). For a fixed value of the damping parameter \(\eta\) and a given time step \(\Delta t\) one needs to choose the stage parameter \(s\) in such a way that \(L_{s,\eta} \geq L \Delta t\). Typically for ODE integration one takes \(\eta\) small (e.g \(\eta=0.05\)) which gives a nearly optimal stability domain size \(L_{s,\eta} \geq (2-\frac{4}{3}\eta)s^{2}\), and a stage parameter \[s=\left\lceil \sqrt{\frac{L \Delta t}{2-\frac{4}{3}\eta}}\right\rceil\] where \(L\) is an upper bound for the largest eigenvalue of the Jacobian of \(g\) evaluated at \(y_{n}\). For larger values of damping parameter \(\eta\), the width \(L_{s,\eta}\) of the stability domain reduces further but still grows quadratically with \(s\) and the number of stages \(s\) is then chosen large enough accordingly, see for instance [18] for details.

2.4 Application to optimization↩︎

In the recent paper [15] the authors utilised [eq:rkc1] to solve efficiently stiff optimization problems by an algorithm called Runge-Kutta Chebyshev descent (RKCD). Here, by stiff optimization problems, we mean the case where the Hessian of the objective function has large eigenvalues distributed in the positive half-plane close to the real axis. In particular, consider the problem \[\label{eq:min} \min_{x\in\;\mathbb{R}^d}f(x),\tag{13}\] where the dimension \(d\) is large and \(f:\mathbb{R}^d\to\mathbb{R}\), is \(\ell\)-strongly convex differentiable function, with \(L\)-Lipschitz derivative, and such that the real parts of the eigenvalues of its Hessian matrix \(\mathcal{H}f(x)\) are strictly positive for all \(x\in\mathbb{R}^d\). The condition number of the Hessian is bounded by \(\kappa=\frac{L}{\ell}\).

A starting point for solving the minimization problem 13 is to consider discretizations of the gradient flow \[\label{eq:gf} \dot{x}(t) = -\nabla f(x(t)), \qquad x(0)=x_0 \in \mathbb{R}^d,\tag{14}\] where \(x_0\) is an initialization. Applying an explicit Euler method with step size² \(h\) to 14 is equivalent to the gradient descent method \[\label{eq:gradient32descent} x_{k+1}=x_k-h\nabla f(x_k),\qquad k=0,1,2,...\tag{15}\] where \(x_k\)³ is the numerical approximation of \(x(t_k)\). In the strongly convex setting described above when (\(L\gg l\)), this method faces a severe restriction on \(h\) to be stable. Furthermore, it is known that the optimal \(h\) in terms of convergence speed towards the minimizer is given by \(h=\frac{2}{L+\ell}\) [19]. This implies that the number of iterations needed to find the minimizer of \(f\) scales linearly with the condition number \(\kappa\). Note that this behaviour is suboptimal since the most efficient optimization methods that make use only of the gradient of \(f\) need \(\mathcal{O}(\sqrt{\kappa})\) iterations to find the minimum [19].

To improve the convergence speed of the explicit Euler method one can apply a more sophisticated discretization of the gradient flow 14 . For example, one can use the first order RKC [eq:rkc1] to integrate 14 . In particular, if one has estimates for the smallest and largest eigenvalues of the Hessian matrix \(\nabla^{2} f\), then by choosing the parameters \(\eta, h,\) and \(s\) in [eq:rkc1] appropriately one can obtain a very efficient optimization method called RKCD [15] for strongly convex problems, see also Algorithm 2. In the case of RKCD, a calculation of \(x_{k+1}\) for \(x_{k}\) requires \(s\) gradient evaluation. Now the following proposition from [15] characterizes the speed of convergence of iterates \(x_{k}\) towards the minimizer \(x_{\star}\) when \(f\) is quadratic in 13 .

A straightforward calculation shows that the learning rate \(h\) in [eq:defn32of32parameters1] satisfies \(h= \tfrac{\eta}{\ell} + \mathcal{O}(\eta^2 + s^{-2}),\) for \(\eta\rightarrow 0, s\rightarrow +\infty\) which makes it mildly dependent on \(s\). Precisely, the following estimates holds.

Using this proposition and the properties of \(\alpha_{s}(\eta)\) it is possible to show that the effective rate of convergence of function iterates towards the minimizer \(x_{*}\) scales like \(1-c(\eta)/\sqrt{\kappa}\) when \(\kappa \gg 1\) [15], exhibiting thus accelerated behavior from an optimization point of view. This, in turn, implies that to calculate \(x_{\star}\) with target accuracy \(\varepsilon\) one needs to do \(\mathcal{O}(\sqrt{\kappa}\log{\varepsilon^{-1}})\) gradient evaluations.

3 A new explicit implementation of SDIRK methods↩︎

As discussed in Section 2.2, implicit methods are usually implemented using Newton’s method or its variants which can be memory demanding and delicate for ill-conditioned problems in particular for high-dimensional problems arising from discretizing PDEs. Here, we propose an alternative way of implementing implicit Runge-Kutta methods 6 7 based on explicit stabilised methods. To do this, one needs to find the solution \(\boldsymbol{Y}_{n+1}\) to the system of non-linear equations 8 given the current state \(y_{n}\). Our starting point for achieving this is to introduce the auxiliary system of ODEs, \[\label{eq:gf1} \dot{\boldsymbol{Y}}(t)=-\boldsymbol{Y}(t) +\mathbb{1}_m \otimes y_n + \Delta t (\mathcal{A}\otimes I_d) g[\boldsymbol{Y}(t)]\tag{17}\] that satisfies under appropriate assumptions \[\lim_{t \rightarrow +\infty} \boldsymbol{Y}(t) = \boldsymbol{Y}_{n+1}.\] We have thus reformulated the solution of an \(m\)-stage implicit Runge-Kutta method at time \(t_{n+1}\) as the steady-state solution of the ODE 17 . In the simplest case of the implicit Euler method and assuming a gradient structure \(g(y)=-\nabla V(y)\) in 2 , it is easy to see that 17 also has a gradient structure 14 with \(f(y)=\Delta t V(y)+\frac{1}{2} \| y - y_n\|_2^2.\) In this simple situation, following the ideas in [15] one can apply RKCD directly to compute \(\boldsymbol{Y}_{n+1}\). However, this gradient structure assumption is not satisfied in general for a diagonally implicit Runge-Kutta method, even when 1 has a gradient structure. In particular, 17 in this case becomes \[\label{eq:gff} \dot{\boldsymbol{Y}}(t)=-\boldsymbol{Y}(t) +\mathbb{1}_m \otimes y_n + \Delta t (\mathcal{A}\otimes I_d)(F_{D}[{\boldsymbol{Y}}(t)]+F_{A}[{\boldsymbol{Y}}(t)]+F_{R}[{\boldsymbol{Y}}(t)])\tag{18}\] As we will see, without this gradient assumption, applying Algorithm 2 to 18 would lead to arbitrarily slow convergence to \(\boldsymbol{Y}_{n+1}\) as the stiffness grows. To tackle this issue, and inspired by the partitioning \(\mathcal{A}=\gamma I_m+(\mathcal{A}-\gamma I_m)\) of the Runge-Kutta coefficients into diagonal part and strictly lower diagonal part proposed in [3], we consider the following partition of 18 , \[\label{eq:gff95part} \begin{align} \dot{{\boldsymbol{Y}}}(t) \;=\;& \mathbb{1}_m \otimes y_n - {\boldsymbol{Y}}(t) \;+\; \Delta t\, \gamma\, I_m \otimes I_d \, F_D[ {\boldsymbol{Y}}(t)] &&=:~G_D( {\boldsymbol{Y}}(t))\\ &\quad +~ \Delta t\, (\mathcal{A}-\gamma I_m)\otimes I_d \, F_D[ {\boldsymbol{Y}}(t)] \;+\; \Delta t\, \mathcal{A}\otimes I_d \, F_A[{\boldsymbol{Y}}(t)] &&=:~G_A( {\boldsymbol{Y}}(t))\\ &\quad +~ \Delta t\,\mathcal{A}\otimes I_d \, F_R[ {\boldsymbol{Y}}(t)] &&=:~G_R( {\boldsymbol{Y}}(t)), \end{align}\tag{19}\] where we split the coefficient matrix \(\mathcal{A}\) of the singly diagonally implicit Runge-Kutta method into its diagonal part \(\gamma I_m\) and its strictly lower triangular part \(\mathcal{A}-\gamma I_m\). The new Algorithm [alg:3] is then defined as follows.

In Algorithm [alg:3], we set \(\ell=1\) for simplicity, as the spectral gap for the Jacobian of \(G_D\). Note that if \(\lambda_{\min}\geq 0\) denotes the smallest eigenvalue of the Jacobian of \(-F_D\) one could use the more precise value \(\ell=1+\Delta t\lambda_{\min}\). Concerning the choice of the damping parameter \(\eta\), while the analysis applies for any value \(\eta>0\), it was shown numerically in [15] that the value \(\eta=1.17\) is an efficient choice for RKCD. For Algorithm [alg:3], we choose \(\eta=4\) which reveals efficient in the considered numerical test problems (see Section 5). Concerning the stopping criteria for the main iteration of Algorithm [alg:3], we use for simplicity the condition \(\|\boldsymbol{x}_k-\boldsymbol{x}_{k-1}\|<tol_{it}\), with \(tol_{it}=10^{-12}\). For the initialisation, we use for simplicity \(\boldsymbol{x}_0=\mathbb{1}_m\otimes y_n\) where \(y_n\) is the current state.

4 Convergence and stability analysis↩︎

In this section, we prove the convergence of Algorithm [alg:3] in two special cases, pure linear diffusion and linear advection-diffusion. Additionally, we provide the rationale behind the partition in equation 19 for a general advection-diffusion-reaction equation.

4.1 The case of pure linear diffusion↩︎

We consider \[\label{eq:heat} \partial_t u=a\Delta u, \quad (t,x)\in[0,+\infty)\times(0,1)^{N}\tag{20}\] with periodic boundary conditions and a smooth initial condition \(u(x,0)=1+e^{-50\left\lVert x-c \right\rVert^2}\) where \(c\) is the center of the spatial domain. We now discretise in space using a finite difference method with a uniform mesh size \(\Delta x\) in each spatial dimension. We thus have \[\label{eq:pure95diffusion} F_{D}(y)= -Dy, \quad F_{A}(y)=0 \quad F_{R}(y)=0\tag{21}\] in 1 with \(D \in \mathbb{R}^{d\times d}\) is the discrete Laplacian with periodic boundary conditions. In this case, we have that the bound for the spectral radius of \(\nabla F_{D}(y)=-D\) is given by \(\lambda_{\text{heat},max}=4 N a /\Delta x^{2}\), for \(N=1,2,3\).

4.1.0.1 Implicit Euler

We consider the solution of 1 by implicit Euler when the different terms are given by 21 . In this case, we have that \(m=1\) and \(\mathcal{A}=\gamma=1\) in 18 implying that \[\label{pure95diffusion1} G_{D}(y)= \mathbb{1}_m\otimes y_{n}-y-\Delta t D y, \quad G_{A}(y)=0 \quad G_{R}(y)=0\tag{22}\] in 19 . We now use Algorithm 3 to calculate the solution of implicit Euler at time \(t_{n+1}\) given the solution \(y_{n}\) at time \(t_{n}\). Note that in this case, because \(G_{A}(y)=G_{R}(y)=0\) and \(m=1\) 18 has a gradient structure and Algorithm [alg:3] coincides with Algorithm 2. One would thus expect from Proposition 3 that given \(y_{n}\), the cost to approximate \(y_{n+1}\) scales proportionally to the \(\sqrt{\kappa}\) where \(\kappa\) is the condition number of \(G_{D}\) given by \(\kappa=1+\gamma \Delta t \lambda_{\text{heat},max}=1+4Na/\Delta x^{2}\). To illustrate this, we now vary \(\kappa\) by changing the size of the mesh refinement \(\Delta x= 1/4,1/8,1/16,1/32,1/64\) while keeping \(\Delta t=5 \times 10^{-2}\). The results of this experiment are presented in Figure 4 for \(N=1\) and \(N=3\) spatial dimensions, and show that the cost in terms of gradient evaluations indeed scales like \(\mathcal{O}(\sqrt{\kappa})\) independently of the value of \(N\). Here we have used a prescribed tolerance \(\varepsilon=10^{-10}\) in the criteria for stopping the iterations.

4.1.0.2 SDIRK4

We now consider the solution of 1 by SDIRK4 when the different terms are given by 21 . In this case \(m=5\) and the Runge-Kutta matrix \(\mathcal{A}\) is given by 10 , hence \(\gamma=1/4\) and \[\label{eq:pure95diffusion95sdirk} \begin{align} G_{D}(\mathbf{Y}) &= \mathbb{1}_{m}\otimes y_{n}-\mathbf{Y}-\frac{1}{4}\Delta t (I_m \otimes D) \mathbf{Y}, \\ G_{A}(\mathbf{Y})&= \Delta t \left[\left(\mathcal{A}-\frac{1}{4} I_{m} \right)\otimes D \right]\mathbf{Y} \\ G_{R}(\mathbf{Y})&=0. \end{align}\tag{23}\]

Unlike the implicit Euler method, we see that in the case of SDIRK4 Algorithm 3 doesn’t coincide with Algorithm 2, hence Proposition [prop:quad95RKC] does not directly apply. Nevertheless, as we observe⁴ in Figure 4 given \(y_{n}\), the cost to approximate \(y_{n+1}\) with Algorithm [alg:3] scales again proportionally to \(\mathcal{O}(\sqrt{\kappa})\) independently of the value of \(N\), where \(\kappa\) is the condition number of \(G_{D}\). We will prove in Theorem 7 that this is indeed a property of Algorithm [alg:3] for the implementation of an \(m\)-stage SDIRK method in the pure linear diffusion case, equation 21 .

4.1.1 Why the partitioning?↩︎

We now discuss in detail why the partitioning ?? is necessary for the construction of Algorithm [alg:3] and why applying Algorithm 2 without the gradient assumption to 18 would fail in general, even in the simplest possible case of pure linear diffusion where \(F_{D}, F_{R}, F_{A}\) are given by 21 . More precisely, then 18 becomes \[\label{eq:naive95sdirk4} \dot{ \boldsymbol{Y}}(t)= \mathbb{1}_{m} \otimes y_{n} +\mathcal{V} \boldsymbol{Y}(t), \quad \mathcal{V}\mathrel{\vcenter{:}}=-I_{md}-\Delta t (\mathcal{A}\otimes D).\tag{24}\] Note that because \(\mathcal{A}\) is lower triangular with all entries in the diagonal equal to \(\gamma\) when \(m \geq 2\) neither \(\mathcal{A}\) nor \(\mathcal{V}\) can be put in diagonal form. In addition, the eigenvalues of \(\mathcal{V}\) all have multiplicity \(m\) and are given by \[\lambda_{j}=-1-\Delta t \gamma \lambda_{\text{heat},j}\] where \(\lambda_{\text{heat},j}\) are the eigenvalues of \(D\), and the corresponding eigenspace associated to \(\lambda_{j}\) has dimension \(1\).

We start by considering the case of Implicit Euler (\(\gamma=1\), \(m=1\)) for \(\Delta t=5\times 10^{-2}, N=1\). Since now \(G_{A}=G_{R}=0\) as discussed above Algorithm [alg:3] coincides with Algorithm 2 and thus Proposition 3 applies. This is indeed the case up to round-off errors as we can see in Figure 5 where we plot the error between the iterations \(x_{k}\) for Algorithm [alg:3] and the exact solution \(\boldsymbol{Y}_{n+1}\) calculated using standard Newton iteration.

If we now consider the case of SDIRK4 (\(\gamma=1/4\) and \(m=5\)), \(G_{A} \neq 0\) and then 18 even though linear in \(\boldsymbol{Y}\) is not of gradient form. As we can see in Figure 5, and we also prove in Theorem 7, the partitioning ?? is essential for ensuring the convergence of iterates, since the iterates calculated with the non-partitioned Algorithm 2 (corresponding to Algorithm [alg:3] with the definition of \(G(\boldsymbol{x})\) replaced by \(G(\boldsymbol{x})=G_D(\boldsymbol{x})+G_A(\boldsymbol{x})\) and with \(G_R(\boldsymbol{x})=0\), or equivalently applying Algorithm 2 to equation 18 with \(F_A=F_R=0\)) suffer from high-round off errors that deteriorate as stiffness increases (see Appendix 7 for details). In contrast, iterates calculated with Algorithm [alg:3] do not suffer from these issues.

We explain this difference in depth in Section 7, but in a nutshell, the issue is caused by the fact that the matrix \(\mathcal{V}\) in 24 is non-symmetric when \(m \geq 2\) which creates issues for the non-partitioned method. In the case of the partitioning ?? , we obtain \[G_{D}(\boldsymbol{Y})=\mathbb{1}_m \otimes y_n +\mathcal{V}_{1}\boldsymbol{Y},\qquad G_{A}(\boldsymbol{Y})=\mathcal{V}_{2}\boldsymbol{Y},\qquad G_{R}(\boldsymbol{Y}) =0,\] where the Jacobian \(\mathcal{V}\) is decomposed in two parts \(\mathcal{V}=\mathcal{V}_{1}+\mathcal{V}_{2}\) as \[\begin{align} \mathcal{V}_{1} &=& -I_{m}\otimes \hat{\mathcal{V}}_{1}, \quad \hat{\mathcal{V}_{1}}= (I_{d}+\gamma \Delta t D), \\ \mathcal{V}_{2} &=&-\left(\mathcal{A}-\gamma I_{m} \right)\otimes \hat{\mathcal{V}}_{2}, \quad \hat{\mathcal{V}}_{2}=-\Delta t D. \end{align}\] By construction, \(\mathcal{V}_{1}\) is symmetric and has the same eigenvalues as \(\mathcal{V}\), while \(\mathcal{V}_{2}\) is non-symmetric and nilpotent.

The iterates \(\boldsymbol{x}_{k}\) for Algorithm [alg:3] satisfy the following recurrence relation \[\begin{align} \boldsymbol{x}_{k+1} & = & R_s(h \mathcal{V}_{1}) \boldsymbol{x}_k + hB_s(h \mathcal{V}_1) (\mathbb{1}_m \otimes y_n + \mathcal{V}_2 \boldsymbol{x}_k) \tag{25} \nonumber \\ &=& \boldsymbol{x}_k+ h B_s(h \mathcal{V}_{1}) ( \mathbb{1}_m \otimes y_n + \mathcal{V}\boldsymbol{x}_{k}) \tag{26} \end{align}\] where \(B_s(z)=(R_s(z)-1)/z\), and we used \(R_s(z)=1+B_s(z)z\). It is important to note that the steady state of 26 coincides with the steady state of 24 as desired. Furthermore, as we can see in 25 , only the symmetric part of the matrix \(\hat{\mathcal{V}}\) is taken as an argument of the stability function \(R_{s}\), while the non-diagonalizable part \(\hat{\mathcal{V}}_2\) is stabilised by \(B_s\). This is, in fact, key for establishing the convergence properties of Algorithm [alg:3] as detailed in the next section.

4.1.2 Theoretical analysis of Algorithm 3↩︎

In the proof of Theorem 7, observe that one can improve the estimate 29 for the quantity \(B(h\hat{\mathcal{V}_1})h \hat{\mathcal{V}}_2\), but this is not useful here because \(\varepsilon\) can be chosen arbitrarily small.

4.2 The case of linear advection-diffusion↩︎

We now consider the following linear advection-diffusion equation \[\label{eq:advdiffexact} \partial_t u=a\Delta u+b \cdot \nabla u, \quad (t,x)\in [0,+\infty)\times (0,1)^{N}\tag{30}\] with periodic boundary conditions and a smooth initial condition \(u(x,0)\), where \(a>0\) and \(b\in \mathbb{R}^N_+\) are fixed. We now discretise in space using a finite difference method with a uniform mesh size \(\Delta x\) in each spatial dimension. We thus have \[\label{eq:adv95diff} F_{D}(y)= -Dy, \quad F_{A}(y)=E y \quad F_{R}(y)=0.\tag{31}\] where \(D\) is the discrete Laplacian with periodic boundary conditions and \(E\) is the discrete centred advection operator. Following the presentation in [10] the eigenvalues \(\lambda_{\text{heat-adv},j}=\lambda_{\text{heat},j}+i\mu_{\text{adv},j}\) of \(F=F_D+F_A\) are located on the following ellipse in the complex plane \[\label{eq:ellipse95pq} \left(\frac{2p}{\alpha}+1\right)^2 + \left(\frac{q}{\beta}\right)^2=1,\tag{32}\] with half-height \(\beta=\left\lVert b \right\rVert_{1}\Delta x^{-1}\) and width \(\alpha=4aN\Delta x^{-2}\), and where \(p=\lambda_{\text{heat},j}\in[-\alpha,0]\) located on the negative real axis are the eigenvalues of \(F_D\) and \(iq=i\mu_{\text{adv},j}\in[-i\beta,i\beta]\) on the imaginary axis are the eigenvalues of \(F_A\).

Before we present a theorem for an \(m\)-stage SDIRK method it is useful to introduce the following function \[\label{eq:stabRIE} \mathcal{R}_{s}(P,Q) = R_s(P) + B_s(P)Q.\tag{33}\] that will play an important role in establishing the stability and convergence properties in the general case. For simplicity, we start by considering the implicit Euler case for which the iteration error for Algorithm [alg:3] satisfies \[\label{eq:iterrIE} x_{k+1}-y_{n+1} = \mathcal{R}_{s}(-hI_d - h\Delta t D,hE)(x_k-y_{n+1})\tag{34}\] and consider the related stability domain \[\label{eq:stab95domain} {\mathcal{S}_\gamma} = \{(p,q) \in \mathbb{R}^2\;;\;R_{\text{stab}}({\gamma} p,q) <1 \}.\tag{35}\] with \(\gamma=1\) for the implicit Euler method and where \[R_{\text{stab}}(p,q)\mathrel{\vcenter{:}}= |R_{s}(-h+h p)|^2 + |B_{s}(-h+h p)|^2|h q|^2\] We first prove the stability property of Algorithm [alg:3] in the simplest case of the implicit Euler method and based on the stability domain 35 . We then show that it also implies the convergence property of Algorithm [alg:3] for an \(m\)-stage SDIRK method.

In proving Proposition 9, we have explicitly taken into account that because of the periodic boundary conditions, the matrices \(D\) and \(E\) are commuting. This will, in general, not hold, for example in the case of Neumann or Dirichlet boundary conditions. It is still possible to establish a stability criterion, as stated in the next proposition.

We now state the general result, that the estimate for implicit Euler generalises for an arbitrary \(m\)-diagonal SDIRK method.

Note that in the absence of advection (\(b=0\)), the estimate of Theorem 11 reduces to the one of Theorem 7 with \(\hat{\rho} \leq \rho=\alpha_s(\eta)\).

We now plot in Figure 6 the stability domain \(\mathcal{S}_\gamma\) given by 35 with \(\gamma=1/4\) for \(s=20\) and \(\eta=4\). This domain covers a large portion of the negative real axis as expected from the properties of the Chebyshev polynomials. On the other hand, the width on the \(q\)-axis near the origin is much narrower. The blue ellipse in Figure 6 indicates the largest ellipse of the form 32 that contains the origin with principal axes parallel to the \(p\)-\(q\) axes and that can be included in the stability domain \(\mathcal{S}_\gamma\). If we denote by \(d_s\) and \(a_s\) respectively the width and the half-height of this ellipse we observe that \(d_{s}\) grows quadratically with the number of stages \(s\) while \(a_{s}\) grows linearly with \(s\) as illustrated in Figure 7. More precisely, we find numerically for this range of parameters that \(d_s\) and \(a_s\) grow like \[\label{eq:ellipsegrowth} d_s = 2s^2,\qquad a_s = 0.62s.\tag{36}\] For comparison, we also plot in in Figure 6 in red dotted lines the largest ellipse of the form 32 that is included in the modified stability domain \[\label{eq:stabdomainmod} \widehat{\mathcal{S}}_\gamma = \{(p,q) \in \mathbb{R}^2\;;\;R_{\text{stab}}(p,q) <1/2 \}\tag{37}\] where the iteration error is at least halved at each iteration using \(\hat{\rho}<1/2\) in ?? . This is different to blue ellipse corresponding to the stability domain \(\mathcal{S}\) in 35 where for \(\hat{\rho}<1\) only the decay of the iteration error is guaranteed. This illustrates the favorable stabilization of the explicit stabilised iteration taking advantage of the diffusion term even for advection terms with relatively large Peclet numbers. Additionally, note that in the case where \({\rho}<1\) in Theorem 11 is an appropriate bound for an SDIRK method, the behaviour of the stability function of the implicit Euler studied in Figures 6-7 dictates the choice of parameters, with the only difference is that the stability domain size is scaled with respect to the diffusion (horizontal \(p\)-axis in Figure 6) by a factor \(\gamma^{-1}\)

5 Numerical tests↩︎

In this section, we present a number of numerical experiments on the Brusselator problem for the SDIRK4 method. More precisely, we study the case of a highly stiff reaction as well as a large advection setting. We denote with exSDIRK4 the proposed explicit stabilised implementation of SDIRK4, with Matlab prototype codes made publicly available at [20], and we compare its performance with the order 2 PIROCK method. Note that established explicit methods of order \(4\) such as ROCK4 [6], are not suitable for the two problems studied here due to the stiff nature of the advection and reaction terms. In all our numerical experiments, the time step \(\Delta t\) is chosen adaptively as described in the remark below.

Throughout the numerical experiments we make the following implementation assumption: the evaluations of \(F_D[\boldsymbol{Y}]\), \(F_A[\boldsymbol{Y}]\), and \(F_R[\boldsymbol{Y}]\) for the \(m\) internal stages, are computed in parallel. This is natural given the block-diagonal structure of the stabilisation operators \(R_s(h\mathcal{V}_1)\) and \(B_s(h\mathcal{V}_1)\), which decouple across stages, and reflects the remark following Theorem 7. For a vector \(v\in R^d\) that represents a discretized function \(f:\mathbb{R}^d\rightarrow\mathbb{R}^d\), we define the discrete \(L^2-\)norm by \[\|v\|_{L^2} = \sqrt{\frac{\sum_{i=1}^{d}v_i^2}{d}}.\]

5.1 1D Brusselator problem with highly stiff reaction↩︎

Consider the following Brusselator problem \[\label{eq:bruss1d} \begin{align} \partial_t u &= \alpha\Delta u + A +u^2v - (B+1)u,\\ \partial_t v &= \alpha\Delta v -u^2v + Bu,\\ u(0,t)&=u(1,t), \quad v(0,t)=v(1,t), \end{align}\tag{38}\] where \(x\in (0,1)\), \(t\in (0,T)\), \(\alpha=0.2\), \(A=1\), \(B=3\times 10^7\), and \(T=1\). We consider the following initial condition \[u(x,0)=1+\sin(2\pi x) \quad \text{and} \quad v(x,0)=3,\] and periodic boundary conditions. We use second order finite differences for the space discretisation with mesh size \(\Delta x =1/200\), and the starting time step \(\Delta t_0=10^{-6}\).

Table 1 and Figure 8 report the cost comparison between exSDIRK4 and PIROCK for the 1D Brusselator problem 38 . For loose tolerances (\(tol = 10^{-1}\) and \(10^{-3}\)), PIROCK requires fewer function evaluations overall; however, exSDIRK4 already delivers a substantially lower error for the same requested tolerance, reflecting its higher order of accuracy. As the tolerance tightens, the balance shifts decisively in favour of exSDIRK4. At \(tol = 10^{-5}\), both methods have a comparable number of \(F_D\) evaluations, but exSDIRK4 requires more than \(4\) times fewer \(F_R\) evaluations and \(6\) times fewer \(F'_R\) evaluations/inversions. This advantage grows rapidly: at \(tol = 10^{-7}\), PIROCK needs approximately \(3\) times more \(F_D\), \(18\) times more \(F_R\), and \(24\) times more \(F'_R\) evaluations/inversions than exSDIRK4. At the tightest tolerance \(tol = 10^{-9}\), PIROCK accumulates 55,888 accepted steps with 2,928 rejections, while exSDIRK4 completes the integration in only 468 accepted steps with 5 rejections—a reduction of more than two orders of magnitude in the number of steps. These trends are clearly visible in Figure 8, which displays the number of \(F_D\), \(F_R\), and \(F_R'\) evaluations as a function of the tolerance. The work-precision curves show that the crossover point beyond which exSDIRK4 is uniformly cheaper occurs around \(tol \approx 10^{-4}\), and the gap widens as the tolerance decreases, which is consistent with the high-order convergence of exSDIRK4.

5.2 2D Brusselator with large advection↩︎

We now consider the following Brusselator model in 2D: \[\label{eq:bruss2d} \begin{align} \partial_t u &= \nu\Delta u + \mu U\cdot \nabla u + A +u^2v - (B+1)u,\\ \partial_t v &= \nu\Delta v + \mu V\cdot \nabla v -u^2v + Bu,\\ \end{align}\tag{39}\] where \(x\in (0,1)^2\), \(t\in (0,T)\), \(\nu=0.1\), \(\mu=2\) \(A=1.3\), \(B=10^7\), \(U=(-0.5,1)^T\), \(V=(0.4,0.7)^T\) and \(T=0.5\). Following [11], we choose \(u(x,0)=22x_2(1-x_2)^{3/2} , \quad v(x,0)=27x_1(1-x_1)^{3/2},\) and periodic boundary conditions. For the space discretisation, we use again second order finite differences with mesh size equals to \(10^{-2}\) in both dimensions and the starting time step \(\Delta t_0=10^{-6}\).

Figure 9 reports analogous results for the 2D Brusselator problem with large advection. The qualitative behaviour is very similar to the 1D case: exSDIRK4 is outperformed by PIROCK for loose tolerances in terms of raw function evaluation counts, but becomes significantly more efficient as the tolerance tightens. The three panels—\(F_D\) evaluations, \(F_A\)/\(F_R\) evaluations, and \(F_R'\) evaluations—all exhibit the same crossover pattern observed in Figure 8, confirming that the efficiency advantage of exSDIRK4 is not an artefact of the one-dimensional setting but carries over to higher-dimensional problems with combined diffusion, advection, and stiff reaction terms.

6 Conclusion↩︎

We introduced an explicit stabilised implementation of singly diagonally implicit Runge-Kutta (SDIRK) methods for advection-diffusion-reaction PDEs, with the fourth-order method SDIRK4 as the primary example. The central idea is to reformulate the implicit Runge-Kutta stage equations 8 as the steady state of an auxiliary ODE 17 , and to compute that steady state by a partitioned Runge-Kutta-Chebyshev iteration (Algorithm [alg:3]) inspired by the RKCD optimisation method of [15]. The partition of the Runge-Kutta coefficient matrix \(\mathcal{A}= \gamma I_m + (\mathcal{A}- \gamma I_m)\) is essential: it isolates a symmetric, block-diagonal stiff part \(G_D\) that governs the stability of the iteration through the stability polynomial \(R_s\), avoiding the non-diagonalisable amplification matrix that would otherwise prevent convergence, as analysed in detail in Section 7. The off-diagonal diffusion and advection contributions enter through the damping polynomial \(B_s\), whose extended stability along the imaginary axis accommodates large advection velocities under a mild CFL-like condition on the iteration step size \(h\) — significantly more permissive than the classical CFL restriction on the physical time step \(\Delta t\). The stiff reaction term is handled implicitly at each inner iteration via a local inversion of \(J = I_{md} - hG'_R(\boldsymbol{x}_k)\), which remains computationally inexpensive because the reaction Jacobian carries no spatial coupling between grid points and is therefore block-diagonal. The result is a method that requires no large-scale Newton solves or global linear algebra, and yet achieves the high-order accuracy and strong stability of a fourth-order \(L\)-stable integrator. The only implicit component is this local reaction Jacobian inversion, in sharp contrast to the dense, high-dimensional Jacobian systems that a standard Newton-based implementation of SDIRK4 would require.

The convergence and stability properties of the new method are rigorously established in two settings. In the pure linear diffusion case (Theorem 7), we prove that Algorithm [alg:3] converges with decay factor \(\alpha_s(\eta)+\varepsilon\) for any \(\varepsilon>0\), independently of the spatial dimension \(d\), and that implementing one step of an \(m\)-stage SDIRK method requires only \(\mathcal{O}(\sqrt{\Delta t}/\Delta x)\) evaluations of \(F_D\) — the same asymptotic cost as a standard first-order explicit stabilised method. In the linear advection-diffusion case (Theorem 11 and Propositions 9–10), we show that convergence is maintained provided the advection-to-diffusion ratio satisfies a mild CFL-like condition, a condition that is significantly more permissive than the one imposed by a classical explicit scheme on the physical time step \(\Delta t\).

Assuming parallel implementation of the function evaluations across the method stages, the numerical experiments on the 1D and 2D Brusselator problems with highly stiff reaction confirms the theoretical predictions and demonstrate a clear advantage over the second-order method PIROCK [11] whenever tight tolerances are required. The method occupies a natural middle ground: more implicit than a pure explicit stabilised method such as PIROCK, but far cheaper than a fully Newton-based implicit solver, and this balance is precisely what drives its efficiency at high accuracy demands.

The present framework relies on the lower-triangular, singly diagonal structure of SDIRK methods to decouple the implicit system stage by stage and to identify the symmetric stabilising part \(G_D\). A natural and challenging next step would be to extend the approach to fully implicit Runge-Kutta methods such as Radau IIA or Gauss-Legendre collocation schemes. A major difficulty for their implementation is their full matrix structure of the Runge-Kutta coefficients \(\mathcal{A}\), and studied for instance in [21], [22] where splitting approaches based on triangular matrices and low rank linear algebra techniques are proposed. In such fully implicit case, the Runge-Kutta matrix \(\mathcal{A}\) is a full, non-triangular matrix with non-real eigenvalues, so the partitioning strategy proposed in this paper and the associated convergence analysis do not apply straightforwardly.

6.0.0.1 Acknowledgments

The authors would like to thank Ernst Hairer and Francesca Mazzia for helpful discussions on the implementation of SDIRK and RADAU methods. G.V. was partially supported by the Swiss National Science Foundation, projects No. 200020_214819, No. 200020_192129 and No. 10009199. I.A. was partially supported by the Swiss National Science Foundation, projects No. P500PT_210968 and No. P5R5PT_222213.

7 Naive implementation↩︎

We now explain in more detail the issues observed in Figure 5 with respect to the naive stabilised explicit implementation of SDIRK4 applied without partitioning (corresponding to Algorithm [alg:3] with the definition of \(G(\boldsymbol{x})\) replaced by \(G(\boldsymbol{x})=G_D(\boldsymbol{x})+G_A(\boldsymbol{x})\) with \(G_R(\boldsymbol{x})=0\), or equivalently applying Algorithm 2 to equation 18 with \(F_A=F_R=0\)). Considering for simplicity only one mode of the discrete linear heat equation, corresponding to the Dahlquist test equation \(\dot{y} =\lambda y\) with \(\lambda=-\Delta t \gamma \lambda_{\text{heat},j}<0\) corresponds to a possibly large parameter \(|\lambda|\gg 1\). We then write one-step recursion induced by Algorithm 1 applied to a SDIRK method with \(m>1\) internal stages, namely \[\mathbf{Y}_{n+1}=R_{s}(-hI_m+ \lambda h\mathcal{A})\mathbf{Y}_{n}\] where \(R_{s}\) given by 11 is the stability function of RKC, and \(h\) is the time step used by the RKC iteration. Consider the Jordan form of the matrix \(\mathcal{A}=P\mathcal{J}P^{-1}\) where \(\mathcal{J}\) is a Jordan block defined in 27 . Using the following formula for any analytic function \(f\) applied to the Jordan matrix \(J\), \[f(\mathcal{J})=\begin{pmatrix}f(\gamma)&f^\prime(\gamma)&\cdots&\frac{f^{(m-1)}(\gamma)}{(m-1)!}\\&f(\gamma)&\ddots&\vdots\\&&\ddots&f^\prime(\gamma)\\0&&&f(\gamma)\end{pmatrix},\] we deduce taking \(f(\mathcal{J})=R(-hI_m+\gamma z \mathcal{J})\) with \(z=h\lambda\) and \(\hat{z}=-h+\gamma z\), \[P^{-1}R_s(-hI_m +\lambda h \mathcal{A})P = \begin{pmatrix} R_s(\hat{z})&R_s^\prime(\hat{z}) \gamma z&\cdots& \frac{R_s^{(m-1)}(\hat{z})(\gamma z)^{m-1}}{(m-1)!} \\&R_s(\hat{z})&\ddots&\vdots\\&&\ddots&R_s^\prime(\hat{z}) \gamma z\\0&&&R_s(\hat{z}) \end{pmatrix}\] and it turns out that the off-diagonal terms \(\frac{R_s^{(j)}(\hat{z})(\gamma z)^{j}}{j!}, 1\leq j\leq m-1\) are arbitrarily large for large \(|z|\) for stiff problems where \(|\lambda| \gg 1\). This is due to the \(z^j\) factors, growing like \(\mathcal{O}(s^{2j})\), and can be shown that the derivatives \(R_s^{(j)}\) of the stability function oscillate with an amplitude of size \(\mathcal{O}(s^{-j})\) which is not sufficient to damp the \(z^j\) terms. We obtain \[\sup_{s\geq 1, \hat{z}\in [-L_{s,\eta},0]} \|R_s(-hI_m+z \mathcal{A})\|=+\infty\] and this makes the convergence factor of the naive implementation without partitioning non-uniformly bounded as the stiffness increases when refining the spatial mesh (\(\Delta x\rightarrow 0\)), as illustrate in Figure 5. Note that in the proof of the main Theorem 7, we used the same Jordan decomposition of the Runge-Kutta matrix \(\mathcal{A}= P\mathcal{J}P^{-1}\), with the main difference in the design of the proposed method with partitioning and its stability analysis that the stability function \(R_s\) of the explicit stabilised Chebyshev method is always applied to a diagonalizable matrix.

References↩︎