# An improved Robin-Robin coupling method for parabolic-parabolic interface problems

Erik Burman, Miguel Angel Fernández, Johnny Guzmán, Sijing Liu

► To cite this version:

Erik Burman, Miguel Angel Fernández, Johnny Guzmán, Sijing Liu. An improved Robin-Robin coupling method for parabolic-parabolic interface problems. 2025. hal-05290627

HAL Id: hal-05290627

https://inria.hal.science/hal-05290627v1

Preprint submitted on 30 Sep 2025

HAL is a multi-disciplinary open access
archive for the deposit and dissemination of sci-
entific research documents, whether they are pub-
lished or not. The documents may come from
teaching and research institutions in France or
abroad, or from public or private research centers.

L'archive ouverte pluridisciplinaire HAL, est
destinée au dépôt et à la diffusion de documents
scientifiques de niveau recherche, publiés ou non,
émanant des établissements d'enseignement et de
recherche français ou étrangers, des laboratoires
publics ou privés.

<figure><img src="image_1.png" /></figure>

Distributed under a Creative Commons Attribution 4.0 International License
# AN IMPROVED ROBIN-ROBIN COUPLING METHOD FOR PARABOLIC-PARABOLIC INTERFACE PROBLEMS

ERIK BURMAN, MIGUEL A. FERNÁNDEZ, JOHNNY GUZMÁN, AND SIJING LIU

**ABSTRACT.** We consider a loosely coupled, non-iterative Robin-Robin coupling method proposed and analyzed in *Numer. Algorithms*, 99:921-948, 2025 for a parabolic-parabolic interface problem. We modify the first step of the scheme so that several error difference quantities remain higher order convergence without requiring additional assumptions. Numerical results are presented to support our findings.

## 1. INTRODUCTION

For the Robin-Robin coupling method considered in [4], we observed in [3] that the following quantities

$$Z^N(\partial_{\Delta t}W, \partial_{\Delta t}U, \partial_{\Delta t}\Lambda) + \sum_{n=1}^{N-1} S^{n+1}(\partial_{\Delta t}W, \partial_{\Delta t}U, \partial_{\Delta t}\Lambda)$$

and

$$Z^N(\partial_{\Delta t}^2 W, \partial_{\Delta t}^2 U, \partial_{\Delta t}^2 \Lambda) + \sum_{n=2}^{N-1} S^{n+1}(\partial_{\Delta t}^2 W, \partial_{\Delta t}^2 U, \partial_{\Delta t}^2 \Lambda)$$

do not converge as $O((\Delta t)^2)$ as expected, where Z and S are defined as follows:

$$Z^{n+1}(\psi, \phi, \theta) := \frac{1}{2} \|\phi^{n+1}\|_{L^2(\Omega_f)}^2 + \frac{1}{2} \|\psi^{n+1}\|_{L^2(\Omega_s)}^2 + \frac{\Delta t \alpha}{2} \|\phi^{n+1}\|_{L^2(\Sigma)}^2 + \frac{\Delta t}{2\alpha} \|\theta^{n+1}\|_{L^2(\Sigma)}^2, \\ S^{n+1}(\psi, \phi, \theta) := \Delta t \left( \nu_f \|\nabla \phi^{n+1}\|_{L^2(\Omega_f)}^2 + \nu_s \|\nabla \psi^{n+1}\|_{L^2(\Omega_s)}^2 \right) + \frac{1}{2} \left( \|\psi^{n+1} - \psi^n\|_{L^2(\Omega_s)}^2 \right. \\ \left. + \|\phi^{n+1} - \phi^n\|_{L^2(\Omega_f)}^2 \right) + \frac{\alpha \Delta t}{2} \|\phi^{n+1} - \phi^n + \frac{1}{\alpha}(\theta^{n+1} - \theta^n)\|_{L^2(\Sigma)}^2.$$

In particular, we observed that the sum terms in $\sum S^{n+1}$ do not converge as $O((\Delta t)^2)$, even though we do observe $O((\Delta t)^2)$ convergence for the terms in $Z^N$ numerically in [4]. However, the $O((\Delta t)^2)$ convergence for the above expressions is a crucial ingredient for the defect correction methods proposed and analyzed in [2]. We examined this issue in [3] and pointed out that, under suitable assumptions on the exact solution of the interface problem, we can obtain $O((\Delta t)^2)$ convergence for both the terms in $Z^N$ and in $\sum S^{n+1}$. Specifically, the additional assumptions state that some of the time derivatives of the exact solutions must vanish during the initial steps.

In this note, we present a modified scheme that differs from that in [1, 4] only at the first time step, but improves the convergence rates of $\sum S^{n+1}$ to $O((\Delta t)^2)$ without requiring any additional assumptions on the exact solutions of the interface problem. However, the first three time steps of the modified scheme are coupled, so it requires solving a large system to obtain the numerical solutions at these time levels. We also briefly discuss the analysis of this modified scheme. The main idea is to modify the first time step so that the initial terms affecting the convergence rates are canceled during the analysis. We follow the notation used in [4] throughout this note.
## 2. THE ORIGINAL ROBIN-ROBIN COUPLING METHOD

We consider the interface problem

$$
\begin{align}
\partial_t u - \nu_f \Delta u &= 0, \quad \text{in } [0, T] \times \Omega_f, \tag{2.1a} \\
u(0, x) &= u_0(x), \quad \text{in } \Omega_f, \nonumber \\
u &= 0, \quad \text{on } [0, T] \times \Gamma_D^f, \nonumber \\
\frac{\partial u}{\partial n} &= 0, \quad \text{on } [0, T] \times \Gamma_{Ne}^f, \nonumber
\end{align}
$$

