Estimates of Discrete Time Derivatives for the Parabolic-Parabolic Robin-Robin Coupling Method


We consider a loosely coupled, non-iterative Robin-Robin coupling method proposed and analyzed in [J. Numer. Math., 31(1):59–77, 2023] for a parabolic-parabolic interface problem and prove estimates for the discrete time derivatives of the scalar field in different norms. When the interface is flat and perpendicular to two of the edges of the domain we prove error estimates in the \(H^2\)-norm. Such estimates are key ingredients to analyze a defect correction method for the parabolic-parabolic interface problem. Numerical results are shown to support our findings.

1 Introduction↩︎

Time splitting methods are popular for fluid-structure interaction (FSI) problems (see, e.g., [1][11])). One of the first stable splitting methods was proposed by Burman and Fernández [2]. Later the same authors [4] developed a related method which was coined the genuine Robin-Robin splitting method. Both these methods however suffered from a coupling between the space and time discretization parameters that reduced the accuracy. This constraint was lifted by eliminating the mesh dependence of the Robin parameter in [12]. The resulting method has been analyzed in [12][14]. A very similar method was developed and analyzed by Bukac and Seboldt [15]. In [13] it was proved for the FSI problems that the method converges as \(\mathcal{O}(\sqrt{\Delta t})\) where \(\Delta t\) is the time step. Numerical evidence suggested that those estimates were not sharp, and nearly first-order accuracy was proved (mod possibly a logarithmic factor) in [14] for the analogue method applied to the parabolic-parabolic and wave-parabolic problems. The analysis was extended to the FSI problem in [16], [17]. For parabolic-parabolic couplings, there is a rich literature on splitting schemes motivated by models of ocean-atmosphere interaction. In these models, friction forces on the interface render the physical coupling dissipative through a Robin-type coupling condition, as discussed in [18]. This aspect has been successfully exploited in the design of splitting methods [19][26]. In our case, the coupling conditions consist of continuity of both the primal variables and the fluxes across the interface. This coupling is conservative, and hence the approach suggested in the above references fails. Instead, the splitting method uses a Robin condition for the computational coupling, which turns out to lead to an unconditionally stable algorithm. An approach for conservative fluid-fluid coupling problems was proposed in [27], using Nitsche or Robin type couplings similar to those introduced in [4]. In a domain decomposition framework, a splitting method based on subcycling, i.e., iterative solution, was proposed and analyzed in [28], [29]. In a similar spirit, but focusing on a multi-timestep approach, a Robin-Robin coupling for time-dependent advection–diffusion was introduced in [30], with numerical investigation of the stability.

It is well known that discrete time differences for time stepping methods (e.g., backward Euler method) superconverge. In particular, when the backward Euler method is applied to a parabolic problem, the first-time difference of the errors converge with order \(\mathcal{O}((\Delta t)^2)\). In this paper, we prove similar results for the splitting method [14] applied to an interface problem. In particular, we will prove second-order convergence for the scalar fields living on the two sub-domains in the \(L^2\)-norm. Moreover, in special configurations (i.e., the interface is flat and perpendicular to two sides of the domain), we will be able to prove second-order convergence in the \(H^2\)-norm. Numerically, the \(H^2\) second-order convergence rates seem to hold on more general configurations. This appears to be the first work where error estimates in stronger norms for splitting methods applied to interface problems have been considered.

Error estimates of derivatives of the solution are, of course, of interest in their own right in many applications. However, our main motivation for this work is the application of these estimates in the analysis of a prediction-correction method [31] that we are concurrently developing. The aim is to improve the first order convergence of [14] to second order convergence in time through a defect-correction procedure [32]. The method in [31] uses a prediction step, which is exactly the splitting method we analyze here. The second-order accuracy of the correction step depends on the second-order accuracy of time differences of the prediction step, which is precisely the subject of this paper.

As for a single parabolic problem, the idea to prove higher convergence for time differences is to use that the time differences satisfy a similar discrete equation with new right-hand sides that have time differences themselves. Then, one uses the error analysis for the original method to proceed. In our case, we will use the analysis provided in [14] to do this. The main difference here is that we consider a problem with Neumann boundary conditions on two of the sides of a square instead of pure Dirichlet boundary conditions, which were considered in [14]. This, in fact, simplifies the analysis slightly, and additionally, one can remove the logarithmic factor that appears in [14]. It should be mentioned that we only consider the time discrete case. The fully discrete case is more involved.

The rest of the paper is organized as follows. In Section 2, we introduce a parabolic-parabolic interface problem and the corresponding Robin-Robin coupling method. In Section 3, we present stability results for a Robin-Robin method. Section 4 is devoted to the error estimates. Finally, we provide some numerical results in Section 5 and end with some concluding remarks in Section 6.

2 The Parabolic-Parabolic interface problem and Robin-Robin coupling method↩︎

Let \(\Omega=(0,1)^2\) and suppose that \(\Omega=\Omega_f \cup \Omega_s \cup \Sigma\). The interface \(\Sigma\) is assumed to be a line segment that intersects \(\Omega\) on the two side edges; see Figure 1. We let \(\Gamma_N\) denote the two side edges of \(\Omega\) and we let \(\Gamma_D\) be the bottom and top edges of \(\Omega\). We let \(\Gamma_N^i= \Gamma_N^i \cap \partial \Omega_i\) for \(i=s,f\).

Figure 1: The domains \(\Omega_f\) and \(\Omega_s\) with interface \(\Sigma\) and Neumann boundaries.

2.1 The Parabolic-parabolic problem↩︎

