Decoupling Runge–Kutta schemes for
elliptic–parabolic problems

Robert Altmann\({}^\dagger\), Abdullah Mujahid\({}^{\star}\), Benjamin Unger\({}^\ddagger\)

January 01, 1970

Abstract

We study the construction and convergence of semi-explicit and iterative decoupling schemes for an elliptic–parabolic problem using higher-order Runge–Kutta methods. For the semi-explicit schemes, which are constructed using a nearby delay system with \(k\) time delays, we establish the convergence of \(k\)th-order Runge–Kutta methods under a weak coupling condition. We develop the convergence analysis by adapting the Fourier stability and perturbation techniques of [Lubich, Ostermann, Math. Comp., 64(210):601–627, 1995]. The key tool is the generating function framework, in which the Runge–Kutta discretization is encoded through an operator-valued function. Stability estimates are then obtained via Parseval’s identity on the unit circle. We further present convergence results for iterative (fixed-stress and undrained-split) higher-order Runge–Kutta schemes. Here, a spectral decomposition of the Schur complement operator is central. Finally, we provide numerical examples to verify the proven convergence results.

Keywords: Runge–Kutta methods, Fourier stability, semi-explicit schemes, iterative decoupling AMS subject classification: 65M12, 65J10

1 Introduction↩︎

This article explores decoupling time integration schemes based on implicit Runge–Kutta (RK) methods for linear elliptic–parabolic problems, including the equations of poroelasticity [1]. Typical applications involve biomechanics, where the human brain and heart are modeled as poroelastic media with multiple fluid networks [2], [3], as well as geomechanics [4].

The well-posedness of the considered elliptic–parabolic problem is studied in [5], its spatial discretization discussed in [6]. To reduce the computational effort of the coupling, several decoupling strategies exist, based on fixed-point iterations combined with an implicit Euler discretization in time [7]–[9]. In [10], fixed-point iterations for operator splitting were combined with higher-order backward differentiation formulae (BDF) up to order \(5\). An alternative approach based on semi-explicit methods, which decouples the system through a time delay approximation, was introduced for the implicit Euler method in [11] and extended to higher-order BDF methods in [12], [13]. Such semi-explicit methods require a certain weak coupling condition, which may be relaxed by the implementation of an additional inner iteration [14], [15].

The convergence analysis for BDF-based semi-explicit schemes presented in [13] relies on the construction of a weighted norm with a symmetric positive definite matrix enabling a telescoping argument (G-stability), cf. [16]. The present article takes a fundamentally different analytical approach. Instead of G-stability, we adapt the Fourier stability and perturbation techniques developed by Lubich and Ostermann in [17] for RK methods of quasi-linear parabolic equations. The core idea is to encode the RK discretization through generating functions and an associated operator-valued resolvent, decompose it via Schur decomposition, and obtain stability estimates through Parseval’s identity on the unit circle.

Extending the framework of [17] from a single parabolic equation to the coupled elliptic–parabolic system with delayed coupling operators is the main analytical challenge tackled in this paper. To be more precise, the coupling introduces a Schur complement operator and the delay approximation modifies the structure of the operator which needs to be inverted. Through a spectral decomposition of the Schur complement and a Rayleigh quotient argument for the coercive diffusion operator, we reduce the invertibility analysis to a scalar condition on the eigenvalues of the Schur complement, bounded by the coupling strength. The analysis then reveals the same critical coupling bounds as in the BDF setting considered in [13].

For iterative RK schemes, we combine the contraction analysis of the iteration with RK consistency estimates. Here again, the spectral decomposition of the Schur complement is crucial to establish the contraction property.

To summarize, the main contributions of this paper are:

  1. Convergence analysis for RK-based semi-explicit decoupling schemes using Fourier stability techniques, establishing convergence of order \(k\) under a weak coupling condition and sufficient spatial regularity.

  2. Unified perspective connecting the BDF analysis of [13] with the RK analysis, showing that both approaches lead to equivalent stability conditions.

  3. Convergence analysis for iterative (fixed-stress and undrained-split) RK decoupling schemes, combining contraction analysis with RK consistency estimates.

The remainder of this article is organized as follows. After this introduction, the abstract model problem is introduced in 2 with the particular example of poroelasticity. The delay approximation, the resulting semi-explicit scheme, and its Fourier stability and convergence analysis are presented in 3. This is followed by the proof of convergence for iterative RK schemes in 4. Finally, we present a numerical study of the convergence results in 5.

Notation↩︎

Throughout the article, we write \(a \lesssim b\) to indicate that there exists a generic constant \(C > 0\), independent of spatial and temporal discretization parameters, such that \(a \leq C b\).

2 Problem Setting and Preliminaries↩︎

We consider a linear elliptic–parabolic system in an \(m\)-dimensional bounded Lipschitz domain \(\Omega \subseteq \mathbb{R}^m\), \(m\in\{2,3\}\), over a time interval \([0,T]\). Let \[\mathcal{V}\vcentcolon=[H_{0}^{1}(\Omega)]^m, \qquad \mathcal{Q}\vcentcolon=H_{0}^{1}(\Omega), \qquad \mathcal{H}_{\scalebox{.5}{\mathcal{V}}}\vcentcolon=[L^{2}(\Omega)]^m, \qquad \mathcal{H}_{\scalebox{.5}{\mathcal{Q}}}\vcentcolon=L^{2}(\Omega)\] and denote by \(\mathcal{V}\hookrightarrow\mathcal{H}_{\scalebox{.5}{\mathcal{V}}}\simeq\mathcal{H}_{\scalebox{.5}{\mathcal{V}}}^{*}\hookrightarrow\mathcal{V}^{*}\) and \(\mathcal{Q}\hookrightarrow\mathcal{H}_{\scalebox{.5}{\mathcal{Q}}}\simeq\mathcal{H}_{\scalebox{.5}{\mathcal{Q}}}^{*}\hookrightarrow\mathcal{Q}^{*}\) the associated Gelfand triples [18]. Given source terms \(f\colon[0,T] \to \mathcal{V}^*\) and \(g\colon[0,T] \to \mathcal{Q}^*\) of sufficient regularity, we seek \(u\colon [0,T]\rightarrow\mathcal{V}\) and \(p\colon [0,T]\rightarrow\mathcal{Q}\) such that \[\tag{1} \begin{align} a(u,v) - d(v, p) &= \langle f, v \rangle, \tag{2} \\ d(\dot{u}, q) + c(\dot{p},q) + b(p,q) &= \langle g, q\rangle \tag{3} \end{align} for all test functions~v\in \mathcal{V}, q \in \mathcal{Q}, and for almost every~t \in (0,T]. The initial data \begin{align} u(0) = u^0 \in \mathcal{V}, \qquad p(0) = p^0 \in \mathcal{H}_{\scalebox{.5}{\mathcal{Q}}}\tag{4} \end{align}\] are assumed to be consistent, i.e., to satisfy the consistency condition \[\begin{align} a(u^0,v) - d(v, p^0) = \langle f(0),v\rangle\qquad\text{ for all~v\in \mathcal{V}}. \end{align}\] Throughout the manuscript, we assume the bilinear forms \(a\colon \mathcal{V}\times\mathcal{V}\to \mathbb{R}\), \(b\colon \mathcal{Q}\times\mathcal{Q}\to\mathbb{R}\), and \(c\colon \mathcal{H}_{\scalebox{.5}{\mathcal{Q}}}\times\mathcal{H}_{\scalebox{.5}{\mathcal{Q}}}\to\mathbb{R}\) to be symmetric, continuous, and elliptic. We write \(c_{\mathfrak{a}}\) and \(C_{\mathfrak{a}}\) for the ellipticity and continuity constants of \(\mathfrak{a}\in\{a, b, c\}\), respectively. The norms \[\|\cdot\|_{b} \vcentcolon= b(\cdot,\cdot)^{1/2} \qquad\text{and}\qquad \|\cdot\|_{c} \vcentcolon= c(\cdot,\cdot)^{1/2},\] induced by the bilinear forms \(b\) and \(c\), satisfy \[\begin{align} \tfrac{1}{C_b}\,\Vert \cdot \Vert^{2}_{b} \le \Vert \cdot \Vert_{\mathcal{Q}}^{2} \le \tfrac{1}{c_b}\,\Vert \cdot \Vert^{2}_{b} \qquad\text{and}\qquad \tfrac{1}{C_c}\,\Vert \cdot \Vert^{2}_{c} \le \Vert \cdot \Vert_{\mathcal{H}_{\scalebox{.5}{\mathcal{Q}}}}^{2} \le \tfrac{1}{c_c}\,\Vert \cdot \Vert^{2}_{c}. \end{align}\] The dual norms on \(\mathcal{Q}^*\) and \(\mathcal{V}^*\) are written as \(\Vert \cdot \Vert_{\mathcal{Q}^*}\) and \(\Vert \cdot \Vert_{\mathcal{V}^*}\), respectively. The coupling form \(d\colon\mathcal{V}\times\mathcal{H}_{\scalebox{.5}{\mathcal{Q}}}\to\mathbb{R}\) is bounded, i.e., there exists a constant \(C_d > 0\) such that \(d(u,p) \leq C_d\, \Vert u \Vert_{\mathcal{V}} \Vert p\Vert_{\mathcal{H}_{\scalebox{.5}{\mathcal{Q}}}}\) for all \(u \in \mathcal{V}\), \(p\in \mathcal{H}_{\scalebox{.5}{\mathcal{Q}}}\). Under these assumptions, well-posedness is established in [5].