$$
\begin{align}
\partial_t w - \nu_s \Delta w &= 0, \quad \text{in } [0, T] \times \Omega_s, \tag{2.1b} \\
w(0, x) &= w_0(x), \quad \text{in } \Omega_s, \nonumber \\
w &= 0, \quad \text{on } [0, T] \times \Gamma_D^s, \nonumber \\
\frac{\partial w}{\partial n} &= 0, \quad \text{on } [0, T] \times \Gamma_{Ne}^s, \nonumber
\end{align}
$$

$$
(2.1c) \qquad w - u = 0, \quad \text{on } [0, T] \times \Sigma,
$$

$$
(2.1d) \qquad \nu_s \nabla w \cdot n_s + \nu_f \nabla u \cdot n_f = 0, \quad \text{on } [0, T] \times \Sigma,
$$

where $n_f$ and $n_s$ are the outward facing normal vectors for $\Omega_f$ and $\Omega_s$, respectively. We assume that the initial data is smooth and that $u$ and $w$ is smooth on $\Omega_f$ and $\Omega_s$ respectively.

The Robin–Robin splitting method is given as in [4]. Find $w^{n+1} \in V_s$, $u^{n+1} \in V_f$, and $\lambda^{n+1} \in V_g$ such that, for $n \ge 0$,

$$
(2.2a) \qquad (\partial_{\Delta t} w^{n+1}, z)_s + \nu_s (\nabla w^{n+1}, \nabla z)_s + \alpha \langle w^{n+1} - u^n, z \rangle + \langle \lambda^n, z \rangle = 0, \quad z \in V_s,
$$

$$
(2.2b) \qquad (\partial_{\Delta t} u^{n+1}, v)_f + \nu_f (\nabla u^{n+1}, \nabla v)_f - \langle \lambda^{n+1}, v \rangle = 0, \quad v \in V_f,
$$

$$
(2.2c) \qquad \langle \alpha(u^{n+1} - w^{n+1}) + (\lambda^{n+1} - \lambda^n), \mu \rangle = 0, \quad \mu \in V_g.
$$

## 3. THE IMPROVED METHOD

The improved scheme is defined as follows. For $n = 0$, find $w^1 \in V_s$, $u^1 \in V_f$, and $\lambda^1 \in V_g$ such that

$$
(3.1a) \qquad \begin{aligned}
& (\partial_{\Delta t} w^2, z)_s + \nu_s (\nabla w^1, \nabla z)_s + \alpha \langle w^1 + u^2 - 2u^1, z \rangle + \langle 2\lambda^1 - \lambda^2, z \rangle \\
& \qquad = (\partial_{\Delta t} w^2 - \partial_{\Delta t} w^1, z)_s + \alpha \Delta t \langle G_1^2, z \rangle - \Delta t \langle G_2^2, z \rangle, \quad z \in V_s,
\end{aligned}
$$

$$
(3.1b) \qquad \begin{aligned}
& (\partial_{\Delta t} u^2, v)_f + \nu_f (\nabla u^1, \nabla v)_f + \langle \alpha(\delta u^3 - \delta u^2) + \lambda^3 - 2\lambda^2, v \rangle \\
& \qquad = (\partial_{\Delta t} u^2 - \partial_{\Delta t} u^1, v)_f + \alpha \Delta t \langle G_1^3, v \rangle + \Delta t \langle G_2^3, v \rangle, \quad v \in V_f,
\end{aligned}
$$

$$
(3.1c) \qquad \begin{aligned}
& \langle \alpha(-u^3 + 2u^2 - w^1), \mu \rangle + \langle \lambda^2 - \lambda^1, \mu \rangle \\
& \qquad = -\alpha \Delta t \langle G_1^3, \mu \rangle + \Delta t \langle G_2^2, \mu \rangle, \quad \mu \in V_g.
\end{aligned}
$$

and for $n \ge 1$, find $w^{n+1} \in V_s$, $u^{n+1} \in V_f$, and $\lambda^{n+1} \in V_g$ such that,

$$
(3.2a) \qquad (\partial_{\Delta t} w^{n+1}, z)_s + \nu_s (\nabla w^{n+1}, \nabla z)_s + \alpha \langle w^{n+1} - u^n, z \rangle + \langle \lambda^n, z \rangle = 0, \quad z \in V_s,
$$

$$
(3.2b) \qquad (\partial_{\Delta t} u^{n+1}, v)_f + \nu_f (\nabla u^{n+1}, \nabla v)_f - \langle \lambda^{n+1}, v \rangle = 0, \quad v \in V_f,
$$

$$
(3.2c) \qquad \langle \alpha(u^{n+1} - w^{n+1}) + (\lambda^{n+1} - \lambda^n), \mu \rangle = 0, \quad \mu \in V_g.
$$

Notice that this modified scheme differs from the original scheme only at the very first step when $n = 0$. However, the first step (3.1) requires information about the numerical solutions at $t_2$ and $t_3$, as well as the exact solutions at $t_1$, $t_2$, and $t_3$. This implies that the first three steps of the method are coupled, hence, they should be solved as a single large system. We also believe that the exact solutions
on the right-hand side of the method can be replaced by high-order numerical approximations. This could potentially be achieved using a higher-order monolithic method for the first few steps.

We now briefly discuss the analysis of the modified scheme (3.1)-(3.2). For $n = 0$, we write the exact solution as follows:

$$ 
\begin{align*}
(\partial_{\Delta t} w^2, z)_s + \nu_s (\nabla w^1, \nabla z)_s + \alpha \langle w^1 + u^2 - 2u^1, z \rangle + \langle 2|^1 - |^2, z \rangle &= -(h_3^2, z)_s + \alpha \langle g_1^2, z \rangle - \langle g_2^2, z \rangle, && z \in V_s, \\
(\partial_{\Delta t} u^2, v)_f + \nu_f (\nabla u^1, \nabla v)_f + \alpha (\delta u^3 - \delta u^2) + |^3 - 2|^2, v \rangle &= \alpha \Delta t \langle G_1^3, v \rangle + \Delta t \langle G_2^3, v \rangle - (h_4^2, v)_f, && v \in V_f, \\
\alpha \langle -u^3 + 2u^2 - w^1, \mu \rangle + \langle |^2 - |^1, \mu \rangle &= - \alpha \Delta t \langle G_1^3, \mu \rangle + \langle g_2^2, \mu \rangle, && \mu \in V_g.
\end{align*}
$$

For $n \ge 1$, we write the exact solution as follows:

$$ (3.3a) \quad (\partial_{\Delta t} w^{n+1}, z)_s + \nu_s (\nabla w^{n+1}, \nabla z)_s + \alpha \langle w^{n+1} - u^n, z \rangle + \langle |^n, z \rangle = L_1^{n+1}(z) - L_3^{n+1}(z), \quad z \in V_s, $$

$$ (3.3b) \quad (\partial_{\Delta t} u^{n+1}, v)_f + \nu_f (\nabla u^{n+1}, \nabla v)_f - \langle |^{n+1}, v \rangle = L_2^{n+1}(v), \quad v \in V_f, $$

$$ (3.3c) \quad \langle \alpha(u^{n+1} - w^{n+1}) + (|^{n+1} - |^n), \mu \rangle = L_3^{n+1}(\mu), \quad \mu \in V_g, $$

where

$$ L_1^{n+1}(z) := \alpha \langle g_1^{n+1}, z \rangle - (h_1^{n+1}, z)_s, $$

$$ L_2^{n+1}(v) := -(h_2^{n+1}, v)_f, $$

$$ L_3^{n+1}(\mu) := \langle g_2^{n+1}, \mu \rangle, $$

and

$$ h_1^{n+1} := \partial_t w^{n+1} - \partial_{\Delta t} w^{n+1}, \quad g_1^{n+1} := u^{n+1} - u^n, $$

$$ h_2^{n+1} := \partial_t u^{n+1} - \partial_{\Delta t} u^{n+1}, \quad g_2^{n+1} := |^{n+1} - |^n, $$

$$ h_3^{n+1} := \partial_t w^n - \partial_{\Delta t} w^{n+1}, \quad h_4^{n+1} := \partial_t u^n - \partial_{\Delta t} u^{n+1}, $$

$$ G_1^{n+1} := \partial_{\Delta t} g_1^{n+1}, \quad G_2^{n+1} := \partial_{\Delta t} g_2^{n+1}. $$

### 3.1. Error Equations. We use the notations

$$ 
\begin{align*}
U^n &= u^n - u^n, & W^n &= w^n - w^n, & \Lambda^n &= |^n - \lambda^n \\
U^{n+1} &= \partial_{\Delta t} U^{n+1}, & W^{n+1} &= \partial_{\Delta t} W^{n+1}, & L^{n+1} &= \partial_{\Delta t} \Lambda^{n+1} \\
U^{n+1} &= \partial_{\Delta t} U^{n+1}, & W^{n+1} &= \partial_{\Delta t} W^{n+1}, & \mathcal{L}^{n+1} &= \partial_{\Delta t} \mathcal{L}^{n+1},
\end{align*}
$$

and

$$ H_j^{n+1} = \partial_{\Delta t} h_j^{n+1}, $$

$$ G_j^{n+1} = \partial_{\Delta t} G_j^{n+1}, $$

$$ \mathcal{H}_j^{n+1} = \partial_{\Delta t} H_j^{n+1}, $$

for $j = 1, 2$.

By subtracting the exact solutions and the numerical solutions, we have the following error equations. For $n = 0$, we have

$$ 
\begin{align}
(3.4a) \quad & (\partial_{\Delta t} W^2, z)_s + \nu_s (\nabla W^1, \nabla z)_s + \alpha \langle W^1 + U^2 - 2U^1, z \rangle + \langle 2\Lambda^1 - \Lambda^2, z \rangle \nonumber \\
& = -(h_1^1, z)_s + \alpha \langle g_1^1, z \rangle - \langle g_2^1, z \rangle, \tag{1}
\end{align}
$$

$$ (3.4b) \quad (\partial_{\Delta t} U^2, z)_f + \nu_f (\nabla U^1, \nabla v)_f + \langle \alpha(\delta U^3 - \delta U^2) + \Lambda^3 - 2\Lambda^2, v \rangle = -(h_2^1, v)_f, $$

$$ (3.4c) \quad \alpha \langle -U^3 + 2U^2 - W^1, \mu \rangle + \langle \Lambda^2 - \Lambda^1, \mu \rangle = \langle g_2^1, \mu \rangle $$
and for $n \ge 1$:

$$
(3.5a) \quad (\partial_{\Delta t} W^{n+1}, z)_s + \nu_s (\nabla W^{n+1}, \nabla z)_s + \alpha \langle (W^{n+1} - U^n), z \rangle + \langle \Lambda^n, z \rangle = L_1^{n+1}(z) - L_3^{n+1}(z),
$$