We consider the interface problem \[\label{eq:ppinterface} \begin{alignat}{2} \partial_{t}\mathsf{u}-\nu_f \Delta \mathsf{u}=&0, \quad && \text{ in } [0, T] \times \Omega_f, \nonumber\\ \mathsf{u}(0, x)=& \mathsf{u}_0(x), \quad && \text{ on } \Omega_f, \\ \mathsf{u}=&0, \quad && \text{ on } [0, T] \times \Gamma^f_D,\nonumber\\ {\partial}_{\boldsymbol{n}} \mathsf{u}=&0, \quad && \text{ on } [0, T] \times \Gamma^f_N,\nonumber\\ \nonumber\\ \partial_{t}\mathsf{w}-\nu_s \Delta \mathsf{w}=&0, \quad && \text{ in } [0, T] \times \Omega_s, \nonumber\\ \mathsf{w}(0, x)=& \mathsf{w}_0(x), \quad && \text{ on } \Omega_s, \\ \mathsf{w}=&0, \quad && \text{ on } [0, T] \times \Gamma^s_D,\nonumber\\ {\partial}_{\boldsymbol{n}} \mathsf{w}=&0, \quad && \text{ on } [0, T] \times \Gamma^s_N,\nonumber\\ \nonumber\\ \mathsf{w}- \mathsf{u}=& 0, \quad && \text{ in } [0,T] \times \Sigma, \\ \nu_s \nabla \mathsf{w}\cdot \boldsymbol{n}_s + \nu_f \nabla \mathsf{u}\cdot \boldsymbol{n}_f =& 0, \quad && \text{ in } [0,T] \times \Sigma, \end{alignat}\tag{1}\] where \(\boldsymbol{n}_f\) and \(\boldsymbol{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 \(\mathsf{u}\) and \(\mathsf{w}\) are smooth on \(\Omega_f\) and \(\Omega_s\), respectively.

2.2 Variational form↩︎

Let \((\cdot, \cdot)_i\) be the \(L^2\)-inner product on \(\Omega_i\) for \(i=f, s\). Moreover, let \(\bigl\langle\cdot, \cdot \bigr\rangle\) be the \(L^2\)-inner product on \(\Sigma\). Let \(N>0\) be an integer, and define \(\Delta t:= \frac{T}{N}\), and let \(\mathsf{u}^{n} := \mathsf{u}(t_n, \cdot)\), where \(t_n := n\Delta t\) for \(n \in \{0, 1, 2, \dots, N\}\). We consider the spaces \[\label{eq:disspaces} \begin{alignat}{1} V_f:=&\{ v \in H^1(\Omega_f): v=0 \text{ on } \Gamma^f_D \}, \\ V_s:= &\{ v \in H^1(\Omega_s): v=0 \text{ on } \Gamma^s_D \}, \\ V_g:= & L^2(\Sigma). \end{alignat}\tag{2}\]

By setting \(\mathsf{l}^{n+1} := \nu_f {\partial}_{\boldsymbol{n}_f}\mathsf{u}^{n+1}\) and assuming that \(\mathsf{l}^{n+1} \in L^2(\Sigma)\) for all \(n\), the solution to 1 also satisfies the following variational formulation, for \(n=0,\ldots,N-1\): \[\tag{3} \begin{alignat}{2} (\partial_t \mathsf{w}^{n+1},z )_s+ \nu_s (\nabla \mathsf{w}^{n+1}, \nabla z)_s+ \bigl\langle\mathsf{l}^{n+1}, z\bigr\rangle=&0 , \quad && z \in V_s,\tag{4} \\ (\partial_t \mathsf{u}^{n+1}, v)_f+\nu_f (\nabla \mathsf{u}^{n+1}, \nabla v)_f- \bigl\langle\mathsf{l}^{n+1}, v\bigr\rangle=&0, \quad &&v \in V_f, \tag{5}\\ \bigl\langle\mathsf{w}^{n+1}-\mathsf{u}^{n+1}, \mu \bigr\rangle=&0, \quad && \mu \in V_g. \tag{6} \end{alignat}\]

2.3 Robin-Robin coupling: time discrete method↩︎

We define the discrete time derivatives: \[\begin{align} {1} \partial_{\Delta t}v^{n+1}=\frac{v^{n+1}-v^n}{\Delta t}, \end{align}\] and \[\begin{align} {1} \partial_{\Delta t}^2 v^{n+1}=\frac{v^{n+1}-2 v^n+ v^{n-1}}{(\Delta t)^2}. \end{align}\]

The Robin-Robin method solves sequentially: \[\label{eq1} \begin{alignat}{2} \partial_{\Delta t}w^{n+1}- \nu_s \Delta w^{n+1}=&0, \quad && \text{ on } \Omega_s, \\ w^{n+1}=&0 , \quad && \text{ on } \Gamma^s_D, \\ {\partial}_{\boldsymbol{n}} w^{n+1}=&0, \quad && \text{ on } \Gamma^s_N,\\ \alpha w^{n+1}+ \nu_s {\partial}_{\boldsymbol{n}_s} w^{n+1} = & \alpha u^{n} - \nu_f {\partial}_{\boldsymbol{n}_f} u^n \quad && \text{ on } \Sigma. \end{alignat}\tag{7}\]

\[\label{eq2} \begin{alignat}{2} \partial_{\Delta t}u^{n+1}- \nu_f\Delta u^{n+1}=&0, \quad && \text{ in } \Omega_f, \\ u^{n+1}=&0 , \quad && \text{ on } \Gamma^f_D, \\ {\partial}_{\boldsymbol{n}} u^{n+1}=&0, \quad && \text{ on } \Gamma^f_N,\\ \alpha u^{n+1}+ \nu_f {\partial}_{\boldsymbol{n}_f} u^{n+1} = & \alpha w^{n+1} + \nu_f {\partial}_{\boldsymbol{n}_f} u^n \quad && \text{ on } \Sigma. \end{alignat}\tag{8}\] We let \(\lambda^{n+1}=\nu_f {\partial}_{\boldsymbol{n}_f}u^{n+1}\). Then, the time semi-discrete solution solves the following: Find \(w^{n+1} \in V_s\), \(u^{n+1} \in V_f\), and \(\lambda^{n+1} \in V_g\) such that, for \(n\geq 0\), \[\tag{9} \begin{alignat}{2} (\partial_{\Delta t}w^{n+1}, z)_s+\nu_s(\nabla w^{n+1}, \nabla z)_s + \alpha \bigl\langle w^{n+1} - u^n, z \bigr\rangle+ \bigl\langle\lambda^n, z\bigr\rangle=&0 , \quad && z \in V_s, \tag{10}\\ (\partial_{\Delta t}u^{n+1}, v)_f+\nu_f(\nabla u^{n+1}, \nabla v)_f- \bigl\langle\lambda^{n+1}, v\bigr\rangle=&0, \quad &&v \in V_f, \tag{11}\\ \bigl\langle\alpha(u^{n+1}-w^{n+1})+(\lambda^{n+1}-\lambda^n), \mu \bigr\rangle=&0, \quad && \mu \in V_g, \tag{12} \end{alignat}\] with \(u^0=\mathsf{u}_0(x)\) and \(w^0=\mathsf{w}_0(x)\).

One can re-write 12 as \[\label{con} \alpha(u^{n+1}-w^{n+1})=\lambda^n-\lambda^{n+1}, \quad \text{ on } \Sigma.\tag{13}\]

3 Stability Result↩︎

In this section we will prove stability results for the Robin-Robin method with a more general right-hand side: Find \(w^{n+1} \in V_s\), \(u^{n+1} \in V_f\), such that, for \(n\geq 0\), \[\label{seq1} \begin{alignat}{2} \partial_{\Delta t}w^{n+1}- \nu_s \Delta w^{n+1}=&b_1^{n+1}, \quad && \text{ on } \Omega_s, \\ w^{n+1}=&0 , \quad && \text{ on } \Gamma^s_D, \\ {\partial}_{\boldsymbol{n}} w^{n+1}=&0, \quad && \text{ on } \Gamma^s_N,\\ \alpha w^{n+1}+ \nu_s {\partial}_{\boldsymbol{n}_s} w^{n+1} = & \alpha u^{n} - \nu_f {\partial}_{\boldsymbol{n}_f} u^n+ \varepsilon_1^{n+1} \quad && \text{ on } \Sigma. \end{alignat}\tag{14}\]

\[\tag{15} \begin{alignat}{2} \partial_{\Delta t}u^{n+1}- \nu_f\Delta u^{n+1}=&b_2^{n+1}, \quad && \text{ in } \Omega_f, \tag{16} \\ u^{n+1}=&0 , \quad && \text{ on } \Gamma^f_D, \\ {\partial}_{\boldsymbol{n}} u^{n+1}=&0, \quad && \text{ on } \Gamma^f_N,\\ \alpha u^{n+1}+ \nu_f {\partial}_{\boldsymbol{n}_f} u^{n+1} = & \alpha w^{n+1} + \nu_f {\partial}_{\boldsymbol{n}_f} u^n + \varepsilon_2^{n+1} \quad && \text{ on } \Sigma, \end{alignat}\] with \[\label{eq:uwinitial} u^0=w^0=b_i^0=\varepsilon_i^0=0\quad i=1,2.\tag{17}\] We assumed zero initial conditions for simplicity and this will be enough for the error analysis below. When analyzing the error of the Robin-Robin method the terms \(\{b_i^n\}\), \(\{\varepsilon_i^n\}\) will be the residual terms.

Remark 1 (Regularity). We assume that \(u^{n+1}\in H^2(\Omega_f)\) and \(w^{n+1}\in H^2(\Omega_s)\) for the remainder of this paper.

Remark 2 (Notations). For any negative superscripts, we set the values of the terms on the right-hand side to be zero, i.e., \[\label{eq:uwnegative} u^j=w^j= b_i^j=\varepsilon_i^j=0\quad j<0,\;i=1,2.\tag{18}\]

We let \(\lambda^{n+1}=\nu_f {\partial}_{\boldsymbol{n}_f}u^{n+1}\) and we see that the above solution satisfies, for \(n\ge0\), \[\tag{19} \begin{alignat}{2} (\partial_{\Delta t}w^{n+1}, z)_s+\nu_s(\nabla w^{n+1}, \nabla z)_s + \alpha \bigl\langle w^{n+1} - u^n, z \bigr\rangle+ \bigl\langle\lambda^n, z\bigr\rangle=& L_1(z) , \quad && z \in V_s, \tag{20}\\ (\partial_{\Delta t}u^{n+1}, v)_f+\nu_f(\nabla u^{n+1}, \nabla v)_f- \bigl\langle\lambda^{n+1}, v\bigr\rangle=&L_2(v), \quad &&v \in V_f, \tag{21}\\ \bigl\langle\alpha(u^{n+1}-w^{n+1})+(\lambda^{n+1}-\lambda^n), \mu \bigr\rangle=&L_3(\mu), \quad && \mu \in V_g, \tag{22} \end{alignat}\] where \[\begin{align} {1} L_1(z) :=& ( b_1^{n+1}, z)_s + \bigl\langle\varepsilon_1^{n+1}, z \bigr\rangle,\\ L_2(v):=&(b_2^{n+1},v)_f,\\ L_3(\mu):=& \bigl\langle\varepsilon_2^{n+1}, \mu \bigr\rangle. \end{align}\] For convenience we rewrite 22 as: \[\label{errorStrong} \alpha(u^{n+1} - w^{n+1} ) = \lambda^n - \lambda^{n+1}+ \varepsilon_2^{n+1}, \quad \text{ on } \Sigma.\tag{23}\]

We will need to define the following quantities for the stability estimates. \[\begin{align} {1} 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(\nu_f \|\nabla \phi^{n+1}\|_{L^2(\Omega_f)}^2+ \nu_s \|\nabla \psi^{n+1}\|_{L^2(\Omega_s)}^2) + \frac{1}{2} (\|\psi^{n+1}-\psi^n\|_{L^2(\Omega_s)}^2 + \|\phi^{n+1}-\phi^n\|_{L^2(\Omega_f)}^2)\\ & + \frac{\alpha \Delta t}{2}\|\phi^{n+1}-\phi^n+ \frac{1}{\alpha} ( \theta^{n+1}-\theta^n)\|_{L^2(\Sigma)}^2. \end{align}\]

We first state a preliminary result.

Lemma 1. Let \(w, u\) solve 14 and 15 then the following identity holds. \[\begin{align} {1}\label{eq:errorlemma1} Z^{n+1}(w, u, \lambda)+S^{n+1}(w, u, \lambda)=Z^{n}(w, u, \lambda) + \Delta tF^{n+1}(w, u) +\frac{\Delta t}{\alpha} \bigl\langle\varepsilon_2^{n+1}, \lambda^{n+1} \bigr\rangle, \end{align}\tag{24}\] where \[\begin{align} {1} F^{n+1}(w, u):= & (b_1^{n+1}, w^{n+1})_s+(b_2^{n+1}, u^{n+1})_f + \bigl\langle\varepsilon_1^{n+1}+\varepsilon_2^{n+1}, w^{n+1} \bigr\rangle+ \bigl\langle u^{n+1}-u^n, \varepsilon_2^{n+1} \bigr\rangle. \end{align}\]

Proof. To begin, we set \(z= \Delta t\,w^{n+1}\) in 20 and \(v = \Delta t\,u^{n+1}\) in 21 to get \[\begin{align} {1} & \frac{1}{2}\|w^{n+1}\|_{L^2(\Omega_s)}^2+ \frac{1}{2} \|u^{n+1}\|_{L^2(\Omega_f)}^2 + \frac{1}{2}\|w^{n+1}-w^n\|_{L^2(\Omega_s)}^2+ \frac{1}{2} \|u^{n+1}-u^n\|_{L^2(\Omega_f)}^2 \nonumber \\ & + \nu_s\Delta t\|\nabla w^{n+1}\|_{L^2(\Omega_s)}^2+ \nu_f \Delta t\|\nabla u^{n+1}\|_{L^2(\Omega_f)}^2 \label{erroraux213} \\ & = \frac{1}{2}\|w^{n}\|_{L^2(\Omega_s)}^2+ \frac{1}{2} \|u^{n}\|_{L^2(\Omega_f)}^2+\Delta t\mathsf{J}^{n+1},\nonumber \end{align}\tag{25}\] where \[\label{eq:Jn1} \begin{align} \mathsf{J}^{n+1}:=& - \alpha \bigl\langle w^{n+1} - u^n, w^{n+1} \bigr\rangle- \bigl\langle\lambda^n, w^{n+1}\bigr\rangle+ \bigl\langle\lambda^{n+1}, u^{n+1}\bigr\rangle+L_1(w^{n+1})+L_2(u^{n+1}). \end{align}\tag{26}\]

Manipulating the first three terms in 26 and using 23 , we obtain \[\begin{gather} \label{eq:pre} - \alpha \bigl\langle w^{n+1} - u^n, w^{n+1} \bigr\rangle- \bigl\langle\lambda^n, w^{n+1}\bigr\rangle+ \bigl\langle\lambda^{n+1}, u^{n+1}\bigr\rangle\\ = \mathbb{J}^{n+1} + \frac{1}{\alpha} \bigl\langle\varepsilon_2^{n+1}, \lambda^{n+1} \bigr\rangle+ \bigl\langle\varepsilon_2^{n+1}, w^{n+1} \bigr\rangle- \bigl\langle u^n-u^{n+1}, \varepsilon_2^{n+1} \bigr\rangle, \end{gather}\tag{27}\] with \[\mathbb{J}^{n+1} := \alpha \bigl\langle u^n-u^{n+1}, u^{n+1} \bigr\rangle+ \frac{1}{\alpha} \bigl\langle\lambda^n-\lambda^{n+1}, \lambda^{n+1} \bigr\rangle- \bigl\langle u^n-u^{n+1}, \lambda^n-\lambda^{n+1} \bigr\rangle.\] One can easily show that \[\begin{align} {1} \mathbb{J}^{n+1}= & \frac{\alpha}{2}(\|u^n\|_{L^2(\Sigma)}^2 - \|u^{n+1}\|_{L^2(\Sigma)}^2) + \frac{1}{2\alpha}(\|\lambda^n\|_{L^2(\Sigma)}^2 - \|\lambda^{n+1}\|_{L^2(\Sigma)}^2) \\ &- \frac{\alpha}{2}\|(u^n-u^{n+1})+ \frac{1}{\alpha}(\lambda^n-\lambda^{n+1}) \|_{L^2(\Sigma)}^2. \end{align}\] By combining this identity with 26 and 27 , we arrive at \[\begin{align} {1} \mathsf{J}^{n+1}=& \frac{\alpha}{2}(\|u^n\|_{L^2(\Sigma)}^2 - \|u^{n+1}\|_{L^2(\Sigma)}^2) + \frac{1}{2\alpha}(\|\lambda^n\|_{L^2(\Sigma)}^2 - \|\lambda^{n+1}\|_{L^2(\Sigma)}^2) \\ &- \frac{\alpha}{2}\|(u^n-u^{n+1})+ \frac{1}{\alpha}(\lambda^n-\lambda^{n+1}) \|_{L^2(\Sigma)}^2\\ & +L_1(w^{n+1}) +L_2(u^{n+1}) + \frac{1}{\alpha} \bigl\langle\varepsilon_2^{n+1}, \lambda^{n+1} \bigr\rangle+ \bigl\langle\varepsilon_2^{n+1}, w^{n+1} \bigr\rangle+\bigl\langle\varepsilon_2^{n+1}, u^{n+1}-u^{n} \bigr\rangle. \end{align}\]

If we plug in these results to 25 we arrive at the identity. ◻

We can now state an identity for the last term in 24 .

Lemma 2. Let \(w, u\) solve 14 and 15 and assuming that \(\varepsilon_2^{n+1} \in V_f\) then the following identity holds \[\begin{align} {1} \frac{\Delta t}{\alpha} \sum_{n=0}^{N-1} \bigl\langle\varepsilon_2^{n+1}, \lambda^{n+1} \bigr\rangle= & - \frac{\Delta t}{\alpha} \sum_{n=1}^{N-1} (u^n, \partial_{\Delta t}\varepsilon_2^{n+1})_f+\frac{\Delta t}{\alpha}\sum_{n=0}^{N-1} \Big(\nu_f (\nabla u^{n+1}, \nabla \varepsilon_2^{n+1})_f- (b_2^{n+1}, \varepsilon_2^{n+1})_f \Big) \\ &+ \frac{1}{\alpha} (u^N, \varepsilon_2^N)_f. \end{align}\]

Proof. We first take \(v = \Delta t \varepsilon_2^{n+1}\) in 21 to obtain \[\begin{align} {1} \Delta t\bigl\langle\varepsilon_2^{n+1}, \lambda^{n+1} \bigr\rangle= (u^{n+1}-u^n, \varepsilon_2^{n+1})_f +\Delta t\nu_f(\nabla u^{n+1}, \nabla \varepsilon_2^{n+1})_f- \Delta t(b_2^{n+1}, \varepsilon_2^{n+1})_f. \end{align}\] If we take the sum over \(n=0,\ldots,N-1\) and use summation by parts, we get \[\begin{align} {1} \Delta t\sum_{n=0}^{N-1} \bigl\langle\varepsilon_2^{n+1}, \lambda^{n+1} \bigr\rangle= & - \Delta t\sum_{n=1}^{N-1} (u^n, \partial_{\Delta t}\varepsilon_2^{n+1})_f+\Delta t\sum_{n=0}^{N-1} \Big( \nu_f (\nabla u^{n+1}, \nabla \varepsilon_2^{n+1})_f- (b_2^{n+1}, \varepsilon_2^{n+1})_f \Big) \\ &+ (u^N, \varepsilon_2^N)_f-(u^0, \varepsilon_2^1). \end{align}\] We conclude the proof by using 17 . ◻

To state the stability estimate we need the next definition. \[\begin{align} {1} \Xi^N(m_1, m_2, s_1, s_2):= &\Delta t\sum_{n=0}^{N-1} \left[\frac{1}{\nu_s}\| m_1^{n+1}\|_{L^2(\Omega_s)}^2+\left(\frac{1}{\nu_f}+1\right)\| m_2^{n+1}\|_{L^2(\Omega_f)}^2\right] \\ &+ \Delta t\sum_{n=0}^{N-1} \Big( \frac{1}{\nu_f\alpha^2} \| \partial_{\Delta t}s_2^{n+1}\|_{L^2(\Omega_f)}^2 + \frac{\nu_f}{\alpha^2} \| \nabla s_2^{n+1}\|_{L^2(\Omega_f)}^2 + \frac{1}{\alpha}\| s_2^{n+1} \|_{L^2(\Omega_f)}^2 \Big) \\ &+ \Delta t\sum_{n=0}^{N-1} \Big(\frac{1}{\nu_s} \| s_1^{n+1}+s_2^{n+1}\|_{L^2(\Sigma)}^2+\frac{1}{\nu_f}\| s_2^{n+1}\|_{L^2(\Sigma)}^2 \Big)+ \frac{1}{\alpha^2}\|s_2^N\|_{L^2(\Omega_f)}^2. \end{align}\]

Theorem 3. Let \(w, u\) solve 14 and 15 and assuming that \(\varepsilon_2^{n+1} \in V_f\) then the following estimate holds \[\begin{align} {1} Z^N(w, u, \lambda)+ \sum_{n=0}^{N-1} S^{n+1}(w, u, \lambda) \le C \Xi^N(b_1, b_2, \varepsilon_1, \varepsilon_2). \end{align}\]

Proof. Using Lemma 1 and taking the sum we get \[\begin{align} {1} Z^N(w, u, \lambda)+ \sum_{n=0}^{N-1} S^{n+1}(w, u, \lambda) = Z^0(w, u, \lambda) + \Delta t\sum_{n=0}^{N-1} F^{n+1}(w, u) + \frac{\Delta t}{\alpha} \sum_{n=0}^{N-1} \bigl\langle\varepsilon_2^{n+1}, \lambda^{n+1} \bigr\rangle. \end{align}\]

Using the Poincaré and trace inequalities, we easily obtain the following bound \[\begin{align} {1} \Delta t\sum_{n=0}^{N-1} F^{n+1}(w, u) \le & \frac{1}{8} \sum_{n=0}^{N-1} S^{n+1}(w, u, \lambda) + C \Delta t\sum_{n=0}^{N-1} (\frac{1}{\nu_s}\| b_1^{n+1}\|_{L^2(\Omega_s)}^2+\frac{1}{\nu_f}\| b_2^{n+1}\|_{L^2(\Omega_f)}^2) \\ & +C \Delta t\sum_{n=0}^{N-1} \Big(\frac{1}{\nu_s} \| \varepsilon_1^{n+1}+\varepsilon_2^{n+1}\|_{L^2(\Sigma)}^2+\frac{1}{\nu_f}\| \varepsilon_2^{n+1}\|_{L^2(\Sigma)}^2\Big). \end{align}\]

Using Lemma 2 and the Poincaré inequality we have \[\begin{align} {1} \frac{\Delta t}{\alpha} \sum_{n=0}^{N-1} \bigl\langle\varepsilon_2^{n+1}, \lambda^{n+1} \bigr\rangle\le & \frac{1}{4} \| u^N\|_{L^2(\Omega_f)}^2+\frac{1}{8} \sum_{n=0}^{N-1} S^{n+1}(w, u, \lambda) \\ &+ C \Delta t\sum_{n=0}^{N-1} \Big( \frac{1}{\nu_f\alpha^2} \| \partial_{\Delta t}\varepsilon_2^{n+1}\|_{L^2(\Omega_f)}^2 + \frac{\nu_f}{\alpha^2} \| \nabla \varepsilon_2^{n+1}\|_{L^2(\Omega_f)}^2 + \frac{1}{\alpha}\|\varepsilon_2^{n+1} \|_{L^2(\Omega_f)}^2 \Big) \\ &+ C \frac{\Delta t}{\alpha} \sum_{n=0}^{N-1} \| b_2^{n+1} \|_{L^2(\Omega_f)}^2 +C \frac{1}{\alpha^2}\|\varepsilon_2^N\|_{L^2(\Omega_f)}^2. \end{align}\] It follows from 17 and the definition of \(Z^0\) that \(Z^0(w, u, \lambda)=0\). We finish the proof by combining the above estimates. ◻

Due to Remark 2, the discrete time derivative of \(u\) and \(w\) solves 14 and 15 with discrete time derivatives of the data as the right-hand sides, we immediately get the following.

Corollary 1. Let \(w, u\) solve 14 and 15 and assuming that \(\partial_{\Delta t}\varepsilon_2^{n+1} \in V_f\) then \[\begin{align} {1} Z^N(\partial_{\Delta t}w, \partial_{\Delta t}u, \partial_{\Delta t}\lambda)+ \sum_{n=0}^{N-1} S^{n+1}(\partial_{\Delta t}w, \partial_{\Delta t}u, \partial_{\Delta t}\lambda) \le & C \Xi^N(\partial_{\Delta t}b_1, \partial_{\Delta t}b_2, \partial_{\Delta t}\varepsilon_1, \partial_{\Delta t}\varepsilon_2) \end{align}\]

Similarly, the second-order discrete time derivative of \(u\) and \(w\) solves 14 and 15 with the second-order discrete time derivatives of the data as the right-hand sides, we also have the following.

Corollary 2. Let \(w, u\) solve 14 and 15 and assuming that \(\partial_{\Delta t}^2\varepsilon_2^{n+1} \in V_f\) then \[\begin{align} {1} Z^N(\partial_{\Delta t}^2 w, \partial_{\Delta t}^2 u, \partial_{\Delta t}^2 \lambda)+ \sum_{n=0}^{N-1} S^{n+1}(\partial_{\Delta t}^2w, \partial_{\Delta t}^2 u, \partial_{\Delta t}^2\lambda) \le & C \Xi^N(\partial_{\Delta t}^2 b_1, \partial_{\Delta t}^2 b_2, \partial_{\Delta t}^2 \varepsilon_1, \partial_{\Delta t}^2 \varepsilon_2) \end{align}\]

Remark 4. We can also analyze the problem 14 -15 but replacing the homogeneous Neumann boundary conditions with homogeneous Dirichlet boundary conditions. In other words, the problem with pure Dirichlet boundary conditions that is \(\Gamma_D^i=\partial \Omega_i \backslash \Sigma\) for \(i=s,f\). In this case exactly the same estimates hold. Now of course, we assume \(\varepsilon_2^{n} \in V_f\) where \(V_f\) has zero Dirichlet boundary conditions on \(\Gamma_D^f=\partial \Omega_f \backslash \Sigma\).

3.1 \(H^2\) stability in a special case↩︎

In this section we prove \(H^2\) estimates for \(u\) in a special configuration. Again, we assume that \(\Omega=(0,1)^2\), and now assume \(\Sigma\) is parallel to the \(x\)-axis (see Figure 2). We take advantage of the fact that the sides composing \(\Gamma_N\) are perpendicular to the \(x\)-axis and, moreover, we also notice that \(\Sigma\) is parallel to the \(x\)-axis.

Figure 2: The domains \(\Omega_f\) and \(\Omega_s\) with horizontal interface \(\Sigma\).

Then, in this particular case we have, for \(n\ge0\), \[\label{seq1x} \begin{alignat}{2} \partial_{\Delta t}{\partial}_x w^{n+1}- \nu_s \Delta \partial_x w^{n+1}=&{\partial}_x b_1^{n+1}, \quad && \text{ on } \Omega_s, \\ {\partial}_x w^{n+1}=&0 , \quad && \text{ on } \Gamma^s_D, \\ {\partial}_x w^{n+1}=&0, \quad && \text{ on } \Gamma^s_N,\\ \alpha{\partial}_xw^{n+1}+ \nu_s {\partial}_{\boldsymbol{n}_f} {\partial}_x w^{n+1} = & \alpha {\partial}_x u^{n} - \nu_f {\partial}_{\boldsymbol{n}_f} {\partial}_x u^n+ {\partial}_x \varepsilon_1^{n+1} \quad && \text{ on } \Sigma. \end{alignat}\tag{28}\]

\[\tag{29} \begin{alignat}{2} \partial_{\Delta t}{\partial}_x u^{n+1}- \nu_f\Delta {\partial}_x u^{n+1}=&{\partial}_x b_2^{n+1}, \quad && \text{ in } \Omega_f, \tag{30} \\ {\partial}_x u^{n+1}=&0 , \quad && \text{ on } \Gamma^s_D, \\ {\partial}_x u^{n+1}=&0, \quad && \text{ on } \Gamma^s_N,\\ \alpha {\partial}_x u^{n+1}+ \nu_f {\partial}_{\boldsymbol{n}_f} {\partial}_x u^{n+1} = & \alpha {\partial}_x w^{n+1} + \nu_f {\partial}_{\boldsymbol{n}_f} {\partial}_x u^n + {\partial}_x \varepsilon_2^{n+1} \quad && \text{ on } \Sigma, \end{alignat}\] with \[{\partial}_x w^0={\partial}_xu^0={\partial}_xb_i^0={\partial}_x\varepsilon_i^0=0\quad i=1,2,\] according to 17 .

We then get an immediate Corollary from Remark 4 and Theorem 3.

Corollary 3. Suppose that \(\Sigma\) is perpendicular to the two sides of \(\Gamma_N\) as in Figure 2. Let \(w, u\) solve 14 and 15 and assuming that \({\partial}_{x} \varepsilon_2^{n} \in \{ v \in H^1(\Omega_f): v=0 \text{ on } {\partial}\Omega_f \backslash \Sigma\}\) on \(\Gamma_N^f\) then the following estimate holds \[\begin{align} {1} Z^N({\partial}_x w, {\partial}_x u, {\partial}_x \lambda)+ \sum_{n=0}^{N-1} S^{n+1}({\partial}_x w, {\partial}_x u, {\partial}_x \lambda) \le C \Xi^N({\partial}_x b_1, {\partial}_x b_2, {\partial}_x \varepsilon_1, {\partial}_x\varepsilon_2). \end{align}\]

Corollary 4. Under the hypothesis of Corollaries 3 and 1 we have \[\begin{align} {1} \Delta t\sum_{n=0}^{N-1} \nu_f \|D^2 u^{n+1}\|_{L^2(\Omega_f)}^2 \le & C \Big( \Xi^N(\partial_{\Delta t}b_1, \partial_{\Delta t}b_2, \partial_{\Delta t}\varepsilon_1, \partial_{\Delta t}\varepsilon_2)+ \Xi^N({\partial}_x b_1, {\partial}_x b_2, {\partial}_x \varepsilon_1, {\partial}_x\varepsilon_2)\Big) \\ &+ C \Delta t\sum_{n=0}^{N-1} \nu_f \|b_1^{n+1}\|_{L^2(\Omega_f)}^2. \end{align}\]

Proof. From Corollary 3 we get \[\begin{align} {1} \Delta t\sum_{n=0}^{N-1} \nu_f \|\nabla {\partial}_x u^{n+1}\|_{L^2(\Omega_f)}^2 \le C \Xi^N({\partial}_x b_1, {\partial}_x b_2, {\partial}_x \varepsilon_1, {\partial}_x\varepsilon_2). \end{align}\] Moreover, using 16 and Corollary 1, we have \[\begin{align} {1} \Delta t\sum_{n=0}^{N-1} \nu_f \|\Delta u^{n+1}\|_{L^2(\Omega_f)}^2 = & \Delta t\sum_{n=0}^{N-1} \nu_f \|\partial_{\Delta t}u^{n+1}-b_1^{n+1}\|_{L^2(\Omega_f)}^2 \\ \le & C \Xi^N(\partial_{\Delta t}b_1, \partial_{\Delta t}b_2, \partial_{\Delta t}\varepsilon_1, \partial_{\Delta t}\varepsilon_2)+ 2 \Delta t\sum_{n=0}^{N-1} \nu_f \|b_1^{n+1}\|_{L^2(\Omega_f)}^2. \end{align}\] Finally, using the following estimate, \[\begin{align} {1} \Delta t\sum_{n=0}^{N-1} \nu_f \|{\partial}_y^2 u^{n+1}\|_{L^2(\Omega_f)}^2 \le 2 \Delta t\sum_{n=0}^{N-1} \nu_f \Big( \|{\partial}_x^2 u^{n+1}\|_{L^2(\Omega_f)}^2+ \|\Delta u^{n+1}\|_{L^2(\Omega_f)}^2 \Big), \end{align}\] and combining the above estimates we obtain the result. ◻

4 Error estimates of the Robin-Robin method↩︎

In this section we apply the stability results of the previous sections to obtain error estimates of the Robin-Robin splitting method 7 -8 applied to 1 . We use the following notation for the errors : \[\begin{align} {1} U^n:=& \mathsf{u}^n-u^n, \quad W^n:=\mathsf{w}^n-w^n, \quad \Lambda^{n}:=\mathsf{l}^n-\lambda^n. \end{align}\] We use the convention \(u^j=\mathsf{u}^j=\mathsf{u}_0\) and \(w^j=\mathsf{w}^j=\mathsf{w}_0\) for \(j<0\).

Then the error equations read, for \(n\ge0\),

\[\label{erroreq1} \begin{alignat}{2} \partial_{\Delta t}W^{n+1}- \nu_s \Delta W^{n+1}=&-h_1^{n+1}, \quad && \text{ on } \Omega_s, \\ W^{n+1}=&0 , \quad && \text{ on } \Gamma^s_D, \\ {\partial}_{\boldsymbol{n}} W^{n+1}=&0, \quad && \text{ on } \Gamma^s_N,\\ \alpha W^{n+1}+ \nu_s {\partial}_{\boldsymbol{n}_s} W^{n+1} = & \alpha U^{n} -\nu_f {\partial}_{\boldsymbol{n}_f} U^n + \alpha g_1^{n+1}-g_2^{n+1} \quad && \text{ on } \Sigma. \end{alignat}\tag{31}\]

\[\label{erroreq2} \begin{alignat}{2} \partial_{\Delta t}U^{n+1}- \nu_f\Delta U^{n+1}=& -h_2^{n+1}, \quad && \text{ in } \Omega_f, \\ U^{n+1}=&0 , \quad && \text{ on } \Gamma^s_D, \\ {\partial}_{\boldsymbol{n}} U^{n+1}=&0, \quad && \text{ on } \Gamma^s_N,\\ \alpha U^{n+1}+ \nu_f {\partial}_{\boldsymbol{n}_f} U^{n+1} = & \alpha W^{n+1} + \nu_f {\partial}_{\boldsymbol{n}_f} U^n+ g_2^{n+1} \quad && \text{ on } \Sigma. \end{alignat}\tag{32}\] where \[\begin{align} {2} h_1^{n+1} &:= \partial_{t}\mathsf{w}^{n+1} - \partial_{\Delta t}\mathsf{w}^{n+1}, \quad && g_1^{n+1} := \mathsf{u}^{n+1}-\mathsf{u}^n, \\ h_2^{n+1} &:= \partial_{t}\mathsf{u}^{n+1} - \partial_{\Delta t}\mathsf{u}^{n+1}, \quad && g_2^{n+1} := \mathsf{l}^{n+1} - \mathsf{l}^n. \end{align}\]

The equations 31 -32 are well-defined and we also have \[W^0=U^0=h_i^0=g_i^0=0,\;i=1,2,\] \[\partial_{\Delta t}W^0=\partial_{\Delta t}U^0=\partial_{\Delta t}h_i^0=\partial_{\Delta t}g_i^0=0,\;i=1,2,\] \[\partial_{\Delta t}^2 W^0=\partial_{\Delta t}^2 U^0=\partial_{\Delta t}^2 h_i^0=\partial_{\Delta t}^2 g_i^0=0,\;i=1,2.\]

We then extend \(\mathsf{l}\) to \(\Omega_f\) in a natural way. We let \(\tilde{\mathsf{l}}= \phi \nu_f\nabla \mathsf{u}\cdot \boldsymbol{n}_f\) where \(\phi\) is a function that is one on \(\Sigma\) and vanishes on \(\Gamma_D^f\). Then, we define \(\tilde{g}_2^{n+1}= \tilde{\mathsf{l}}^{n+1} - \tilde{\mathsf{l}}^n\). By construction \(\tilde{g}_2^{n+1} \in V_f\) and \(\tilde{g}_2^n= g_2^{n}\) on \(\Sigma\). We immediately get the following result if we apply Theorem 3 and Corollary 1.

Corollary 5. Let \(\mathsf{u}, \mathsf{w}\) solve 1 and \(w, u\) solve 7 and 8 then \[\begin{align} {1}\label{eq:WUerror} Z^N(W, U, \Lambda)+ \sum_{n=0}^{N-1} S^{n+1}(W, U, \Lambda) \le C \Xi^N(h_1, h_2, \alpha g_1-g_2,\tilde{g}_2), \end{align}\tag{33}\] \[\label{eq:WUdterror} \begin{align} & Z^N(\partial_{\Delta t}W, \partial_{\Delta t}U, \partial_{\Delta t}\Lambda)+ \sum_{n=0}^{N-1} S^{n+1}(\partial_{\Delta t}W, \partial_{\Delta t}U, \partial_{\Delta t}\Lambda) \\ &\le C \Xi^N(\partial_{\Delta t}h_1, \partial_{\Delta t}h_2, \alpha \partial_{\Delta t}g_1-\partial_{\Delta t}g_2, \partial_{\Delta t}\tilde{g}_2), \end{align}\tag{34}\] and \[\label{eq:WUdt2error} \begin{align} & Z^N(\partial_{\Delta t}^2 W, \partial_{\Delta t}^2 U, \partial_{\Delta t}^2 \Lambda)+ \sum_{n=0}^{N-1} S^{n+1}(\partial_{\Delta t}^2 W, \partial_{\Delta t}^2 U, \partial_{\Delta t}^2 \Lambda) \\ &\le C \Xi^N(\partial_{\Delta t}^2 h_1, \partial_{\Delta t}^2 h_2, \alpha \partial_{\Delta t}^2 g_1-\partial_{\Delta t}^2 g_2, \partial_{\Delta t}^2 \tilde{g}_2). \end{align}\tag{35}\]

Then, it is quite straightforward to get a convergence rate by estimating the right-hand sides. The proof of the following Corollary can be found in Appendix 7.

Corollary 6. Let \(\mathsf{u}, \mathsf{w}\) solve 1 and \(w, u\) solve 7 and 8 then \[\begin{align} {1}\label{eq:errorW} Z^N(W, U, \Lambda)+ \sum_{n=0}^{N-1} S^{n+1}(W, U, \Lambda) \le C (\Delta t)^2 \mathsf{Y}, \end{align}\tag{36}\] \[\begin{align} {1}\label{eq:errordW} & Z^N(\partial_{\Delta t}W, \partial_{\Delta t}U, \partial_{\Delta t}\Lambda)+ \sum_{n=0}^{N-1} S^{n+1}(\partial_{\Delta t}W, \partial_{\Delta t}U, \partial_{\Delta t}\Lambda) \le C (\Delta t)^2 \mathcal{Y}, \end{align}\tag{37}\] and \[\begin{align} {1}\label{eq:errorddW} & Z^N(\partial_{\Delta t}^2 W, \partial_{\Delta t}^2 U, \partial_{\Delta t}^2 \Lambda)+ \sum_{n=0}^{N-1} S^{n+1}(\partial_{\Delta t}^2 W, \partial_{\Delta t}^2 U, \partial_{\Delta t}^2 \Lambda) \le C (\Delta t)^2 \mathfrak{Y}, \end{align}\tag{38}\] where \(\mathsf{Y}\) is defined as \[\begin{align} \mathsf{Y}:=& \frac{1}{\nu_s}\|\partial_{t}^2 \mathsf{w}\|_{L^2(0,T;L^2(\Omega_s))}^2 + (\frac{1}{\nu_f}+1)\|\partial_{t}^2 \mathsf{u}\|_{L^2(0,T;L^2(\Omega_f))}^2) \\ & + (\frac{\nu_f}{\alpha^2}+\frac{\nu_f^2}{\alpha}) \|\partial_{t}\mathsf{u}\|_{L^2(0,T;H^1(\Omega_f))}^2+ \frac{(\nu_f)^3}{\alpha^2} \|\partial_{t}\mathsf{u}\|_{L^2(0,T;H^2(\Omega_f))}^2\\ & +\frac{\alpha^2}{\nu_s}\|\partial_{t}\mathsf{u}\|_{L^2(0,T;L^2(\Sigma))}^2 + \frac{1}{\nu_f}\|\partial_{t}\mathsf{l}\|_{L^2(0,T;L^2(\Sigma))}^2+\frac{\nu_f^2}{\alpha^2} \|\partial_{t}\mathsf{u}\|_{L^\infty(0,T;H^1(\Omega_f))}^2, \end{align}\] \(\mathcal{Y}\) is defined as \[\label{eq:caly} \begin{align} \mathcal{Y}:=& \frac{1}{\nu_s}\|\partial_{t}^3 \mathsf{w}\|_{L^2(0,T;L^2(\Omega_s))}^2 + (\frac{1}{\nu_f}+1)\|\partial_{t}^3 \mathsf{u}\|_{L^2(0,T;L^2(\Omega_f))}^2) \\ & + (\frac{\nu_f}{\alpha^2}+\frac{\nu_f^2}{\alpha}) \|\partial_{t}^2 \mathsf{u}\|_{L^2(0,T;H^1(\Omega_f))}^2+ \frac{(\nu_f)^3}{\alpha^2} \|\partial_{t}^2 \mathsf{u}\|_{L^2(0,T;H^2(\Omega_f))}^2\\ & +\frac{\alpha^2}{\nu_s}\|\partial_{t}^2 \mathsf{u}\|_{L^2(0,T;L^2(\Sigma))}^2 + \frac{1}{\nu_f}\|\partial_{t}^2 \mathsf{l}\|_{L^2(0,T;L^2(\Sigma))}^2+\frac{\nu_f^2}{\alpha^2} \|\partial_{t}^2 \mathsf{u}\|_{L^\infty(0,T;H^1(\Omega_f))}^2, \end{align}\tag{39}\] and \(\mathfrak{Y}\) is defined as \[\label{eq:frakY} \begin{align} \mathfrak{Y}:=& \frac{1}{\nu_s}\|\partial_{t}^4 \mathsf{w}\|_{L^2(0,T;L^2(\Omega_s))}^2 + (\frac{1}{\nu_f}+1)\|\partial_{t}^4 \mathsf{u}\|_{L^2(0,T;L^2(\Omega_f))}^2) \\ & + (\frac{\nu_f}{\alpha^2}+\frac{\nu_f^2}{\alpha}) \|\partial_{t}^3 \mathsf{u}\|_{L^2(0,T;H^1(\Omega_f))}^2+ \frac{(\nu_f)^3}{\alpha^2} \|\partial_{t}^3 \mathsf{u}\|_{L^2(0,T;H^2(\Omega_f))}^2\\ & +\frac{\alpha^2}{\nu_s}\|\partial_{t}^3 \mathsf{u}\|_{L^2(0,T;L^2(\Sigma))}^2 + \frac{1}{\nu_f}\|\partial_{t}^3 \mathsf{l}\|_{L^2(0,T;L^2(\Sigma))}^2+\frac{\nu_f^2}{\alpha^2} \|\partial_{t}^3 \mathsf{u}\|_{L^\infty(0,T;H^1(\Omega_f))}^2. \end{align}\tag{40}\]