Example 1 (linear poroelasticity). The quasi-static Biot model [1] of linear poroelasticity with homogeneous Dirichlet boundary conditions fits into the framework 1 . Here, the unknowns are the displacement \(u\colon [0,T]\times\Omega\rightarrow\mathbb{R}^{m}\) and the pore pressure \(p\colon [0,T]\times\Omega\rightarrow\mathbb{R}\), satisfying \[\label{eq:pdes} \begin{align} - \nabla\cdot\sigma(u) + \alpha \nabla p &= {\hat{f}} \qquad\text{in }(0,T]\times\Omega, \label{eq:pdes:a}\\ \partial_{t} \big(\alpha \nabla\cdot u + \tfrac{1}{M} p\big) - \nabla\cdot(\kappa\nabla p) &= {\hat{g}} \qquad\text{in } (0,T]\times\Omega. \label{eq:pdes:b} \end{align}\] {#eq: sublabel=eq:eq:pdes,eq:eq:pdes:a,eq:eq:pdes:b} Therein, \(\sigma(u) = {\mu}\, \big(\nabla u + (\nabla u)^\mathsf{T}\big) + {\lambda}\, (\nabla \cdot u) \mathop{\mathrm{id}}\) is the stress tensor with Lamé coefficients \({\lambda}\) and \({\mu}\), \(\kappa\) denotes the permeability, \(\alpha\) the Biot–Willis coupling coefficient, and \(M\) the Biot modulus. The identification of the bilinear forms with the abstract setting 1 is standard; see, e.g., [6].

We define the coupling strength \[\label{eqn:coupling:strength} \omega\vcentcolon=\frac{C_d^{2}}{c_a c_c},\tag{5}\] which governs the convergence of all decoupling schemes considered in this article, and plays a central role in the upcoming analysis.

We further introduce the operators \(\mathcal{A}\colon\mathcal{V}\rightarrow\mathcal{V}^*\), \(\mathcal{B}\colon\mathcal{Q}\rightarrow\mathcal{Q}^*\), \(\mathcal{C}\colon\mathcal{H}_{\scalebox{.5}{\mathcal{Q}}}\rightarrow\mathcal{H}_{\scalebox{.5}{\mathcal{Q}}}^{*}\), and \(\mathcal{D}\colon\mathcal{V}\rightarrow\mathcal{H}_{\scalebox{.5}{\mathcal{Q}}}^{*}\) associated with \(a\), \(b\), \(c\), and \(d\), respectively. In operator notation, system 1 becomes \[\tag{6} \begin{align} \mathcal{A}u - \mathcal{D}^{*} p &= f \qquad \text{in } \mathcal{V}^{*}, \tag{7}\\ \mathcal{D}{\dot{u}} + \mathcal{C}{\dot{p}} + \mathcal{B}p &= g\text{in } \mathcal{Q}^{*}.\tag{8} \end{align}\] Introducing the vectors and operator matrices \[\label{eq:dae:operators} y \vcentcolon= \begin{bmatrix} u \\ p \end{bmatrix}, \qquad \mathcal{E} \vcentcolon= \begin{bmatrix} 0 & 0 \\ \mathcal{D}& \mathcal{C}\end{bmatrix}, \qquad \mathcal{F} \vcentcolon= \begin{bmatrix} -\mathcal{A}& \phantom{-}\mathcal{D}^* \\ \phantom{-}0 & -\mathcal{B}\end{bmatrix}, \qquad h \vcentcolon= \begin{bmatrix} f \\ g \end{bmatrix},\tag{9}\] we can rewrite 6 in the form \[\label{eq:dae} \mathcal{E}\dot{y} = \mathcal{F}y + h.\tag{10}\] Since \(\mathcal{A}\) is invertible by the ellipticity of \(a\), the displacement can be eliminated from 7 , reducing the system to the parabolic equation \[\begin{align} \label{eq:ppde} (\mathcal{M}+ \mathcal{C})\dot{p} + \mathcal{B}p &= r, \end{align}\tag{11}\] where \(r \vcentcolon= g - \mathcal{D}\mathcal{A}^{-1}\dot{f}\) and \(\mathcal{M}\vcentcolon= \mathcal{D}\mathcal{A}^{-1}\mathcal{D}^{*}\) is the self-adjoint, non-negative Schur complement operator.

2.1 Implicit Runge–Kutta methods↩︎

We now describe the RK time discretization of 10 . The delay approximation that yields the semi-explicit scheme is recalled in 3 below. To construct a numerical approximation of the solution \(y\) of 10 on the time interval \([0,T]\), we rely on \(s\)-stage implicit RK methods, i.e., for a given invertible matrix \(\mathbb{A}\in\mathbb{R}^{s\times s}\) and a vector \(\mathbf{\beta}\in\mathbb{R}^s\), the RK method is given by the Butcher tableau \[\renewcommand{\arraystretch}{1.2} \begin{array}{c|c} \mathbf{\mathbb{\chi}}& \mathbb{A}\\\hline & \mathbf{\beta}^\mathsf{T} \end{array}\qquad \text{with \mathbf{\mathbb{\chi}}\vcentcolon= \mathbb{A}\mathbb{1}},\] where \(\mathbb{1} = [1,\ldots,1]^\mathsf{T}\in \mathbb{R}^s\). We thus use the short notation \((\mathbb{A},\mathbf{\beta})\) to denote a specific RK method. In more detail, we consider RK methods applied to operator equations of the form 10 . Let us consider a time grid \(t^n = n\tau\) with time step \(\tau>0\). Given an approximation \(y^{n-1}\) to \(y(t^{n-1})\), the RK approximation \(y^n\) to the solution of system 10 at time point \(t^n\) is computed in two steps (cf. [19]): In step one, approximations \(\dot{Y}^n_\ell\) of the stage derivatives \(\dot{y}(t^n_\ell)\) at the intermediate stage points \(t^n_\ell \vcentcolon= t^{n-1} + \mathbf{\mathbb{\chi}}_\ell\tau\), \(\ell \in \{1, \ldots, s\}\), are computed from \[\label{eq:rk:stage} \mathcal{E}\dot{Y}^n_\ell = \mathcal{F}Y^n_\ell + h(t^{n}_\ell), \qquad \text{where} \quad Y^n_\ell = y^{n-1} + \tau \sum_{j=1}^s \mathbb{A}_{\ell,j} \dot{Y}^n_j,\tag{12}\] where \(\mathcal{E}\) and \(\mathcal{F}\) are defined in 9 . Then, in the second step, we set \[\label{eq:rk:update} y^n = y^{n-1} + \tau \sum_{\ell=1}^{s} \mathbf{\beta}_\ell \dot{Y}^n_\ell.\tag{13}\] Introducing the compact notation \[\begin{align} Y^n \vcentcolon= \begin{bmatrix} Y^n_1\\ \vdots\\ Y^n_s \end{bmatrix} \qquad\text{and}\qquad \dot{Y}^n \vcentcolon= \begin{bmatrix} \dot{Y}^n_1\\ \vdots\\ \dot{Y}^n_s \end{bmatrix}, \end{align}\] we see that the RK stage derivative approximations satisfy the identity \[\label{eq:rk:stage:deriv} \dot{Y}^n = \tfrac{1}{\tau}\, \mathbb{A}^{-1} \big(Y^n - \mathbb{1} y^{n-1}\big) \qquad\text{with }\; \mathbb{1}y^{n-1} = \begin{bmatrix}y^{n-1}\\\vdots\\y^{n-1}\end{bmatrix}.\tag{14}\] Using the stability function defined as \[\label{eq:stabfun} R(z) = 1 + z\mathbf{\beta}^\mathsf{T}(I_s - z\mathbb{A})^{-1}\mathbb{1},\tag{15}\] we can thus use the vector notation to write the update formula 13 as \[\label{eq:rk:update:stabilityFunction} y^n = R(\infty)y^{n-1} + \mathbf{\beta}^\mathsf{T}\mathbb{A}^{-1}Y^n,\tag{16}\] where \(R(\infty) = 1 - \mathbf{\beta}^\mathsf{T}\mathbb{A}^{-1}\mathbb{1}\). The analysis throughout the paper requires the following stability properties of the RK method.

Assumption 1 (RK method). The \(s\)-stage RK method \((\mathbb{A}, \mathbf{\beta})\) is A-stable, i.e., A(\(\theta\))-stable with \(\theta \geq \pi/2\) (so \(|R(z)| \leq 1\) for \(\mathop{\mathrm{Re}}z \leq 0\)), with \(\mathbb{A}\in \mathbb{R}^{s\times s}\) invertible and \(|R(\infty)| < 1\).

Remark 1. A canonical family satisfying 1 is given by the Radau IIA methods [20]. With \(s = 1, 2, 3\) stages they have classical orders \(k = 1, 3, 5\) and stage orders \(q = 1, 2, 3\), respectively. They are also L-stable (\(R(\infty) = 0\)) and stiffly accurate, so that the update 16 reduces to \(y^n = \mathbf{\beta}^\mathsf{T}\mathbb{A}^{-1} Y^n\). The strict bound \(|R(\infty)| < 1\) in 1 excludes the Gauss–Legendre methods (for which \(|R(\infty)| = 1\)), which we do not analyze in the present work.

Assumption 2 (Resolvent smoothing). Let \(\sigma> 0\) denote the resolvent-smoothing exponent associated with the elliptic operator \(\mathcal{B}\). Its value depends on the operator itself, the boundary conditions, and the spatial dimension. cf. [17]. The solution under consideration is assumed to possess the additional spatial regularity required for the resolvent-smoothing estimates of [17] to apply.

Example 2. For a second-order strongly elliptic operator on a smooth bounded domain with homogeneous Dirichlet boundary conditions and a smooth solution, \(\sigma= 3/4 - \varepsilon\) for arbitrary \(\varepsilon > 0\) in two and three spatial dimensions [17].

Recall that the values produced by the RK method are stage vectors. For any normed space \((\mathcal{X}, \Vert\cdot\Vert_\mathcal{X})\) and stage vector \(V = (V_\ell)_{\ell=1}^s \in \mathcal{X}^s\), we define the stage-product norm via \[\label{eqn:norms} \Vert V \Vert_{\mathcal{X}^s}^2 \vcentcolon= \sum_{\ell=1}^s \Vert V_\ell \Vert_\mathcal{X}^2\tag{17}\] together with the stage duality pairing \(\langle F, V\rangle_s \vcentcolon= \sum_{\ell=1}^s \langle F_\ell, V_\ell\rangle\) between \(F \in (\mathcal{X}^*)^s\) and \(V \in \mathcal{X}^s\). This convention covers \(\mathcal{X}= \mathcal{Q}, \mathcal{V}\) and their duals \(\mathcal{Q}^*, \mathcal{V}^*\) used below. Based on the norms induced by the bilinear forms \(a\) and \(c\), we likewise write \[\Vert V \Vert_{c,s}^2 \vcentcolon= \sum_{\ell=1}^s \Vert V_\ell \Vert_c^2, \qquad \Vert V \Vert_{a,s}^2 \vcentcolon= \sum_{\ell=1}^s \Vert V_\ell \Vert_a^2.\]

3 Semi-explicit Runge–Kutta Decoupling↩︎

This section introduces the semi-explicit RK scheme, based on a delay approximation of the elliptic variable ([subsec:delay,subsec:semiexplicit]), and proves its stability and convergence using a generating-function framework adapted from [17]. The analysis proceeds in four steps:

  1. transform the scheme using generating functions (3.3),

  2. eliminate the elliptic variable and derive a closed equation for \(P(\zeta)\) (3.4),

  3. analyze the resulting operator \(\mathcal{L}(\zeta)\) via spectral arguments (3.4),

  4. transfer stability to the time domain using Parseval’s identity (3.5).

3.1 Delay approximation↩︎

The decoupling strategy of [12], [13] replaces the pressure \(p\) in the elliptic equation 7 by a Lagrange interpolation polynomial of degree \((k-1)\) based on values at the \(k\) preceding time levels, i.e., \[\label{eq:p:approx:lagpol} {\widehat p}(t;\tau) = \sum_{\delta=1}^{k}c_{k,\delta}\,{p}(t - \delta\tau), \qquad c_{k,\delta} = (-1)^{(\delta-1)}\binom{k}{\delta}.\tag{18}\] The central idea for the forthcoming decoupling time integration scheme is to use the time delay \(\tau>0\) in 18 as the time step for the RK method. Substituting 18 into 6 yields the delay system \[\tag{19} \begin{align} \mathcal{A}{\bar u} - \mathcal{D}^{*} \bigg( \sum_{\delta=1}^{k} c_{k,\delta}\, \bar{p}(t - \delta\tau) \bigg) &= f \qquad \text{in } \mathcal{V}^{*}, \tag{20}\\ \mathcal{D}{\dot{\bar u}} + \mathcal{C}{\dot{\bar p}} + \mathcal{B}{\bar p} &= g\text{in } \mathcal{Q}^{*},\tag{21} \end{align}\] in which the elliptic and parabolic equations are decoupled in the following sense. Given the solution \(\bar{p}\) until some time \(t>0\), we can solve 20 for \(\bar{u}\) on the interval \([t,t+\tau]\) independently of the current pressure. This solution can then be used to compute \(\bar{p}\) on the interval \([t,t+\tau]\) from 21 . Eliminating \(\bar{u}\) as before gives the inherent delay parabolic equation \[\label{eqn:par:opt:delay} \mathcal{C}\dot{\bar{p}} + \mathcal{M}\bigg( \sum_{\delta=1}^{k} c_{k,\delta}\, \dot{\bar{p}}(t - \delta\tau) \bigg) + \mathcal{B}\bar{p} = r,\tag{22}\] where \(\mathcal{M}\) is the Schur complement operator introduced in 11 . The approximation error introduced by the delay is controlled by the following result.

Proposition 2 ([13]). Under sufficient smoothness assumptions, the solutions of 6 and 19 satisfy for almost every \(t \in [0,T]\) the estimate \[\Vert {\bar u}(t) - u(t)\Vert_{\mathcal{V}}^{2} + \Vert {\bar p}(t) - p(t)\Vert_{\mathcal{Q}}^{2}\, \lesssim\, \tau^{2k}.\]

3.2 Semi-explicit schemes↩︎

To construct a decoupling time-integration scheme for 1 , respectively the operator formulation 6 , we apply an implicit \(s\)-stage RK method \((\mathbb{A},\mathbf{\beta})\) to the time delay approximation 19 . To this end, we interpret 19 as a system for \(y = (u,p)\), where the elliptic equation 20 acts as constraint, while 21 contains the time derivative. At time step \(n\), we denote the stage values by \[\begin{align} U^n = [U_1^n,\ldots,U_s^n]^\mathsf{T}\in \mathcal{V}^s, \qquad P^n = [P_1^n,\ldots,P_s^n]^\mathsf{T}\in \mathcal{Q}^s \end{align}\] with corresponding stage derivatives \(\dot{U}^n\), \(\dot{P}^n\). Applying the RK method to 19 yields a decoupled scheme at time step \(t^n\). In vector form, this reads \[\tag{23} \begin{align} (I_s \otimes \mathcal{A}) {U}^{n} - (I_s\otimes \mathcal{D}^{*}) \sum_{\delta=1}^{k}c_{k,\delta}{P}^{n-\delta} &= {F}^{n}, \tag{24}\\ (I_s \otimes \mathcal{D}) {\dot{U}}^{n} + (I_s \otimes \mathcal{C}) {\dot{P}}^{n} + (I_s \otimes \mathcal{B}) {P}^{n} &= {G}^{n}. \tag{25} \end{align} Here, \otimes denotes the Kronecker product combining the s-dimensional stage structure with the function-space operators. This means, in particular, that (I_s \otimes \mathcal{A}) acts as \mathcal{A} on each stage independently. Moreover, the right-hand sides consists of F_\ell^n = f(t^{n-1} + \mathbf{\mathbb{\chi}}_\ell\tau) and G_\ell^n = g(t^{n-1} + \mathbf{\mathbb{\chi}}_\ell\tau). Following~\eqref{eq:rk:stage:deriv}, the stage derivatives are related to the stage values by \begin{equation} \tag{26} \dot{U}^n = \tfrac{1}{\tau}\, \mathbb{A}^{-1}\big(U^n - \mathbb{1} u^{n-1}\big), \qquad \dot{P}^n = \tfrac{1}{\tau}\, \mathbb{A}^{-1}\big(P^n - \mathbb{1} p^{n-1}\big). \end{equation}\] The key observation is that the delay approximation eliminates any dependence of 24 on the current stage values \(P^n\). Hence, the stage values \(U^n\) can be computed solely from previous time steps, and system 23 becomes semi-explicit: first solve 24 for \(U^n\), then compute \(P^n\) from 25 and 26 . Following 16 , the update formulae are \[\begin{align} \label{eq:rk:trans} u^{n} = R(\infty) u^{n-1} + \mathbf{\beta}^\mathsf{T}\mathbb{A}^{-1} {U}^{n} \qquad\text{and}\qquad p^{n} = R(\infty) p^{n-1} + \mathbf{\beta}^\mathsf{T}\mathbb{A}^{-1} {P}^{n}. \end{align}\tag{27}\]

Remark 3 (Commutativity of the Schur complement construction and the RK discretization). The semi-explicit scheme 23 is independent of the order in which the RK discretization and elimination of \(u\) variables to construct the Schur complement are performed.

3.3 Generating functions and the \(\Delta(\zeta)\) operator↩︎

Following the (formal) generating power series framework for RK methods introduced in [21], [22], we define \[\begin{gather} u(\zeta) \vcentcolon= \sum_{n=1}^{\infty} u^{n} \zeta^{n}, \qquad p(\zeta) \vcentcolon= \sum_{n=1}^{\infty} p^{n} \zeta^{n} \end{gather}\] and, analogously, \(U\), \(P\), \(F\), and \(G\). The stage-product norm and pairing conventions from 2.1 extend to these generating-function-valued quantities termwise in \(\zeta\). Note that, compared to [17], the summation starts at \(n = 1\) (rather than \(n = 0\)) so that the initial data \(u^0, p^0\) do not appear in the generating functions and are treated separately. The sequences \(u^n, p^n\), etc., are defined by the scheme for \(n = 1,\ldots, N\) with \(N = T/\tau\). For the generating function analysis, we extend the scheme to all \(n > N\) by setting the data to zero, i.e., \(F^n = G^n = 0\) for \(n > N\). This uniquely determines \(P^n, U^n\), etc., for all \(n \geq 1\) and ensures that the algebraic manipulations below hold as identities of formal power series. The RK update formula 27 thus yields \[\begin{align} \label{eq:tran:2} u(\zeta) = \frac{R(\infty)\zeta}{1 - R(\infty)\zeta}\, u^{0} + \frac{\mathbf{\beta}^\mathsf{T}\mathbb{A}^{-1}}{1 - R(\infty)\zeta}\, U(\zeta) \end{align}\tag{28}\] and analogously for \(p(\zeta)\). Following [17], we define the \(\Delta\)-operator \[\label{eq:t:rk:delta} \Delta(\zeta) \vcentcolon= \Big(\mathbb{A}+ \frac{\zeta}{1 - \zeta}\mathbb{1}\mathbf{\beta}^\mathsf{T}\Big)^{-1}\tag{29}\] which encodes the RK structure in a single matrix-valued function of \(\zeta\). As indicated in [17], it satisfies the identities \[\begin{align} \label{eq:delta:identities} \Delta(\zeta) = \mathbb{A}^{-1} - \frac{\zeta\mathbb{A}^{-1}\mathbb{1}\mathbf{\beta}^\mathsf{T}\mathbb{A}^{-1}}{1 - R(\infty)\zeta}, \qquad \frac{\Delta(\zeta)\mathbb{1}}{1-\zeta} = \frac{\mathbb{A}^{-1}\mathbb{1}}{1 - R(\infty)\zeta}. \end{align}\tag{30}\]

Lemma 1 (Spectral property of \(\Delta(\zeta)\)). Assume that the RK method \((\mathbb{A},\mathbf{\beta})\) is \(A(\theta)\)-stable with \(\theta > 0\). Then, for \(|\zeta|\le 1\), all eigenvalues \(\lambda\) of \(\Delta(\zeta)\) satisfy \[\label{eq:DeltaSpec} |\arg \lambda| \le \pi - \theta.\qquad{(1)}\] In particular, for A-stable methods, i.e., \(\theta \ge \pi/2\), all eigenvalues satisfy \(\mathop{\mathrm{Re}}(\lambda) \ge 0\).

Proof. By [17], the eigenvalues of \(\Delta(\zeta)\) are either eigenvalues of \(\mathbb{A}^{-1}\) or satisfy \(R(\lambda) = 1/\zeta\). By the assumed A(\(\theta\))-stability, both classes lie in the sector \(|\arg \lambda| \leq \pi - \theta\), which gives ?? . For \(\theta \ge \pi/2\), the sector \(|\arg \lambda| \le \pi/2\) is contained in the closed right half-plane. ◻

3.4 The operator \(\mathcal{L}(\zeta)\) and its structure↩︎

We derive the transformed system by passing to the generating-function representation of the scheme 23 and eliminating \(U(\zeta)\).

Lemma 2. The generating functions for the semi-explicit RK scheme 23 satisfy \[\label{eq:L:zeta} \mathcal{L}(\zeta) P(\zeta) = \mathcal{R}(\zeta),\qquad{(2)}\] where the operator \(\mathcal{L}(\zeta)\colon \mathcal{Q}^s \to (\mathcal{Q}^*)^s\) is defined by \[\label{eq:L:zeta:def} \mathcal{L}(\zeta) \vcentcolon= (I_s\otimes \mathcal{B}) + \frac{\Delta(\zeta)}{\tau} \otimes \Big(\mathcal{C}+ \Psi_{k}(\zeta)\, \mathcal{M}\Big)\qquad{(3)}\] with \(\Psi_{k}(\zeta) = \sum_{\delta=1}^{k}c_{k,\delta}\, \zeta^{\delta}\) and the right-hand side \(\mathcal{R}(\zeta) \in (\mathcal{Q}^*)^s\) is given by \[\begin{align} \label{eq:R:zeta:def} \mathcal{R}(\zeta) &\vcentcolon= G(\zeta) - \frac{1}{\tau}\Big(\Delta(\zeta)\otimes \mathcal{D}\mathcal{A}^{-1}\Big)F(\zeta) \nonumber\\ &\qquad - \frac{1}{\tau}\Big(\Delta(\zeta)\otimes \mathcal{M}\Big) \sum_{\delta=1}^{k}c_{k,\delta}\sum_{n=1}^{\delta} P^{n-\delta}\zeta^{n} + \frac{1}{\tau}\frac{\mathbb{A}^{-1}\mathbb{1}\zeta}{1 - R(\infty)\zeta}\otimes \big(\mathcal{D}u^{0} + \mathcal{C}p^{0}\big). \end{align}\qquad{(4)}\]

Proof. Using 25 and the stage derivative identity 26 yields \[\begin{align} \label{eq:z:transform:parabolic} \begin{aligned} G(\zeta) &= \frac{1}{\tau}\left( \big(I_s\otimes \mathcal{D}\big) \mathbb{A}^{-1}\right)\left(U(\zeta) - \zeta \mathbb{1}\big(u(\zeta) + u^0\big)\right)\\ &\qquad + \frac{1}{\tau}\left( \big(I_s\otimes \mathcal{C}\big) \mathbb{A}^{-1}\right)\left(P(\zeta) - \zeta \mathbb{1}\big(p(\zeta) + p^0\big)\right) + \big(I_s \otimes \mathcal{B}\big) P(\zeta). \end{aligned} \end{align}\tag{31}\] Substituting the update formula 28 for \(u(\zeta)\) together with \(\Delta(\zeta)\) with the representation given in 30 , we get \[\begin{align} \mathbb{A}^{-1}\left(U(\zeta) - \zeta \mathbb{1}\big(u(\zeta) + u^0\big)\right) &= \Delta(\zeta)U(\zeta) - \frac{\mathbb{A}^{-1} \mathbb{1}\zeta}{1-R(\infty)\zeta}\, u^0 \end{align}\] and, similarly, for terms related to \(P(\zeta)\) and \(p(\zeta)\). Substituting these expressions into 31 yields \[\begin{align} \label{eq:z:transform:parabolic:b} \begin{aligned} G(\zeta) &= \frac{1}{\tau}\big(\Delta(\zeta)\otimes \mathcal{D}\big) U(\zeta) + \frac{1}{\tau}\big(\Delta(\zeta)\otimes \mathcal{C}\big) P(\zeta) + \big(I_s \otimes \mathcal{B}\big) P(\zeta) \\ &\qquad - \frac{1}{\tau}\frac{\mathbb{A}^{-1}\mathbb{1}\zeta}{1-R(\infty)\zeta}\otimes\big(\mathcal{D}u^0 + \mathcal{C}p^0\big), \end{aligned} \end{align}\tag{32}\] where we have used \((I_s\otimes \mathcal{D})(\Delta(\zeta)\otimes \mathrm{Id}) = \Delta(\zeta)\otimes \mathcal{D}\). Next, we eliminate \(U(\zeta)\) by observing that 24 yields \[U(\zeta) = (I_s \otimes \mathcal{A}^{-1})F(\zeta) + (I_s \otimes \mathcal{A}^{-1}\mathcal{D}^{*}) \sum_{\delta=1}^{k}c_{k,\delta}\Big( P(\zeta)\zeta^{\delta} + \sum_{n=1}^{\delta} P^{n-\delta} \zeta^{n}\Big).\] Inserting this into 32 and using \(\mathcal{M}= \mathcal{D}\mathcal{A}^{-1}\mathcal{D}^{*}\), we collect all terms involving \(P(\zeta)\) on the left-hand side, leading to \[\begin{align} \bigg((I_s\otimes \mathcal{B}) &+ \frac{1}{\tau}\big(\Delta(\zeta)\otimes \mathcal{C}\big) + \frac{1}{\tau}\big(\Delta(\zeta)\otimes \mathcal{M}\big)\sum_{\delta=1}^k c_{k,\delta} \zeta^{\delta}\bigg) P(\zeta) \\ &= G(\zeta) - \frac{1}{\tau}\big(\Delta(\zeta)\otimes \mathcal{D}\mathcal{A}^{-1}\big)F(\zeta) - \frac{1}{\tau}\big(\Delta(\zeta)\otimes \mathcal{M}\big) \sum_{\delta=1}^{k}c_{k,\delta}\sum_{n=1}^{\delta} P^{n-\delta} \zeta^{n} \\ &\qquad + \frac{1}{\tau}\frac{\mathbb{A}^{-1}\mathbb{1}\zeta}{1-R(\infty)\zeta}\otimes\big(\mathcal{D}u^0 + \mathcal{C}p^0\big), \end{align}\] which completes the proof. ◻

Compared to the classical case [17], [23], the operator \(\mathcal{L}(\zeta)\) contains the additional perturbation \(\Psi_{k}(\zeta)\mathcal{M}\) originating from the delay approximation. Controlling this term is the key difficulty in the forthcoming analysis.

Remark 4 (Convergence of the generating functions). The identity ?? holds as a formal power series by construction (cf.3.3). Whether the series converges on the unit circle depends on the particular RK method. For L-stable methods (\(R(\infty) = 0\)), identity 30 gives \(\Delta(\zeta) = \mathbb{A}^{-1} - \zeta\,\mathbb{A}^{-1}\mathbb{1}\mathbf{\beta}^\mathsf{T}\mathbb{A}^{-1}\), i.e., a polynomial of degree one. Hence, \(\mathcal{L}(\zeta)\) and \(\mathcal{R}(\zeta)\) are matrix polynomials in \(\zeta\) and convergence is immediate. For \(|R(\infty)| < 1\), the factor \(\zeta/(1-R(\infty)\zeta)\) in ?? has its pole at \(|\zeta| = 1/|R(\infty)| > 1\), so \(\mathcal{R}(\zeta)\) is analytic on the closed unit disc and convergence follows from the uniform bound on \(\mathcal{L}(\zeta)^{-1}\) established in 3 below.

To show the invertibility of \(\mathcal{L}(\zeta)\), we need the following observation regarding the coupling strength \(\omega\) defined in 5 .

Proposition 5 (Coupling bounds). Consider the delay operator \(\Psi_{k}\) from 2. Then \(\mathop{\mathrm{Re}}\bigl(1 + \mu \Psi_{k}(\zeta)\bigr) > 0\) for all \(|\zeta| \le 1\), \(\mu \in [0,\omega]\) if and only if the strict weak coupling condition \[\begin{align} \label{eqn:weakCoupling} \omega < \frac{1}{2^k - 1} \end{align}\qquad{(5)}\] holds. Moreover, \(\mathop{\mathrm{Re}}\bigl(1 + \mu \Psi_{k}(\zeta)\bigr) = 0\) if and only if \(\omega = \mu = 1/(2^k-1)\) and \(\zeta = -1\).

Proof. Recall from 18 that the delay coefficients satisfy \(c_{k,\delta} = (-1)^{\delta-1}\binom{k}{\delta}\) such that the binomial theorem \(\sum_{\delta=0}^{k}(-1)^{\delta}\binom{k}{\delta}\zeta^\delta = (1-\zeta)^k\) implies \[\label{eq:delay:poly:identity} \Psi_{k}(\zeta) = \sum_{\delta=1}^{k} (-1)^{\delta-1}\binom{k}{\delta}\zeta^\delta = 1 - (1 - \zeta)^k.\tag{33}\] Consequently, \(\Psi_{k}\) is analytic and, hence, \(\mathop{\mathrm{Re}}\bigl(1 + \mu \Psi_{k}(\zeta)\bigr) = 1+ \mu\mathop{\mathrm{Re}}\bigl(\Psi_{k}(\zeta)\bigr)\) is harmonic on the unit disc, showing that its minimal value is attained on the boundary of the unit disc. Let \(\theta\in [0,2\pi]\) and define \(\varphi = \tfrac{\theta-\pi}{2}\). Then \(\zeta = \mathrm{e}^{i\theta} = -\mathrm{e}^{2i\varphi}\) and, hence, \[\begin{align} 1 - \mathrm{e}^{i\theta} = 1 + \mathrm{e}^{2i\varphi} = \mathrm{e}^{i\varphi}(\mathrm{e}^{-i\varphi} + \mathrm{e}^{i\varphi}) = 2\cos(\varphi)\, \mathrm{e}^{i\varphi}. \end{align}\] Thus, \(\mathop{\mathrm{Re}}\big((1 - \mathrm{e}^{i\theta})^k\big) = 2^k \cos^k\!(\varphi)\cos(k\varphi) \leq 2^k\) and the maximum is attained at \(\varphi = 0\), translating to \(\mathop{\mathrm{Re}}\big((1 - \zeta)^k\big) = \mathop{\mathrm{Re}}\big((1 - \mathrm{e}^{i\theta})^k\big) = 2^k\) if and only if \(\theta = \pi\). We conclude \[\begin{align} 1+ \mu\mathop{\mathrm{Re}}\bigl(\Psi_{k}(\zeta)\bigr) \geq 1 + \mu\,(1 - 2^k). \end{align}\] The expression is minimized for \(\mu = \omega\), which concludes the proof. ◻

Remark 6 (Connection to the energy-based analysis). The same threshold \(\omega \leq 1/(2^k-1)\) also arises in the G-stability analysis, where it is the condition for the energy identity underlying the summation lemma approach; see [13] for the analogous BDF result. This confirms that the bound is intrinsic to the delay approximation 18 and independent of the proof technique.

Lemma 3 (Invertibility of \(\mathcal{L}(\zeta)\)). Consider the notation from 2 and let the RK method be A-stable. Assume that the weak coupling condition ?? holds such that \[\label{eq:delayCond} \mathop{\mathrm{Re}}\bigl(1 + \mu \Psi_{k}(\zeta)\bigr) > 0 \quad \text{for all } |\zeta| \le 1, \;\mu \in [0,\omega].\qquad{(6)}\] Then the operator \(\mathcal{L}(\zeta)\) from ?? is invertible for all \(|\zeta| \leq 1\) and there exists a constant \(C>0\) independent of \(\tau\) such that \[\sup_{|\zeta| \leq 1} \|\mathcal{L}(\zeta)^{-1}\|_{(\mathcal{Q}^*)^s \to \mathcal{Q}^s} \leq C.\]

Proof. Fix \(\zeta\) with \(|\zeta| \le 1\). We prove uniform invertibility by establishing injectivity and applying a Fredholm argument. Let \(\Delta(\zeta) = V(\zeta)T(\zeta)V(\zeta)^*\) denote a Schur decomposition of \(\Delta(\zeta)\) with unitary \(V(\zeta) \in \mathbb{C}^{s \times s}\) and upper triangular matrix \(T(\zeta) \in \mathbb{C}^{s \times s}\). We study the transformed operator \[\begin{align} \tilde{\mathcal{L}}(\zeta) \vcentcolon= (V(\zeta)^* \otimes \mathrm{Id})\mathcal{L}(\zeta)(V(\zeta) \otimes \mathrm{Id}) = (I_s\otimes \mathcal{B}) + \frac{T(\zeta)}{\tau}\otimes \big(\mathcal{C}+\Psi_{k}(\zeta)\mathcal{M}\big). \end{align}\] Since \(V(\zeta)\) is unitary, \(V(\zeta) \otimes \mathrm{Id}\) is an isometry on \(\mathcal{Q}^s\) (with the product norm). Hence, \(\tilde{\mathcal{L}}(\zeta)\) is injective if and only if \(\mathcal{L}(\zeta)\) is and \(\|\tilde{\mathcal{L}}(\zeta)\|_{\mathcal{Q}^s \to (\mathcal{Q}^*)^s} = \|\mathcal{L}(\zeta)\|_{\mathcal{Q}^s \to (\mathcal{Q}^*)^s}\). Since \(\tilde{\mathcal{L}}(\zeta)\) is block upper triangular on the stage structure, injectivity reduces to injectivity of the diagonal blocks \[\begin{align} \mathcal{Y}(\zeta) \vcentcolon= \mathcal{B}+ \frac{\lambda(\zeta)}{\tau}\big(\mathcal{C}+\Psi_{k}(\zeta)\mathcal{M}\big) \colon \mathcal{Q}\to \mathcal{Q}^*, \end{align}\] where \(\lambda(\zeta)\) is an eigenvalue of \(\Delta(\zeta)\). Let \(q\in\mathcal{Q}\) and assume \(\langle\mathcal{Y}(\zeta)q,q\rangle = 0\). Define the Rayleigh quotient \[\begin{align} \mu(q) \vcentcolon= \frac{\langle \mathcal{M}q,q\rangle}{\|q\|_c^2}\in[0,\omega], \end{align}\] where we exploit that \(\mathcal{M}\) is self-adjoint and non-negative. We thus obtain \[\begin{align} 0 = \langle \mathcal{Y}(\zeta)q,q\rangle = b(q,q) + \frac{\lambda(\zeta)}{\tau}\, \big(1 + \mu(q)\Psi_{k}(\zeta)\big)\,\|q\|_c^2 =\vcentcolon b(q,q) + \frac{\lambda(\zeta)}{\tau}\, w\,\|q\|_c^2. \end{align}\] Since \(b(q,q)\) and \(\|q\|_c^2\) are real, separating real and imaginary parts gives: \[\begin{align} {2} &\text{Im:}\quad & \mathop{\mathrm{Im}}(\lambda w)\,\|q\|_c^2 &= 0, \\ &\text{Re:}\quad & b(q,q) + \tfrac{1}{\tau}\mathop{\mathrm{Re}}(\lambda w)\,\|q\|_c^2 &= 0. \end{align}\] From the imaginary part, either \(q = 0\) (done) or \(\mathop{\mathrm{Im}}(\lambda w) = 0\), i.e., \(\lambda w \in \mathbb{R}\). In the latter case, note that \(\mathop{\mathrm{Re}}(\lambda) \ge 0\) by 1 and \(\mathop{\mathrm{Re}}(w) > 0\) by the strict coupling condition ?? . We distinguish three cases:

  • Case \(\lambda = 0\): The real part equation reduces to \(b(q,q) = 0\), which by coercivity forces \(q = 0\).

  • Case \(\mathop{\mathrm{Re}}(\lambda) > 0\): From \(\mathop{\mathrm{Im}}(\lambda w) = \mathop{\mathrm{Re}}(\lambda)\mathop{\mathrm{Im}}(w) + \mathop{\mathrm{Im}}(\lambda)\mathop{\mathrm{Re}}(w) = 0\), we obtain \(\mathop{\mathrm{Im}}(w) = -\mathop{\mathrm{Im}}(\lambda)\mathop{\mathrm{Re}}(w)/\mathop{\mathrm{Re}}(\lambda)\), so \[\mathop{\mathrm{Re}}(\lambda w) = \mathop{\mathrm{Re}}(\lambda)\mathop{\mathrm{Re}}(w) + \frac{\mathop{\mathrm{Im}}(\lambda)^2\mathop{\mathrm{Re}}(w)}{\mathop{\mathrm{Re}}(\lambda)} = \frac{\mathop{\mathrm{Re}}(w)\,|\lambda|^2}{\mathop{\mathrm{Re}}(\lambda)} \ge 0.\] Together with \(b(q,q) \ge c_b\|q\|_\mathcal{Q}^2 > 0\) for \(q \neq 0\), the real part equation forces \(q = 0\).

  • Case \(\mathop{\mathrm{Re}}(\lambda) = 0\) with \(\lambda \neq 0\): Since \(\mathop{\mathrm{Re}}(w) > 0\), we have \(\mathop{\mathrm{Im}}(\lambda w) = \mathop{\mathrm{Im}}(\lambda)\mathop{\mathrm{Re}}(w) = 0\), which forces \(\mathop{\mathrm{Im}}(\lambda) = 0\), contradicting \(\lambda \neq 0\).

Hence, \(\mathcal{L}(\zeta)\) is injective. Towards surjectivity, we factor \[\begin{align} \mathcal{L}(\zeta) = (I_s\otimes \mathcal{B})\big((I_s\otimes \mathrm{Id}) + \mathcal{K}(\zeta)\big) \end{align}\] with \(\mathcal{K}(\zeta) = \tfrac{1}{\tau}(I_s \otimes \mathcal{B}^{-1})\big(\Delta(\zeta)\otimes (\mathcal{C}+ \Psi_{k}(\zeta)\mathcal{M})\big)\). Since \((\mathcal{C}+ \Psi_{k}(\zeta)\mathcal{M})\colon \mathcal{H}_{\scalebox{.5}{\mathcal{Q}}}\to \mathcal{H}_{\scalebox{.5}{\mathcal{Q}}}^* \subset \mathcal{Q}^*\) and \(\mathcal{B}^{-1}\colon \mathcal{Q}^* \to \mathcal{Q}\), the operator \(\mathcal{K}(\zeta)\) maps \(\mathcal{Q}^s \to \mathcal{Q}^s\) and factors through the compact embedding \(\mathcal{Q}\hookrightarrow\mathcal{H}_{\scalebox{.5}{\mathcal{Q}}}\). Hence \(\mathcal{K}(\zeta)\) is compact on \(\mathcal{Q}^s\). By the Fredholm alternative [24], \((I_s\otimes \mathrm{Id}) + \mathcal{K}(\zeta)\) has index zero on \(\mathcal{Q}^s\), so injectivity implies bijectivity. Since \(\mathcal{B}\) is an isomorphism, \(\mathcal{L}(\zeta)\colon \mathcal{Q}^s \to (\mathcal{Q}^*)^s\) is boundedly invertible.

To conclude the proof, we observe that the mapping \(\zeta \mapsto \mathcal{L}(\zeta)\) is continuous in the operator norm on the compact set \(\{|\zeta|\le 1\}\) and that the inverse exists everywhere. This implies the claimed boundedness independent of \(\tau\). ◻

Remark 7 (Sharpness of the coupling bound). The strict coupling condition ?? is essential for the injectivity argument. At the boundary case \(\omega = 1/(2^k-1)\), by 5, \(\mathop{\mathrm{Re}}(1 + \mu\Psi_{k}(\zeta)) = 0\) occurs at \(\zeta = -1\) and \(\mu = \omega\). In this case, \(\mathop{\mathrm{Re}}(w) = 0\) and the injectivity proof breaks down, as the case \(\mathop{\mathrm{Re}}(\lambda) = 0\) with \(\lambda \neq 0\) can no longer be excluded.

3.5 Stability and convergence↩︎

With the invertibility of \(\mathcal{L}(\zeta)\) established in the previous subsection, we can now derive stability estimates using Parseval’s identity.

Theorem 8 (Stability). Let 1 hold and assume that the weak coupling condition ?? is satisfied. Then the semi-explicit RK scheme 23 satisfies \[\begin{gather} \label{eq:stability} \tau^2 \sum_{n=1}^{N} \|P^n\|_{\mathcal{Q}^s}^2 + \tau\sum_{n=1}^{N} \|P^n\|_{c,s}^2 \\ \leq C \bigg(\tau^2 \sum_{n=1}^{N} \|G^n\|_{(\mathcal{Q}^*)^s}^2 + \sum_{n=1}^{N}\|F^n\|_{(\mathcal{V}^*)^s}^2 + \|u^0\|_{\mathcal{V}}^2 + \|p^0\|_{\mathcal{H}_{\scalebox{.5}{\mathcal{Q}}}}^2 + \sum_{n=-k}^{0} \|P^n\|_{c,s}^2\bigg), \end{gather}\qquad{(7)}\] where \(C\) depends on the method, the coupling parameter, and \(1/(1-|R(\infty)|)\), but is independent of \(\tau\) and \(N\).

Proof. We extend the scheme 23 to all \(n > N\) by setting \(F^n = G^n = 0\), as described in 3.3. By 2, the generating functions then satisfy \[\label{eq:proof:identity} \mathcal{L}(\zeta)\, P(\zeta) = \mathcal{R}(\zeta)\tag{34}\] as an identity of formal power series. Due to the scheme extension, the data generating functions \[G(\zeta) = \sum_{n=1}^{N}G^n\zeta^n, \qquad F(\zeta) = \sum_{n=1}^{N}F^n\zeta^n\] and the initial-value delay term \(\sum_{\delta=1}^{k}c_{k,\delta}\sum_{n=1}^{\delta} P^{n-\delta}\zeta^{n}\) are polynomials. The only singularities in \(\mathcal{R}(\zeta)\) and \(\mathcal{L}(\zeta)\) arise from the rational functions \[\label{eq:proof:rational:terms} \frac{\zeta}{1 - R(\infty)\zeta} \qquad\text{and}\qquad \Delta(\zeta) = \mathbb{A}^{-1} - \frac{\zeta\,\mathbb{A}^{-1}\mathbb{1}\mathbf{\beta}^\mathsf{T}\mathbb{A}^{-1}}{1 - R(\infty)\zeta},\tag{35}\] which have a pole at \(\zeta = 1/R(\infty)\). Since \(|R(\infty)| < 1\) by assumption, this pole lies at \[|\zeta| = \frac{1}{|R(\infty)|} > 1,\] so \(\mathcal{R}(\zeta)\) is analytic on a neighbourhood of the closed unit disc. By 3, \(P(\zeta) = \mathcal{L}(\zeta)^{-1}\mathcal{R}(\zeta)\) is analytic on the closed unit disc and \(\sum_{n \geq 1}\|P^n\|_{c,s}^2 < \infty\).

Testing 34 with \(P(\zeta)\) and taking the real part, the coercivity of \(b\) and the coupling condition ?? yield (cf. the proof of 3) \[\label{eq:proof:coercivity} c_b\,\|P(\zeta)\|_{\mathcal{Q}^s}^2 + \frac{c_0}{\tau}\,\|P(\zeta)\|_{c,s}^2 \leq \mathop{\mathrm{Re}}\langle\mathcal{R}(\zeta), P(\zeta)\rangle_s \qquad\text{for all } |\zeta| \leq 1,\tag{36}\] where \(c_0 > 0\) depends on the method and coupling parameter. Since \(P\) and \(\mathcal{R}\) are analytic on the closed unit disc, integrating 36 over \(|\zeta| = 1\) and applying Parseval’s identity to the left-hand side gives \[\label{eq:proof:integrated} c_b\sum_{n=1}^{\infty}\|P^n\|_{\mathcal{Q}^s}^2 + \frac{c_0}{\tau}\sum_{n=1}^{\infty}\|P^n\|_{c,s}^2 \leq \frac{1}{2\pi}\int_0^{2\pi}\mathop{\mathrm{Re}}\langle\mathcal{R}(\mathrm{e}^{i\theta}), P(\mathrm{e}^{i\theta})\rangle_s\,\mathrm{d}\theta.\tag{37}\] We estimate the right-hand side while keeping \(\mathcal{R}\) in the transform domain. The Cauchy–Schwarz inequality in the \((\mathcal{Q}^*)^s\)–\(\mathcal{Q}^s\) duality applied pointwise on the unit circle yields \[\mathop{\mathrm{Re}}\langle\mathcal{R}(\mathrm{e}^{i\theta}), P(\mathrm{e}^{i\theta})\rangle_s \leq \|\mathcal{R}(\mathrm{e}^{i\theta})\|_{(\mathcal{Q}^*)^s}\,\|P(\mathrm{e}^{i\theta})\|_{\mathcal{Q}^s}.\] Integration, followed by the Cauchy–Schwarz inequality for the integral, then gives \[\frac{1}{2\pi}\int_0^{2\pi}\mathop{\mathrm{Re}}\langle\mathcal{R}(\mathrm{e}^{i\theta}), P(\mathrm{e}^{i\theta})\rangle_s\,\mathrm{d}\theta \leq \bigg(\frac{1}{2\pi}\int_0^{2\pi}\|\mathcal{R}(\mathrm{e}^{i\theta})\|_{(\mathcal{Q}^*)^s}^2\,\mathrm{d}\theta\bigg)^{\!1/2} \bigg(\frac{1}{2\pi}\int_0^{2\pi}\|P(\mathrm{e}^{i\theta})\|_{\mathcal{Q}^s}^2\,\mathrm{d}\theta\bigg)^{\!1/2}.\] Applying Parseval’s identity to the \(P\)-factor and the weighted Young inequality, this is bounded by \[\frac{1}{2c_b}\,\frac{1}{2\pi}\int_0^{2\pi}\|\mathcal{R}(\mathrm{e}^{i\theta})\|_{(\mathcal{Q}^*)^s}^2\,\mathrm{d}\theta + \frac{c_b}{2}\sum_{n=1}^{\infty}\|P^n\|_{\mathcal{Q}^s}^2.\] Substituting into 37 and absorbing the term with constant \(c_b\) into the left-hand side, we obtain \[\frac{c_b}{2}\sum_{n=1}^{\infty}\|P^n\|_{\mathcal{Q}^s}^2 + \frac{c_0}{\tau}\sum_{n=1}^{\infty}\|P^n\|_{c,s}^2 \leq \frac{1}{2c_b}\,\frac{1}{2\pi}\int_0^{2\pi}\|\mathcal{R}(\mathrm{e}^{i\theta})\|_{(\mathcal{Q}^*)^s}^2\,\mathrm{d}\theta,\] which after multiplication by \(\tau\) becomes \[\label{eq:proof:stability:pre} \frac{c_b\tau}{2}\sum_{n=1}^{\infty}\|P^n\|_{\mathcal{Q}^s}^2 + c_0\sum_{n=1}^{\infty}\|P^n\|_{c,s}^2 \leq \frac{\tau}{2c_b}\,\frac{1}{2\pi}\int_0^{2\pi}\|\mathcal{R}(\mathrm{e}^{i\theta})\|_{(\mathcal{Q}^*)^s}^2\,\mathrm{d}\theta.\tag{38}\]

It remains to estimate the right-hand side of 38 . Since \(\mathcal{R}(\zeta)\) in ?? is a sum of four terms, the inequality \(\|a_1+\cdots+a_4\|^2 \leq 4\, (\|a_1\|^2+\cdots+\|a_4\|^2)\) gives \[\label{eq:proof:Rsplit} \frac{1}{2\pi}\int_0^{2\pi}\|\mathcal{R}(\mathrm{e}^{i\theta})\|_{(\mathcal{Q}^*)^s}^2\,\mathrm{d}\theta \leq 4\,\big(T_1 + T_2 + T_3 + T_4\big),\tag{39}\] where each \(T_k\) is defined as a contour integral over the unit circle \(|\zeta|=1\), namely \[\begin{align} T_1 &\vcentcolon= \frac{1}{2\pi}\int_0^{2\pi}\|G(\mathrm{e}^{i\theta})\|_{(\mathcal{Q}^*)^s}^2\,\mathrm{d}\theta, \nonumber\\[4pt] T_2 &\vcentcolon= \frac{1}{2\pi}\int_0^{2\pi}\bigg\|\frac{\Delta(\mathrm{e}^{i\theta})\otimes\mathcal{D}\mathcal{A}^{-1}}{\tau}\,F(\mathrm{e}^{i\theta})\bigg\|_{(\mathcal{Q}^*)^s}^2\,\mathrm{d}\theta, \nonumber\\[4pt] T_3 &\vcentcolon= \frac{1}{2\pi}\int_0^{2\pi}\bigg\|\frac{\Delta(\mathrm{e}^{i\theta})\otimes\mathcal{M}}{\tau}\,\varphi(\mathrm{e}^{i\theta})\bigg\|_{(\mathcal{Q}^*)^s}^2\,\mathrm{d}\theta, \nonumber\\[4pt] T_4 &\vcentcolon= \frac{1}{2\pi}\int_0^{2\pi}\bigg\|\frac{\mathbb{A}^{-1}\mathbb{1}\,\mathrm{e}^{i\theta}}{\tau(1-R(\infty)\mathrm{e}^{i\theta})}\otimes(\mathcal{D}u^0+\mathcal{C}p^0)\bigg\|_{(\mathcal{Q}^*)^s}^2\,\mathrm{d}\theta, \end{align}\] with the delay initial-value polynomial \[\varphi(\zeta) \vcentcolon= \sum_{\delta=1}^{k}c_{k,\delta}\sum_{m=1}^{\delta} P^{m-\delta}\zeta^{m}.\] Exchanging the order of summation (the inner sum contributes to the coefficient of \(\zeta^j\) when \(m = j\) and \(\delta \geq j\)), this polynomial takes the form \[\label{eq:proof:phi:expanded} \varphi(\zeta) = \varphi^1\zeta + \varphi^2\zeta^2 + \cdots + \varphi^k\zeta^k,\tag{40}\] where the \(j\)-th coefficient is \[\label{eq:proof:phi:coeff:early} \varphi^j = \sum_{\delta=j}^{k} c_{k,\delta}\, P^{j-\delta}, \qquad j = 1,\ldots,k.\tag{41}\] In particular, the first and last coefficients read \[\begin{align} \varphi^1 = c_{k,1}\,P^{0} + c_{k,2}\,P^{-1} + \cdots + c_{k,k}\,P^{1-k}, \qquad \varphi^k = c_{k,k}\,P^{0}. \end{align}\] Since \(j - \delta \leq 0\) for every term, each \(P^{j-\delta}\) is an initial value with time index in \(\{-k+1,\ldots,0\}\), and \(\varphi^j = 0\) for \(j > k\).

We now apply Parseval’s identity to each \(T_k\) independently. Since \(G(\zeta)\) is a polynomial of degree \(N\), Parseval’s identity gives \[\label{eq:proof:T1} T_1 = \sum_{n=1}^{N}\|G^n\|_{(\mathcal{Q}^*)^s}^2.\tag{42}\]

For \(T_2\), we first collect the needed operator bounds. From 30 and the triangle inequality, \[\label{eq:proof:Delta:bound} \sup_{\theta\in[0,2\pi]}\|\Delta(\mathrm{e}^{i\theta})\| \leq \|\mathbb{A}^{-1}\| + \frac{\|\mathbb{A}^{-1}\mathbb{1}\|\,\|\mathbf{\beta}^\mathsf{T}\mathbb{A}^{-1}\|}{1-|R(\infty)|} \eqqcolon \frac{C_\Delta}{1-|R(\infty)|},\tag{43}\] where we used \(|1-R(\infty)\mathrm{e}^{i\theta}| \geq 1 - |R(\infty)|\), and \(C_\Delta > 0\) depends only on the RK method. From the coercivity of \(a\) and the continuity of \(d\), \[\label{eq:proof:DA:bound} \|\mathcal{D}\mathcal{A}^{-1}\|_{\mathcal{V}^*\to\mathcal{Q}^*} \,\leq\, \frac{C_d}{c_a}.\tag{44}\] Using the submultiplicativity \(\|(\Delta\otimes\mathcal{D}\mathcal{A}^{-1})\,v\|_{(\mathcal{Q}^*)^s} \leq \|\Delta\|\,\|\mathcal{D}\mathcal{A}^{-1}\|_{\mathcal{V}^*\to\mathcal{Q}^*}\,\|v\|_{(\mathcal{V}^*)^s}\) in the integrand and applying Parseval’s identity to the polynomial \(F\), we obtain \[\begin{align} T_2 &= \frac{1}{\tau^2}\,\frac{1}{2\pi}\int_0^{2\pi}\big\|\big(\Delta(\mathrm{e}^{i\theta})\otimes\mathcal{D}\mathcal{A}^{-1}\big)\,F(\mathrm{e}^{i\theta})\big\|_{(\mathcal{Q}^*)^s}^2\,\mathrm{d}\theta\nonumber\\ &\leq \frac{C_\Delta^2}{\tau^2(1-|R(\infty)|)^2}\,\bigg(\frac{C_d}{c_a}\bigg)^{\!2}\sum_{n=1}^{N}\|F^n\|_{(\mathcal{V}^*)^s}^2, \end{align}\] using 43 , 44 , and Parseval’s identity for \(F\).

For \(T_3\), recall from 40 –41 that \(\varphi(\zeta)\) is a polynomial of degree at most \(k\) with coefficients \(\varphi^j\) depending only on the initial values \(P^{-k+1},\ldots,P^0\). Applying the same operator-bound argument to \(\mathcal{M}= \mathcal{D}\mathcal{A}^{-1}\mathcal{D}^*\) and using \(\omega = C_d^2/(c_a c_c)\), \[\|\mathcal{M}q\|_{\mathcal{Q}^*} \,\leq\, c_c\,\omega\,\|q\|_c \qquad\text{for all } q \in \mathcal{H}_{\scalebox{.5}{\mathcal{Q}}}.\] Using the submultiplicativity \(\|(\Delta\otimes\mathcal{M})\,v\|_{(\mathcal{Q}^*)^s} \leq \|\Delta\|\,c_c\omega\,\|v\|_{c,s}\) in the integrand together with 43 and Parseval’s identity for the polynomial \(\varphi\), we obtain \[\begin{align} T_3 &= \frac{1}{\tau^2}\,\frac{1}{2\pi}\int_0^{2\pi}\big\|\big(\Delta(\mathrm{e}^{i\theta})\otimes\mathcal{M}\big)\,\varphi(\mathrm{e}^{i\theta})\big\|_{(\mathcal{Q}^*)^s}^2\,\mathrm{d}\theta\nonumber\\ &\leq \frac{C_\Delta^2\,(c_c\omega)^2}{\tau^2(1-|R(\infty)|)^2}\sum_{j=1}^{k}\|\varphi^j\|_{c,s}^2. \end{align}\] It remains to estimate the coefficients. By the triangle inequality applied to 41 , \[\|\varphi^j\|_{c,s} \leq \sum_{\delta=j}^{k}|c_{k,\delta}|\,\|P^{j-\delta}\|_{c,s}.\] Squaring with \(\bigl(\sum_{i=1}^{r} a_i\bigr)^2 \leq r\sum_{i=1}^{r} a_i^2\) for \(r = k-j+1\) terms and summing over \(j = 1,\ldots,k\) yields \[\sum_{j=1}^{k}\|\varphi^j\|_{c,s}^2 \leq \sum_{j=1}^{k}(k-j+1)\sum_{\delta=j}^{k}|c_{k,\delta}|^2\,\|P^{j-\delta}\|_{c,s}^2 \leq k\sum_{j=1}^{k}\sum_{\delta=j}^{k}|c_{k,\delta}|^2\,\|P^{j-\delta}\|_{c,s}^2.\] After the change of indices \(\ell = j - \delta \in \{-k+1,\ldots,0\}\), each initial value \(P^\ell\) appears at most \(k\) times in the double sum, with coefficient bounded by \(\max_\delta|c_{k,\delta}|^2\), so \[T_3 \leq \frac{C_\Delta^2\,(c_c\omega)^2\,k^2\max_\delta|c_{k,\delta}|^2}{\tau^2(1-|R(\infty)|)^2}\sum_{\ell=-k+1}^{0}\|P^\ell\|_{c,s}^2.\]

For \(T_4\), since \(\mathcal{D}u^0 + \mathcal{C}p^0 \in \mathcal{Q}^*\) is a fixed vector, the integrand factorises as \[\bigg\|\frac{\mathbb{A}^{-1}\mathbb{1}\,\mathrm{e}^{i\theta}}{\tau(1-R(\infty)\mathrm{e}^{i\theta})}\otimes(\mathcal{D}u^0+\mathcal{C}p^0)\bigg\|_{(\mathcal{Q}^*)^s}^2 = \frac{\|\mathbb{A}^{-1}\mathbb{1}\|^2}{\tau^2\,|1-R(\infty)\mathrm{e}^{i\theta}|^2}\;\|\mathcal{D}u^0+\mathcal{C}p^0\|_{\mathcal{Q}^*}^2,\] where we used \(|\mathrm{e}^{i\theta}| = 1\). Expanding \[\frac{1}{1-R(\infty)\zeta} = \sum_{n=0}^{\infty}R(\infty)^{n}\zeta^{n}\] and applying Parseval’s identity to compute the scalar integral, we obtain \[\label{eq:proof:T4:bound} T_4 = \frac{\|\mathbb{A}^{-1}\mathbb{1}\|^2}{\tau^2}\;\|\mathcal{D}u^0+\mathcal{C}p^0\|_{\mathcal{Q}^*}^2\,\frac{1}{2\pi}\int_0^{2\pi}\frac{d\theta}{|1-R(\infty)\mathrm{e}^{i\theta}|^2} = \frac{\|\mathbb{A}^{-1}\mathbb{1}\|^2}{\tau^2(1-|R(\infty)|^2)}\;\|\mathcal{D}u^0+\mathcal{C}p^0\|_{\mathcal{Q}^*}^2.\tag{45}\] Substituting 42 –45 into 39 and then into 38 , and multiplying by \(\tau\), we arrive at \[\begin{gather} \frac{c_b\tau^2}{2}\sum_{n=1}^{\infty}\|P^n\|_{\mathcal{Q}^s}^2 + c_0\tau\sum_{n=1}^{\infty}\|P^n\|_{c,s}^2 \\ \leq C\bigg(\tau^2\sum_{n=1}^{N}\|G^n\|_{(\mathcal{Q}^*)^s}^2 + \sum_{j=1}^{N}\|F^j\|_{(\mathcal{V}^*)^s}^2 + \sum_{n=-k}^{0}\|P^n\|_{c,s}^2 + \|\mathcal{D}u^0+\mathcal{C}p^0\|_{\mathcal{Q}^*}^2\bigg), \end{gather}\] where \(C > 0\) depends on the method, the coupling parameter, and \(1/(1-|R(\infty)|)\). Restricting the left-hand side to \(n = 1,\ldots,N\) gives ?? . ◻

We now combine the stability estimates with consistency to obtain convergence.

Theorem 9 (Convergence of RK stages). Consider the solution \({\bar p}\) of the delay equation 22 for sufficiently smooth right-hand sides. Let \(P^n\) be the RK stage approximation given by 23 and \(\bar{P}^n\) the exact stage values of the delay solution. Let \(q\) denote the stage order of the RK method (cf. 1). Then we have under the assumptions of 8, \[\label{eq:convergence:rk} \tau^2\sum_{n=1}^{N}\|P^n - \bar{P}^n\|_{\mathcal{Q}^s}^2 + \tau\sum_{n=1}^{N}\|P^n - \bar{P}^n\|_{c,s}^2 \lesssim \tau^{2r} + \sum_{\delta=-k}^{0}\|P^{\delta} - \bar{P}^{\delta}\|_{c,s}^2,\qquad{(8)}\] where the exponent \(r\) is determined by

  1. \(r = \min(k,\,q+1)\) in the general case and

  2. \(r = \min(k,\,q+1+\sigma)\) under 2.

Proof. Let \(E^n \vcentcolon= P^n - \bar{P}^n\) denote the stage error and \(e^n \vcentcolon= p^n - \bar{p}(t^n)\) the grid error. Inserting the exact stage values \(\bar{P}^n\) into the scheme 23 produces a defect \(D^n \in (\mathcal{Q}^*)^s\). For a method of stage order \(q\), the defect satisfies (cf. [17]) \[\|D^n\|_{(\mathcal{Q}^*)^s} \lesssim \tau^{q+1}.\] The error satisfies the operator equation ?? with the defect as right-hand side, i.e., \[\label{eq:proof:error:eq} \mathcal{L}(\zeta)\,E(\zeta) = D(\zeta) + \text{(delay initial-value and initial-data errors)}.\tag{46}\] Here, the defect \(D\) plays the role of the data term \(G\) in ?? , while the elliptic equation 24 is satisfied exactly at each stage, so there is no contribution from the \(F\)-term.

Part (i). Applying 8 to 46 yields \[\begin{align} \label{eq:proof:convergence:bound} \tau^2\sum_{n=1}^{N}\|E^n\|_{\mathcal{Q}^s}^2 + \tau\sum_{n=1}^{N}\|E^n\|_{c,s}^2 &\lesssim \tau^2\sum_{n=1}^{N}\|D^n\|_{(\mathcal{Q}^*)^s}^2 + \sum_{\delta=-k}^{0}\|E^{\delta}\|_{c,s}^2 \nonumber\\ &\lesssim \tau^2 \, N \, \tau^{2(q+1)} + \sum_{\delta=-k}^{0}\|E^{\delta}\|_{c,s}^2 \lesssim \tau^{2q+2}, \end{align}\tag{47}\] where the last step uses \(N\tau = T\) and absorbs initial errors of order \(\tau^{2q}\) or better. Capping the stage-order bound \(\tau^{q+1}\) by the classical order \(\tau^k\) gives ?? .

Part (ii). The basic bound 47 estimates the defect \(D\) as generic data in \((\mathcal{Q}^*)^s\). Under 2, however, \(D\) is itself smoothed by the parabolic component of \(\mathcal{L}(\zeta)^{-1}\), gaining a factor \(\tau^{2\sigma}\) in 47 . This is the resolvent-smoothing mechanism of [17], which closes the gap between the stage-order rate and the full classical order \(k\). In our setting, the elliptic constraint 24 enters through \(\mathcal{D}\) without altering the exponent, which is determined by \(\mathcal{B}\). This yields the improved exponent \(r = \min(k,\,q+1+\sigma)\). ◻

Remark 10 (Stiffly accurate methods and grid values). For stiffly accurate methods, we have \(R(\infty) = 0\) and \(p^n = \mathbf{\beta}^\mathsf{T}\mathbb{A}^{-1}P^n\) by 27 , so the grid error satisfies \[\|p^n - \bar{p}(t^n)\|_c \leq \|\mathbf{\beta}^\mathsf{T}\mathbb{A}^{-1}\|\,\|P^n - \bar{P}^n\|_{c,s}.\]

Corollary 1 (Convergence of semi-explicit RK scheme). For stiffly accurate methods, combining 9 with 2 and 10 via the triangle inequality yields the total error \[\label{eq:total:error:rk} \Vert u^{n} - u(t^{n})\Vert^{2}_{\mathcal{V}} + \Vert p^{n} - p(t^{n})\Vert^{2}_{\mathcal{H}_{\scalebox{.5}{\mathcal{Q}}}} \lesssim \tau^{2r}\qquad{(9)}\] with \(r\) defined in 9.

Remark 11 (Convergence rates for Radau IIA). For Radau IIA-\(s\), \(q = s\) and \(k = 2s-1\), so 1 predicts the baseline orders \(\min(k,\,q+1) = 1,\,3,\,4\) for \(s = 1, 2, 3\). With \(\sigma= 3/4 - \varepsilon\) from 2, the sharpened bound predicts orders \(1\), \(3\), and \(\approx 4.75\) for \(s = 1, 2, 3\); the third rate falls short of the classical order \(k = 5\) by \(1/4\), consistent with the numerically observed rates \(\approx 4.5\)–\(4.8\) in 5.

4 Iterative Runge–Kutta Decoupling↩︎

This section is devoted to iterative decoupling schemes for the elliptic–parabolic system 6 using fixed-stress and undrained-split strategies. Similar to the semi-explicit schemes studied in 3, iterative schemes decouple the fully coupled system 6 by solving the two equations alternatingly: Given approximations \(u^{(i-1)}\), \(p^{(i-1)}\) from the previous iteration, the next iterates \(u^{(i)}\), \(p^{(i)}\) are computed by solving a sequence of two subproblems. The advantage of the iterative decoupling methods is that they use a stabilization parameter to avoid a restriction on the coupling condition as for semi-explicit schemes.

Assumption 3 (RK method for iterative decoupling). In addition to 1, the \(s\)-stage RK method \((\mathbb{A}, \mathbf{\beta})\) is algebraically stable, i.e., \(\mathop{\mathrm{diag}}(\mathbf{\beta})\,\mathbb{A}+ \mathbb{A}^\mathsf{T}\!\mathop{\mathrm{diag}}(\mathbf{\beta}) - \mathbf{\beta}\mathbf{\beta}^\mathsf{T}\succeq 0\), where \(\succeq\) denotes positive semidefiniteness, with weights \(\mathbf{\beta}_\ell > 0\) for all \(\ell\); see [20].

Remark 12. As already discussed in 1, 3 is satisfied by the Radau IIA methods with \(s = 1, 2, 3\) stages.

In the following, we apply an \(s\)-stage RK method \((\mathbb{A}, \mathbf{\beta})\) to an iterative scheme in order to obtain a fully discrete method. On each time interval \([t^{n-1}, t^n]\) of size \(\tau\), we denote the RK stage values at iteration \(i\) by \[U^{n,i} = [U_1^{n,i}, \ldots, U_s^{n,i}]^\mathsf{T}\in \mathcal{V}^s \quad\text{and}\quad P^{n,i} = [P_1^{n,i}, \ldots, P_s^{n,i}]^\mathsf{T}\in \mathcal{Q}^s\] with corresponding stage derivatives \(\dot{U}^{n,i}\), \(\dot{P}^{n,i}\) given by the identity 26 . The iteration at each time step is initialized by \(U^{n,0} = \mathbb{1}\, u^{n-1}\) and \(P^{n,0} = \mathbb{1}\, p^{n-1}\), i.e., all stages are set to the solution from the previous time step.

Since \(\mathcal{C}^{-1}\mathcal{M}\) is self-adjoint, non-negative on \(\mathcal{H}_{\scalebox{.5}{\mathcal{Q}}}\), and compact, the spectral theorem (see, e.g., [24]) yields eigenpairs \((\mu_j, \phi_j)\) with \[\label{eq:spectral:M} \mathcal{M}\phi_{j} = \mu_j\mathcal{C}\phi_{j}, \qquad c(\phi_j,\phi_{j'}) = \langle \mathcal{C}\phi_j,\phi_{j'}\rangle = \delta_{jj'}, \qquad 0 \leq \mu_j \leq \omega.\tag{48}\] This provides an orthonormal basis of \(\mathcal{H}_{\scalebox{.5}{\mathcal{Q}}}\) with respect to the \(c\)-norm. In addition to the unweighted stage-product norms introduced in 2.1, the iterative analysis uses the \(\mathbf{\beta}\)-weighted variants \[\begin{gather} \Vert V \Vert_{a,\mathbf{\beta}}^2 \vcentcolon= \sum_{\ell=1}^s \mathbf{\beta}_\ell\,\Vert V_\ell\Vert_a^2, \qquad \Vert Q \Vert_{c,\mathbf{\beta}}^2 \vcentcolon= \sum_{\ell=1}^s \mathbf{\beta}_\ell\,\Vert Q_\ell\Vert_c^2, \\[2pt] \Vert \theta \Vert_\mathbf{\beta}^2 \vcentcolon= \sum_{\ell=1}^s \mathbf{\beta}_\ell\, \theta_\ell^2, \qquad \langle\theta, \theta'\rangle_\mathbf{\beta}\vcentcolon= \sum_{\ell=1}^s \mathbf{\beta}_\ell\, \theta_\ell\,\theta'_\ell \qquad \text{for } \theta, \theta' \in \mathbb{R}^s. \end{gather}\] The same notation \(\langle f, g\rangle_\mathbf{\beta}\vcentcolon= \sum_{\ell=1}^s \mathbf{\beta}_\ell\,\langle f_\ell, g_\ell\rangle\) is used for the stage-tensor duality pairing on \((\mathcal{X}^*)^s\times\mathcal{X}^s\) for \(\mathcal{X}\in \{\mathcal{V}, \mathcal{Q}\}\). Since \(\mathbf{\beta}_\ell > 0\) by 3, the \(\mathbf{\beta}\)-weighted norms are equivalent to the unweighted stage-product versions with constants \(\min_\ell \mathbf{\beta}_\ell\) and \(\max_\ell \mathbf{\beta}_\ell\).

In the following, we introduce two splitting strategies which yield a contraction, namely fixed-stress (4.1) and undrained-split (4.2). Afterwards, we establish convergence of the two iterative schemes (4.3).

4.1 Fixed-stress splitting↩︎

The fixed-stress splitting decouples the fully coupled system by adding a stabilization term to the flow equation; see [7] for the underlying idea. In each iteration step, the flow equation is solved first for the pressure, followed by the mechanics equation to update the displacement. In the continuous setting, at iteration \(i\), the scheme reads \[\tag{49} \begin{align} \mathcal{A}\, u^{(i)} - \mathcal{D}^{*}\, p^{(i)} &= f, \tag{50}\\ \mathcal{D}\, \dot{u}^{(i-1)} + \mathcal{C}\, \dot{p}^{(i)} + \mathcal{B}\, p^{(i)} + L\, \mathcal{C}\, (\dot{p}^{(i)} - \dot{p}^{(i-1)}) &= g, \tag{51} \end{align}\] where \(L > 0\) is a stabilization parameter. Applying the RK method to 49 , the fixed-stress RK scheme at time step \(n\) and iteration \(i\) reads \[\tag{52} \begin{align} (I_s \otimes \mathcal{A})\, U^{n,i} - (I_s \otimes \mathcal{D}^{*})\, P^{n,i} &= F^{n} \tag{53}\\ (I_s \otimes \mathcal{D})\, {\dot{U}}^{n,i-1} + (I_s \otimes \mathcal{C})\, {\dot{P}}^{n,i} + (I_s \otimes \mathcal{B})\, P^{n,i} + L\, (I_s \otimes \mathcal{C})\, ({\dot{P}}^{n,i} - {\dot{P}}^{n,i-1}) &= G^{n} \tag{54} \end{align}\] in \((\mathcal{V}^{*})^s\) and \((\mathcal{Q}^{*})^s\), respectively. Note that the two equations are decoupled since one can first solve for \(P^{n,i}\) with the second equation.

Theorem 13 (Contraction of fixed-stress RK iteration). Under 3, the fixed-stress iteration 52 satisfies \[\label{eq:rk:contraction:fs} \Vert P^{n,i} - P^{n,i-1} \Vert_{c,\mathbf{\beta}} \leq \rho_{\mathrm{FS}} \, \Vert P^{n,i-1} - P^{n,i-2} \Vert_{c,\mathbf{\beta}}\qquad{(10)}\] with contraction rate \[\label{eq:rk:contraction:fs:rate} \rho_{\mathrm{FS}} = \max_{\mu \in [0,\omega]} \frac{|L - \mu|}{1 + L}.\qquad{(11)}\] In particular, \(\rho_{\mathrm{FS}} < 1\) for all \(L > \max\big(0, \tfrac{\omega-1}{2}\big)\) and the optimal choice \(L = \omega/2\) yields \(\rho_{\mathrm{FS}} = \omega/(2 + \omega)\).

Proof. Define the iterate differences \[\Theta_{u}^{n,i} \vcentcolon= U^{n,i} - U^{n,i-1} \in \mathcal{V}^s \quad\text{and}\quad \Theta_{p}^{n,i} \vcentcolon= P^{n,i} - P^{n,i-1} \in \mathcal{Q}^s\] as well as \({\dot{\Theta}}^{n,i}_{u}\) and \({\dot{\Theta}}^{n,i}_{p}\) accordingly. Subtracting 52 for consecutive iterates yields \[\tag{55} \begin{align} (I_s \otimes \mathcal{A})\, \Theta^{n,i}_{u} - (I_s \otimes \mathcal{D}^{*})\, \Theta^{n,i}_{p} &= 0, \tag{56}\\ (I_s \otimes \mathcal{D})\, {\dot{\Theta}}^{n,i-1}_{u} + (I_s \otimes \mathcal{C})\, {\dot{\Theta}}^{n,i}_{p} + (I_s \otimes \mathcal{B})\, \Theta^{n,i}_{p} + L\, (I_s \otimes \mathcal{C})\, ({\dot{\Theta}}^{n,i}_{p} - {\dot{\Theta}}^{n,i-1}_{p}) &= 0. \tag{57} \end{align}\] Eliminating \(\Theta_u^{n,i}\) via 56 and using the stage derivative formula 14 , in which the \(\mathbb{1}\otimes p^{n-1}\) contribution cancels under the iterate subtraction, we obtain \[\label{eq:rk:fp:compact} \bigg(I_s \otimes \mathcal{B}+ \frac{1+L}{\tau}\,\mathbb{A}^{-1} \otimes \mathcal{C}\bigg)\Theta^{n,i}_{p} = \frac{\mathbb{A}^{-1}}{\tau} \otimes (L\mathcal{C}- \mathcal{M})\,\Theta^{n,i-1}_{p}.\tag{58}\] We expand \(\Theta^{n,i}_p = \sum_j \theta_j^{n,i} \otimes \phi_j\) with coefficients \(\theta_j^{n,i} \in \mathbb{R}^s\) using the eigenpairs \((\mu_j, \phi_j)\) from 48 . Note that this expansion diagonalizes the \(\mathcal{C}\)- and \(\mathcal{M}\)-terms (by orthonormality and \(\mathcal{M}\phi_j = \mu_j\mathcal{C}\phi_j\)), but not the \(\mathcal{B}\)-term.

Multiplying 58 by \(\mathbb{A}\) on the stage structure and testing with \((\mathop{\mathrm{diag}}(\mathbf{\beta})\otimes \mathrm{Id})\,\Theta^{n,i}_p\) in the stage-tensor \(\mathcal{Q}\)-duality pairing gives \[\label{eq:fs:test:scalar} \underbrace{\big\langle (\mathbb{A}\otimes\mathcal{B})\,\Theta^{n,i}_p,\, \Theta^{n,i}_p\big\rangle_\mathbf{\beta}}_{\eqqcolon\, T_\mathcal{B}} + \frac{1+L}{\tau}\,\Vert\Theta^{n,i}_p\Vert_{c,\mathbf{\beta}}^2 = \frac{1}{\tau}\,\big\langle \big(I_s\otimes(L\mathcal{C}- \mathcal{M})\big)\,\Theta^{n,i-1}_p,\, \Theta^{n,i}_p\big\rangle_\mathbf{\beta}.\tag{59}\] Writing out the stage indices, we have \[T_\mathcal{B} = \sum_{\ell,\ell'=1}^s \mathbf{\beta}_{\ell'}\mathbb{A}_{\ell'\ell}\,b(\Theta^{n,i}_{p,\ell}, \Theta^{n,i}_{p,\ell'}) = \mathop{\mathrm{tr}}(\mathop{\mathrm{diag}}(\mathbf{\beta})\mathbb{A}\,B),\] where \(B\) is the Gram matrix of the stage errors in the \(b\)-inner product, i.e., \(B_{\ell\ell'} = b(\Theta^{n,i}_{p,\ell}, \Theta^{n,i}_{p,\ell'})\). Since \(\mathop{\mathrm{tr}}(\mathop{\mathrm{diag}}(\mathbf{\beta})\mathbb{A}\, B) = \mathop{\mathrm{tr}}(\mathrm{sym}(\mathop{\mathrm{diag}}(\mathbf{\beta})\mathbb{A})\, B)\) and \(B \succeq 0\) by the ellipticity of \(b\), the assumed algebraic stability implies \(\mathrm{sym}(\mathop{\mathrm{diag}}(\mathbf{\beta})\mathbb{A}) \succeq \tfrac{1}{2}\mathbf{\beta}\mathbf{\beta}^\mathsf{T}\succeq 0\) and, hence, \(T_\mathcal{B}\geq 0\). Dropping \(T_\mathcal{B}\) and substituting the spectral expansion on both sides yields \[\frac{1+L}{\tau}\, \sum_j \|\theta_j^{n,i}\|_{\mathbf{\beta}}^2 \;\leq\;\frac{1}{\tau}\, \sum_j (L - \mu_j)\, \big\langle\theta_j^{n,i-1}, \theta_j^{n,i}\big\rangle_{\mathbf{\beta}}.\] The Cauchy–Schwarz inequality in the \(\mathbf{\beta}\)-inner product (first on each \(\langle\theta_j^{n,i-1}, \theta_j^{n,i}\rangle_{\mathbf{\beta}}\), then on the sum over \(j\)) gives \[\label{eq:rk:fp:percomp} \Big(\sum_j \|\theta_j^{n,i}\|_{\mathbf{\beta}}^2\Big)^{1/2} \leq\, \frac{\max_{\mu \in [0,\omega]} |L - \mu|}{1 + L}\, \Big(\sum_j \|\theta_j^{n,i-1}\|_{\mathbf{\beta}}^2\Big)^{1/2}.\tag{60}\]

Introducing the contraction rate \(\rho_{\mathrm{FS}}\) as in ?? and noting that \(\Vert \Theta^{n,i}_p \Vert_{c,\mathbf{\beta}}^2 = \sum_j \|\theta_j^{n,i}\|_{\mathbf{\beta}}^2\) via the \(c\)-orthonormality of \(\{\phi_j\}\), estimate 60 becomes \[\label{eq:rk:contraction:fs:final} \Vert \Theta^{n,i}_{p} \Vert_{c,\mathbf{\beta}} \leq \rho_{\mathrm{FS}} \, \Vert \Theta^{n,i-1}_{p} \Vert_{c,\mathbf{\beta}},\tag{61}\] which is ?? . This is a contraction provided \(\rho_{\mathrm{FS}} < 1\), which imposes conditions on the stabilization parameter \(L\):

  • For \(L \geq \omega/2\): \(\max_{\mu \in [0,\omega]} |L - \mu| = L\), giving \(\rho_{\mathrm{FS}} = L/(1+L) < 1\) for all \(L \geq \omega/2\).

  • For \(L < \omega/2\): \(\max_{\mu \in [0,\omega]} |L - \mu| = \omega - L\), giving \(\rho_{\mathrm{FS}} = (\omega - L)/(1+L) < 1\) if and only if \(L > (\omega - 1)/2\).

  • The minimum of \(\rho_{\mathrm{FS}}\) over \(L \geq 0\) is attained at \(L = \omega/2\), yielding \(\rho_{\mathrm{FS}} = \omega/(2 + \omega)\).

In particular, \(\rho_{\mathrm{FS}} < 1\) for all \(L > \max\big(0, \tfrac{\omega-1}{2}\big)\). ◻

Remark 14 (Case \(L = 0\)). Without stabilization, i.e., for \(L=0\), we obtain \(\rho_{\mathrm{FS}} = \omega\). Hence, the iteration is contractive only for \(\omega < 1\).

The proven contraction of the pressure differences in ?? translates – via the elliptic coupling and a stopping criterion argument – into convergence of the displacement as well as the pressure iterates to the fully coupled monolithic RK solution at the optimal order in \(\tau\). The precise statement and proof are given in 4.3 below.

4.2 Undrained-split decoupling↩︎

In contrast to the previous approach, the undrained-split stabilizes the mechanics equation; see [7], [8] for the underlying idea. For this, we define the operator \[\label{eq:us:tildeM} \widetilde{\mathcal{M}} \vcentcolon= \mathcal{D}^*\mathcal{C}^{-1}\mathcal{D}\colon \mathcal{V}\to \mathcal{V}^*,\tag{62}\] where \(\mathcal{C}^{-1}\colon \mathcal{H}_{\scalebox{.5}{\mathcal{Q}}}^* \to \mathcal{H}_{\scalebox{.5}{\mathcal{Q}}}\) is well-defined due to the ellipticity of \(c\). Moreover, we introduce the associated semi-norm \[\label{eq:us:seminorm} \vert v \vert_{\widetilde{\mathcal{M}}}^2 \vcentcolon= \big\langle \widetilde{\mathcal{M}}\, v, v \big\rangle = \big\langle \mathcal{C}^{-1}\mathcal{D}\, v, \mathcal{D}v \big\rangle, \qquad v \in \mathcal{V},\tag{63}\] together with its (unweighted and \(\mathbf{\beta}\)-weighted) stage-product extensions \[\vert V \vert_{\widetilde{\mathcal{M}},s}^2 \vcentcolon= \sum_{\ell=1}^s \vert V_\ell\vert_{\widetilde{\mathcal{M}}}^2, \qquad \vert V \vert_{\widetilde{\mathcal{M}},\mathbf{\beta}}^2 \vcentcolon= \sum_{\ell=1}^s \mathbf{\beta}_\ell\,\vert V_\ell\vert_{\widetilde{\mathcal{M}}}^2 \qquad \text{for } V \in \mathcal{V}^s.\] In the continuous setting, the undrained-split scheme at iteration \(i\) reads \[\begin{align} \mathcal{A}\, u^{(i)} - \mathcal{D}^{*}\, p^{(i-1)} + L\,\widetilde{\mathcal{M}}\, (u^{(i)} - u^{(i-1)}) &= f, \\ \mathcal{D}\, \dot{u}^{(i)} + \mathcal{C}\, \dot{p}^{(i)} + \mathcal{B}\, p^{(i)} &= g, \end{align}\] where \(L > 0\) is again a stabilization parameter. Applying a RK method, the undrained-split RK scheme at time step \(n\) and iteration \(i\) reads \[\tag{64} \begin{align} (I_s \otimes \mathcal{A})\, U^{n,i} - (I_s \otimes \mathcal{D}^{*})\, P^{n,i-1} + L\,(I_s \otimes \widetilde{\mathcal{M}})\, (U^{n,i} - U^{n,i-1}) &= F^{n} & &\text{in } (\mathcal{V}^{*})^s, \tag{65}\\ (I_s \otimes \mathcal{D})\, {\dot{U}}^{n,i} + (I_s \otimes \mathcal{C})\, {\dot{P}}^{n,i} + (I_s \otimes \mathcal{B})\, P^{n,i} &= G^{n} & &\text{in } (\mathcal{Q}^{*})^s. \tag{66} \end{align}\] The contraction property is subject of the following theorem.

Theorem 15 (Contraction of undrained-split RK iteration). Given 3, the undrained-split iteration 64 satisfies \[\label{eq:rk:contraction:us} \big\vert U^{n,i} - U^{n,i-1}\big\vert_{\widetilde{\mathcal{M}},\mathbf{\beta}} \leq \rho_{\mathrm{US}}\, \big\vert U^{n,i-1} - U^{n,i-2}\big\vert_{\widetilde{\mathcal{M}},\mathbf{\beta}}\qquad{(12)}\] with contraction rate \[\label{eq:rk:contraction:us:rate} \rho_{\mathrm{US}} = \frac{\omega\,\max(L,\,1-L)}{1 + L\omega}.\qquad{(13)}\] In particular, \(\rho_{\mathrm{US}} < 1\) for all \(L \geq \tfrac{1}{2}\) (unconditionally) and for \(0 \leq L < \tfrac{1}{2}\) provided \(\omega < 1/(1 - 2L)\). The optimal choice \(L = \tfrac{1}{2}\) yields \(\rho_{\mathrm{US}} = \omega/(2 + \omega)\).

Proof. Set \[\Theta_u^{n,i} \vcentcolon= U^{n,i} - U^{n,i-1} \in \mathcal{V}^s, \qquad \Theta_p^{n,i} \vcentcolon= P^{n,i} - P^{n,i-1} \in \mathcal{Q}^s,\] and define \({\dot{\Theta}}^{n,i}_u, {\dot{\Theta}}^{n,i}_p\) accordingly. Subtracting 64 for two consecutive iterates yields \[\tag{67} \begin{align} (I_s \otimes \mathcal{A}+ L\,I_s \otimes \widetilde{\mathcal{M}})\,\Theta_u^{n,i} &= (I_s \otimes \mathcal{D}^*)\,\Theta_p^{n,i-1} + L\,(I_s \otimes \widetilde{\mathcal{M}})\,\Theta_u^{n,i-1}, \tag{68}\\ (I_s \otimes \mathcal{D})\,{\dot{\Theta}}^{n,i}_u + (I_s \otimes \mathcal{C})\,{\dot{\Theta}}^{n,i}_p + (I_s \otimes \mathcal{B})\,\Theta_p^{n,i} &= 0. \tag{69} \end{align}\] Using the stage-derivative formula 14 in 69 , where the constant contributions \(\mathbb{1} \otimes u^{n-1}\) and \(\mathbb{1} \otimes p^{n-1}\) cancel, and multiplying by \(\tau\mathbb{A}\) on the stage structure gives the compact form \[\label{eq:us:flow:compact} \mathcal{T}\,\Theta_p^{n,i} = -\,(I_s \otimes \mathcal{D})\,\Theta_u^{n,i}, \qquad \mathcal{T}\vcentcolon= I_s\otimes\mathcal{C}+ \tau\mathbb{A}\otimes\mathcal{B}.\tag{70}\] Evaluating 70 at iterate \(i{-}1\) and substituting the resulting expression for \(\Theta_p^{n,i-1}\) into 68 , the mechanics equation becomes a closed equation in \(\Theta_u^{n,i}\), namely \[\label{eq:us:closed} \big(I_s \otimes \mathcal{A}+ L\,I_s \otimes \widetilde{\mathcal{M}}\big)\,\Theta_u^{n,i} = L\,(I_s \otimes \widetilde{\mathcal{M}})\,\Theta_u^{n,i-1} - (I_s\otimes\mathcal{D}^*)\,\mathcal{T}^{-1}\,(I_s\otimes\mathcal{D})\,\Theta_u^{n,i-1}.\tag{71}\] Testing this equation with \((\mathop{\mathrm{diag}}(\mathbf{\beta}) \otimes \mathrm{Id})\,\Theta_u^{n,i}\) in the stage-tensor \(\mathcal{V}\)-duality pairing \(\langle\cdot,\cdot\rangle_\mathbf{\beta}\) and applying \(\langle\widetilde{\mathcal{M}}u, v\rangle = \langle\mathcal{C}^{-1}\mathcal{D}u, \mathcal{D}v\rangle\) on the right-hand side, we obtain \[\label{eq:us:test:scalar} \Vert\Theta_u^{n,i}\Vert_{a,\mathbf{\beta}}^2 + L\,\big\vert\Theta_u^{n,i}\big\vert_{\widetilde{\mathcal{M}},\mathbf{\beta}}^2 = \big\langle \big[L\,(I_s\otimes\mathcal{C}^{-1}) - \mathcal{T}^{-1}\big]\,(I_s\otimes\mathcal{D})\,\Theta_u^{n,i-1},\;(I_s\otimes\mathcal{D})\,\Theta_u^{n,i}\big\rangle_\mathbf{\beta}.\tag{72}\]

Since \(\mathcal{C}^{-1}\mathcal{B}\colon \mathcal{H}_{\scalebox{.5}{\mathcal{Q}}}\to \mathcal{H}_{\scalebox{.5}{\mathcal{Q}}}\) is self-adjoint and positive on the \(c\)-inner product (by the ellipticity of \(b\)), the spectral theorem yields eigenpairs \((\nu_j, \tilde{\phi}_j)\) with \[\label{eq:us:spectral:B} \mathcal{B}\,\tilde{\phi}_j = \nu_j\,\mathcal{C}\,\tilde{\phi}_j, \qquad c(\tilde{\phi}_j, \tilde{\phi}_k) = \delta_{jk}, \qquad \nu_j \geq c_b/C_c > 0.\tag{73}\] We expand \((I_s\otimes\mathcal{D})\Theta_u^{n,i-1}\) and \((I_s\otimes\mathcal{D})\Theta_u^{n,i}\) in the basis \(\{\mathcal{C}\tilde{\phi}_j\}\) of \(\mathcal{Q}^*\) stage-wise, leading to \[(I_s\otimes\mathcal{D})\Theta_u^{n,i-1} = \sum_j \xi_j \otimes \mathcal{C}\tilde{\phi}_j, \qquad (I_s\otimes\mathcal{D})\Theta_u^{n,i} = \sum_j \eta_j \otimes \mathcal{C}\tilde{\phi}_j\] with coefficient vectors \(\xi_j, \eta_j \in \mathbb{R}^s\) whose components are \(\xi_{j,\ell} = \langle\mathcal{D}\Theta^{n,i-1}_{u,\ell}, \tilde{\phi}_j\rangle\) and \(\eta_{j,\ell} = \langle\mathcal{D}\Theta^{n,i}_{u,\ell}, \tilde{\phi}_j\rangle\). Applying \(\mathcal{T}^{-1}\) and \(I_s\otimes\mathcal{C}^{-1}\), respectively, we obtain by 73 \[\mathcal{T}^{-1}\,(I_s\otimes\mathcal{D})\Theta_u^{n,i-1} = \sum_j (I_s + \tau\nu_j\mathbb{A})^{-1}\xi_j \otimes \tilde{\phi}_j, \quad (I_s\otimes\mathcal{C}^{-1})(I_s\otimes\mathcal{D})\Theta_u^{n,i-1} = \sum_j \xi_j \otimes \tilde{\phi}_j.\] The right-hand side of 72 therefore becomes the diagonal sum \[\label{eq:us:rhs:spectral} \big\langle \big[L\,(I_s\otimes\mathcal{C}^{-1}) - \mathcal{T}^{-1}\big]\,(I_s\otimes\mathcal{D})\Theta_u^{n,i-1},\;(I_s\otimes\mathcal{D})\Theta_u^{n,i}\big\rangle_\mathbf{\beta} = \sum_j \big\langle K_j \xi_j,\;\eta_j\big\rangle_\mathbf{\beta},\tag{74}\] where \[K_j \vcentcolon= L I_s - S_j \;\in\; \mathbb{R}^{s\times s}, \qquad \text{with } S_j \vcentcolon= (I_s + \tau\nu_j\mathbb{A})^{-1}.\] We claim that \[\label{eq:us:amp:bound} \Vert K_j\Vert_\mathbf{\beta}\leq \max(L,\,1-L) \qquad \text{for every } j,\tag{75}\] where \(\Vert\cdot\Vert_\mathbf{\beta}\) denotes the matrix norm induced by the \(\mathbf{\beta}\)-weighted inner product on \(\mathbb{R}^s\). The assumed algebraic stability from 3 implies \(\mathrm{sym}(\mathop{\mathrm{diag}}(\mathbf{\beta})\mathbb{A}) \succeq 0\) and, hence, \[\label{eq:us:accretive} x^\mathsf{T}\!\mathrm{sym}(\mathop{\mathrm{diag}}(\mathbf{\beta})\,(I_s + \tau\nu_j\mathbb{A}))\,x = \Vert x\Vert_\mathbf{\beta}^2 + \tau\nu_j\,x^\mathsf{T}\!\mathrm{sym}(\mathop{\mathrm{diag}}(\mathbf{\beta})\mathbb{A})\,x \geq \Vert x\Vert_\mathbf{\beta}^2\tag{76}\] for all \(x\in\mathbb{R}^s\). Inserting \(x = S_j y\) in 76 and using \((I_s + \tau\nu_j\mathbb{A})\,S_j = I_s\) yields \(\langle S_j y,\,y\rangle_\mathbf{\beta}\geq \Vert S_j y\Vert_\mathbf{\beta}^2\). An application of Cauchy–Schwarz then yields \(\Vert S_j y\Vert_\mathbf{\beta}\leq \Vert y\Vert_\mathbf{\beta}\) for every \(y\in\mathbb{R}^s\). Taking the supremum yields the contraction property of the resolvent, namely \[\Vert S_j\Vert_\mathbf{\beta}\leq 1.\] An expansion of the squared \(\mathbf{\beta}\)-norm yields \[\Vert(L I_s - S_j)\,y\Vert_\mathbf{\beta}^2 = L^2 \Vert y\Vert_\mathbf{\beta}^2 - (2L-1)\,\Vert S_j y\Vert_\mathbf{\beta}^2 - 2L\,\big(\langle y, S_j y\rangle_\mathbf{\beta}- \Vert S_j y\Vert_\mathbf{\beta}^2\big),\] where the last term is non-negative by 76 . Dropping it yields \[\label{eq:us:LI-S:bound} \Vert K_j y \Vert_\mathbf{\beta}^2 = \Vert(L I_s - S_j)\,y\Vert_\mathbf{\beta}^2 \leq L^2\Vert y\Vert_\mathbf{\beta}^2 - (2L-1)\,\Vert S_j y\Vert_\mathbf{\beta}^2.\tag{77}\] A case distinction on \(L\) closes 75 ,

  • \(L \geq \tfrac{1}{2}\): \((2L-1) \geq 0\) and 77 yields \(\Vert K_j y \Vert_\mathbf{\beta}\leq L\, \Vert y\Vert_\mathbf{\beta}\);

  • \(L \leq \tfrac{1}{2}\): using \(\Vert S_j y\Vert_\mathbf{\beta}\leq \Vert y\Vert_\mathbf{\beta}\) in 77 yields \(\Vert K_j y\Vert_\mathbf{\beta}\leq (1-L)\Vert y\Vert_\mathbf{\beta}\).

By the \(c\)-orthonormality of \(\{\tilde{\phi}_j\}\), Parseval’s identity gives \[\sum_j \Vert\xi_j\Vert_\mathbf{\beta}^2 = \big\vert\Theta_u^{n,i-1}\big\vert_{\widetilde{\mathcal{M}},\mathbf{\beta}}^2, \qquad \sum_j \Vert\eta_j\Vert_\mathbf{\beta}^2 = \big\vert\Theta_u^{n,i}\big\vert_{\widetilde{\mathcal{M}},\mathbf{\beta}}^2.\] Applying Cauchy–Schwarz to each summand of 74 as well as to the sum over the modes, we conclude \[\label{eq:us:rhs:bound} \bigg|\sum_j \langle K_j \xi_j,\,\eta_j\rangle_\mathbf{\beta}\bigg| \leq \max(L,\,1-L)\,\big\vert\Theta_u^{n,i-1}\big\vert_{\widetilde{\mathcal{M}},\mathbf{\beta}}\, \big\vert\Theta_u^{n,i}\big\vert_{\widetilde{\mathcal{M}},\mathbf{\beta}}.\tag{78}\] With the coupling strength \(\omega\) from 5 , we have the Rayleigh inequality \[\big\vert v\big\vert_{\widetilde{\mathcal{M}}}^2 = \big\Vert \mathcal{C}^{-1/2}\mathcal{D}\,v\big\Vert^2 \leq \omega\,\Vert v\Vert_a^2.\] Weighted stage-wise by \(\mathbf{\beta}\), this lifts to \[\label{eq:us:rayleigh:stage} \vert\Theta_u^{n,i}\vert_{\widetilde{\mathcal{M}},\mathbf{\beta}}^2 \,\leq\, \omega\,\Vert\Theta_u^{n,i}\Vert_{a,\mathbf{\beta}}^2.\tag{79}\] Now, combining 72 , 74 , 78 , and 79 gives \[\frac{1+L\omega}{\omega}\,\big\vert\Theta_u^{n,i}\big\vert_{\widetilde{\mathcal{M}},\mathbf{\beta}} \,\leq\, \max(L,\,1-L)\,\big\vert\Theta_u^{n,i-1}\big\vert_{\widetilde{\mathcal{M}},\mathbf{\beta}},\] which is ?? with rate \(\rho_{\mathrm{US}}\) as defined in ?? . The conditions for \(\rho_{\mathrm{US}} < 1\) follow from a case analysis of the stabilization parameter \(L\):

  • \(L \geq \tfrac{1}{2}\): \(\rho_{\mathrm{US}} = L\omega/(1+L\omega) < 1\) unconditionally;

  • \(0 \leq L < \tfrac{1}{2}\): \(\rho_{\mathrm{US}} = (1-L)\omega/(1+L\omega) < 1\) iff \(\omega < 1/(1-2L)\).

The minimum is attained at \(L = \tfrac{1}{2}\), yielding \(\rho_{\mathrm{US}} = \omega/(2+\omega) < 1\). ◻

Remark 16 (Case \(L = 0\)). Without stabilization, i.e., for \(L=0\), the contraction constant satisfies \(\rho_{\mathrm{US}}\big|_{L=0} = \omega\). Hence, the iteration is contractive only for \(\omega < 1\).

Remark 17 (Contraction in a genuine norm). The form \(\vert\cdot\vert_{\widetilde{\mathcal{M}}}\) vanishes on \(\ker(\mathcal{D})\), so \(\vert\cdot\vert_{\widetilde{\mathcal{M}},\mathbf{\beta}}\) is, in general, only a seminorm on \(\mathcal{V}^s\). One can show, however, that the iterates \(\Theta_u^{n,i}\) live in a subspace on which \(\vert\cdot\vert_{\widetilde{\mathcal{M}},\mathbf{\beta}}\) is a norm equivalent to \(\Vert\cdot\Vert_{a,\mathbf{\beta}}\), so that ?? is indeed a genuine norm contraction.

4.3 Convergence of iterative RK splittings↩︎

As shown in the previous subsections, the fixed-stress (13) as well as the undrained-split (15) iterations provide a contraction with rate \(\rho < 1\) in their respective norms if the stabilization parameter is chosen appropriately. The convergence argument is identical for both splittings. We write it generically using \(\Vert\cdot\Vert_*\) acting on \(\Theta^{n,i} \vcentcolon= X^{n,i} - X^{n,i-1}\) with

  • \(X = P\) and \(\Vert\cdot\Vert_* = \Vert\cdot\Vert_{c,\mathbf{\beta}}\), \(\rho = \rho_{\mathrm{FS}}\) for fixed-stress,

  • \(X = U\) with \(\Vert\cdot\Vert_* = \vert\cdot\vert_{\widetilde{\mathcal{M}},\mathbf{\beta}}\), \(\rho = \rho_{\mathrm{US}}\) for undrained-split.

Since \(\mathbf{\beta}_\ell > 0\) for all \(\ell\) by 3, the weighted norms \(\Vert\cdot\Vert_{c,\mathbf{\beta}}\) and \(\Vert\cdot\Vert_{a,\mathbf{\beta}}\) are equivalent to the unweighted stage norms \(\Vert\cdot\Vert_c\) and \(\Vert\cdot\Vert_a\), respectively. The same equivalence applies to \(\vert\cdot\vert_{\widetilde{\mathcal{M}},\mathbf{\beta}}\) versus \(\vert\cdot\vert_{\widetilde{\mathcal{M}},s}\).

Theorem 18 (Convergence of iterative RK splittings). Given [ass:iterative:rk,ass:smoothing], let \(u, p\) be the solutions of 6 with sufficient temporal regularity, and consider an \(s\)-stage RK method of stage order \(q\) and classical order \(k\). Let \((u^{n,J_n}, p^{n,J_n})\) denote the iterative solution (fixed-stress or undrained-split) after \(J_n\) iterations satisfying the stopping criterion \[\label{eq:iter:stopping} \big\Vert \Theta^{n,J_n} \big\Vert_* \leq \mathop{\mathrm{TOL}}.\qquad{(14)}\] Then, \[\label{eq:conv:rk} \big\Vert u^{n,J_n} - u(t^n) \big\Vert^2_{\mathcal{V}} + \big\Vert p^{n,J_n} - p(t^n) \big\Vert^2_{\mathcal{H}_{\scalebox{.5}{\mathcal{Q}}}}\, \lesssim\, \frac{\mathop{\mathrm{TOL}}^2}{\tau^3} + \tau^{2\min(k,\,q+1+\sigma)},\qquad{(15)}\] where the hidden constant contains an exponential factor of the form \(\mathrm{e}^{Ct_n}\). Setting \(\mathop{\mathrm{TOL}}= \tau^{\min(k,\,q+1+\sigma)+3/2}\) balances both terms, yielding an overall error of order \(\min(k,\,q+1+\sigma)\).

Proof. Let \((u^n, p^n)\) denote the exact fully coupled RK solution at time \(t^n\) and \(X^n\) the corresponding stage vector (\(X = P\) for fixed-stress, \(X = U\) for undrained-split), so that \[\label{eq:iter:err:decomp} \big\Vert p^{n,J_n} - p(t^n) \big\Vert_{\mathcal{H}_{\scalebox{.5}{\mathcal{Q}}}} \le \big\Vert p^{n,J_n} - p^n \big\Vert_{\mathcal{H}_{\scalebox{.5}{\mathcal{Q}}}} + \big\Vert p^n - p(t^n) \big\Vert_{\mathcal{H}_{\scalebox{.5}{\mathcal{Q}}}}\tag{80}\] and analogously for \(u\). By the contraction results ([thm:rk:contraction:fs,thm:rk:contraction:us]), we have after \(J_n\) iterations \[\Vert X^{n,J_n} - X^{n,J_n-1} \Vert_* \leq \rho^{J_n-1}\, \Vert X^{n,1} - X^{n,0} \Vert_*.\] Summing the geometric series, the stopping criterion ?? yields \(\Vert X^{n,J_n} - X^n \Vert_* \,\lesssim\, \mathop{\mathrm{TOL}}\). Since \(p^{n-1}\) (resp.\(u^{n-1}\)) is fixed across iterates and stiffly-accurate RK methods satisfy \(p^n = \mathbf{\beta}^\mathsf{T}\mathbb{A}^{-1} P^n\) (resp.\(u^n = \mathbf{\beta}^\mathsf{T}\mathbb{A}^{-1} U^n\)) by 27 , the time-step iteration error of the contracted component is bounded by the stage iteration error, \[\label{eq:iter:solution:bound} \big\Vert x^{n,J_n} - x^n \big\Vert \leq \|\mathbf{\beta}^\mathsf{T}\mathbb{A}^{-1}\|\, \Vert X^{n,J_n} - X^n \Vert_* \lesssim \mathop{\mathrm{TOL}},\tag{81}\] where \(x \in \{p,u\}\) matches the contracted component \(X\) and the LHS norm is \(\Vert\cdot\Vert_{\mathcal{H}_{\scalebox{.5}{\mathcal{Q}}}}\) for fixed-stress and \(\Vert\cdot\Vert_\mathcal{V}\) for undrained-split. The error of the other component is then controlled by the elliptic coupling \(\mathcal{A}u = \mathcal{D}^* p + f\) and the stability of \(\mathcal{A}^{-1}\). For the discretization error of the fully coupled implicit RK scheme for 6 , we do not re-derive a Fourier stability estimate. Instead, we invoke the resolvent-smoothing analysis of [17], which, by 2, yields \[\Vert p^n - p(t^n) \Vert_{\mathcal{H}_{\scalebox{.5}{\mathcal{Q}}}} \lesssim \tau^{\min(k,\,q+1+\sigma)}.\] The per-step iteration error 81 of order \(\mathop{\mathrm{TOL}}\) propagates through the parabolic structure, since the stage derivatives involve a factor \(1/\tau\) (from 14 ). A discrete Gronwall argument over \(N = T/\tau\) steps hence gives \[\max_{1 \leq n \leq N} \big\Vert p^{n,J_n} - p^n \big\Vert_{\mathcal{H}_{\scalebox{.5}{\mathcal{Q}}}}^2 \lesssim \frac{\mathop{\mathrm{TOL}}^2}{\tau^3}.\] Combining this with the discretization bound via 80 yields \[\big\Vert p^{n,J_n} - p(t^n) \big\Vert_{\mathcal{H}_{\scalebox{.5}{\mathcal{Q}}}}^2 \lesssim \frac{\mathop{\mathrm{TOL}}^2}{\tau^3} + \tau^{2\min(k,\,q+1+\sigma)}.\] The choice \(\mathop{\mathrm{TOL}}= \tau^{\min(k,\,q+1+\sigma) + 3/2}\) balances both terms, yielding the overall error \(\mathcal{O}(\tau^{2\min(k,\,q+1+\sigma)})\). The bound for \(u\) follows from the elliptic equation \(\mathcal{A}u = \mathcal{D}^* p + f\) and the stability of \(\mathcal{A}^{-1}\). ◻

To summarize, both iterative approaches reach rate \(\rho = \omega/(2 + \omega)\) if the optimal stabilization parameter (\(L = \omega/2\) for fixed-stress, \(L = 1/2\) for undrained-split) is chosen. Hence, both methods converge rapidly for small \(\omega\).

5 Numerical Experiments↩︎

We verify the theoretical convergence rates obtained in [sec:semiexplicit,sec:iterative] using a manufactured solution for linear poroelasticity on the unit square \(\Omega = (0,1)^2\) with \(T = 1\). We prescribe the exact solution as \[\begin{align} \label{eq:mfg:solution} u(t,x,y) &= -\mathrm{e}^{-At} \begin{bmatrix}\sin(\pi x)\sin(\pi y)\\\sin(\pi x)\sin(\pi y)\end{bmatrix}, & p(t,x,y) &= \mathrm{e}^{-At}\,\sin(\pi x)\sin(\pi y). \end{align}\tag{82}\] The forcing terms are chosen accordingly. The decay rate \[\label{eq:decay:rate} A = \frac{2\pi^2\,\kappa/\nu}{\alpha + 1/M}\tag{83}\] depends on all material parameters of the Biot system, which are chosen as

\(\lambda\) \(\mu\) \(\tfrac{\kappa}{\nu}\) \(M\) \(\alpha\)
\(1\) \(0.5\) \(0.1\) \(1\) \(0.1\)

The coupling strength \(\omega\) entering the weak coupling condition ?? is estimated at the continuous level by \(\omega \approx \alpha^2 M / (2\mu) = 10^{-2}\), which is well below \(1/(2^5-1) = 1/31\), the bound for delay order \(k = 5\) (Radau IIA-3, the highest-order method tested). With these parameters, the decay rate 83 evaluates to \(A \approx 1.79\), giving \(\mathrm{e}^{-AT} \approx 0.17\) at \(T = 1\).

Spatial discretization is performed using Taylor–Hood finite elements for \((u, p)\) on a uniform triangular mesh with \(h = 2^{-6}\) (\(64\times 64\) cells). To ensure that temporal errors dominate over spatial discretization effects across the tested time-step range, the polynomial degree is increased with the RK order: \[\begin{align} \text{Radau~IIA-1}&: (P_4,P_3), & \text{Radau~IIA-2}&: (P_6, P_5), & \text{Radau~IIA-3}&: (P_7, P_6). \end{align}\] Each method is paired with \(k = 2s{-}1\) delays matching its classical order. Recall that Radau IIA methods are stiffly accurate, which ensures full classical order convergence for both the algebraic variable \(u\) and the differential variable \(p\); cf. [20]. All errors are reported in the \(L^\infty(0,T)\) norm, i.e., the maximum error over all time steps. The convergence results are summarized in [fig:conv:semiexpl,fig:conv:iter].

Both figures include gray dotted reference lines of slopes \(1\), \(3\), and \(5\), corresponding to the classical orders of the Radau IIA family. For benchmarking, the monolithic implicit RK errors (dashed red) are plotted in both figures. For the implicit scheme, Radau IIA-\(1\) and IIA-\(2\) attain their classical orders \(1\) and \(3\) in both the \(\mathcal{V}\)- and \(\mathcal{H}_{\scalebox{.5}{\mathcal{Q}}}\)-norms. Radau IIA-\(3\) attains a rate of \(\approx 4.2\) for the pressure in agreement with the bound \(r = \min(k,\,q+1+\sigma)\) of 9. For the displacement, this rate can be seen only for coarse time steps before it saturates for \(\tau \leq 2^{-6}\) at the spatial discretization floor. The semi-explicit pressure errors (1) essentially coincide with the implicit ones for all three methods. The semi-explicit displacement errors, however, lie clearly above the implicit ones across all tested \(\tau\), while sharing the same asymptotic rate.

Figure 1: Comparison of implicit (dashed red) and semi-explicit (solid blue) RK discretization.

For the iterative schemes (2), with optimal stabilizations \(L=\omega/2\) (fixed-stress, cf. 13) and \(L=1/2\) (undrained-split, cf. 15), the converged pressure and displacement errors essentially coincide with the monolithic ones, confirming convergence to the fully coupled solution. For Radau IIA-\(3\), the iterative \(u\)-error plateaus instead of tracking the monolithic errors.

Finally, the average iteration counts presented in 1 only show a mild growth with decreasing \(\tau\), using \(\mathop{\mathrm{TOL}}= \tau^{k+3/2}\).

Table 1: Average number of inner iterations per time step for the fixed-stress (\(L=\omega/2\)) and undrained-split (\(L=1/2\)) schemes, with stopping criterion [eq:iter:stopping] and \(\tol = \tau^{k+3/2}\).
Method \(\tau = 2^{-4}\) \(\tau = 2^{-5}\) \(\tau = 2^{-6}\) \(\tau = 2^{-7}\)
fixed-stress, Radau IIA-\(1\) 2.44 2.69 2.92 3.00
fixed-stress, Radau IIA-\(2\) 3.00 3.00 4.00 4.00
fixed-stress, Radau IIA-\(3\) 4.00 4.00 5.00 6.00
undrained-split, Radau IIA-\(1\) 2.38 2.97 3.00 3.00
undrained-split, Radau IIA-\(2\) 3.00 3.00 4.00 4.00
undrained-split, Radau IIA-\(3\) 4.00 4.00 5.00 6.00
Figure 2: Comparicon of iterative RK splittings schemes (fixed-stress in blue, undrained-split in green) with the monolithic implicit RK method (dashed red).

6 Conclusions↩︎

We have presented a convergence analysis for decoupling RK schemes applied to elliptic–parabolic problems. For the semi-explicit schemes based on a delay approximation, we adapted the Fourier stability framework of [17] and established convergence of order \(k\) under weak coupling conditions for each delay order \(k\), matching the bounds obtained for BDF methods in [13]. The consistency of results across BDF and RK time integrators suggests that the coupling bounds are sharp and intrinsic to the delay approximation structure. For the iterative schemes (fixed-stress and undrained-split), we combined contraction analysis with RK consistency estimates, using a spectral decomposition of the Schur complement operator to establish the contraction property. Future directions include the extension to nonlinear problems and the treatment of variable time steps.

Acknowledgments↩︎

This project is funded by the Deutsche Forschungsgemeinschaft (DFG, German Research Foundation) – 467107679. BU further acknowledges funding by the Deutsche Forschungsgemeinschaft (DFG, German Research Foundation) – Project-ID 258734477 – SFB 1173. AM and BU acknowledge support by the Stuttgart Center for Simulation Science (SimTech).

References↩︎

[1]
M. A. Biot. General theory of three-dimensional consolidation. J. Appl. Phys, 12(2):155–164, 1941.
[2]
J. C. Vardakis, D. Chou, B. J. Tully, C. C. Hung, T. H. Lee, P. H. Tsui, and Y. Ventikos. Investigating cerebral oedema using poroelasticity. Med. Eng. Phys., 38(1):48–57, 2016.
[3]
E. Eliseussen, M. E. Rognes, and T. B. Thompson. A posteriori error estimation and adaptivity for multiple-network poroelasticity. ESAIM: Math. Model. Numer. Anal., 57(4):1921–1952, 2023.
[4]
M. D. Zoback. Reservoir geomechanics. Cambridge University Press, Cambridge, 2010.
[5]
R. E. Showalter. Diffusion in poro-elastic media. J. Math. Anal. Appl., 251(1):310–340, 2000.
[6]
A. Ern and S. Meunier. A posteriori error analysis of Euler-Galerkin approximations to coupled elliptic-parabolic problems. ESAIM: Math. Model. Numer. Anal., 43(2):353–375, 2009.
[7]
A. Mikelić and M. F. Wheeler. Convergence of iterative coupling for coupled flow and geomechanics. Comput. Geosci., 17(3):455–461, 2013.
[8]
J. Kim, H.A. Tchelepi, and R. Juanes. Stability and convergence of sequential methods for coupled flow and geomechanics: Drained and undrained splits. Comput. Meth. Appl. Mech. Eng., 200(23):2094–2116, 2011.
[9]
J. Kim, H.A. Tchelepi, and R. Juanes. Stability and convergence of sequential methods for coupled flow and geomechanics: Fixed-stress and fixed-strain splits. Comput. Meth. Appl. Mech. Eng., 200(13):1591–1606, 2011.
[10]
R. Altmann, A. Mujahid, and B. Unger. Higher-order iterative decoupling for poroelasticity. Adv. Comput. Math., 50:11, 2024.
[11]
R. Altmann, R. Maier, and B. Unger. Semi-explicit discretization schemes for weakly-coupled elliptic-parabolic problems. Math. Comp., 90:1089–1118, 2021.
[12]
R. Altmann, R. Maier, and B. Unger. Semi-explicit integration of second order for weakly coupled poroelasticity. BIT Numer. Math., 64:20, 2024.
[13]
R. Altmann, A. Mujahid, and B. Unger. Decoupling multistep schemes for elliptic–parabolic problems. SMAI J. Comput. Math., 2026. accepted for publication.
[14]
R. Altmann and M. Deiml. A novel iterative time integration scheme for linear poroelasticity. Electron. Trans. Numer. Anal., 60:256–275, 2024.
[15]
R. Altmann and M. Deiml. A second-order iterative time integration scheme for linear poroelasticity. SIAM J. Sci. Comput., 47(4):B875–B898, 2025.
[16]
O. Nevanlinna and F. Odeh. Multiplier techniques for linear multistep methods. Numer. Funct. Anal. Optim., 3(4):377–423, 1981.
[17]
C. Lubich and A. Ostermann. approximation of quasi-linear parabolic equations. Math. Comp., 64(210):601–627, 1995.
[18]
E. Zeidler. Nonlinear functional analysis and its applications. 2A: Linear monotone operators. Springer, New York, Berlin, Heidelberg, 1990.
[19]
P. Kunkel and V. Mehrmann. Differential-Algebraic Equations. Analysis and Numerical Solution. European Mathematical Society, Zürich, 2006.
[20]
E. Hairer and G. Wanner. Solving Ordinary Differential Equations II: Stiff and Differential–Algebraic Problems. Springer Berlin Heidelberg, Berlin, Heidelberg, 1996.
[21]
C. Lubich. Convolution quadrature and discretized operational calculus. I. Numer. Math., 52(2):129–145, 1988.
[22]
C. Lubich. Convolution quadrature and discretized operational calculus. II. Numer. Math., 52(4):413–425, 1988.
[23]
C. Lubich and A. Ostermann. methods for parabolic equations and convolution quadrature. Math. Comp., 60(201):105–131, 1993.
[24]
H. Brezis. Functional analysis, Sobolev spaces and partial differential equations. Universitext. Springer, New York, London, 2011.

Subjects

Updated on Academus

2026/06/22 07:33 UTC