$$
(3.5b) \quad (\partial_{\Delta t} U^{n+1}, v)_f + \nu_f (\nabla U^{n+1}, \nabla v)_f - \langle \Lambda^{n+1}, v \rangle = L_2^{n+1}(v),
$$

$$
(3.5c) \quad \alpha \langle U^{n+1} - W^{n+1}, \mu \rangle + \langle \Lambda^{n+1} - \Lambda^n, \mu \rangle = L_3^{n+1}(\mu).
$$

3.1.1. *First-order time difference.* Furthermore, using (3.4) and (3.5), we have, for $n = 1$

$$
(3.6a) \quad \nu_s (\nabla W^2, \nabla z)_s + \alpha \langle W^2 - U^2, z \rangle + \langle L^2, z \rangle = R_1^2(z) - R_3^2(z),
$$

$$
(3.6b) \quad \nu_f (\nabla U^2, \nabla v)_f - \alpha \langle U^3 - U^2, v \rangle - \langle L^3, v \rangle = R_2^2(v),
$$

$$
(3.6c) \qquad \alpha \langle U^3 - W^2, \mu \rangle = R_3^2(\mu),
$$

and $n \ge 2$:

$$
(3.7a) \quad (\partial_{\Delta t} \mathcal{W}^{n+1}, z)_s + \nu_s (\nabla \mathcal{W}^{n+1}, \nabla z)_s + \alpha \langle (\mathcal{W}^{n+1} - \mathcal{U}^n), z \rangle + \langle \mathcal{L}^n, z \rangle = R_1^{n+1}(z) - R_3^{n+1}(z),
$$

$$
(3.7b) \quad (\partial_{\Delta t} \mathcal{U}^{n+1}, v)_f + \nu_f (\nabla \mathcal{U}^{n+1}, \nabla v)_f - \langle \mathcal{L}^{n+1}, v \rangle = R_2^{n+1}(v),
$$

$$
(3.7c) \qquad \alpha \langle \mathcal{U}^{n+1} - \mathcal{W}^{n+1}, \mu \rangle + \langle \mathcal{L}^{n+1} - \mathcal{L}^n, \mu \rangle = R_3^{n+1}(\mu),
$$

where

$$
R_1(z) := \alpha \langle G_1^{n+1}, z \rangle - (H_1^{n+1}, z)_s,
$$

$$
R_3(v) := -(H_2^{n+1}, v)_f,
$$

$$
R_4(\mu) := \langle G_2^{n+1}, \mu \rangle.
$$

3.1.2. *Second-order time difference.* Subtract (3.6) and (3.7) and divide them by $\Delta t$, we have, for $n=2$,

$$
(3.8a) \quad \frac{1}{\Delta t} (\mathfrak{W}^3, z)_s + \nu_s (\nabla \mathfrak{W}^3, \nabla z)_s + \alpha \langle \mathfrak{W}^3, z \rangle = \mathcal{R}_1^3(z) - \mathcal{R}_3^3(z),
$$

$$
(3.8b) \quad \frac{1}{\Delta t} (\mathfrak{U}^3, v)_f + \nu_f (\nabla \mathfrak{U}^3, \nabla v)_f + \alpha \langle \mathfrak{U}^3, v \rangle = \mathcal{R}_2^3(v),
$$

$$
(3.8c) \qquad \langle \mathcal{L}^3 - \alpha \mathfrak{W}^3, \mu \rangle = \mathcal{R}_3^3(\mu),
$$

and $n \ge 3$,

$$
(3.9a) \quad (\partial_{\Delta t} \mathfrak{W}^{n+1}, z)_s + \nu_s (\nabla \mathfrak{W}^{n+1}, \nabla z)_s + \alpha \langle (\mathfrak{W}^{n+1} - \mathcal{U}^n), z \rangle + \langle \mathcal{L}^n, z \rangle = \mathcal{R}_1^{n+1}(z) - \mathcal{R}_3^{n+1}(z),
$$

$$
(3.9b) \quad (\partial_{\Delta t} \mathcal{U}^{n+1}, v)_f + \nu_f (\nabla \mathcal{U}^{n+1}, \nabla v)_f - \langle \mathcal{L}^{n+1}, v \rangle = \mathcal{R}_2^{n+1}(v),
$$

$$
(3.9c) \qquad \alpha \langle \mathcal{U}^{n+1} - \mathfrak{W}^{n+1}, \mu \rangle + \langle \mathcal{L}^{n+1} - \mathcal{L}^n, \mu \rangle = \mathcal{R}_3^{n+1}(\mu),
$$

where

$$
\mathcal{R}_1^{n+1}(z) := \alpha \langle S_1^{n+1}, z \rangle - (\mathcal{H}_1^{n+1}, z)_s,
$$

$$
\mathcal{R}_2^{n+1}(v) := -(\mathcal{H}_2^{n+1}, v)_f,
$$

$$
\mathcal{R}_3^{n+1}(\mu) := \langle S_2^{n+1}, \mu \rangle.
$$

3.2. **Sketch of the Analysis.** Take $z = \mathfrak{W}^{n+1}$ in (3.9a) and $v = \mathcal{U}^{n+1}$ in (3.9b) and use (3.9c), we have (cf. [4, Lemma 3.3]),

**Lemma 3.1.** It holds, for $n \ge 3$,

$$
Z^{n+1}(\mathfrak{U}, \mathfrak{W}, \mathcal{L}) + S^{n+1}(\mathfrak{U}, \mathfrak{W}, \mathcal{L}) = Z^n(\mathfrak{U}, \mathfrak{W}, \mathcal{L}) + \Delta t \mathfrak{F}^{n+1} + \frac{\Delta t}{\alpha} \langle S_2^{n+1}, \mathcal{L}^{n+1} \rangle,
$$