4.1 \(H^2\) error estimates in a special case↩︎

Here we assume that we have the configuration as in Figure 1. Then, we see that \(\partial_x \tilde{g_2}^{n+1}\) vanishes on \(\Gamma_N^f\) and \(\Gamma_D^f\) and hence belongs to \(V_f\). Hence, Corollary 4 gives the following corollary.

Corollary 7. Suppose that \(\Sigma\) is perpendicular to the two sides of \(\Gamma_N\) as in Figure 2. Let \(\mathsf{u}, \mathsf{w}\) solve 1 and \(w, u\) solve 7 and 8 then \[\begin{align} {1} \Delta t\sum_{n=0}^{N-1} \nu_f \|D^2 (U^{n+1})\|_{L^2(\Omega_f)}^2 \le & C \Xi^N(\partial_{\Delta t}h_1, \partial_{\Delta t}h_2, \partial_{\Delta t}(\alpha g_1-g_2), \partial_{\Delta t}\tilde{g}_2) \\ & + C \Xi^N({\partial}_x h_1, {\partial}_x h_2, {\partial}_x (\alpha g_1-g_2), \partial_{\Delta t}{\partial}_x \tilde{g}_2) \\ &+ C \Delta t\sum_{n=0}^{N-1} \nu_f \| h_1^{n+1}\|_{L^2(\Omega_f)}^2. \end{align}\]

