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\).

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.

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.

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

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}\]

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.

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.

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.

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.

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.

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.

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.

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

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.

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.

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

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.

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.

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.

Acknowledgement↩︎

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↩︎

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}\]

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}\]

8 Sketch of Proof: Corollary 9↩︎

References↩︎