where

$$
\mathfrak{F}^{n+1} := -(\mathcal{H}_1^{n+1}, \mathfrak{W}^{n+1})_s - (\mathcal{H}_2^{n+1}, \mathcal{U}^{n+1})_f + \alpha \langle S_1^{n+1}, \mathfrak{W}^{n+1} \rangle + \langle \mathcal{U}^{n+1} - \mathcal{U}^n, S_2^{n+1} \rangle.
$$

Hence, by taking the sum from $n=3$ to $n=N-1$,
**Lemma 3.2.** We have,

$$ (3.10) \quad Z^N(\mathfrak{u}, \mathfrak{w}, \mathcal{L}) + \sum_{n=3}^{N-1} S^{n+1}(\mathfrak{u}, \mathfrak{w}, \mathcal{L}) = Z^3(\mathfrak{u}, \mathfrak{w}, \mathcal{L}) + \Delta t \sum_{n=3}^{N-1} \mathfrak{F}^{n+1} + \frac{\Delta t}{\alpha} \sum_{n=3}^{N-1} \langle \mathfrak{S}_2^{n+1}, \mathfrak{L}^{n+1} \rangle. $$

Using Lemma 3.2, we obtain

**Theorem 3.3.** We have,

$$ (3.11) \quad Z^N(\mathfrak{u}, \mathfrak{w}, \mathcal{L}) + \sum_{n=3}^{N-1} S^{n+1}(\mathfrak{u}, \mathfrak{w}, \mathcal{L}) + \frac{1}{2} \|\mathfrak{w}^3\|_{L^2(\Omega_s)}^2 + \frac{1}{2} \|\mathfrak{u}^3\|_{L^2(\Omega_f)}^2 \\ + \Delta t \nu_s \|\nabla \mathfrak{w}^3\|_{L^2(\Omega_s)}^2 + \Delta t \nu_f \|\nabla \mathfrak{u}^3\|_{L^2(\Omega_f)}^2 + \frac{\Delta t}{2\alpha} \|\mathfrak{L}^3\|_{L^2(\Sigma)}^2 + \frac{\alpha \Delta t}{2} \|\mathfrak{u}^3\|_{L^2(\Sigma)}^2 \le C(\Delta t)^2. $$

*Sketch of the proof.* Take $z = \mathfrak{w}^3$ in (3.8a) and $v = \mathfrak{u}^3$ in (3.8b) and use (3.8c), we have

$$ (3.12) \quad \begin{aligned} & \|\mathfrak{w}^3\|_{L^2(\Omega_s)}^2 + \|\mathfrak{u}^3\|_{L^2(\Omega_f)}^2 + \Delta t \nu_s \|\nabla \mathfrak{w}^3\|_{L^2(\Omega_s)}^2 + \Delta t \nu_f \|\nabla \mathfrak{u}^3\|_{L^2(\Omega_f)}^2 \\ & \quad + \frac{\Delta t}{\alpha} \|\mathfrak{L}^3\|_{L^2(\Sigma)}^2 + \alpha \Delta t \|\mathfrak{u}^3\|_{L^2(\Sigma)}^2 \\ & = \Delta t \left[ \mathfrak{R}_1^3(\mathfrak{w}^3) - \mathfrak{R}_3^3(\mathfrak{w}^3) + \mathfrak{R}_2^3(\mathfrak{u}^3) + \langle \mathfrak{S}_2^3, \mathfrak{w}^3 \rangle + \frac{1}{\alpha} \langle \mathfrak{S}_2^3, \mathfrak{L}^3 \rangle \right] := \Delta t \mathfrak{J}^3 \end{aligned} $$

Adding (3.12) and (3.10), we obtain

$$ (3.13) \quad \begin{aligned} & Z^N(\mathfrak{u}, \mathfrak{w}, \mathcal{L}) + \sum_{n=3}^{N-1} S^{n+1}(\mathfrak{u}, \mathfrak{w}, \mathcal{L}) + \frac{1}{2} \|\mathfrak{w}^3\|_{L^2(\Omega_s)}^2 + \frac{1}{2} \|\mathfrak{u}^3\|_{L^2(\Omega_f)}^2 \\ & \quad + \Delta t \nu_s \|\nabla \mathfrak{w}^3\|_{L^2(\Omega_s)}^2 + \Delta t \nu_f \|\nabla \mathfrak{u}^3\|_{L^2(\Omega_f)}^2 + \frac{\Delta t}{2\alpha} \|\mathfrak{L}^3\|_{L^2(\Sigma)}^2 + \frac{\alpha \Delta t}{2} \|\mathfrak{u}^3\|_{L^2(\Sigma)}^2 \\ & = \Delta t \mathfrak{J}^3 + \Delta t \sum_{n=3}^{N-1} \mathfrak{F}^{n+1} + \frac{\Delta t}{\alpha} \sum_{n=3}^{N-1} \langle \mathfrak{S}_2^{n+1}, \mathfrak{L}^{n+1} \rangle. \end{aligned} $$

The terms involving $\mathfrak{J}^3$, $\mathfrak{F}^{n+1}$ and the last term in (3.13) can be estimated similarly as in [4] and then we obtain the results immediately. $\square$