Since \(\partial_{\Delta t}U\) and \(\partial_{\Delta t}W\) satisfy the same equations as \(U, W\) with time difference right-hand sides, we have

Corollary 8. Suppose that \(\Sigma\) is perpendicular to the two sides of \(\Gamma_N\) as in Figure 2. Let \(\mathsf{u}, \mathsf{w}\) solve 1 and \(w, u\) solve 7 and 8 , then \[\label{eq:h2duesti} \begin{align} \Delta t\sum_{n=0}^{N-1} \nu_f \|D^2 (\partial_{\Delta t}U^{n+1})\|_{L^2(\Omega_f)}^2 \le & C \Xi^N(\partial_{\Delta t}^2 h_1, \partial_{\Delta t}^2 h_2, \partial_{\Delta t}^2 (\alpha g_1-g_2), \partial_{\Delta t}^2 \tilde{g}_2) \\ & + C \Xi^N(\partial_{\Delta t}{\partial}_x h_1, \partial_{\Delta t}{\partial}_x h_2, \partial_{\Delta t}{\partial}_x (\alpha g_1-g_2), \partial_{\Delta t}{\partial}_x \tilde{g}_2) \\ &+ C \Delta t\sum_{n=0}^{N-1} \nu_f \| \partial_{\Delta t}h_1^{n+1}\|_{L^2(\Omega_f)}^2. \end{align}\tag{41}\]