*Remark 3.4.* Theorem 3.3 indicates that the term $(\Delta t \sum_{n=2}^{N-1} \nu_f \|\nabla(U^{n+1} - 2U^n + U^{n-1})\|_{L^2(\Omega_f)}^2)^{\frac{1}{2}}$ is of order $O((\Delta t)^3)$, which is a crucial ingredient for a defect correction scheme in [2]. The idea of the modified scheme is to construct the first time step (3.1) such that the first step of the error equations of the second-order time difference (3.8) have this particular form, so that when performing the analysis, the $Z^3$ term on the right-hand side of (3.10) is canceled. We also observe $O((\Delta t)^2)$ convergence rates for the first difference of the errors and we also believe this idea can be used to prove the $H^2$ estimate in [4]. We omit the details here.

## 4. NUMERICAL RESULTS

In this section, we present three numerical examples to support our findings. We begin by re-computing Example 5.1 from [4] to demonstrate that the improved scheme indeed achieves better convergence rates for the terms in $\Sigma_{n=2}S^{n+1}$. We then include two additional examples with different exact solutions. We define the following notations:

$$ \begin{aligned} e_u &= \|U^N\|_{L^2(\Omega_f)}, & e_{du} &= \|U^N - U^{N-1}\|_{L^2(\Omega_f)}, \\ e_{dw} &= \|W^N - W^{N-1}\|_{L^2(\Omega_s)}, & e_{gdu} &= \|\nabla(U^N - U^{N-1})\|_{L^2(\Omega_f)}. \end{aligned} $$
and

$$e_{gdus} = (\Delta t \sum_{n=1}^{N-1} \|\nabla(U^{n+1} - U^n)\|_{L^2(\Omega_f)}^2)^{\frac{1}{2}}, \quad e_{gdws} = (\Delta t \sum_{n=1}^{N-1} \|\nabla(W^{n+1} - W^n)\|_{L^2(\Omega_s)}^2)^{\frac{1}{2}},$$

$$e_{dls} = (\Delta t \sum_{n=1}^{N-1} \|\Lambda^{n+1} - \Lambda^n\|_{L^2(\Sigma)}^2)^{\frac{1}{2}}, \quad e_{gdu2s} = (\Delta t \sum_{n=2}^{N-1} \|\nabla(U^{n+1} - 2U^n + U^{n-1})\|_{L^2(\Omega_f)}^2)^{\frac{1}{2}},$$

$$e_{ggdus} = (\Delta t \sum_{n=1}^{N-1} \|D^2(U^{n+1} - U^n)\|_{L^2(\Omega_f)}^2)^{\frac{1}{2}}.$$

*Example 4.1. We take*

$$u = w = e^{-2\pi^2 t} \cos(\pi x_1) \sin(\pi x_2),$$

$h = \Delta t$ and $\alpha = 4$ in this example as in [4, Example 5.1]. We set $\Omega_f = [0, 1] \times [0, 0.75]$ and $\Omega_s = [0, 1] \times [0.75, 1]$. Here at level $k$ we set $\Delta t = (\frac{1}{2})^{k+1}$.

We report the convergence rates in Tables 1-2. We observe the same convergence rates for quantities at final time $N$ as in [4, Example 5.1]. In contrast to [3], we observe $O((\Delta t)^2)$ for the sum of the first difference of the errors and $O((\Delta t)^3)$ for the sum of the second difference of the errors. These convergence rates coincide with our assumptions in [2].

<table><caption>TABLE 1. Convergence results for Example 4.3 at final time with linear finite element ($h = \Delta t$).</caption><thead><tr><th>k</th><th>e<sub>u</sub></th><th>Order</th><th>e<sub>du</sub></th><th>Order</th><th>e<sub>dw</sub></th><th>Order</th><th>e<sub>gdu</sub></th><th>Order</th></tr></thead><tbody><tr><td>1</td><td>4.65e-01</td><td>-</td><td>2.05e-02</td><td>-</td><td>1.76e-02</td><td>-</td><td>8.18e-02</td><td>-</td></tr><tr><td>2</td><td>4.34e-01</td><td>0.10</td><td>6.09e-02</td><td>-1.57</td><td>3.67e-02</td><td>-1.06</td><td>2.33e-01</td><td>-1.51</td></tr><tr><td>3</td><td>4.25e-02</td><td>3.35</td><td>3.77e-02</td><td>0.69</td><td>2.35e-02</td><td>0.64</td><td>1.45e-01</td><td>0.69</td></tr><tr><td>4</td><td>1.50e-02</td><td>1.50</td><td>7.59e-03</td><td>2.31</td><td>4.02e-03</td><td>2.55</td><td>2.95e-02</td><td>2.30</td></tr><tr><td>5</td><td>5.77e-03</td><td>1.38</td><td>1.52e-03</td><td>2.32</td><td>7.17e-04</td><td>2.49</td><td>5.97e-03</td><td>2.30</td></tr><tr><td>6</td><td>2.46e-03</td><td>1.23</td><td>3.24e-04</td><td>2.23</td><td>1.43e-04</td><td>2.33</td><td>1.28e-03</td><td>2.22</td></tr><tr><td>7</td><td>1.12e-03</td><td>1.14</td><td>7.35e-05</td><td>2.14</td><td>3.12e-05</td><td>2.19</td><td>2.93e-04</td><td>2.14</td></tr><tr><td>8</td><td>5.26e-04</td><td>1.09</td><td>1.72e-05</td><td>2.09</td><td>7.20e-06</td><td>2.12</td><td>6.87e-05</td><td>2.09</td></tr></tbody></table>

<table><caption>TABLE 2. Convergence results for Example 4.3 of the sum with linear finite element ($h = \Delta t$).</caption><thead><tr><th>k</th><th>e<sub>gdus</sub></th><th>Order</th><th>e<sub>gdws</sub></th><th>Order</th><th>e<sub>gdu2s</sub></th><th>Order</th></tr></thead><tbody><tr><td>3</td><td>7.37e-02</td><td>-</td><td>9.42e-02</td><td>-</td><td>4.58e-02</td><td>-</td></tr><tr><td>4</td><td>4.07e-02</td><td>0.86</td><td>5.20e-02</td><td>0.86</td><td>1.18e-02</td><td>1.96</td></tr><tr><td>5</td><td>1.53e-02</td><td>1.41</td><td>1.92e-02</td><td>1.44</td><td>2.09e-03</td><td>2.50</td></tr><tr><td>6</td><td>4.64e-03</td><td>1.72</td><td>5.77e-03</td><td>1.73</td><td>3.08e-04</td><td>2.76</td></tr><tr><td>7</td><td>1.25e-03</td><td>1.89</td><td>1.56e-03</td><td>1.89</td><td>4.09e-05</td><td>2.91</td></tr><tr><td>8</td><td>3.08e-04</td><td>2.02</td><td>3.96e-04</td><td>1.98</td><td>5.03e-06</td><td>3.02</td></tr></tbody></table>

*Example 4.2. We take*

$$u = w = (t^3 + 1) \cos(\pi x_1) \sin(\pi x_2),$$
in this example. Other parameters are identical as those in Example 4.1.

We report the convergence rates in Tables 3–4. In this example, we use quadratic finite elements to eliminate any impact of the spatial discretization on the convergence. Similar convergence rates are observed as in Example 4.1.

<table><caption>TABLE 3. Convergence results for Example 4.2 at final time with quadratic finite element ($h = \Delta t$).</caption><thead><tr><th>k</th><th>e<sub>u</sub></th><th>Order</th><th>e<sub>du</sub></th><th>Order</th><th>e<sub>dw</sub></th><th>Order</th><th>e<sub>gdu</sub></th><th>Order</th></tr></thead><tbody><tr><td>1</td><td>2.10e-01</td><td>-</td><td>1.45e-01</td><td>-</td><td>2.53e-02</td><td>-</td><td>5.87e-01</td><td>-</td></tr><tr><td>2</td><td>2.77e-02</td><td>2.93</td><td>1.88e-02</td><td>2.95</td><td>3.45e-03</td><td>2.87</td><td>7.59e-02</td><td>2.95</td></tr><tr><td>3</td><td>8.34e-04</td><td>5.05</td><td>2.87e-04</td><td>6.03</td><td>3.73e-04</td><td>3.21</td><td>1.15e-03</td><td>6.04</td></tr><tr><td>4</td><td>4.46e-04</td><td>0.90</td><td>7.33e-05</td><td>1.97</td><td>1.14e-04</td><td>1.71</td><td>2.93e-04</td><td>1.98</td></tr><tr><td>5</td><td>2.31e-04</td><td>0.95</td><td>1.83e-05</td><td>2.00</td><td>3.09e-05</td><td>1.88</td><td>7.32e-05</td><td>2.00</td></tr><tr><td>6</td><td>1.18e-04</td><td>0.98</td><td>4.55e-06</td><td>2.01</td><td>8.06e-06</td><td>1.94</td><td>1.83e-05</td><td>2.00</td></tr><tr><td>7</td><td>5.93e-05</td><td>0.99</td><td>1.14e-06</td><td>2.00</td><td>2.05e-06</td><td>1.97</td><td>4.56e-06</td><td>2.00</td></tr><tr><td>8</td><td>2.98e-05</td><td>0.99</td><td>2.83e-07</td><td>2.00</td><td>5.19e-07</td><td>1.99</td><td>1.14e-06</td><td>2.00</td></tr></tbody></table>

<table><caption>TABLE 4. Convergence results for Example 4.2 of the sum with quadratic finite element ($h = \Delta t$).</caption><thead><tr><th>k</th><th>e<sub>gdus</sub></th><th>Order</th><th>e<sub>gdws</sub></th><th>Order</th><th>e<sub>gdu2s</sub></th><th>Order</th></tr></thead><tbody><tr><td>3</td><td>4.07e-04</td><td>-</td><td>8.89e-04</td><td>-</td><td>4.81e-06</td><td>-</td></tr><tr><td>4</td><td>1.26e-04</td><td>1.69</td><td>2.58e-04</td><td>1.78</td><td>1.65e-06</td><td>1.55</td></tr><tr><td>5</td><td>3.43e-05</td><td>1.88</td><td>6.86e-05</td><td>1.91</td><td>3.77e-07</td><td>2.13</td></tr><tr><td>6</td><td>8.90e-06</td><td>1.95</td><td>1.76e-05</td><td>1.96</td><td>6.70e-08</td><td>2.49</td></tr><tr><td>7</td><td>2.26e-06</td><td>1.98</td><td>4.47e-06</td><td>1.98</td><td>1.06e-08</td><td>2.66</td></tr><tr><td>8</td><td>5.70e-07</td><td>1.99</td><td>1.13e-06</td><td>1.99</td><td>1.55e-09</td><td>2.78</td></tr></tbody></table>

Example 4.3. We take

$$u = w = e^t \cos(\pi x_1) \sin(\pi x_2),$$

in this example. Other parameters are identical as those in Example 4.1.