Then, it is straight-forward to obtain the convergence rate by estimate the right-hand sides. See Appendix 8 for a proof for the following Corollary.

Corollary 9. Suppose that \(\Sigma\) is perpendicular to the two sides of \(\Gamma_N\) as in Figure 2. Let \(\mathsf{u}, \mathsf{w}\) solve 1 and \(w, u\) solve 7 and 8 then \[\begin{align} {1} \Delta t\sum_{n=0}^{N-1} \nu_f \|D^2 (\partial_{\Delta t}U^{n+1})\|_{L^2(\Omega_f)}^2 \le C(\Delta t)^2 (\mathfrak{Y}+\mathbb{Y}+\nu_f\|\partial_{t}^3\mathsf{w}\|^2_{L^2((0,T),L^2(\Omega_s))}) \end{align}\] where \(\mathfrak{Y}\) is defined in Corollary 6 and \(\mathbb{Y}\) is defined as, \[\label{eq:bbY} \begin{align} \mathbb{Y}:=& \frac{1}{\nu_s}\|{\partial}_x\partial_{t}^3 \mathsf{w}\|_{L^2(0,T;L^2(\Omega_s))}^2 + (\frac{1}{\nu_f}+1)\|{\partial}_x\partial_{t}^3 \mathsf{u}\|_{L^2(0,T;L^2(\Omega_f))}^2) \\ & + (\frac{\nu_f}{\alpha^2}+\frac{\nu_f^2}{\alpha}) \|{\partial}_x\partial_{t}^2 \mathsf{u}\|_{L^2(0,T;H^1(\Omega_f))}^2+ \frac{(\nu_f)^3}{\alpha^2} \|{\partial}_x\partial_{t}^2 \mathsf{u}\|_{L^2(0,T;H^2(\Omega_f))}^2\\ & +\frac{\alpha^2}{\nu_s}\|{\partial}_x\partial_{t}^2 \mathsf{u}\|_{L^2(0,T;L^2(\Sigma))}^2 + \frac{1}{\nu_f}\|{\partial}_x\partial_{t}^2 \mathsf{l}\|_{L^2(0,T;L^2(\Sigma))}^2+\frac{\nu_f^2}{\alpha^2} \|{\partial}_x\partial_{t}^2 \mathsf{u}\|_{L^\infty(0,T;H^1(\Omega_f))}^2. \end{align}\tag{42}\]

5 Numerical Experiments↩︎

In this section we provide numerical experiments that agree with our theoretical results. We let \[\begin{align} {3} e_u&= \|U^N\|_{L^2(\Omega_f)},\quad e_{1,u}= \|U^N-U^{N-1}\|_{L^2(\Omega_f)},\\ e_{2,u}&= \|(U^N-U^{N-1})-(U^{N-1}-U^{N-2})\|_{L^2(\Omega_f)},\\ e_{\lambda}&=\|\Lambda^N\|_{L^2(\Sigma)},\quad e_{1,\lambda}=\|\Lambda^N-\Lambda^{N-1}\|_{L^2(\Sigma)}\\ e_{1,u,2}&=\|U^N-U^{N-1}\|_{H^2(\Omega_f)}. \end{align}\]

All numerical experiments are performed using FEniCS and multiphenics [33], [34]. Although we only analyze the semi-discrete method, here we present the results for a fully discrete method where we use the piecewise linear finite element method for the spatial discretization, except for the computation of \(e_{1,u,2}\), where we use the piecewise quadratic finite element method because piecewise linear function do not approximation a function well in \(H^2\). In addition, we also present convergence rates for the Lagrange multiplier.

Example 1. We consider the domain \(\Omega=(0,1)^2\), \(\Omega_f =(0,1)\times (0,.75)\) and \(\Omega_s=(0,1) \times (.75,1)\). See Figure 2 for an illustration. We take \(\nu_f=1=\nu_s\) and take the solution of 3 to be \[\mathsf{w}=\mathsf{u}= e^{-2\pi^2 t} \cos(\pi x_1) \sin(\pi x_2).\]

We take \(h=\Delta t\), \(T=0.25\) and \(\alpha=4\) where \(h\) is the mesh size of the triangulation.

Table 1: Errors and convergence rates of \(U^N\) for Example 1
\(\Delta t\) \(e_u\) rates \(e_{1,u}\) rates \(e_{2,u}\) rates
\((1/2)^2\) 7.65e-02 - 7.73e-02 - 6.57e+01 -
\((1/2)^3\) 3.72e-02 1.04 5.87e-02 0.40 1.55e-01 8.73
\((1/2)^4\) 1.74e-02 1.10 1.47e-02 1.99 7.91e-03 4.29
\((1/2)^5\) 7.95e-03 1.13 3.58e-03 2.04 1.27e-03 2.64
\((1/2)^6\) 3.52e-03 1.17 8.41e-04 2.09 1.75e-04 2.85
\((1/2)^7\) 1.62e-03 1.12 1.96e-04 2.10 2.15e-05 3.03
\((1/2)^8\) 7.70e-04 1.07 4.69e-05 2.06 2.62e-06 3.04
\((1/2)^9\) 3.75e-04 1.04 1.14e-05 2.04 3.22e-07 3.02
Table 2: Error and convergence rates of \(\Lambda^N\) for Example 1
\(\Delta t\) \(e_{\lambda}\) rates \(e_{1,\lambda}\) rates \(e_{1,u,2}\) rates
\((1/2)^2\) 2.55e-01 - 2.55e-01 - 1.33e+01 -
\((1/2)^3\) 9.73e-02 1.39 2.41e-01 0.08 1.11e+01 0.26
\((1/2)^4\) 3.06e-02 1.67 4.41e-02 2.45 2.61e-01 2.08
\((1/2)^5\) 1.62e-02 0.92 6.69e-03 2.72 5.76e-02 2.18
\((1/2)^6\) 9.20e-03 0.82 2.06e-03 1.70 1.40e-02 2.04
\((1/2)^7\) 4.55e-03 1.01 5.28e-04 1.97 3.32e-03 2.07
\((1/2)^8\) 2.23e-03 1.03 1.31e-04 2.01 8.03e-04 2.05
\((1/2)^9\) 1.10e-03 1.02 3.26e-05 2.01 1.97e-04 2.03

As we can see from Tables 1-2, the \(L^2\) error at the final step, \(e_u\) is of order \((\Delta t)\) whereas the difference of two consecutive errors, \(e_{1,u}\) is \((\Delta t)^2\) and the second difference, \(e_{2,u}\) is \((\Delta t)^3\). The \(L^2\) error of the Lagrange multiplier at the final time, \(e_{\lambda}\) is of order \(\Delta t\) and the difference \(e_{1,\lambda}\), is of order \((\Delta t)^2\). It is also clear that the \(H^2\) error of the difference of \(U^N\), \(e_{1,u,2}\) is of order \((\Delta t)^2\) as we proved in Corollary 9.

Example 2. In this example, we test our algorithm for a non-horizontal interface problem. We consider the domain \(\Omega=(0,1)^2\) and we let \(\Sigma\) be defined as the straight line connecting \((0,0.25)\) and \((1,0.75)\). We then define \(\Omega_s\) as the region above \(\Sigma\) and \(\Omega_f\) as the region below \(\Sigma\). We take \(\nu_f=1=\nu_s\) and take the solution of 3 to be \[\mathsf{w}=\mathsf{u}= e^{-2\pi^2 t} \cos(\pi x_1) \sin(\pi x_2).\] Other parameters are identical to those of Example 1.

We report the convergence results in Tables 3-4. We again observe expected convergence rates for both \(U^N\) and \(\Lambda^N\). It indicates that our methods also work for a more general interface problem.