We report the convergence rates in Tables 5–6. In this example, we use quadratic finite elements to eliminate any impact of the spatial discretization on the convergence. Similar convergence rates are observed as in Example 4.1.
<table><caption>TABLE 5. Convergence results for Example 4.3 at final time with quadratic finite element ($h = \Delta t$).</caption><thead><tr><th>k</th><th>e<sub>u</sub></th><th>Order</th><th>e<sub>du</sub></th><th>Order</th><th>e<sub>dw</sub></th><th>Order</th><th>e<sub>gdu</sub></th><th>Order</th></tr></thead><tbody><tr><td>1</td><td>5.34e-01</td><td>-</td><td>2.23e-01</td><td>-</td><td>5.90e-02</td><td>-</td><td>9.03e-01</td><td>-</td></tr><tr><td>2</td><td>2.18e-01</td><td>1.29</td><td>8.17e-02</td><td>1.45</td><td>2.36e-02</td><td>1.32</td><td>3.30e-01</td><td>1.45</td></tr><tr><td>3</td><td>6.58e-04</td><td>8.37</td><td>3.83e-05</td><td>11.06</td><td>5.10e-04</td><td>5.53</td><td>2.05e-04</td><td>10.65</td></tr><tr><td>4</td><td>3.29e-04</td><td>1.00</td><td>1.00e-05</td><td>1.93</td><td>1.31e-04</td><td>1.96</td><td>5.08e-05</td><td>2.01</td></tr><tr><td>5</td><td>1.65e-04</td><td>1.00</td><td>2.55e-06</td><td>1.98</td><td>3.33e-05</td><td>1.98</td><td>1.29e-05</td><td>1.98</td></tr><tr><td>6</td><td>8.24e-05</td><td>1.00</td><td>6.40e-07</td><td>1.99</td><td>8.40e-06</td><td>1.99</td><td>3.24e-06</td><td>1.99</td></tr><tr><td>7</td><td>4.12e-05</td><td>1.00</td><td>1.60e-07</td><td>2.00</td><td>2.11e-06</td><td>1.99</td><td>8.13e-07</td><td>2.00</td></tr><tr><td>8</td><td>2.06e-05</td><td>1.00</td><td>4.01e-08</td><td>2.00</td><td>5.28e-07</td><td>2.00</td><td>2.04e-07</td><td>2.00</td></tr></tbody></table>

<table><caption>TABLE 6. Convergence results for Example 4.3 of the sum with quadratic finite element ($h = \Delta t$).</caption><thead><tr><th>k</th><th>e<sub>gdus</sub></th><th>Order</th><th>e<sub>gdws</sub></th><th>Order</th><th>e<sub>gdu2s</sub></th><th>Order</th></tr></thead><tbody><tr><td>3</td><td>7.13e-05</td><td>-</td><td>1.39e-03</td><td>-</td><td>3.14e-06</td><td>-</td></tr><tr><td>4</td><td>2.05e-05</td><td>1.80</td><td>4.19e-04</td><td>1.73</td><td>6.45e-07</td><td>2.28</td></tr><tr><td>5</td><td>5.45e-06</td><td>1.91</td><td>1.12e-04</td><td>1.90</td><td>1.04e-07</td><td>2.63</td></tr><tr><td>6</td><td>1.40e-06</td><td>1.96</td><td>2.90e-05</td><td>1.96</td><td>1.56e-08</td><td>2.74</td></tr><tr><td>7</td><td>3.54e-07</td><td>1.98</td><td>7.35e-06</td><td>1.98</td><td>2.25e-09</td><td>2.79</td></tr><tr><td>8</td><td>8.91e-08</td><td>1.99</td><td>1.85e-06</td><td>1.99</td><td>3.13e-10</td><td>2.85</td></tr></tbody></table>

REFERENCES

[1] E. Burman, R. Durst, M. Fernández, and J. Guzmán. Loosely coupled, non-iterative time-splitting scheme based on Robin–Robin coupling: Unified analysis for parabolic/parabolic and parabolic/hyperbolic problems. *Journal of Numerical Mathematics*, 31(1):59–77, 2023.

[2] E. Burman, R. Durst, M. A. Fernández, J. Guzmán, and S. Liu. A second-order correction method for loosely coupled discretizations applied to parabolic-parabolic interface problems. *IMA Journal of Numerical Analysis*, page drae075, 2024.

[3] E. Burman, R. Durst, M. A. Fernández, J. Guzmán, and S. Liu. Erratum: Estimates of discrete time derivatives for the parabolic-parabolic Robin-Robin coupling method. 2025.

[4] E. Burman, R. Durst, M. A. Fernández, J. Guzmán, and S. Liu. Estimates of discrete time derivatives for the parabolic-parabolic Robin-Robin coupling method. *Numerical Algorithms*, 99(2):921–948, 2025.

<footnote n="1">¹DEPARTMENT OF MATHEMATICS, UNIVERSITY COLLEGE LONDON, LONDON, UK–WC1E 6BT, UNITED KINGDOM
*Email address: e.burman@ucl.ac.uk*</footnote>

<footnote n="2">²SORBONNE UNIVERSITÉ & CNRS, UMR 7598 LJLL, 75005 PARIS FRANCE – INRIA PARIS, 75012 PARIS, FRANCE
*Email address: miguel.fernandez@inria.fr*</footnote>

<footnote n="3">³DIVISION OF APPLIED MATHEMATICS, BROWN UNIVERSITY, 182 GEORGE STREET, PROVIDENCE, RI 02912, USA
*Email address: johnny_guzman@brown.edu*</footnote>

<footnote n="4">⁴DEPARTMENT OF MATHEMATICAL SCIENCES, WORCESTER POLYTECHNIC INSTITUTE, 100 INSTITUTE ROAD, WORCESTER, MA 01609, USA
*Email address: sliu13@wpi.edu*</footnote>