Table 3: Errors and convergence rates of \(U^N\) for Example 2
\(\Delta t\) \(e_u\) rates \(e_{1,u}\) rates \(e_{2,u}\) rates
\((1/2)^2\) 8.05e-02 - 9.64e-02 - 4.87e+01 -
\((1/2)^3\) 5.10e-02 0.66 4.23e-02 1.19 1.36e-01 8.48
\((1/2)^4\) 2.46e-02 1.05 1.57e-02 1.43 4.94e-03 4.78
\((1/2)^5\) 9.20e-03 1.42 4.07e-03 1.95 1.47e-03 1.74
\((1/2)^6\) 3.52e-03 1.38 8.46e-04 2.27 1.82e-04 3.02
\((1/2)^7\) 1.52e-03 1.22 1.86e-04 2.18 2.11e-05 3.11
\((1/2)^8\) 7.00e-04 1.11 4.33e-05 2.10 2.49e-06 3.08
\((1/2)^9\) 3.36e-04 1.06 1.04e-05 2.06 3.02e-07 3.04
Table 4: Error and convergence rates of \(\Lambda^N\) for Example 2
\(\Delta t\) \(e_{\lambda}\) rates \(e_{1,\lambda}\) rates \(e_{1,u,2}\) rates
\((1/2)^2\) 8.51e-01 - 8.54e-01 - 2.04e+01 -
\((1/2)^3\) 5.01e-01 0.76 3.81e-01 1.16 9.67e-01 1.08
\((1/2)^4\) 1.94e-01 1.37 1.45e-01 1.39 3.66e-01 1.40
\((1/2)^5\) 5.34e-02 1.86 2.44e-02 2.58 9.08e-02 2.01
\((1/2)^6\) 1.84e-02 1.54 4.38e-03 2.48 1.86e-02 2.29
\((1/2)^7\) 7.49e-03 1.30 9.10e-04 2.27 4.04e-03 2.20
\((1/2)^8\) 3.39e-03 1.14 2.07e-04 2.13 9.31e-04 2.11
\((1/2)^9\) 1.61e-03 1.07 4.95e-05 2.07 2.23e-04 2.06

6 Concluding Remarks↩︎

We analyzed the Robin-Robin coupling methods [14] for parabolic\parabolic interface problems and proved higher convergence rates in time for the first-order and second-order discrete time derivatives. We also prove \(H^2\) estimates of the discrete time derivatives in a special case. All the estimates in this work are key ingredients in proving that a prediction correction method [31] produces a \(O((\Delta t)^2)\) convergence rate.


This material is based upon work supported by the National Science Foundation under Grant No. DMS-1929284 while the last author was in residence at the Institute for Computational and Experimental Research in Mathematics in Providence, RI, during the "Numerical PDEs: Analysis, Algorithms, and Data Challenges" program. EB was supported by EPSRC grants EP/J002313/2 and EP/W007460/1. EB and MF were supported by the IMFIBIO Inria-UCL associated team.

7 Sketch of Proof: Corollary 6↩︎

Proof of 36 in Corollary 6. To prove 36 , it is suffice to bound the term \(\Xi^N(h_1, h_2, \alpha g_1-g_2,\tilde{g}_2)\). Note that since \(g_2^{n+1}=\tilde{g}_2^{n+1}\) on \(\Sigma\), we have \[\begin{align} {1} \Xi^N(h_1, h_2, \alpha g_1-g_2,\tilde{g}_2):= &\Delta t\sum_{n=0}^{N-1} (\frac{1}{\nu_s}\| h_1^{n+1}\|_{L^2(\Omega_s)}^2+(\frac{1}{\nu_f}+1)\| h_2^{n+1}\|_{L^2(\Omega_f)}^2) \\ &+ \Delta t\sum_{n=0}^{N-1} \Big( \frac{1}{\nu_f\alpha^2} \| \partial_{\Delta t}\tilde{g}_2^{n+1}\|_{L^2(\Omega_f)}^2 + \frac{\nu_f}{\alpha^2} \| \nabla \tilde{g}_2^{n+1}\|_{L^2(\Omega_f)}^2 + \frac{1}{\alpha}\| \tilde{g}_2^{n+1} \|_{L^2(\Omega_f)}^2 \Big) \\ &+ \Delta t\sum_{n=0}^{N-1} \Big(\frac{1}{\nu_s} \| \alpha g_1\|_{L^2(\Sigma)}^2+\frac{1}{\nu_f}\| \tilde{g}_2^{n+1}\|_{L^2(\Sigma)}^2 \Big) + \frac{1}{\alpha^2}\|\tilde{g}_2^N\|_{L^2(\Omega_f)}^2 \\ &=T_1+T_2+\ldots+T_8. \end{align}\] All the terms in \(\Xi^N(h_1, h_2, \alpha g_1-g_2,\tilde{g}_2)\) can be easily estimated by 44 except the \(T_3\) and \(T_4\). For \(T_3\), it follows from 48 that, \[\begin{align} \Delta t\sum_{n=0}^{N-1} \frac{1}{\nu_f\alpha^2} \| \partial_{\Delta t}\tilde{g}_2^{n+1}\|_{L^2(\Omega_f)}^2&=(\Delta t)^3 \sum_{n=0}^{N-1} \frac{1}{\nu_f\alpha^2} \| \partial_{\Delta t}^2 \tilde{\mathsf{l}}^{n+1}\|_{L^2(\Omega_f)}^2\\ &\le C\frac{(\Delta t)^2}{\nu_f\alpha^2}\|\partial_{t}^2 \tilde{\mathsf{l}}\|_{L^2(0,T;L^2(\Omega_f))}^2\\ &\le C(\Delta t)^2\frac{\nu_f}{\alpha^2}\|\partial_{t}^2 \mathsf{u}\|_{L^2(0,T;H^1(\Omega_f))}^2. \end{align}\] The term \(T_4\) can be bounded as follows, \[\Delta t\frac{\nu_f}{\alpha^2} \sum_{n=0}^{N-1}\|\nabla \tilde{g}_2^{n+1}\|_{L^2(\Omega_f)}^2 \le C \Delta t\frac{(\nu_f)^3}{\alpha^2} \sum_{n=0}^{N-1} \Big( \|\nabla(\mathsf{u}^{n+1}-u^n)\|_{L^2(\Omega_f)}^2 + \|D^2 (\mathsf{u}^{n+1}-\mathsf{u}^n)\|_{L^2(\Omega_f)}^2 \Big).\] Here \(D^2 \mathsf{u}\) denotes the Hessian of \(u\). Using 43 we get \[\begin{align} {1} \Delta t\frac{\nu_f}{\alpha^2} \sum_{n=0}^{N-1}\|\nabla \tilde{g}_2^{n+1}\|_{L^2(\Omega_f)}^2 \le & C (\Delta t)^2\frac{(\nu_f)^3}{\alpha^2} \|\partial_{t}\mathsf{u}\|_{L^2(0,T;H^2(\Omega_f))}^2. \end{align}\] The estimates above imply 36 . ◻

Before we prove the estimate 37 , we need the following preliminary results. We use the Bochner norms \(\|v\|_{L^2(a,b; X)}=\Big(\int_a^b \|v(\cdot, s)\|_X^2 ds\Big)^{1/2}\) and \(\|v\|_{L^\infty(a,b; X)}=\text{ess sup}_{ a \le s \le b} \|v(\cdot, s)\|_X\). For a Sobolev space \(X\), it is well known that \[\begin{alignat}{1}\tag{43} \| v^{n+1}-v^n\|_X^2 \le& C \Delta t\int_{t_n}^{t_{n+1}} \|\partial_{t}v(\cdot,s)\|_{X}^2 ds, \\ \| \partial_{\Delta t}v^{n+1}-\partial_{t}v^{n+1}\|_X^2 \le& C \Delta t\int_{t_n}^{t_{n+1}} \|\partial_{t}^2 v(\cdot,s)\|_{X}^2 ds, \tag{44}\\ \int_{a}^{b} \|v(\cdot,s)\|_X^2 ds \leq & (b-a) \|v\|_{L^\infty(a,b;X)}^2.\tag{45} \end{alignat}\] The following identities can easily be shown \[\label{eq:secordfd} \partial_{\Delta t}^2 v^n= \frac{1}{(\Delta t)^2} \int_{-\Delta t}^{\Delta t} (\Delta t-|s|)\partial_{t}^2 v(\cdot, t_n+s) ds,\tag{46}\] and \[\label{eq:cg23} \begin{align} \frac{\partial_{\Delta t}^2 v^{n}- \partial_{\Delta t}^2 v^{n-1}}{\Delta t} =& \frac{1}{(\Delta t)^3} \int_{-\Delta t}^{\Delta t} (\Delta t-|s|)\big(\partial_{t}^2 v(\cdot, t_n+s)- \partial_{t}^2 v(\cdot, t_{n-1}+s)\big) ds\\ =& \frac{1}{(\Delta t)^3} \int_{-\Delta t}^{\Delta t} (\Delta t-|s|)\int_{t_{n-1}+s}^{t_{n}+s} \partial_{t}^3 v(\cdot, r) dr ds. \end{align}\tag{47}\] From these we can show that \[\label{621} \| \partial_{\Delta t}^2 v^n\|_{X}^2 \le \frac{C}{\Delta t} \int_{t_{n-1}}^{t_{n+1}} \|\partial_{t}^2 v(\cdot,s)\|_{X}^2 ds,\tag{48}\] and \[\begin{align} {1} \|\frac{\partial_{\Delta t}^2 v^{n}- \partial_{\Delta t}^2 v^{n-1}}{\Delta t}\|_X^2 \le \frac{C}{\Delta t} \int_{t_{n-2}}^{t_{n+1}} \|\partial_{t}^3 v(\cdot, r)\|_X^2 dr. \label{v3d} \end{align}\tag{49}\]

Proof of 37 in Corollary 6. It is similar to prove 37 . Indeed, let us define, for \(j=1,2\) \[G_j^{n+1}=\partial_{\Delta t}g_j^{n+1},\quad\tilde{G}_2^{n+1}= \partial_{\Delta t}\tilde{g}_2^{n+1}, \quad H_j^{n+1}=\partial_{\Delta t}h_j^{n+1}.\] We then need to bound the term \[\begin{align} {1} \Xi^N&(H_1, H_2, \alpha G_1-G_2, \tilde{G}_2)\\ =& \Delta t\sum_{n=0}^{N-1} (\frac{1}{\nu_s}\| H_1^{n+1}\|_{L^2(\Omega_s)}^2+(\frac{1}{\nu_f}+1)\| H_2^{n+1}\|_{L^2(\Omega_f)}^2) \\ &+ \Delta t\sum_{n=0}^{N-1} \Big( \frac{1}{\nu_f\alpha^2} \| \partial_{\Delta t}\tilde{G}_2^{n+1}\|_{L^2(\Omega_f)}^2 + \frac{\nu_f}{\alpha^2} \| \nabla \tilde{G}_2^{n+1}\|_{L^2(\Omega_f)}^2 + \frac{1}{\alpha}\| \tilde{G}_2^{n+1} \|_{L^2(\Omega_f)}^2 \Big) \\ &+ \Delta t\sum_{n=0}^{N-1} \Big(\frac{1}{\nu_s} \| \alpha G_1\|_{L^2(\Sigma)}^2+\frac{1}{\nu_f}\| \tilde{G}_2^{n+1}\|_{L^2(\Sigma)}^2 \Big) +\frac{1}{\alpha^2}\|\tilde{G}_2^N\|_{L^2(\Omega_f)}^2 \\ =&R_1+R_2+\ldots+R_8. \end{align}\] For \(R_1\), it follows from the definitions of \(H_1^{n+1}\) and \(h_1^{n+1}\) that, \[\label{eq:h1esti} \begin{align} \|H_1^{n+1}\|^2_{L^2(\Omega_s)}&=\int_{\Omega_s}\left(\frac{h_1^{n+1}-h_1^n}{\Delta t}\right)^2\;\!dx\\ &=\int_{\Omega_s}\left(\frac{\partial_{t}\mathsf{w}^{n+1} - \partial_{\Delta t}\mathsf{w}^{n+1}-\partial_{t}\mathsf{w}^{n} + \partial_{\Delta t}\mathsf{w}^{n}}{\Delta t}\right)^2\;\!dx\\ &=\int_{\Omega_s}\left(\partial_{\Delta t}(\partial_t\mathsf{w})^{n+1}-\partial_{\Delta t}^2\mathsf{w}^{n}\right)^2\;\!dx\\ &=\int_{\Omega_s}\left(\partial_{\Delta t}(\partial_t\mathsf{w})^{n+1}-\partial_{t}^2\mathsf{w}^n+\partial_{t}^2\mathsf{w}^n-\partial_{\Delta t}^2\mathsf{w}^{n+1}\right)^2\;\!dx\\ &\le C\int_{\Omega_s}\left(\frac{1}{\Delta t}\int_{t_n}^{t_{n+1}}(\partial_{t}^2\mathsf{w}(t_n)-\partial_{t}^2\mathsf{w}(s))\;\!ds\right)^2\;\!dx\\ &\quad +C \int_{\Omega_s}\left(\frac{1}{(\Delta t)^2}\int_{t_{n-1}}^{t_{n+1}}\big(\Delta t-|s-t_n|\big)\big(\partial_{t}^2\mathsf{w}(s)-\partial_{t}^2\mathsf{w}(t_n)\big)\;\!ds\right)^2\;\!dx\\ &\le C\Delta t\|\partial_{t}^3\mathsf{w}\|^2_{L^2((t_{n-1},t_{n+1}),L^2(\Omega_s))}. \end{align}\tag{50}\] Therefore, we obtain \[\Delta t\sum_{n=0}^{N-1}\frac{1}{\nu_s}\|H_1^{n+1}\|^2_{L^2(\Omega_s)}\le C(\Delta t)^2\frac{1}{\nu_s}\|\partial_{t}^3\mathsf{w}\|^2_{L^2((0,T),L^2(\Omega_s))}\] The estimate of \(R_2\) is similar to that of \(R_1\). The term \(R_3\) can be estimated as follow by 49 , \[\begin{align} \Delta t\sum_{n=0}^{N-1} \frac{1}{\nu_f\alpha^2} \| \partial_{\Delta t}\tilde{G}_2^{n+1}\|_{L^2(\Omega_f)}^2&=(\Delta t)^3\sum_{n=2}^{N-1} \frac{1}{\nu_f\alpha^2}\| \frac{\partial_{\Delta t}^2 \tilde{\mathsf{l}}^{n}- \partial_{\Delta t}^2 \tilde{\mathsf{l}}^{n-1}}{\Delta t}\|_{L^2(\Omega_f)}^2\\ &\le C\frac{(\Delta t)^2}{\nu_f\alpha^2}\|\partial_{t}^3 \tilde{\mathsf{l}}\|_{L^2(0,T;L^2(\Omega_f))}^2\\ &\le C(\Delta t)^2\frac{\nu_f}{\alpha^2}\|\partial_{t}^3 \mathsf{u}\|_{L^2(0,T;H^1(\Omega_f))}^2. \end{align}\] For \(R_4\), it follows from the definition of \(\tilde{\mathsf{l}}\) and the fact \(\|\nabla\phi\|_{L^2(\Omega_f)}\le C\) that, \[\label{eq:ng2} \begin{align} \|\nabla\tilde{G}_2^{n+1}\|_{L^2(\Omega_f)}^2&=\int_{\Omega_f}\left(\frac{\nabla\tilde{\mathsf{l}}^{n+1}-2\nabla\tilde{\mathsf{l}}^{n}+\nabla\tilde{\mathsf{l}}^{n-1}}{\Delta t}\right)^2\;\!dx\\ &=\int_{\Omega_f}\left(\frac{1}{\Delta t}\int_{t_{n-1}}^{t_{n+1}}\big(\Delta t-|s-t_n|\big)\partial_{t}^2\nabla\tilde{\mathsf{l}}(s)\;\!ds\right)^2\;\!dx\\ &\le C\int_{\Omega_f}\left(\int_{t_{n-1}}^{t_{n+1}}|\partial_{t}^2\nabla\tilde{\mathsf{l}}(s)|\;\!ds\right)^2\;\!dx\\ &\le C\Delta t\int_{\Omega_f}\int_{t_{n-1}}^{t_{n+1}}(\partial_{t}^2\nabla\tilde{\mathsf{l}}(s))^2\;\!ds\;\!dx\\ &\le C \nu_f^2\Delta t\|\partial_{t}^2\mathsf{u}\|_{L^2(t_{n-1},t_{n+1};H^2(\Omega_f))}^2, \end{align}\tag{51}\] and hence, \[\label{eq:R4} \Delta t\sum_{n=0}^{N-1}\frac{\nu_f}{\alpha^2}\| \nabla \tilde{G}_2^{n+1}\|_{L^2(\Omega_f)}^2\le C(\Delta t)^2\frac{(\nu_f^3)}{\alpha^2}\|\partial_{t}^2\mathsf{u}\|_{L^2(0,T;H^2(\Omega_f))}^2.\tag{52}\] The remaining terms \(R_5\) to \(R_8\) can be estimated similarly, we give the estimate of \(R_6\) here: \[\label{eq:g1s} \begin{align} \|G_1^{n+1}\|^2_{L^2(\Sigma)}&=\int_{\Sigma}\left(\frac{g_1^{n+1}-g_1^n}{\Delta t}\right)^2\;\!d\sigma\\ &=\int_{\Sigma}\left(\frac{\mathsf{u}^{n+1}-2\mathsf{u}^{n}+\mathsf{u}^{n-1}}{\Delta t}\right)^2\;\!d\sigma\\ &=\int_{\Sigma}\left(\frac{1}{\Delta t}\int_{t_{n-1}}^{t_{n+1}}\big(\Delta t-|s-t_n|\big)\partial_{t}^2\mathsf{u}(s)\;\!ds\right)^2\;\!d\sigma\\ &\le C\Delta t\|\partial_{t}^2\mathsf{u}\|^2_{L^2((t_{n-1},t_{n+1}),L^2(\Sigma))}, \end{align}\tag{53}\] and therefore, we obtain \[\label{eq:R6} \Delta t\sum_{n=0}^{N-1} \frac{1}{\nu_s} \| \alpha G_1\|_{L^2(\Sigma)}^2\le C\frac{\alpha^2}{\nu_s}\|\partial_{t}^2\mathsf{u}\|^2_{L^2((0,T),L^2(\Sigma))}.\tag{54}\] We finish the proof by combining all the estimates above. ◻

In order to prove 38 , we define the following quantities for \(j=1,2\), \[\label{eq:curlygh} \mathcal{G}_j^{n+1}=\partial_{\Delta t}G_j^{n+1},\quad \mathcal{H}_j^{n+1}=\partial_{\Delta t}H_j^{n+1}, \tilde{\mathcal{G}}_2^{n+1}= \partial_{\Delta t}\tilde{G}_2^{n+1}.\tag{55}\] Note that we also have the following \[\begin{align}\label{eq:cg21dd} \tilde{\mathcal{G}}_2^{n+1}=\frac{\tilde{g}_2^{n+1}-2\tilde{g}_2^{n}+\tilde{g}_2^{n-1}}{(\Delta t)^2}=\frac{\tilde{\mathsf{l}}^{n+1}-3\tilde{\mathsf{l}}^{n}+3\tilde{\mathsf{l}}^{n-1}-\tilde{\mathsf{l}}^{n-2}}{(\Delta t)^2}. \end{align}\tag{56}\] Denote \[\label{eq:dfq3} \partial_{\Delta t}^3v^{n+1}=\frac{v^{n+1}-3v^n+3v^{n-1}-v^{n-2}}{(\Delta t)^3},\tag{57}\] thus we see that \[\begin{align}\label{eq:cg22dd} \frac{\tilde{\mathcal{G}}_2^{n+1}-\tilde{\mathcal{G}}_2^{n}}{(\Delta t)^2}=\frac{\partial_{\Delta t}^3\tilde{\mathsf{l}}^{n+1}-\partial_{\Delta t}^3\tilde{\mathsf{l}}^{n}}{\Delta t}. \end{align}\tag{58}\] Moreover, we see that \[\begin{align}\label{eq:cg23dd} \frac{\partial_{\Delta t}^3 v^{n}- \partial_{\Delta t}^3 v^{n-1}}{\Delta t} =& \frac{1}{(\Delta t)^4} \Bigg(\int_{-\Delta t}^{\Delta t} (\Delta t-|s|)\big(\partial_{t}^2 v(\cdot, t_n+s)- \partial_{t}^2 v(\cdot, t_{n-1}+s)\big) ds\\ &-\int_{-\Delta t}^{\Delta t} (\Delta t-|s|)\big(\partial_{t}^2 v(\cdot, t_{n-1}+s)- \partial_{t}^2 v(\cdot, t_{n-2}+s)\big) ds\Bigg)\\ =&\frac{1}{(\Delta t)^4} \Bigg(\int_{-\Delta t}^{\Delta t} (\Delta t-|s|)\big(\partial_{t}^2 v(\cdot, t_n+s)- 2\partial_{t}^2 v(\cdot, t_{n-1}+s)+\partial_{t}^2 v(\cdot, t_{n-2}+s)\big) ds\Bigg)\\ =& \frac{1}{(\Delta t)^4} \int_{-\Delta t}^{\Delta t} (\Delta t-|s|)\int_{t_{n-2}+s}^{t_{n}+s} (\Delta t-|r-t_{n-1}|)\partial_{t}^4 v(\cdot, r) dr ds. \end{align}\tag{59}\]

Therefore, we obtain the following estimate \[\begin{align} {1} \|\frac{\partial_{\Delta t}^3 v^{n}- \partial_{\Delta t}^3 v^{n-1}}{\Delta t}\|_X^2 \le \frac{C}{\Delta t} \int_{t_{n-3}}^{t_{n+1}} \|\partial_{t}^4 v(\cdot, r)\|_X^2 dr. \label{v3dd} \end{align}\tag{60}\]

Proof of 38 in Corollary 6. We now bound the term \(\Xi^N(\partial_{\Delta t}^2 h_1, \partial_{\Delta t}^2 h_2, \partial_{\Delta t}^2 (\alpha g_1-g_2), \partial_{\Delta t}^2 \tilde{g}_2)\) which is \(\Xi^N(\mathcal{H}_1, \mathcal{H}_2, (\alpha \mathcal{G}_1-\mathcal{G}_2), \tilde{\mathcal{G}}_2)\) according to 55 . The techniques are similar to those mentioned above. We present the key estimates involving \(\mathcal{H}_1\) and \(\tilde{\mathcal{G}}_2\) first. It follows from the definition of \(\mathcal{H}_1^{n+1}\) that, \[\label{eq:h1estid} \begin{align} \|\mathcal{H}_1^{n+1}\|^2_{L^2(\Omega_s)}&=\int_{\Omega_s}\left(\frac{h_1^{n+1}-2h_1^n+h_1^{n-1}}{(\Delta t)^2}\right)^2\;\!dx\\ &=\int_{\Omega_s}\left(\frac{\partial_{\Delta t}(\partial_{t}\mathsf{w})^{n+1}-\partial_{\Delta t}^2\mathsf{w}^{n+1}-(\partial_{\Delta t}(\partial_{t}\mathsf{w})^{n}-\partial_{\Delta t}^2\mathsf{w}^{n})}{\Delta t}\right)^2\;\!dx\\ &=\int_{\Omega_s}\left(\partial_{\Delta t}^2(\partial_{t}\mathsf{w})^{n+1}-\partial_{\Delta t}^3\mathsf{w}^{n+1}\right)^2\;\!dx\\ &=\int_{\Omega_s}\left(\partial_{\Delta t}^2(\partial_{t}\mathsf{w})^{n+1}-\partial_{t}^3\mathsf{w}^n+\partial_{t}^3\mathsf{w}^n-\partial_{\Delta t}^3\mathsf{w}^{n+1}\right)^2\;\!dx\\ &\le C\int_{\Omega_s}\left(\frac{1}{(\Delta t)^2}\int_{t_{n-1}}^{t_{n+1}}\big(\Delta t-|s-t_n|\big)\big(\partial_{t}^3\mathsf{w}(s)-\partial_{t}^3\mathsf{w}(t_n)\big)\;\!ds\right)^2\;\!dx\\ &\quad+ C\int_{\Omega_s}\left(\frac{1}{(\Delta t)^3}\int_{-\Delta t}^{\Delta t}(\Delta t-|s|)\int_{t_{n-1}+s}^{t_n+s}\partial_{t}^3\mathsf{w}(t_n)-\partial_{t}^3\mathsf{w}(\cdot,r)\;\!ds\right)^2\;\!dx\\ &\le C\Delta t\|\partial_{t}^4\mathsf{w}\|^2_{L^2((t_{n-2},t_{n+1}),L^2(\Omega_s))}. \end{align}\tag{61}\] We also have, according to the definition of \(\tilde{\mathsf{l}}\) and the fact \(\|\nabla\phi\|_{L^2(\Omega_f)}\le C\) that, \[\begin{align} \|\nabla\tilde{\mathcal{G}}_2^{n+1}\|_{L^2(\Omega_f)}^2&=\int_{\Omega_f}\left(\frac{\frac{\nabla\tilde{\mathsf{l}}^{n+1}-2\nabla\tilde{\mathsf{l}}^{n}+\nabla\tilde{\mathsf{l}}^{n-1}}{\Delta t}-\frac{\nabla\tilde{\mathsf{l}}^{n}-2\nabla\tilde{\mathsf{l}}^{n-1}+\nabla\tilde{\mathsf{l}}^{n-2}}{\Delta t}}{\Delta t}\right)^2\;\!dx\\ &=\int_{\Omega_f}\left(\frac{1}{(\Delta t)^2}\int_{-\Delta t}^{\Delta t} (\Delta t-|s|)\big(\partial_{t}^2 \nabla\tilde{\mathsf{l}}(\cdot, t_n+s)- \partial_{t}^2 \nabla\tilde{\mathsf{l}}(\cdot, t_{n-1}+s)\big) ds\right)^2\;\!dx\\ &\le C\int_{\Omega_f}\left(\frac{1}{\Delta t}\int_{-\Delta t}^{\Delta t}|\int_{t_{n-1}+s}^{t_{n}+s} \partial_{t}^3 \nabla\tilde{\mathsf{l}}(\cdot, r) dr |\;\!ds\right)^2\;\!dx\\ &\le C\Delta t\int_{\Omega_f}\int_{t_{n-2}}^{t_{n+1}}(\partial_{t}^3\nabla\tilde{\mathsf{l}}(\cdot,r))^2\;\!dr\;\!dx\\ &\le C \nu_f^2\Delta t\|\partial_{t}^3\mathsf{u}\|_{L^2(t_{n-2},t_{n+1};H^2(\Omega_f))}^2. \end{align}\] It follows from 60 that, \[\begin{align} \Delta t\sum_{n=0}^{N-1} \frac{1}{\nu_f\alpha^2}\| \partial_{\Delta t}\tilde{\mathcal{G}_2}^{n+1}\|_{L^2(\Omega_f)}^2&=(\Delta t)^3 \sum_{n=0}^{N-1} \frac{1}{\nu_f\alpha^2} \|\frac{\partial_{\Delta t}^3 \tilde{\mathsf{l}}^{n}- \partial_{\Delta t}^3 \tilde{\mathsf{l}}^{n-1}}{\Delta t}\|_{L^2(\Omega_f)}^2\\ &\le C\frac{(\Delta t)^2}{\nu_f\alpha^2}\|\partial_{t}^4 \tilde{\mathsf{l}}\|_{L^2(0,T;L^2(\Omega_f))}^2\\ &\le C(\Delta t)^2\frac{\nu_f}{\alpha^2}\|\partial_{t}^4 \mathsf{u}\|_{L^2(0,T;H^1(\Omega_f))}^2. \end{align}\] Then we follow the same idea as before and obtain the following bound: \[\Xi^N(\mathcal{H}_1, \mathcal{H}_2, (\alpha \mathcal{G}_1-\mathcal{G}_2), \tilde{\mathcal{G}}_2)\le C(\Delta t)^2\mathfrak{Y},\] where \(\mathfrak{Y}\) is defined in 40 . ◻

8 Sketch of Proof: Corollary 9↩︎

Proof. According to Corollary 8, we need to bound the terms on the right-hand side of the inequality 41 . The first term in 41 is bounded in Corollary 6 while the last term in 41 can be easily bounded according to 50 . At last, the analysis to bound the term \(\Xi^N(\partial_{\Delta t}{\partial}_x h_1, \partial_{\Delta t}{\partial}_x h_2, \partial_{\Delta t}{\partial}_x (\alpha g_1-g_2), \partial_{\Delta t}{\partial}_x \tilde{g}_2)\) is almost identical to that of the term \(\Xi^N(H_1, H_2, \alpha G_1-G_2, \tilde{G}_2)\) except the additional partial derivative which does not affect the techniques. Therefore, we obtain the bound \[\Xi^N(\partial_{\Delta t}{\partial}_x h_1, \partial_{\Delta t}{\partial}_x h_2, \partial_{\Delta t}{\partial}_x (\alpha g_1-g_2), \partial_{\Delta t}{\partial}_x \tilde{g}_2)\le C(\Delta t)^2\mathbb{Y},\] where \(\mathbb{Y}\) is defined in 42 . This finishes the proof. ◻


S. Badia, A. Quaini, and A. Quarteroni. Splitting methods based on algebraic factorization for fluid-structure interaction. SIAM Journal on Scientific Computing, 30(4):1778–1805, 2008.
E. Burman and M. A. Fernández. Stabilization of explicit coupling in fluid–structure interaction involving fluid incompressibility. Computer Methods in Applied Mechanics and Engineering, 198(5):766–784, 2009.
J. W. Banks, W. D. Henshaw, and D. W. Schwendeman. An analysis of a new stable partitioned algorithm for FSI problems. Part I: Incompressible flow and elastic solids. J. Comput. Phys., 269:108–137, 2014.
E. Burman and M. A. Fernández. Explicit strategies for incompressible fluid-structure interaction problems: Nitsche type mortaring versus RobinRobin coupling. International Journal for Numerical Methods in Engineering, 97(10):739–758, 2014.
M. Bukač, S. Čanić, R. Glowinski, B. Muha, and A. Quaini. A modular, operator-splitting scheme for fluid-structure interaction problems with thick structures. Internat. J. Numer. Methods Fluids, 74(8):577–604, 2014.
M. Bukač. An extension of explicit coupling for fluid–structure interaction problems. Mathematics, 9(15), 2021.
G. Gigante and C. Vergara. On the choice of interface parameters in Robin–Robin loosely coupled schemes for fluid–structure interaction. Fluids, 6(6), 2021.
G. Gigante and C. Vergara. On the stability of a loosely-coupled scheme based on a Robin interface condition for fluid-structure interaction. Comput. Math. Appl., 96:109–119, 2021.
D. A. Serino, J. W. Banks, W. D. Henshaw, and D. W. Schwendeman. A stable added-mass partitioned (amp) algorithm for elastic solids and incompressible flow: Model problem analysis. SIAM Journal on Scientific Computing, 41(4):A2464–A2484, 2019.
P. Minev and R. Usubov. Splitting schemes for the stress formulation of fluid–structure interaction problems. Applications in Engineering Science, 9:100082, 2022.
M. Bucelli, M. Geraint Gabriel, A. Quarteroni, G. Gigante, and C. Vergara. A stable loosely-coupled scheme for cardiac electro-fluid-structure interaction. Journal of Computational Physics, 490:112326, 2023.
E. Burman, R. Durst, and J. Guzmán. Stability and error analysis of a splitting method using Robin–Robin coupling applied to a fluid–structure interaction problem. Numerical Methods for Partial Differential Equations, 38(5):1396–1406, 2022.
E. Burman, R. Durst, M. A. Fernández, and J. Guzmán. Fully discrete loosely coupled Robin-Robin scheme for incompressible fluid–structure interaction: stability and error analysis. Numerische Mathematik, 151(4):807–840, 2022.
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.
A. Seboldt and M. Bukač. A non-iterative domain decomposition method for the interaction between a fluid and a thick structure. Numerical Methods for Partial Differential Equations, 37(4):2803–2832, 2021.
R. Durst. Recent Advances in Splitting Methods Based on Robin-Robin Coupling Conditions. PhD thesis, Brown University, 2022.
E. Burman, R. Durst, M. A. Fernández, J. Guzmán, and O. Ruz. obin-Robin loose coupling for incompressible fluid-structure interaction: non-linear setting and nearly-optimal error analysis. 2023. .
W. Layton and A. Takhirov. Energy stability of a first order partitioned method for systems with general coupling. Int. J. Numer. Anal. Model. Ser. B, 4(3):203–214, 2013.
J. M. Connors, J. S. Howell, and W. J. Layton. Partitioned time stepping for a parabolic two domain problem. SIAM Journal on Numerical Analysis, 47(5):3526–3549, 2009.
J. M. Connors, J. S. Howell, and W. J. Layton. Decoupled time stepping methods for fluid-fluid interaction. SIAM J. Numer. Anal., 50(3):1297–1319, 2012.
J. M. Connors and J. S. Howell. A fluid-fluid interaction method using decoupled subproblems and differing time steps. Numer. Methods Partial Differential Equations, 28(4):1283–1308, 2012.
Y. Zhang, Y. Hou, and L. Shan. Error estimates of a decoupled algorithm for a fluid-fluid interaction problem. J. Comput. Appl. Math., 333:266–291, 2018.
H. Zhang, Z. Liu, E. Constantinescu, and R. Jacob. Stability analysis of interface conditions for ocean-atmosphere coupling. J. Sci. Comput., 84(3):Paper No. 44, 25, 2020.
Y. Zhang, L. Shan, and Y. Hou. New approach to prove the stability of a decoupled algorithm for a fluid-fluid interaction problem. J. Comput. Appl. Math., 371:112695, 19, 2020.
K. C. Sockwell, P. Bochev, K. Peterson, and P. Kuberry. Interface flux recovery framework for constructing partitioned heterogeneous time-integration methods. Numer. Methods Partial Differential Equations, 39(5):3572–3593, 2023.
W. Li, P. Huang, and Y. He. An unconditionally energy stable finite element scheme for a nonlinear fluid-fluid interaction model. IMA J. Numer. Anal., 44(1):157–191, 2024.
M. A. Fernández, J.-F. Gerbeau, and S. Smaldone. Explicit coupling schemes for a fluid-fluid interaction problem arising in hemodynamics. SIAM J. Sci. Comput., 36(6):A2557–A2583, 2014.
M. Beneš, A. Nekvinda, and M. K. Yadav. Multi-time-step domain decomposition method with non-matching grids for parabolic problems. Appl. Math. Comput., 267:571–582, 2015.
M. Beneš. Convergence and stability analysis of heterogeneous time step coupling schemes for parabolic problems. Appl. Numer. Math., 121:198–222, 2017.
C. Canuto and A. Lo Giudice. A multi-timestep Robin-Robin domain decomposition method for time dependent advection-diffusion problems. Appl. Math. Comput., 363:124596, 14, 2019.
E. Burman, R. Durst, M. Fernández, J. Guzmán, and S. Liu. A second-order correction method for loosely coupled discretizations applied to parabolic-parabolic interface problems. Preprint, 2024.
K. Böhmer, P. W. Hemker, and H. J. Stetter. The defect correction approach. Defect Correction Methods: Theory and Applications, pages 1–32, 1984.
M. Alnæs, J. Blechta, J. Hake, A. Johansson, B. Kehlet, A. Logg, C. Richardson, J. Ring, M. E. Rognes, and G. N. Wells. The fenics project version 1.5. Archive of numerical software, 3(100), 2015.
F. Ballarin and G. Rozza., 2016.