remarkRemark hypothesisHypothesis claimClaim Dirichlet-Neumann Learning
Algorithm for Interface ProblemsQ. Sun, X. Xu, and H. Yi
[][nocite]ex_supplement

# Dirichlet-Neumann learning algorithm for solving elliptic interface
problems††thanks: Submitted to the editors DATE. Qi Sun is supported in part
by the National Natural Science Foundation of China under Grant 12201465.

Qi Sun School of Mathematical Sciences, Tongji University, Shanghai 200092,
China<EMAIL_ADDRESS><EMAIL_ADDRESS>Xuejun
Xu22footnotemark: 2$\,\,{}^{,}$ Institute of Computational Mathematics, AMSS,
Chinese Academy of Sciences, Beijing 100190, China<EMAIL_ADDRESS>Haotian Yi22footnotemark: 2

###### Abstract

Non-overlapping domain decomposition methods are natural for solving interface
problems arising from various disciplines, however, the numerical simulation
requires technical analysis and is often available only with the use of high-
quality grids, thereby impeding their use in more complicated situations. To
remove the burden of mesh generation and to effectively tackle with the
interface jump conditions, a novel mesh-free scheme, i.e., Dirichlet-Neumann
learning algorithm, is proposed in this work to solve the benchmark elliptic
interface problem with high-contrast coefficients as well as irregular
interfaces. By resorting to the variational principle, we carry out a rigorous
error analysis to evaluate the discrepancy caused by the boundary penalty
treatment for each decomposed subproblem, which paves the way for realizing
the Dirichlet-Neumann algorithm using neural network extension operators. The
effectiveness and robustness of our proposed methods are demonstrated
experimentally through a series of elliptic interface problems, achieving
better performance over other alternatives especially in the presence of
erroneous flux prediction at interface.

###### keywords:

Elliptic interface problem, discontinuous coefficients, compensated deep Ritz
method, artificial neural networks

65M55, 65Nxx, 92B20, 49S05

##  1 Introduction

Many problems in science and engineering are carried out with domains
separated by curves or surfaces, e.g., the abrupt changes in material
properties between adjacent regions, from which the interface problems
naturally arise. A widely studied benchmark example is the elliptic interface
problem with high-contrast coefficients [39, 6, 35], whose solution lies in
the Sobolev space $H^{1+\epsilon}(\Omega)$ with $\epsilon>0$ possibly close to
zero [41]. Due to the low regularity of solution at the interface, classical
numerical methods, such as finite difference and finite element methods [33,
4], require the generation of an interface-fitted mesh in the discretization
of the computation domain [6], which could be technically involved and time
consuming especially when the geometry of interface gets complicated. To ease
the burden of manual mesh construction, numerical methods based on unfitted
meshes, e.g., the immersed interface method [35] and many others, have emerged
as an attractive alternative for solving the elliptic interface problems [39].
However, using unfitted mesh, e.g, a uniform Cartesian mesh, often makes it
difficult to enforce the jump conditions across subdomain interfaces
accurately [5, 19]. Although both types of methods have shown to be effective
to some extent, their practical implementation is not an easy task due to the
complex geometry of interface and the discontinuity of solution.

Figure 1: Computational graph of our proposed Dirichlet-Neumann learning
algorithm.

Note that the entire computation domain has already been separated as the
union of multiple non-overlapping subregions, each of which corresponds to a
local boundary value problem after endowing the subregion interface with an
appropriate boundary condition [47]. As a result, the elliptic interface
problem being considered can also be tackled by a non-overlapping Dirichlet-
Neumann algorithm at the continuous level [60, 61, 47], where the decomposed
subproblems are typically solved using the mesh-based numerical methods [33,
4, 34]. However, the complex geometry of subdomain interfaces remains a major
concern during the meshing process. Fortunately, the domain decomposition
methods [54, 47, 9] are essentially continuous schemes that greatly differ
from the aforementioned interface-fitted or -unfitted methods, making it
computationally feasible to adopt the mesh-free deep learning techniques [17,
28, 62] as the local problem solvers [20]. Thanks to the rapid development of
artificial intelligence science, much attention has recently been paid to
combining deep learning with insight from the domain decomposition methods.
The physics-informed neural networks [48, 28, 31, 30], among others [62, 63,
51], has been utilized to discretize and solve the Dirichlet and Neumann
subproblems within the classical Dirichlet-Neumann algorithm framework [38],
which is named “DeepDDM” and applied to a simple interface problem as a proof
of concept. To further enhance its scalability properties, the DeepDDM method
is extended with the aid of coarse space correction [42, 47]. However, these
works fail to address the impact of erroneous flux prediction [8, 1] at
interface, which can degrade the overall performance [52]. Note that in the
degenerate case of homogeneous jump conditions, the continuity of averaged
solution between neighbouring subregions, as well as its first and higher-
order derivatives, are explicitly enforced through additional penalty terms in
a series of papers [25, 26, 50, 59, 23], which also suffer from the issue of
erroneous flux predictions. Worse still, these strategies can not be applied
to resolving discontinuous solutions, where the higher-order derivatives of
network solutions at interface could be meaningless due to the lack of global
regularity [41]. Designing specific network architectures is another way of
dealing with the complex geometry and jump condition of elliptic interface
problems [32, 21, 55, 22, 10], e.g., augmenting an additional coordinate
variable as a feature input of the solution ansatz [32], replacing neural
network structures with extreme learning machines [10, 12] or graph neural
networks [53], while the theoretical study is still in need to gain a deeper
understanding.

In contrast to the existing methods that require an explicit computation of
the flux data along subregion interfaces, a novel Dirichlet-Neumann learning
algorithm using neural network extension operators (see Figure 1) is proposed
in this work for solving the elliptic interface problems with high-contrast
coefficients as well as irregular interfaces. By resorting to the variational
principle and applying the Green’s formula to Dirichlet-Neumann algorithm [47,
60, 61], a rigorous error analysis is first established to estimate the
discrepancies caused by the penalty treatment111The boundary conditions are
included as “soft” penalty terms in the training loss function. of boundary
conditions [27, 11, 44], which paves the way for constructing the training
loss functions for each decomposed subproblem. As a direct result, the flux
transmission between neighbouring subdomains is realized without explicitly
computing the derivatives of network solution at interface (see also Figure 1)
and thereby alleviating the deterioration of outer iterations in the presence
of erroneous flux predictions. Moreover, a comparison study is presented to
demonstrate the robustness and effectiveness of our proposed methods over the
DeepDDM scheme [38], followed by a series of numerical examples to validate
our statements.

The rest of this paper is organized as follows. In section 2, we begin by
recalling the Dirichlet-Neumann algorithm for solving elliptic interface
problems in the continuous level, then the numerical methods, i.e., mesh-based
and mesh-free solvers, are briefly reviewed and compared. To combine domain
decomposition methods with ideas from machine learning in a consistent manner,
a rigorous error analysis of boundary penalty method for both the Dirichlet
and Neumann subproblems is presented in section 3, followed by the
implementation details of our proposed Dirichlet-Neumann learning algorithm.
Then, the numerical results on a series of benchmark problems are reported in
section 4. Finally, we summarize our work in section 5.

##  2 Preliminaries

###  2.1 Elliptic Interface Problem with High-Contrast Coefficients

Figure 2: A Lipschitz domain $\Omega\subset\mathbb{R}^{2}$ that is decomposed
into several non-overlapping subdregions.

Let $\Omega\subset\mathbb{R}^{d}$ be a bounded domain with Lipschitz boundary
$\partial\Omega$, which is first assumed to be partitioned into two non-
overlapping subregions as depicted in Figure 2, that is,

$\overline{\Omega}=\overline{\Omega_{1}\cap\Omega_{2}},\ \ \
\Omega_{1}\cap\Omega_{2}=\emptyset,\ \ \
\Gamma=\partial\Omega_{1}\cap\partial\Omega_{2}.$

Then, we consider the elliptic interface problem with high-contrast
coefficients and natural jump conditions, which is often formally written as
[39]

(1) $\begin{array}[]{cl}-\nabla\cdot\left(c(x)\nabla
u(x)\right)+u(x)=f(x)&\text{in}\ \Omega,\\\ u(x)=0&\text{on}\
\partial\Omega,\\\ \llbracket u(x)\rrbracket=0\ \ \text{and}\ \ \llbracket
c(x)\nabla u(x)\cdot\bm{n}\rrbracket=q(x)&\text{on}\ \Gamma,\end{array}$

where $f(x)$ represents a given function of $L^{2}(\Omega)$,
$\bm{n}=\bm{n}_{2}$ ($\bm{n}_{1}$) the unit outer normal vector for subdomain
$\Omega_{2}$ ($\Omega_{1}$) and notation $\llbracket\cdot\rrbracket$ the
difference of quantity across the interface, namely, for any point
$X\in\Gamma$,

(2) $\begin{array}[]{c}\lim\limits_{x\rightarrow
X,\,x\in\Omega_{1}}u(x)=\lim\limits_{x\rightarrow X,\,x\in\Omega_{2}}u(x),\\\
-\lim\limits_{x\rightarrow X,\,x\in\Omega_{1}}c(x)\nabla
u(x)\cdot\bm{n}_{1}-\lim\limits_{x\rightarrow X,\,x\in\Omega_{2}}c(x)\nabla
u(x)\cdot\bm{n}_{2}=q(X).\end{array}$

Notably, $c(x)$ is a piecewise constant function that has a finite jump of
function value across the interface $\Gamma$, that is,

$c(x)=\left\\{\begin{array}[]{cl}c_{1}>0&\text{in}\ \Omega_{1},\\\ c_{2}\gg
c_{1}&\text{in}\ \Omega_{2}.\end{array}\right.$

which is typically caused by the abrupt changes in material properties or the
interaction of fluid dynamics [39]. More broadly, the regular coefficients
$c_{1}\approx c_{2}>0$ or $c_{1}=c_{2}$ can be regarded as a degenerate case
of problem (1).

By setting $V_{i}=\big{\\{}v_{i}\in
H^{1}(\Omega_{i})\,\big{|}\,v_{i}|_{\partial\Omega\cap\partial\Omega_{i}}=0\big{\\}}$,
$V_{i}^{0}=H_{0}^{1}(\Omega_{i})$, and defining

$b_{i}(u_{i},v_{i})=\int_{\Omega_{i}}c_{i}\nabla u_{i}\cdot\nabla
v_{i}+u_{i}v_{i}\,dx,\ \ (f,v_{i})_{i}=\int_{\Omega_{i}}fv_{i}\,dx,\ \
(q,v_{2})_{L^{2}(\Gamma)}=\int_{\Gamma}qv\,ds,$

for $i=1$, 2, the Green’s formula implies that the weak formulation of (1)
reads: find $u_{1}\in V_{1}$ and $u_{2}\in V_{2}$ such that

(3) $\begin{array}[]{cl}b_{1}(u_{1},v_{1})=(f,v_{1})_{1}&\text{for any}\
v_{1}\in V_{1}^{0},\\\ u_{1}=u_{2}&\text{on}\ \Gamma,\\\
b_{2}(u_{2},v_{2})\\!=\\!(f,v_{2})_{2}\\!+\\!(f,R_{1}\gamma_{0}v_{2})_{1}\\!-\\!b_{1}(u_{1},R_{1}\gamma_{0}v_{2})\\!-\\!(q,v_{2})_{L^{2}(\Gamma)}&\text{for
any}\ v_{2}\in V_{2},\end{array}$

where $\gamma_{0}v=v|_{\Gamma}$ is the restriction of $v\in H^{1}(\Omega_{i})$
on $\Gamma$ and $R_{i}:H_{00}^{\frac{1}{2}}(\Gamma)\to V_{i}$ any
differentiable extension operator [47].

As a direct result, by employing a suitable relaxation parameter
$\rho\in(0,\rho_{\textnormal{max}})$ [45, 15], an iterative scheme (also known
as the Dirichlet-Neumann algorithm [60, 47]) can be developed for solving the
elliptic interface problem (3): given the initial guess of the unknown
solution value at interface $u_{\Gamma}^{[0]}\in
H_{00}^{\frac{1}{2}}(\Gamma)$, then solve for $k\geq 0$,

  * 1)

$\displaystyle u_{1}^{[k]}=\operatorname*{arg\,min}_{u_{1}\in
V_{1},\,u_{1}|_{\Gamma}=u_{\Gamma}^{[k]}}\frac{1}{2}b_{1}(u_{1},u_{1})-(f,u_{1})_{1}$
(4)

  * 2)

$\displaystyle u_{2}^{[k]}=\operatorname*{arg\,min}_{u_{2}\in
V_{2}}\frac{1}{2}b_{2}(u_{2},u_{2})-(f,u_{2})_{2}+b_{1}(u_{1}^{[k]},R_{1}\gamma_{0}u_{2})-(f,R_{1}\gamma_{0}u_{2})_{1}+(q,u_{2})_{L^{2}(\Gamma)}$
(5)

  * 3)

$\displaystyle u_{\Gamma}^{[k+1]}=\rho u_{2}^{[k]}+(1-\rho)u_{\Gamma}^{[k]}\ \
\ \text{on}\ \Gamma$ (6)

until certain stopping criteria are met [54]. It is noteworthy that in
contrast to the strong form (1), the jump conditions are enforced without
explicitly computing and exchanging the interface flux of local solutions
[52].

######  Remark 2.1.

Although the solution $u(x)$ of interface problem (1) lies in the space
$H^{1+\epsilon}(\Omega)$ with $\epsilon>0$ possibly close to zero [41], the
decomposed solutions $u_{1}^{[k]}$, $u_{2}^{[k]}$ of subproblems (1)), (2))
are regular under mild geometric assumptions on the boundary [13, 16]. More
precisely, it is assumed that the geometries of subdomains $\Omega_{1}$,
$\Omega_{2}$ are smooth enough to ensure the regularity $u_{1}^{[k]}\in
H^{2}(\Omega_{1})$, $u_{2}^{[k]}\in H^{2}(\Omega_{2})$ of the decomposed
solutions for error estimates in section 3.

###  2.2 Related Work

Traditional numerical methods for solving the elliptic interface problems (1)
can be roughly categorized into two mainstreams by using either an interface-
fitted or -unfitted mesh in the discretization of the domain. Provided a mesh
that aligns exactly with the interface, the interface jump conditions can be
naturally absorbed into a finite element formulation [4], achieving high
accuracy with nearly optimal error bounds [6]. However, a major drawback is
the requirement of an interface-fitted mesh generator, which could be time
consuming and technically involved for complicated interface geometries in two
and higher dimensions. As such, the second type of methods use interface-
unfitted meshes (e.g., the Cartesian mesh) and have attracted enormous
attention due to the easiness of mesh generation. A rich literature can be
found in this direction including, but not limited to immersed interface
methods [35], extended finite element methods [14], and many others. However,
there can be significant difficulties in enforcing jump conditions across the
complicated interfaces [5], as well as a rigorous error analysis that is
typically an ad hoc.

On the other hand, the Dirichlet-Neumann algorithm (1), 2), 3)) illustrated in
section 2.1 is essentially a continuous method for solving the elliptic
interface problem (1), thereby making it feasible to integrate with techniques
from the deep learning community [17, 28]. Consequently, one of the most
straightforward learning approach [38] is based on the Dirichlet-Neumann
algorithm written in terms of differential operators [54], namely, for $k\geq
1$,

(7) $\begin{array}[]{cl}\textnormal{Dirichlet
subproblem}\\!\\!&\\!\\!\left\\{\begin{array}[]{cl}-\nabla\cdot(c_{1}\nabla
u_{1}^{[k]})+u_{1}^{[k]}=f&\text{in}\ \Omega_{1},\\\ u_{1}^{[k]}=0&\text{on}\
\partial\Omega_{1}\setminus\Gamma,\\\ \\!\\!u_{1}^{[k]}=\rho
u_{2}^{[k-1]}+(1-\rho)u_{1}^{[k-1]}&\text{on}\ \Gamma,\end{array}\right.\\\
\textnormal{Neumann
subproblem}\\!\\!&\\!\\!\left\\{\begin{array}[]{cl}-\nabla\cdot(c_{2}\nabla
u_{2}^{[k+1]})+u_{2}^{[k+1]}=f&\text{in}\ \Omega_{2},\\\
u_{2}^{[k+1]}=0&\text{on}\ \partial\Omega_{2}\setminus\Gamma,\\\
\\!\\!c_{2}\nabla u_{2}^{[k+1]}\cdot\bm{n}_{2}=-q-c_{1}\nabla
u_{1}^{[k]}\cdot\bm{n}_{1}&\text{on}\ \Gamma,\end{array}\right.\end{array}$

where all the decomposed subproblems in (7) are approximately solved using the
physics-informed neural networks [48] or the deep Ritz method [62]. However,
though the training loss tends to decrease as the iteration proceeds, the
trained network solutions are prone to returning erroneous flux prediction
along the interfaces [8, 56, 1], thereby hampering the convergence of outer
iteration due to the incorrect transmission conditions being posed (see remark
3.5 or section 4). It is also interesting to note that similar idea has been
applied in areas of overlapping domain decomposition [36, 38, 49, 37], which
achieves empirical success to a certain extent since the updated interface
conditions are taken from the interior of neighbouring subregions rather than
its boundaries. Unfortunately, these overlapping methods face difficulties in
handling problems with large jumps in the coefficient or discontinuous
solutions.

Besides the meshless method mentioned above, recent advances in scientific
machine learning [2] have also led to an increasing interest in solving
elliptic interface problems (1) using various neural network models. For
instance, by removing the inhomogeneous boundary conditions, a deep Ritz-type
approach [62] is developed in [57] for solving the interface problem using one
single neural network. Based on the observation that the solution of interface
problem is typically piecewise-continuous, an improved learning method is to
employ a piecewise neural network to approximate the solution of (1) on
multiple subregions, where the training loss function could be built from
either the least square principle [19] or the variational formulation [18] of
the underlying system (1). In addition, setting appropriate penalty weights
among various loss terms is a crucial but painstaking process, which can be
mitigated by assigning adaptive weights to different loss terms [58, 56]. Note
that in the case of homogeneous interface jump conditions, the continuity of
averaged solution between neighbouring subdomains, as well as its first and
higher-order derivatives [25, 26, 50, 59, 23], can be explicitly enforced
through additional penalty terms posed on the interfaces. This family of
methods are quite general and parallelizable, however, high-order derivatives
at interface are meaningless for low-regularity solutions and no domain
knowledge of the underlying equation is used during the training mode, leading
to large errors near the interface that degrade the overall performance [24].
Unfortunately, all these algorithms fail to address the issue of erroneous
flux estimation at interface, which is the focus of this work. Another way of
dealing with the jump conditions is to design specific network architectures
that can capture the jump discontinuity across subdomain interfaces [32, 21,
55, 22], while the theoretical work will be needed to obtain a deeper
understanding.

##  3 Method

In this section, we start with the convergence analysis of our proposed
Dirichlet-Neumann learning algorithm for solving the elliptic interface
problem (1) in the “continuum-level”, which offers a novel theoretical
understanding of how to combine deep learning models [48, 62, 51] with ideas
from domain decomposition methods [47, 54]. Then, the implementation details
for realizing our proposed method is illustrated and summarized (see also
Figure 1).

Note that the convergence results of the iterative sequence of subregion
solutions (1), 2), 3)) have been well-studied at the continuous level [60], it
is then sufficient to show that these iterative solutions could be accurately
approximated using modern deep learning techniques.

###  3.1 Error Estimates

To realize the iterative schemes (1) – 3)) using neural networks [17], the
essential boundary conditions are treated in a “soft” manner by modifying the
energy functional with boundary penalty terms [62, 28, 51], namely,

(8) $\hat{u}_{1}^{[k]}=\operatorname*{arg\,min}_{\hat{u}_{1}\in
H^{1}(\Omega_{1})}\frac{1}{2}b_{1}(\hat{u}_{1},\hat{u}_{1})-(f,\hat{u}_{1})_{1}+\frac{\beta_{D}}{2}\left(\lVert\hat{u}_{1}\rVert_{L^{2}(\partial\Omega_{1}\cap\partial\Omega)}^{2}+\lVert\hat{u}_{1}-u_{\Gamma}^{[k]}\rVert_{L^{2}(\Gamma)}^{2}\right)$

for the Dirichlet subproblem (1)) at the $k$-th outer iteration, where
$\beta_{D}>0$ is a user-defined penalty coefficient [56]. On the other hand,
by extending the local solution $u_{2}\in V_{2}$ of decomposed problem (2)) to
its neighboring subregion (not relabelled),

(9) $R_{1}\gamma_{0}u_{2}(x)=u_{2}(x)\in V_{1}$

the modified loss functional associated with the Neumann subproblem (2)) gives

(10) $\hat{u}_{2}^{[k]}=\operatorname*{arg\,min}_{\hat{u}_{2}\in
H^{1}(\Omega)}\frac{1}{2}b_{2}(\hat{u}_{2},\hat{u}_{2})-(f,\hat{u}_{2})_{2}+b_{1}(\hat{u}_{1}^{[k]},\hat{u}_{2})-(f,\hat{u}_{2})_{1}+(q,\hat{u}_{2})_{L^{2}(\Gamma)}+\frac{\beta_{N}}{2}\lVert\hat{u}_{2}\rVert_{L^{2}(\partial\Omega)}^{2}$

where $\beta_{N}>0$ denotes another penalty coefficient. It is noteworthy that
the minimizer of functional (10) is now defined globally over the whole
domain, which differs from the traditional mesh-based treatment [54].

Before introducing the neural network parametrization of the unknown solutions
in (8) and (10) , the error estimations induced by the relaxation from exact
boundary conditions to a penalization-based approach are established in what
follows.

######  Theorem 3.1.

Let $u_{1}^{[k]}$ be the solution of problem (1)) and $\hat{u}_{1}^{[k]}$ the
solution of problem (8), then there holds

(11) $\lVert\hat{u}_{1}^{[k]}-u_{1}^{[k]}\rVert_{H^{1}(\Omega_{1})}\leq
C(\Omega_{1},u_{1}^{[k]})\frac{c_{1}}{\beta_{D}}\sqrt{\frac{\hat{c}_{1}}{\check{c}_{1}}}$

where $\hat{c}_{1}=\max\\{c_{1},1\\}$, $\check{c}_{1}=\min\\{c_{1},1\\}$, and
$C(\Omega_{1},u_{1}^{[k]})$ represents a generic constant that depends on the
subdomain $\Omega_{1}$ and the solution $u_{1}^{[k]}$ of Dirichlet subproblem
(1)).

######  Proof 3.2.

Step 1) We first denote by $L_{1}(\hat{u}_{1})$ the loss function of problem
(8), that is, a functional on $H^{1}(\Omega_{1})$

(12)
$\mathcal{L}_{1}(\hat{u}_{1})=\frac{1}{2}b_{1}(\hat{u}_{1},\hat{u}_{1})-(f,\hat{u}_{1})_{1}+\frac{\beta_{D}}{2}\left(\lVert\hat{u}_{1}\rVert_{L^{2}(\partial\Omega_{1}\cap\partial\Omega)}^{2}+\lVert\hat{u}_{1}-u_{\Gamma}^{[k]}\rVert_{L^{2}(\Gamma)}^{2}\right),$

and then derive the optimality conditions that are satisfied by the unique
global minimizer of (12). To be precise, the function $\hat{u}_{1}\in
H^{1}(\Omega_{1})$ is decomposed as a sum of two local functions, i.e.,
$\hat{u}_{1}=\hat{u}_{1}^{[k]}+g$ with $\hat{u}_{1}^{[k]}\in
H^{1}(\Omega_{1})$ satisfying

(13)
$\left\\{\begin{array}[]{cl}\displaystyle-\nabla\cdot(c_{1}\nabla\hat{u}_{1}^{[k]})+\hat{u}_{1}^{[k]}=f&\text{in}\
\Omega_{1},\\\
\displaystyle\hat{u}_{1}^{[k]}+c_{1}\beta_{D}^{-1}\nabla\hat{u}_{1}^{[k]}\cdot\bm{n}_{1}=0&\text{on}\
\partial\Omega_{1}\cap\partial\Omega,\\\
\displaystyle\hat{u}_{1}^{[k]}+c_{1}\beta_{D}^{-1}\nabla\hat{u}_{1}^{[k]}\cdot\bm{n}_{1}=u_{\Gamma}^{[k]}&\text{on}\
\Gamma,\\\ \end{array}\right.$

in the sense of distributions. Then, by applying the Green’s formula to (13),
a direct calculation of (12) implies that222More details can be found in the
technical Appendix A. for any $\hat{u}_{1}\in H^{1}(\Omega_{1})$,

$\mathcal{L}_{1}(\hat{u}_{1})=\mathcal{L}_{1}(\hat{u}_{1}^{[k]})+\int_{\Omega_{1}}\left(\frac{c_{1}}{2}|\nabla
g|^{2}+\frac{1}{2}|g|^{2}\right)dx+\frac{\beta_{D}}{2}\int_{\partial\Omega_{1}}|g|^{2}ds\geq\mathcal{L}_{1}(\hat{u}_{1}^{[k]}).$

Or, equivalently, the unique weak solution of problem (13) is the global
minimizer of functional (12). Notably, when comparing (13) with the original
Dirichlet subproblem (1)) (written in terms of differential operators), i.e.,

(14) $\left\\{\begin{array}[]{cl}\displaystyle-\nabla\cdot(c_{1}\nabla
u_{1}^{[k]})+u_{1}^{[k]}=f&\text{in}\ \Omega_{1},\\\ \displaystyle
u_{1}^{[k]}=0&\text{on}\ \partial\Omega_{1}\cap\partial\Omega,\\\
\displaystyle u_{1}^{[k]}=u_{\Gamma}^{[k]}&\text{on}\ \Gamma,\\\
\end{array}\right.$

the Dirichlet boundary condition is modified to be of a Robin type owing to
the boundary penalty treatment in (8).

Step 2) Now we are ready to quantitatively estimate the error induced from the
“soft” boundary enforcement, that is, the distance between the weak solutions
of (14) and (13). To deal with the inhomogeneous boundary conditions in (14)
and (13), let us write $u_{1}^{[k]}=w_{1}+g_{1}$ with an extension $g_{1}\in
V_{1}$ of $u_{\Gamma}^{[k]}$ into $\Omega_{1}$ [16], namely,

(15) $\left\\{\begin{array}[]{cl}\displaystyle-\nabla\cdot(c_{1}\nabla
w_{1})+w_{1}=f&\\!\\!\text{in}\ \Omega_{1},\\\ \displaystyle
w_{1}=0&\\!\\!\text{on}\ \partial\Omega_{1},\\\ \end{array}\right.\
\left\\{\begin{array}[]{cl}\displaystyle-\nabla\cdot(c_{1}\nabla
g_{1})+g_{1}=0&\\!\\!\text{in}\ \Omega_{1},\\\ \displaystyle
g_{1}=0&\\!\\!\text{on}\ \partial\Omega_{1}\cap\partial\Omega,\\\
\displaystyle g_{1}=u_{\Gamma}^{[k]}&\\!\\!\text{on}\ \Gamma,\\\
\end{array}\right.$

and $\hat{u}_{1}^{[k]}=\hat{w}_{1}+\hat{g}_{1}$ with another extension
$\hat{g}_{1}\in H^{1}(\Omega_{1})$ of $u_{\Gamma}^{[k]}$ into $\Omega_{1}$,
that is,

(16)
$\left\\{\begin{array}[]{cl}\displaystyle-\nabla\cdot(c_{1}\nabla\hat{w}_{1})+\hat{w}_{1}=f&\\!\\!\text{in}\
\Omega_{1},\\\
\displaystyle\hat{w}_{1}+c_{1}\beta_{D}^{-1}\nabla\hat{w}_{1}\cdot\bm{n}_{1}=0&\\!\\!\text{on}\
\partial\Omega_{1},\\\ \end{array}\right.\
\left\\{\begin{array}[]{cl}\displaystyle-\nabla\cdot(c_{1}\nabla\hat{g}_{1})+\hat{g}_{1}=0&\\!\\!\text{in}\
\Omega_{1},\\\
\displaystyle\hat{g}_{1}+c_{1}\beta_{D}^{-1}\nabla\hat{g}_{1}\cdot\bm{n}_{1}=0&\\!\\!\text{on}\
\partial\Omega_{1}\cap\partial\Omega,\\\
\displaystyle\hat{g}_{1}+c_{1}\beta_{D}^{-1}\nabla\hat{g}_{1}\cdot\bm{n}_{1}=u_{\Gamma}^{[k]}&\\!\\!\text{on}\
\Gamma,\\\ \end{array}\right.$

then it immediately follows from the triangle inequality that

$\lVert\hat{u}_{1}^{[k]}-u_{1}^{[k]}\rVert_{H^{1}(\Omega_{1})}=\lVert(\hat{w}_{1}+\hat{g}_{1})-(w_{1}+g_{1})\rVert_{H^{1}(\Omega_{1})}\leq\lVert\hat{w}_{1}-w_{1}\rVert_{H^{1}(\Omega_{1})}+\lVert\hat{g}_{1}-g_{1}\rVert_{H^{1}(\Omega_{1})}.$

Step 3) Based on the variational form, the extension function $g_{1}$ in (15)
satisfies

$b_{1}(g_{1},\hat{g}_{1})\\!=\\!\int_{\Omega_{1}}(c_{1}\nabla
g_{1}\cdot\nabla\hat{g}_{1}+g_{1}\hat{g}_{1})\,dx\\!=\\!\int_{\partial\Omega_{1}}(c_{1}\nabla
g_{1}\cdot\bm{n_{1}})\hat{g}_{1}\,ds\\!=\\!(c_{1}\nabla
g_{1}\cdot\bm{n_{1}},\hat{g}_{1})_{L^{2}(\partial\Omega_{1})}$

with $\hat{g}_{1}\in H^{1}(\Omega_{1})$ being used as the test function, while
the extension function $\hat{g}_{1}$ of (16) is the minimizer of energy
functional333Here and in what follows, we do not distinguish between the
general function and the optimal solution for notational simplicity.

$\begin{array}[]{cl}\mathcal{F}_{1}(\hat{g}_{1})\\!\\!\\!\\!\\!&=\displaystyle\frac{1}{2}b_{1}(\hat{g}_{1},\hat{g}_{1})+\frac{\beta_{D}}{2}\left(\lVert\hat{g}_{1}\rVert^{2}_{L^{2}(\partial\Omega_{1}\cap\partial\Omega)}+(\hat{g}_{1}-2u_{\Gamma}^{[k]},\hat{g}_{1})_{L^{2}(\Gamma)}\right)-b_{1}(g_{1},\hat{g}_{1})+(c_{1}\nabla
g_{1}\cdot\bm{n_{1}},\hat{g}_{1})_{L^{2}(\partial\Omega_{1})}\\\
&=\displaystyle\frac{1}{2}b_{1}(\hat{g}_{1}-g_{1},\hat{g}_{1}-g_{1})+\frac{\beta_{D}}{2}\left(\lVert\hat{g}_{1}+c_{1}\beta_{D}^{-1}\nabla
g_{1}\cdot\bm{n}_{1}\rVert^{2}_{L^{2}(\partial\Omega_{1}\cap\partial\Omega)}+\lVert\hat{g}_{1}+c_{1}\beta_{D}^{-1}\nabla
g_{1}\cdot\bm{n}_{1}-u_{\Gamma}^{[k]}\rVert^{2}_{L^{2}(\Gamma)}\right)\\\ &\ \
\ \displaystyle-\frac{1}{2}b_{1}(g_{1},g_{1})-\frac{\beta_{D}}{2}\left(\lVert
c_{1}\beta_{D}^{-1}\nabla
g_{1}\cdot\bm{n}_{1}\rVert^{2}_{L^{2}(\partial\Omega_{1}\cap\partial\Omega)}+\lVert
c_{1}\beta_{D}^{-1}\nabla
g_{1}\cdot\bm{n}_{1}-u_{\Gamma}^{[k]}\rVert^{2}_{L^{2}(\Gamma)}\right)\end{array}$

and from which we can conclude that the function $\hat{g}_{1}$ is also the
minimizer of

$\mathcal{G}_{1}(\hat{g}_{1})=\frac{1}{2}b_{1}(\hat{g}_{1}-g_{1},\hat{g}_{1}-g_{1})+\frac{\beta_{D}}{2}\left(\lVert\hat{g}_{1}+c_{1}\beta_{D}^{-1}\nabla
g_{1}\cdot\bm{n}_{1}\rVert^{2}_{L^{2}(\partial\Omega_{1}\cap\partial\Omega)}+\lVert\hat{g}_{1}+c_{1}\beta_{D}^{-1}\nabla
g_{1}\cdot\bm{n}_{1}-u_{\Gamma}^{[k]}\rVert^{2}_{L^{2}(\Gamma)}\right).$

On the one hand, by defining $\check{c}_{1}=\min\\{c_{1},1\\}$, it is obvious
that

(17)
$\mathcal{G}_{1}(\hat{g}_{1})\geq\int_{\Omega_{1}}\left(\frac{c_{1}}{2}|\nabla(\hat{g}_{1}-g_{1})|^{2}+\frac{1}{2}|\hat{g}_{1}-g_{1}|^{2}\right)dx\geq\frac{\check{c}_{1}}{2}\lVert\hat{g}_{1}-g_{1}\rVert_{H^{1}(\Omega_{1})}^{2}$

On the other hand, due to the fact that $g_{1}\in H^{2}(\Omega_{1})$ under
mild assumptions [13], we have by using the trace theorem [40] that $(\nabla
g_{1}\cdot\bm{n}_{1})|_{\partial\Omega_{1}}\in
H^{\frac{1}{2}}(\partial\Omega_{1})$ and therefore there exists a function
$\phi\in H^{1}(\Omega_{1})$ such that $\phi|_{\partial\Omega_{1}}=-(\nabla
g_{1}\cdot\bm{n}_{1})|_{\partial\Omega_{1}}$.

Then, by choosing $\bar{g}=c_{1}\beta_{D}^{-1}\phi+g_{1}$ and using the
boundary conditions of $g_{1}$ (15), the optimality of $\hat{g}_{1}$ among all
functions in $H^{1}(\Omega_{1})$ implies that

$\mathcal{G}_{1}(\hat{g}_{1})\leq\mathcal{G}_{1}(\bar{g})=c_{1}^{2}\beta_{D}^{-2}\int_{\Omega_{1}}b_{1}(\phi,\phi)+\frac{\beta_{D}}{2}\left(\lVert
g_{1}\rVert^{2}_{L^{2}(\partial\Omega_{1}\cap\partial\Omega)}+\lVert
g_{1}-u_{\Gamma}^{[k]}\rVert^{2}_{L^{2}(\Gamma)}\right)\leq\frac{\hat{c}_{1}}{2}c_{1}^{2}\beta_{D}^{-2}\lVert\phi\rVert_{H^{1}(\Omega_{1})}^{2}$

where $\hat{c}_{1}=\max\\{c_{1},1\\}$. As a direct result, we have by (17)
that

$\lVert\hat{g}_{1}-g_{1}\rVert_{H^{1}(\Omega_{1})}\leq\frac{c_{1}}{\beta_{D}}\sqrt{\frac{\hat{c}_{1}}{\check{c}_{1}}}\lVert\phi\rVert_{H^{1}(\Omega_{1})}.$

Step 4) It remains to prove that the solution $\hat{w}_{1}\in
H^{1}(\Omega_{1})$ of the Robin problem (16) can converge to the solution
$w_{1}\in H^{1}_{0}(\Omega_{1})$ of the Dirichlet problem (15) as
$\beta_{D}\to\infty$ [27, 11]. Similar as before, by employing $\hat{w}_{1}\in
H^{1}(\Omega_{1})$ as the test function in (15) and resorting to variational
form of (16), the function $\hat{w}_{1}\in H^{1}(\Omega_{1})$ minimizes

$\begin{array}[]{cl}\mathcal{I}_{1}(\hat{w}_{1})\\!\\!\\!\\!\\!&=\displaystyle\frac{1}{2}b_{1}(\hat{w}_{1},\hat{w}_{1})-(f,\hat{w}_{1})_{1}+\frac{\beta_{D}}{2}\lVert\hat{w}_{1}\rVert^{2}_{L^{2}(\partial\Omega_{1})}-b_{1}(w_{1},\hat{w}_{1})+(f,\hat{w}_{1})+(c_{1}\nabla
w_{1}\cdot\bm{n_{1}},\hat{w}_{1})_{L^{2}(\partial\Omega_{1})}\\\
&=\displaystyle\frac{1}{2}b_{1}(\hat{w}_{1}-w_{1},\hat{w}_{1}-w_{1})+\frac{\beta_{D}}{2}\lVert\hat{w}_{1}+c_{1}\beta_{D}^{-1}\nabla
w_{1}\cdot\bm{n}_{1}\rVert^{2}_{L^{2}(\partial\Omega_{1})}-\frac{1}{2}b_{1}(w_{1},w_{1})-\frac{\beta_{D}}{2}\lVert
c_{1}\beta_{D}^{-1}\nabla
w_{1}\cdot\bm{n}_{1}\rVert^{2}_{L^{2}(\partial\Omega_{1})}\end{array}$

and therefore is also the minimizer of the energy functional

$\mathcal{J}_{1}(\hat{w}_{1})=\frac{1}{2}b_{1}(\hat{w}_{1}-w_{1},\hat{w}_{1}-w_{1})+\frac{\beta_{D}}{2}\lVert\hat{w}_{1}+c_{1}\beta_{D}^{-1}\nabla
w_{1}\cdot\bm{n}_{1}\rVert^{2}_{L^{2}(\partial\Omega_{1})}.$

Note that $w_{1}\in H^{2}(\Omega_{1})$ under mild assumptions [13], the trace
theorem implies that $(\nabla w_{1}\cdot\bm{n}_{1})|_{\partial\Omega_{1}}\in
H^{\frac{1}{2}}(\partial\Omega_{1})$ and therefore there exists a function
$\varphi\in H^{1}(\Omega_{1})$ such that
$\varphi|_{\partial\Omega_{1}}=-(\nabla
w_{1}\cdot\bm{n}_{1})|_{\partial\Omega_{1}}$ [40]. As a consequence, by
employing a particular function $\bar{w}=c_{1}\beta_{D}^{-1}\varphi+w_{1}\in
H^{1}(\Omega_{1})$ and using the boundary condition of $w_{1}$ in (15),

$\mathcal{J}_{1}(\hat{w}_{1})\leq\mathcal{J}_{1}(\bar{w})=c_{1}^{2}\beta_{D}^{-2}\int_{\Omega_{1}}\left(\frac{c_{1}}{2}|\nabla\varphi|^{2}+\frac{1}{2}|\varphi|^{2}\right)dx+\frac{\beta_{D}}{2}\lVert
w_{1}\rVert^{2}_{L^{2}(\partial\Omega_{1})}\leq\frac{\hat{c}_{1}}{2}c_{1}^{2}\beta_{D}^{-2}\lVert\varphi\rVert_{H^{1}(\Omega_{1})}^{2}.$

On the other hand, it is obvious that

$\mathcal{J}_{1}(\hat{w}_{1})\geq\frac{1}{2}b_{1}(\hat{w}_{1}-w_{1},\hat{w}_{1}-w_{1})=\int_{\Omega_{1}}\left(\frac{c_{1}}{2}|\nabla(\hat{w}_{1}-w_{1})|^{2}+\frac{1}{2}|\hat{w}_{1}-w_{1}|^{2}\right)dx\geq\frac{\check{c}_{1}}{2}\lVert\hat{w}_{1}-w_{1}\rVert_{H^{1}(\Omega_{1})}^{2}$

which leads to the error estimation that completes the proof

$\lVert\hat{w}_{1}-w_{1}\rVert_{H^{1}(\Omega_{1})}\leq\frac{c_{1}}{\beta_{D}}\sqrt{\frac{\hat{c}_{1}}{\check{c}_{1}}}\lVert\varphi\rVert_{H^{1}(\Omega_{1})}.$

######  Theorem 3.3.

Assume that $\beta_{D}\to\infty$ in (8) (or Theorem 3.1), let $u_{2}^{[k]}$ be
the solution of problem (2)) and $\hat{u}_{2}^{[k]}$ the solution of problem
(10), then there holds444The minimizer $\hat{u}_{2}^{[k]}\in
H^{1}(\Omega_{1})$ of (10) is restricted on the subregion $\Omega_{2}$
(denoted by $\hat{u}_{2}^{[k]}|_{\Omega_{2}}$).

(18)
$\lVert\hat{u}_{2}^{[k]}|_{\Omega_{2}}-u_{2}^{[k]}\rVert_{H^{1}(\Omega_{2})}\leq
C(\Omega_{2},u_{2}^{[k]})\frac{c_{2}}{\beta_{N}}\sqrt{\frac{\hat{c}_{2}}{\check{c}_{2}}}$

where $\hat{c}_{2}=\max\\{c_{2},1\\}$, $\check{c}_{2}=\min\\{c_{2},1\\}$, and
$C(\Omega_{2},u_{2}^{[k]})$ represents a generic constant that depends on the
subdomain $\Omega_{2}$ and the solution $u_{2}^{[k]}$ of Neumann subproblem
(2)).

######  Proof 3.4.

Step 1) We first denote by $L_{2}(\hat{u}_{2})$ the loss function of problem
(10),

(19)
$\mathcal{L}_{2}(\hat{u}_{2})=\frac{1}{2}b_{2}(\hat{u}_{2},\hat{u}_{2})-(f,\hat{u}_{2})_{2}+b_{1}(\hat{u}_{1}^{[k]},\hat{u}_{2})-(f,\hat{u}_{2})_{1}+(q,\hat{u}_{2})_{L^{2}(\Gamma)}+\frac{\beta_{N}}{2}\lVert\hat{u}_{2}\rVert_{L^{2}(\partial\Omega)}^{2},$

i.e., a functional on $H^{1}(\Omega)$, and then derive the optimality
conditions that are satisfied by the global minimizer of (19). It’s of
particular noteworthy that the function $\hat{u}_{2}\in H^{1}(\Omega)$ is
defined over the entire domain, which greatly differs from the standard
Neumann subproblem (2)) that only depends on the subdomain $\Omega_{2}$.

As such, we decompose $\hat{u}_{2}\in H^{1}(\Omega)$ as a sum of two global
functions, namely, $\hat{u}_{2}=\hat{u}_{2}^{[k]}+g$, where the restriction of
$\hat{u}_{2}^{[k]}\in H^{1}(\Omega)$ on subregion $\Omega_{2}$ (not
relabelled) is required to satisfy the equations

(20)
$\left\\{\begin{array}[]{cl}\displaystyle-\nabla\cdot(c_{2}\nabla\hat{u}_{2}^{[k]})+\hat{u}_{2}^{[k]}=f&\text{in}\
\Omega_{2},\\\
\displaystyle\hat{u}_{2}^{[k]}+c_{2}\beta_{N}^{-1}\nabla\hat{u}_{2}^{[k]}\cdot\bm{n}_{2}=0&\text{on}\
\partial\Omega_{2}\cap\partial\Omega,\\\ \displaystyle
c_{2}\nabla\hat{u}_{2}^{[k]}\cdot\bm{n}_{2}=-q-c_{1}\nabla\hat{u}_{1}^{[k]}\cdot\bm{n}_{1}&\text{on}\
\Gamma,\\\ \end{array}\right.$

in the sense of distributions. Furthermore, the extension of function
$\hat{u}_{2}^{[k]}|_{\Omega_{2}}$ to the other subregion $\Omega_{1}$ is
required to be weakly differentiable and to satisfy the Robin boundary
condition

(21)
$\hat{u}_{2}^{[k]}+c_{1}\beta_{N}^{-1}\nabla\hat{u}_{1}^{[k]}\cdot\bm{n}_{1}=0\
\ \ \text{on}\ \partial\Omega_{1}\cap\partial\Omega,$

in the weak sense. Then, by applying the Green’s formula to (13), (20) and
using the jump condition (1), it can be deduced directly from (19) that555More
details can be found in the technical Appendix B. for any $\hat{u}_{2}\in
H^{1}(\Omega)$,

$\mathcal{L}_{2}(\hat{u}_{2})=\mathcal{L}_{2}(\hat{u}_{2}^{[k]})+\int_{\Omega_{2}}\left(\frac{c_{2}}{2}|\nabla
g|^{2}+\frac{1}{2}|g|^{2}\right)dx+\frac{\beta_{N}}{2}\int_{\partial\Omega}|g|^{2}ds\geq\mathcal{L}_{2}(\hat{u}_{2}^{[k]}),$

namely, the global minimizer of (19) can be characterized by the function
$\hat{u}_{2}^{[k]}\in H^{1}(\Omega)$ that satisfies (20) and (21). It is noted
that only the restricted solution $\hat{u}_{2}^{[k]}|_{\Omega_{2}}\in
H^{1}(\Omega_{2})$, or equivalently, the weak solution of subproblem (20) is
of interest for error estimation, which is still denoted by
$\hat{u}_{2}^{[k]}$ for short in the remaining of this proof. It is also
noteworthy that when compared to the original Neumann subproblem (2)) (written
in terms of differential operators), that is,

(22) $\left\\{\begin{array}[]{cl}\displaystyle-\nabla\cdot(c_{2}\nabla
u_{2}^{[k]})+u_{2}^{[k]}=f&\text{in}\ \Omega_{2},\\\ \displaystyle
u_{2}^{[k]}=0&\text{on}\ \partial\Omega_{2}\cap\partial\Omega,\\\
\displaystyle c_{2}\nabla u_{2}^{[k]}\cdot\bm{n}_{2}=-q-c_{1}\nabla
u_{1}^{[k]}\cdot\bm{n}_{1}&\text{on}\ \Gamma,\\\ \end{array}\right.$

a Robin boundary condition is imposed on
$\partial\Omega_{2}\cap\partial\Omega$ instead of the Dirichlet type, while
the interface jump condition can be maintained exactly as long as the penalty
coefficient $\beta_{D}$ in (1)) goes to infinity (see Theorem 3.1).

Step 2) To simplify the error analysis, it is assumed that the penalty
coefficient $\beta_{D}\to\infty$ in (8), i.e., $\hat{u}_{1}^{[k]}=u_{1}^{[k]}$
in (20) and (22). With the solution $\hat{u}_{2}^{[k]}\in H^{1}(\Omega_{2})$
being used as the test function, integration by parts for (22) implies that

$b_{2}(u_{2}^{[k]},\hat{u}_{2}^{[k]})+(q+c_{1}\nabla
u_{1}^{[k]}\cdot\bm{n}_{1},\hat{u}_{2}^{[k]})_{L^{2}(\Gamma)}-(c_{2}\nabla
u_{2}^{[k]}\cdot\bm{n}_{2},\hat{u}_{2}^{[k]})_{L^{2}(\partial\Omega_{2}\cap\partial\Omega)}=(f,\hat{u}_{2}^{[k]})_{2}$

which can be employed to reformulate the energy functional666We do not
distinguish between the general function and the optimal solution for
simplicity. of subproblem (20) as

$\begin{array}[]{cl}\mathcal{F}_{2}(\hat{u}_{2}^{[k]})\\!\\!\\!\\!\\!&=\displaystyle\frac{1}{2}b_{2}(\hat{u}_{2}^{[k]},\hat{u}_{2}^{[k]})-(f,\hat{u}_{2}^{[k]})_{2}+(q+c_{1}\nabla\hat{u}_{1}^{[k]}\cdot\bm{n}_{1},\hat{u}_{2}^{[k]})_{L^{2}(\Gamma)}+\frac{\beta_{N}}{2}\lVert\hat{u}_{2}^{[k]}\rVert^{2}_{L^{2}(\partial\Omega_{2}\cap\partial\Omega)}\\\
&=\displaystyle\frac{1}{2}b_{2}(\hat{u}_{2}^{[k]},\hat{u}_{2}^{[k]})-b_{2}(u_{2}^{[k]},\hat{u}_{2}^{[k]})+\frac{\beta_{N}}{2}\lVert\hat{u}_{2}^{[k]}\rVert^{2}_{L^{2}(\partial\Omega_{1}\cap\partial\Omega)}+(c_{2}\nabla
u_{2}^{[k]}\cdot\bm{n}_{2},\hat{u}_{2}^{[k]})_{L^{2}(\partial\Omega_{2}\cap\partial\Omega)}\\\
&=\displaystyle\frac{1}{2}b_{2}(\hat{u}_{2}^{[k]}-u_{2}^{[k]},\hat{u}_{2}^{[k]}-u_{2}^{[k]})+\frac{\beta_{N}}{2}\lVert\hat{u}_{2}^{[k]}+c_{2}\beta_{N}^{-1}\nabla
u_{2}^{[k]}\cdot\bm{n}_{2}\rVert^{2}_{L^{2}(\partial\Omega_{2}\cap\partial\Omega)}\\\
&\ \ \
\displaystyle-\frac{1}{2}b_{2}(u_{2}^{[k]},u_{2}^{[k]})-\frac{\beta_{N}}{2}\lVert
c_{2}\beta_{N}^{-1}\nabla
u_{2}^{[k]}\cdot\bm{n}_{2}\rVert^{2}_{L^{2}(\partial\Omega_{2}\cap\partial\Omega)}.\end{array}$

As a result, we can conclude that the weak solution $\hat{u}_{2}^{[k]}\in
H^{1}(\Omega_{2})$ of (20) is also the minimizer of functional

$\mathcal{G}_{2}(\hat{u}_{2}^{[k]})=\frac{1}{2}b_{2}(\hat{u}_{2}^{[k]}-u_{2}^{[k]},\hat{u}_{2}^{[k]}-u_{2}^{[k]})+\frac{\beta_{N}}{2}\lVert\hat{u}_{2}^{[k]}+c_{2}\beta_{N}^{-1}\nabla
u_{2}^{[k]}\cdot\bm{n}_{2}\rVert^{2}_{L^{2}(\partial\Omega_{2}\cap\partial\Omega)}$

Clearly, be defining $\check{c}_{2}=\min\\{c_{2},1\\}$, it is obvious that

(23)
$\mathcal{G}_{2}(\hat{u}_{2}^{[k]})\geq\frac{1}{2}b_{2}(\hat{u}_{2}^{[k]}-u_{2}^{[k]},\hat{u}_{2}^{[k]}-u_{2}^{[k]})=\int_{\Omega_{2}}\left(\frac{c_{2}}{2}|\nabla(\hat{u}_{2}^{[k]}-u_{2}^{[k]})|^{2}+\frac{1}{2}|\hat{u}_{2}^{[k]}-u_{2}^{[k]}|^{2}\right)dx\geq\frac{\check{c}_{2}}{2}\lVert\hat{u}_{2}^{[k]}-u_{2}^{[k]}\rVert_{H^{1}(\Omega_{2})}^{2}.$

On the other hand, note that $u_{2}^{[k]}\in H^{2}(\Omega_{2})$ under mild
assumptions [13, 16], we have by using the trace theorem [40] that $(\nabla
u_{2}^{[k]}\cdot\bm{n}_{2})|_{\partial\Omega_{2}}\in
H^{\frac{1}{2}}(\partial\Omega_{2})$ and therefore there exists a function
$\zeta\in H^{1}(\Omega_{2})$ such that $\zeta|_{\partial\Omega_{2}}=-(\nabla
u_{2}^{[k]}\cdot\bm{n}_{2})|_{\partial\Omega_{2}}$. In particular, the trace
operator has a continuous linear right inverse and therefore there holds

$\lVert\zeta\rVert_{H^{1}(\Omega_{2})}\leq C\lVert\nabla
u_{2}^{[k]}\cdot\bm{n}_{2}\rVert_{H^{\frac{1}{2}}(\partial\Omega_{2})}\leq
C\lVert u_{2}^{[k]}\rVert_{H^{2}(\Omega_{2})}$

where $C>0$ is a generic constant that only depends on the subregion
$\Omega_{1}$ [3].

Then by setting $\bar{u}_{2}=c_{2}\beta_{N}^{-1}\zeta+u_{2}^{[k]}$ and using
the homogeneous boundary condition of $u_{2}^{[k]}$ in (22), the optimality of
$\hat{u}_{2}^{[k]}$ among all functions in $H^{1}(\Omega_{2})$ implies

$\mathcal{G}_{2}(\hat{u}_{2}^{[k]})\leq\mathcal{G}_{2}(\bar{u}_{2})=c_{2}^{2}\beta_{N}^{-2}\int_{\Omega_{1}}\left(\frac{c_{1}}{2}|\nabla\zeta|^{2}+\frac{1}{2}|\zeta|^{2}\right)dx+\frac{\beta_{N}}{2}\lVert
u_{2}^{[k]}\rVert^{2}_{L^{2}(\partial\Omega_{1}\cap\partial\Omega)}\leq\frac{\hat{c}_{2}}{2}c_{2}^{2}\beta_{N}^{-2}\lVert\zeta\rVert_{H^{1}(\Omega_{1})}^{2}$

where $\hat{c}_{2}=\max\\{c_{2},1\\}$. Consequently, we have by (23) that

$\lVert\hat{u}_{2}^{[k]}-u_{2}^{[k]}\rVert_{H^{1}(\Omega_{1})}\leq\frac{c_{2}}{\beta_{N}}\sqrt{\frac{\hat{c}_{2}}{\check{c}_{2}}}\lVert\zeta\rVert_{H^{1}(\Omega_{1})}.$

To put it differently, by sending the penalty coefficients
$\beta_{D}\to\infty$ in (8) and $\beta_{N}\to\infty$ in (10), the minimizers
of our relaxed optimization problems (8, 10) could converge to that of the
classical Dirichlet-Neumann algorithm (1), 2)), which paves the way for
constructing the Dirichlet-Neumann learning algorithm in the practical
scenario. Moreover, our analysis (20) also sheds light on the tuning of
penalty coefficients, that is, $\beta_{N}$ should be increased as $c_{2}$
increases, otherwise the trained network solution would fail to capture the
Dirichlet boundary condition on $\partial\Omega_{2}\cap\partial\Omega$.

###  3.2 Dirichlet-Neumann Learning Algorithm

Next, the unknown solutions in (8) and (10) are parametrized using artificial
neural networks [17], that is,

$\hat{u}_{1}^{[k]}(x)=\hat{u}_{1}(x;\theta_{1}^{[k]})\ \ \ \text{and}\ \ \
\hat{u}_{2}^{[k]}(x)=\hat{u}_{2}(x;\theta_{2}^{[k]})$

where $\theta_{i}^{[k]}$ denotes the collection of trainable parameters at the
$k$-th outer iteration for $i=1$, 2. More specifically, the fully connected
neural network [17] or other kinds of architectures can be deployed to
construct the solution ansatz (more details about the model setup are included
in Appendix C). Moreover, thanks to the mesh-free characteristic of artificial
neural networks, the extension operator (9) can be realized in a very
straightforward way, that is,

(24) $R_{1}\gamma_{0}\hat{u}_{2}(x,\theta_{2})=\hat{u}_{2}(x,\theta_{2})$

which is required to be differentiable within $\Omega_{1}$ and to satisfy the
zero boundary values on $\partial\Omega_{1}\cap\partial\Omega$ through an
additional penalty term in (10). One can also employ a piecewise neural
network [19] to realize the extension operation (24).

Accordingly, to discretize the energy functionals (8) and (10), the routine
way of generating training points inside each subregion and at its boundary is
to use the Monte Carlo method or its variants [43], namely,

$X_{\Omega_{i}}=\big{\\{}x_{n}^{\Omega_{i}}\big{\\}}_{n=1}^{N_{\Omega_{i}}},\
\ \ X_{D_{i}}=\big{\\{}x_{n}^{D_{i}}\big{\\}}_{n=1}^{N_{D_{i}}},\ \ \
\text{and}\ \ \
X_{\Gamma}=\big{\\{}x_{n}^{\Gamma}\big{\\}}_{n=1}^{N_{\Gamma}},$

where $D_{i}:=\partial\Omega_{i}\cap\partial\Omega$, $N_{\Omega_{i}}$,
$N_{D_{i}}$ and $N_{\Gamma}$ denote the sample size of training datasets
$X_{\Omega_{i}}$, $X_{D_{i}}$ and $X_{\Gamma}$, respectively.

As a result, by defining the following empirical loss functions for $1\leq
i,j\leq 2$,

$\begin{array}[]{c}\displaystyle
L_{\Omega_{i}}\big{(}\hat{u}_{i}(x;\theta_{i})\big{)}=\frac{1}{N_{\Omega_{i}}}\sum_{n=1}^{N_{\Omega_{i}}}\left(\frac{c_{i}}{2}|\nabla\hat{u}_{i}(x_{n}^{\Omega_{i}};\theta_{i})|^{2}-f(x_{n}^{\Omega_{i}})\hat{u}_{i}(x_{n}^{\Omega_{i}};\theta_{i})\right),\\\
\displaystyle
L_{D_{i}}\big{(}\hat{u}_{j}(x;\theta_{j})\big{)}\\!=\\!\frac{1}{N_{D_{i}}}\sum_{n=1}^{N_{D_{i}}}|\hat{u}_{j}(x_{n}^{D_{i}};\theta_{j})|^{2},\
L_{\Gamma_{N}}\big{(}\hat{u}_{2}(x;\theta_{2})\big{)}\\!=\\!\frac{1}{N_{\Gamma}}\sum_{n=1}^{N_{\Gamma}}q(x_{n}^{\Gamma})\hat{u}_{2}(x_{n}^{\Gamma};\theta_{2}),\\\
\displaystyle
L_{\Gamma_{D}}\big{(}\hat{u}_{1}(x;\theta_{1}),u_{\Gamma}^{[k]}(x)\big{)}=\frac{1}{N_{\Gamma}}\sum_{n=1}^{N_{\Gamma}}|\hat{u}_{1}(x_{n}^{\Gamma};\theta_{1})-u_{\Gamma}^{[k]}(x_{n}^{\Gamma})|^{2},\end{array}$
$\begin{array}[]{c}\displaystyle
L_{N}\big{(}\hat{u}_{2}(x;\theta_{2}),\hat{u}_{1}(x;\theta_{1}^{[k]})\big{)}=\frac{1}{N_{\Omega_{1}}}\sum_{n=1}^{N_{\Omega_{1}}}\left(c_{1}\nabla\hat{u}_{1}(x_{n}^{\Omega_{1}};\theta_{1}^{[k]})\cdot\nabla\hat{u}_{2}(x_{n}^{\Omega_{1}};\theta_{2})-f(x_{n}^{\Omega_{1}})\hat{u}_{2}(x_{n}^{\Omega_{1}};\theta_{2})\right),\end{array}$

the learning task associated with the Dirichlet subproblem (8) is now given by

(25)
$\theta_{1}^{[k]}=\operatorname*{arg\,min}_{\theta_{1}}L_{\Omega_{1}}\big{(}\hat{u}_{1}(x;\theta_{1})\big{)}+\frac{\beta_{D}}{2}\left(L_{D_{1}}\big{(}\hat{u}_{1}(x;\theta_{1})\big{)}+L_{\Gamma_{D}}\big{(}\hat{u}_{1}(x;\theta_{1}),u_{\Gamma}^{[k]}(x)\big{)}\right),$

while that of the Neumann subproblem (10) takes on the form (referred to as
Compensated Deep Ritz Method)

(26)
$\theta_{2}^{[k]}=\operatorname*{arg\,min}_{\theta_{2}}L_{\Omega_{2}}\big{(}\hat{u}_{2}(x;\theta_{2})\big{)}+L_{N}\big{(}\hat{u}_{2}(x;\theta_{2}),\hat{u}_{1}(x;\theta_{1}^{[k]})\big{)}+L_{\Gamma_{N}}\big{(}\hat{u}_{2}(x;\theta_{2})\big{)}+\frac{\beta_{N}}{2}\left(L_{D_{1}}\big{(}\hat{u}_{2}(x;\theta_{2})\big{)}+L_{D_{2}}\big{(}\hat{u}_{2}(x;\theta_{2})\big{)}\right).$

It is noteworthy that although the loss value of (25) continues to decrease as
the training proceeds, the trained model is often observed to possess highly
fluctuating errors at and near the boundaries [56, 8, 1], which eventually
leads to a failure to approximate the interface flux and therefore hampers the
subsequent amendment (22) [52]. Fortunately, by employing our algorithm (26),
the Neumann subproblem can now be solved without explicitly enforcing the flux
jump condition (1), thereby ensuring the convergence of outer iteration in the
presence of inaccurate flux prediction.

Algorithm 1 Dirichlet-Neumann Learning Algorithm

% Initialization

– specify the network architecture $\hat{u}_{i}(x;\theta_{i})$ $(i=1$, 2) for
each subproblem;

– generate the Monte Carlo sampling points $X_{\Gamma}$, $X_{\Omega_{i}}$, and
$X_{D_{i}}$ for $i=1$, 2;

% Outer Iteration Loop

Start with the initial guess $u_{\Gamma}^{[0]}$ of unknown solution values at
interface $\Gamma$;

for $k\leftarrow 0$ to $K$ (maximum number of outer iterations) do

while stopping criteria are not satisfied do

% Dirichlet Subproblem-Solving via Deep Ritz Method

$\theta_{1}^{[k]}=\operatorname*{arg\,min}_{\theta_{1}}L_{\Omega_{1}}(\hat{u}_{1})+\frac{\beta_{D}}{2}\left(L_{D_{1}}(\hat{u}_{1})+L_{\Gamma}(\hat{u}_{1},u_{\Gamma}^{[k]})\right)$

% Neumann Subproblem-Solving via Compensated Deep Ritz Method

$\theta_{2}^{[k]}=\operatorname*{arg\,min}_{\theta_{2}}L_{\Omega_{2}}(\hat{u}_{2})+L_{N}(\hat{u}_{2},\hat{u}_{1}^{[k]})+L_{\Gamma_{N}}(\hat{u}_{2})+\frac{\beta_{N}}{2}\left(L_{D_{1}}(\hat{u}_{2})+L_{D_{2}}(\hat{u}_{2})\right)$

% Update of Unknown Solution Values at Interface

$u_{\Gamma}^{[k+1]}(x_{n}^{\Gamma})=\rho\hat{u}_{2}(x_{n}^{\Gamma};\theta_{2}^{[k]})+(1-\rho)u_{\Gamma}^{[k]}(x_{n}^{\Gamma}),\
\ i=1,\cdots,N_{\Gamma},$

end while

end for

More specifically, the detailed iterative scheme is summarized in Algorithm 1,
where the stopping criteria can be constructed by measuring the difference in
solutions between two consecutive iterations [38, 19]. We also note that as an
alternative to the deep Ritz method (25) [62], the Dirichlet subproblem (8)
can be solved using the physics-informed neural networks (abbreviated as PINNs
in what follows) [48], which is known to empirically work better for problems
with sufficient smooth solutions. To be precise, by incorporating the residual
of the equation (1) into the loss function

$L^{\textnormal{PINNs}}_{\Omega_{1}}\big{(}\hat{u}_{1}(x;\theta_{1})\big{)}=\frac{1}{N_{\Omega_{1}}}\sum_{n=1}^{N_{\Omega_{1}}}\left|-\nabla\cdot\big{(}c_{1}\nabla\hat{u}_{1}(x_{n}^{\Omega_{1}};\theta_{1})\big{)}+\hat{u}_{1}(x_{n}^{\Omega_{1}};\theta_{1})-f(x_{n}^{\Omega_{1}})\right|^{2}$

the learning task of Dirichlet subproblem (8) can alternatively be formulated
as

$\theta_{1}^{[k]}=\operatorname*{arg\,min}_{\theta_{1}}L^{\textnormal{PINNs}}_{\Omega_{1}}\big{(}\hat{u}_{1}(x;\theta_{1})\big{)}+\frac{\beta_{D}}{2}\left(L_{D_{1}}\big{(}\hat{u}_{1}(x;\theta_{1})\big{)}+L_{\Gamma_{D}}\big{(}\hat{u}_{1}(x;\theta_{1})\big{)}\right).$

Similar in spirit, the first term in the training loss function of Neumann
subproblem (26) could also be replaced with

(27)
$L^{\textnormal{PINNs}}_{\Omega_{2}}\big{(}\hat{u}_{2}(x;\theta_{2})\big{)}=\frac{1}{N_{\Omega_{2}}}\sum_{n=1}^{N_{\Omega_{2}}}\left|-\nabla\cdot\big{(}c_{2}\nabla\hat{u}_{2}(x_{n}^{\Omega_{2}};\theta_{2})\big{)}+\hat{u}_{2}(x_{n}^{\Omega_{2}};\theta_{2})-f(x_{n}^{\Omega_{2}})\right|^{2}$

which helps prevent the training process (26) from getting trapped in trivial
solutions especially when $c_{2}\gg 1$.

Moreover, our Dirichlet-Neumann learning algorithm can be easily modified to
adapt to a parallel computing environment, which is based on the theoretical
work [60] (see also Figure 1) and is left for future investigation.

######  Remark 3.5.

Note that when solving the Neumann subproblem, the flux transmission condition
in (22), i.e.,

$c_{2}\nabla u_{2}^{[k]}\cdot\bm{n}_{2}=-q-c_{1}\nabla
u_{1}^{[k]}\cdot\bm{n}_{1}\ \ \ \textnormal{on}\ \Gamma$

can be enforced using either our proposed algorithm or the DeepDDM scheme
[38], where both the neural networks are trained to approximately satisfy this
jump condition and eventually ended up with

$c_{2}\nabla\hat{u}_{2}^{[k]}\cdot\bm{n}_{2}\approx-q-
c_{1}\nabla\hat{u}_{1}^{[k]}\cdot\bm{n}_{1}\ \ \ \textnormal{on}\ \Gamma.$

As such, we immediately arrive at the error estimation

(28) $\nabla\hat{u}_{2}^{[k]}-\nabla
u_{2}^{[k]}\approx\frac{c_{1}}{c_{2}}\big{(}\nabla\hat{u}_{1}^{[k]}-\nabla
u_{1}^{[k]}\big{)}\ \ \ \textnormal{on}\ \Gamma$

which shows that the error incurred by an inaccurate right-hand-side term in
(28) may propagate to nearby interior points during training, and therefore
eventually leads to failure models of the trained network solution [7].

Unfortunately, the DeepDDM method [38] and other related work [19, 25, 26]
rely on on a direct computation of the Dirichlet subproblem (14) using neural
networks, followed by a straightforward evaluation of the trained network
solutions at interface to execute the subsequent operations. However, the
approximation errors for Dirichlet subproblem are typically found to
concentrate and fluctuate at and near the boundary [8, 1, 56, 52], hence the
representation of $\nabla u_{1}^{[k]}$ using $\nabla\hat{u}_{1}^{[k]}$ along
the interface may be of low accuracy. Consequently, the iterative solutions
would fail to converge to the exact solution once the error
$\nabla\hat{u}_{1}^{[k]}-\nabla u_{1}^{[k]}$ could not be eliminated by the
factor $c_{1}/c_{2}$, e.g., $c_{1}\approx c_{2}$ or $c_{1}=c_{2}$ (see section
4 for more experimental evidences).

On the contrary, our proposed scheme (26) is based on the variational
formulation of Dirichlet subproblem (14), which has a satisfactory accuracy in
representing $\nabla u_{1}^{[k]}|_{\Gamma}$ in the weak sense and therefore
substantially differs from other related methods. In other words, the right-
hand-side term of (28) is small whenever the coefficients are set to be
$c_{1}\ll c_{2}$ or $c_{1}\approx c_{2}$ (see section 4 for more details).

##  4 Numerical Experiments

In this section, numerical experiments on a series of elliptic interface
problems (1) are carried out to demonstrate the effectiveness of our
Dirichlet-Neumann learning algorithm (abbreviated as DNLA in what follows,
with text in bracket indicating the type of deep learning solver adopted for
solving the Dirichlet subproblem (14)). In all experiments, the network
architecture deployed for local problem is a fully connected neural network,
while the design of more sophisticated structures [24, 21] is left for future
investigation. Moreover, due to the smoothness requirement of extension
operator (24) (see also (9)), the hyperbolic tangent activation function is
used rather than the ReLU activation function [17].

As a comparison to our method, we also conduct experiments using the most
straightforward learning algorithm, that is, employing the PINNs for solving
all the decomposed subproblems (7) as proposed in [38], which is referred to
as DeepDDM. For a fair comparison, the same relaxation parameter, i.e., $\rho$
in (7) and (3)), is set for both methods, and the mean value and standard
deviation of the relative $L^{2}$ error evaluated on $100\times 100$ uniform
grid points (testing dataset)

$\epsilon=\frac{\lVert u-\hat{u}^{[k]}\rVert_{L^{2}(\Omega)}}{\lVert
u\rVert_{L^{2}(\Omega)}}\ \ \ \text{where}\ \ \
\hat{u}^{[k]}(x;\theta)=\left\\{\begin{array}[]{cl}\hat{u}_{1}(x;\theta_{1}^{[k]})&\text{if}\
x\in\Omega_{1},\\\ \hat{u}_{2}(x;\theta_{2}^{[k]})&\text{if}\
x\in\Omega_{2},\end{array}\right.$

are reported over 5 independent runs along the outer iterations, followed by a
comprehensive study that validates our discussion in remark 3.5. When the
maximum number of outer iteration is reached or the stopping criteria are met,
the resulting network solution is denoted by $\hat{u}(x;\theta)$ instead of
$\hat{u}^{[K]}(x;\theta)$ or $\hat{u}(x;\theta^{[K]})$ for simplicity.

More detailed experimental setups, e.g., the depth and width of the network
solution, the training and testing datasets, are briefly summarized in
Appendix C. Unless otherwise specified, we shall use the following
hyperparameter configurations. The AdamW optimizer [29] is adopted with an
initial learning rate set to be $0.1$, which is then multiplied by $0.1$ when
the training process reaches $60\%$ and $80\%$. In each outer iteration, the
Dirichlet or Neumann subproblem-solving, namely, (25) or (26) is terminated
after $3k$ and $1k$ epochs, respectively, and the model solution with the
minimum training loss is saved to execute the subsequent operations. All
experiments are implemented using PyTorch [46] on Nvidia GeForce RTX 3090
graphic cards. It is of particular noteworthy that the penalty coefficient
$\beta_{N}$ should be increased as the coefficient $c_{2}$ increases (see
Appendix C), which can be easily concluded from our analysis (20) but has not
been addressed in the literature.

###  4.1 Circle Interface in Two Dimension

Figure 3: Geometry of interface, exact solution with
$(c_{1},c_{2})=(1,10^{3})$ and $(1,1)$ for (29). Table 1: Error profiles
$\hat{u}(x,y;\theta)-u(x,y)$ of different methods for example (29) (see also
Appendix C).

Coefficients | DeepDDM | DNLM (PINNs) | DNLM (deep Ritz)
---|---|---|---
$(1,10^{3})$ with $\rho=1$ |  $K=2$ |  $K=2$ |  $K=2$
$(1,1)$ with $\rho=0.5$ |  $K=15$ |  $K=13$ |  $K=13$

Table 2: Comparison of different methods in terms of relative $L^{2}$ errors
for example (29).

Coefficients Outer Iteration | 2 | 4 | 10
---|---|---|---
$\\!(1,10^{3})\\!$ | DeepDDM | 0.013 $\pm$ 0.000 | 0.014 $\pm$ 0.001 | -
DNLA (PINNs) | 0.004 $\pm$ 0.001 | 0.004 $\pm$ 0.000 | -
 DNLA (deep Ritz) | 0.014 $\pm$ 0.001 | 0.014 $\pm$ 0.001 | -
$(1,1)$ | DeepDDM | 11.126 $\pm$ 0.096 | 2.997 $\pm$ 0.177 | 0.910 $\pm$ 0.356
DNLA (PINNs) | 11.110 $\pm$ 0.245 | 2.720 $\pm$ 0.067 | 0.095 $\pm$ 0.011
 DNLA (deep Ritz) | 9.456 $\pm$ 4.506 | 2.322 $\pm$ 1.092 | 0.122 $\pm$ 0.023

Table 3: Error profiles in terms of the pointwise $\ell^{2}$-norm of
$\nabla\hat{u}(x,y;\theta)-\nabla u(x,y)$ for example (29).

Coefficients | DeepDDM | DNLM (PINNs) | DNLM (deep Ritz)
---|---|---|---
$(1,1)$ with $\rho=0.5$ |  $K=15$ |  $K=13$ |  $K=13$

To begin with, we consider an elliptic interface problem (1) whose solution is
continuous [39], namely,

(29) $\begin{array}[]{cl}-\nabla\cdot\left(c_{i}\nabla
u_{i}(x,y)\right)=f(x,y)&\text{in}\ \Omega_{i},\\\
u_{i}(x,y)=g(x,y)&\text{on}\ \partial\Omega_{i}\cap\partial\Omega,\\\
u_{1}(x,y)=u_{2}(x,y)\ \text{and}\ \llbracket c(x)\nabla
u(x)\cdot\bm{n}\rrbracket=q(x,y)&\text{on}\ \Gamma,\end{array}$

where the entire computation domain is $\Omega=(-1,1)\times(-1,1)$ and the
interface $\Gamma=\big{\\{}(x,y)\,|\,x^{2}+y^{2}=0.25\big{\\}}$ is a circle
centered at point $(0,0)$ with radius $R=0.5$. The source term $f(x,y)$ and
the boundary data $g(x,y)$ on each subregion, i.e.,
$\Omega_{1}=\big{\\{}(x,y)\,|\,x^{2}+y^{2}<0.25\big{\\}}$ and
$\Omega_{2}=\Omega\setminus\Omega_{1}$, are calculated from the exact solution

$u(x,y)=\left\\{\begin{array}[]{cl}c_{1}^{-1}e^{10(r^{2}-R^{2})}&\text{if}\
r<R\ \ \ \ (\text{or in}\ \Omega_{1}),\\\
c_{2}^{-1}e^{10(R^{2}-r^{2})}+\big{(}c_{1}^{-1}-c_{2}^{-1}\big{)}e^{r^{2}-R^{2}}&\text{otherwise}\
(\text{or in}\ \Omega_{2}),\end{array}\right.$

where $r=\sqrt{x^{2}+y^{2}}$. More specifically, Figure 3 depicts the exact
solution $u(x,y)$ of (29) when the coefficients are given by
$(c_{1},c_{2})=(1,10^{3})$ and $(1,1)$ respectively, both of which should be
solved reasonably well so as to meet the robustness requirement with respect
to the varying coefficients. The initial guess of the unknown solution’s value
at interface is set to be $u_{\Gamma}^{[0]}=-1000x(x-1)y(y-1)+1$.

Table 3 shows the error profiles $|\hat{u}(x,y;\theta)-u(x,y)|$ of different
methods in a typical simulation, while the relative $L^{2}$ errors (mean $\pm$
standard deviation over 5 independent runs) are displayed in Table 3. Clearly,
the elliptic interface problem (29) with high-contrast coefficient
$(c_{1},c_{2})=(1,10^{3})$ can be solved with desired accuracy through the use
of the DeepDDM scheme [38], however, it fails to be effective when the
coefficient degenerates to the case of $(c_{1},c_{2})=$ $(1,1)$. On the other
hand, our proposed methods can achieve promising performance no matter
$c_{1}\ll c_{2}$ or $c_{1}=c_{2}$.

More specifically, when solving the Dirichlet subproblem using neural
networks, the trained model is often found to satisfy the underlying equations
but for overfitted boundary condition (see Table 3). As such, a direct
computation of $\nabla\hat{u}(x,y;\theta)$ at interface would lead to an
erroneous right-hand-side term in (28), which explains the failure of DeepDDM
when $c_{1}=c_{2}$. On the contrary, our methods is based on the variational
formulation of (13) and (14) with vanishing terms $\hat{u}_{1}^{[k]}$ and
$u_{1}^{[k]}$, i.e.,

(30) $\int_{\Omega_{1}}\big{(}\nabla\hat{u}_{1}^{[k]}-\nabla
u_{1}^{[k]}\big{)}\cdot\nabla
v_{1}\,dx=\int_{\Gamma}\big{(}\nabla\hat{u}_{1}^{[k]}-\nabla
u_{1}^{[k]}\big{)}\cdot\bm{n}_{1}v_{1}\,ds\ \ \ \textnormal{for any}\ v_{1}\in
V_{1},$

which allows the flux transmission without evaluating
$\nabla\hat{u}^{[k]}_{1}$ at interface, and therefore substantially differs
from other deep learning-based methods (see section 2.2). The experimental
results in Table 3 further demonstrate that our proposed algorithms are still
effective in the presence of erroneous flux data. Furthermore, it can be
inferred from (30) that the Neumann subproblem solver could benefit from a
good approximation of $\nabla u_{1}^{[k]}$ inside the subregion $\Omega_{1}$,
and hence DNLA (PINNs) empirically performs better than DNLA (deep Ritz) (see
Table 3).

###  4.2 Zigzag Interface in Two Dimension

Figure 4: Geometry of interface, exact solution with
$(c_{1},c_{2})=(1,10^{3})$ and $(1,1)$ for (31). Table 4: Error profiles
$|\hat{u}(x,y;\theta)-u(x,y)|$ of different methods for example (31) (see also
Appendix C).

Coefficients | DeepDDM | DNLA (PINNs) | DNLA (deep Ritz)
---|---|---|---
$(1,10^{3})$ with $\rho=1.0$ |  $(K=2)$ |  $(K=2)$ |  $(K=2)$
$(c_{1},c_{2})=(1,1)$ $\rho=0.5$ |  $(K=15)$ |  $(K=8)$ |  $(K=12)$

Table 5: Comparison of different methods in terms of relative $L^{2}$ errors
for example (31).

Coefficients Outer Iteration | 2 | 4 | 8
---|---|---|---
$\\!\\!\\!(1,10^{3})\\!\\!\\!$
| DeepDDM |  0.054 $\pm$ 0.023  |  0.041 $\pm$ 0.025  | -
DNLA (PINNs) |  0.059 $\pm$ 0.026  |  0.061 $\pm$ 0.019  | -
 DNLA (deep Ritz) |  0.041 $\pm$ 0.006  |  0.040 $\pm$ 0.009  | -
$\\!\\!\\!(1,1)\\!\\!\\!$
| DeepDDM |  8.453 $\pm$ 0.372  |  9.103 $\pm$ 0.715  |  10.569 $\pm$ 0.682
DNLA (PINNs) |  2.132 $\pm$ 0.772  |  0.157 $\pm$ 0.036  |  0.132 $\pm$ 0.033
 DNLA (deep Ritz) |  21.022 $\pm$ 4.628  |  5.183 $\pm$ 1.469  |  0.814 $\pm$ 0.030

Table 6: Error profiles in terms of the pointwise $\ell^{2}$-norm of
$\nabla\hat{u}(x,y;\theta)-\nabla u(x,y)$ for example (31).

Coefficients | DeepDDM | DNLA (PINNs) | DNLA (deep Ritz)
---|---|---|---
$(1,1)$ with $\rho=0.5$ |  $(K=15)$ |  $(K=8)$ |  $(K=12)$

We then consider a more complicated interface (i.e., the zigzag curve shown in
Figure 4) with non-homogeneous jump conditions to demonstrate the meshless
advantage of our proposed algorithm over the traditional mesh-based numerical
methods. More specifically, the boundary value problem (1) now takes on the
form, for $i=1$, 2,

(31) $\begin{array}[]{cl}-\nabla\cdot\left(c_{i}\nabla
u_{i}(x,y)\right)+u_{i}(x,y)=f(x,y)&\text{in}\ \Omega_{i},\\\
u_{i}(x,y)=g(x,y)&\text{on}\ \partial\Omega_{i}\cap\partial\Omega,\\\
u_{1}(x,y)-u_{2}(x,y)=p(x,y)\ \text{and}\ \llbracket c(x)\nabla
u(x)\cdot\bm{n}\rrbracket=q(x,y)&\text{on}\ \Gamma,\end{array}$

where $\Omega=(0,1)\times(0,1)$ and the interface is characterized by a zigzag
function

$x=Z_{3}(Z_{1}(20y-\text{floor}(20y))+Z_{2})+0.5$

with $Z_{1}=0.05(-1+2\text{mod}(\text{floor}(20y),2))$,
$Z_{2}=-0.05\text{mod}(\text{floor}(20x),2)$ and,
$Z_{3}=-2\text{mod}(\text{floor}(10x),2)+1$. Here, the source term $f(x,y)$,
boundary value $g(x,y)$, jump functions $p(x,y)$ and $q(x,y)$ are derived from
the true solution

$u(x,y)=\left\\{\begin{array}[]{cl}c_{1}^{-1}\sin(2\pi x)(\cos(2\pi
y)-1)&\text{in}\ \Omega_{1},\\\ c_{2}^{-1}\sin(2\pi x)(\cos(2\pi
y)-1)&\text{in}\ \Omega_{2}.\end{array}\right.$

Similar as before, the elliptic coefficient is chosen as
$(c_{1},c_{2})=(1,10^{3})$ and $(1,1)$ to evaluate the robustness against
varying coefficients, where the exact solutions are depicted in Figure 4.
Moreover, the initial guess of the unknown solution’s value at interface is
set to be $u_{\Gamma}^{[0]}(x,y)=\sin(2\pi x)(\cos(2\pi
y)-1)-1000x(x-1)y(y-1)$.

In a typical simulation, the error profiles and relative $L^{2}$ errors for
different methods are reported in Table 6, Table 6, and Table 6. It can be
concluded that the DeepDDM approach performs well in the case of high-contrast
coefficients, but fails to converge when the coefficients degenerate to
$c_{1}=c_{2}=1$, which reveals the fact that DeepDDM breaks down rapidly in
face of inaccurate flux predictions.

On the contrary, the iterative solutions using our algorithms can match well
with the exact solution whenever the coefficient satisfies $c_{1}\ll c_{2}$ or
$c_{1}=c_{2}$ (see Table 6, Table 6). Note that the relaxation parameter is
set to be $\rho=1$ for the case of high-contrast coefficients
$(c_{1},c_{2})=(1,10^{3})$ and therefore the optimal convergence [45, 15] is
nearly obtained in two iterations. Accordingly, we then choose $\rho=0.5$ for
the case of $(c_{1},c_{2})=(1,1)$, where additional iterations are required
due to the numerical error caused by the deep learning solvers. Notably, when
solving the Dirichlet subproblem through PINNs, the second-order derivatives
of $\hat{u}_{1}(x;\theta_{1})$ are explicitly incorporated into the training
loss function, leading to a satisfactory approximation of
$\nabla\hat{u}_{1}(x;\theta_{1})$ inside the subdomain $\Omega_{1}$. As such,
DNLA (PINNs) shows better empirical performance than DNLA (deep Ritz) as
reported in Table 6.

It is also noteworthy that even though the outer iteration converges by using
our methods, the erroneous gradient approximation along the subregion
interfaces remains inevitable as depicted in Table 6. Fortunately, our flux
transmission does not rely on the direct computation of
$\nabla\hat{u}_{1}(x;\theta_{1})$ at interface, which plays a key role.

###  4.3 Checkerboard Lattice in Two Dimension

Figure 5: Geometry of interface, exact solution with
$(c_{1},c_{2})=(1,10^{3})$ and $(1,1)$ for (32). Table 7: Error profiles
$|\hat{u}(x,y;\theta)-u(x,y)|$ of different methods for example (32) (see also
Appendix C).

Coefficients | DeepDDM | DNLA (PINNs) | DNLA (deep Ritz)
---|---|---|---
$(1,10^{3})$ with $\rho=1.0$ |  $(K=2)$ |  $(K=2)$ |  $(K=2)$
$(1,1)$ with $\rho=0.5$ |  $(K=15)$ |  $(K=9)$ |  $(K=10)$

Table 8: Comparison of different methods in terms of relative $L^{2}$ errors
for example (32).

Coefficients Outer Iterations | 1 | 2 | 10
---|---|---|---
$\\!\\!\\!(1,10^{3})\\!\\!\\!$
| DeepDDM | 0.448 $\pm$ 0.002 | 0.015 $\pm$ 0.006 | -
DNLA (PINNs) | 0.447 $\pm$ 0.002 | 0.020 $\pm$ 0.009 | -
 DNLA (deep Ritz) | 0.457 $\pm$ 0.009 | 0.053 $\pm$ 0.001 | -
$\\!\\!\\!(1,1)\\!\\!\\!$
| DeepDDM | 1.820 $\pm$ 0.125 | 8.981 $\pm$ 0.486 | 459.855 $\pm$ 352.586
DNLA (PINNs) | 1.052 $\pm$ 0.244 | 0.644 $\pm$ 0.069 | 0.134 $\pm$ 0.040
 DNLA (deep Ritz) | 0.944 $\pm$ 0.302 | 0.720 $\pm$ 0.428 | 0.365 $\pm$ 0.166

Next, we consider the interface problem with domain partitioned in a
checkerboard fashion (see Figure 5),

(32) $\begin{array}[]{cl}-\nabla\cdot\left(c_{i}\nabla
u_{i}(x,y)\right)+u_{i}(x,y)=f(x,y)&\text{in}\ \Omega_{i},\\\
u_{i}(x,y)=g(x,y)&\text{on}\ \partial\Omega_{i}\cap\partial\Omega,\\\
u_{1}(x,y)-u_{2}(x,y)=p(x,y)\ \text{and}\ \llbracket c(x)\nabla
u(x)\cdot\bm{n}\rrbracket=q(x,y)&\text{on}\ \Gamma,\end{array}$

for $i=1$, 2, where the interior cross-point is well-known to require extra
treatments when using mesh-based methods [54, 47]. Similar as before, the
source term, boundary and interface conditions of our underlying problem are
derived from the exact solution with a large jump across interface (see Figure
5)

$u(x,y)=\left\\{\begin{array}[]{cl}c_{1}^{-1}\sin(4\pi y)\sin(4\pi
x)&\text{in}\ \Omega_{1},\\\ c_{2}^{-1}4x(x-1)y(y-1)&\text{in}\
\Omega_{2}.\end{array}\right.$

As can be seen from Table 8 and Table 8, the DeepDDM, DNLA (PINNs), and DNLA
(deep Ritz) methods are able to provide satisfactory performance for interface
problem (32) with high-contrast coefficients $(c_{1},c_{2})=(1,10^{3})$.
However, the DeepDDM scheme diverges as the ratio $c_{1}/c_{2}$ fails to
eliminate the gradient approximation error (28), while our methods still work
reasonably well. Here, the initial guess takes on the form
$u_{\Gamma}^{[0]}(x,y)=\sin(4\pi y)\sin(4\pi x)+100x(x-1)^{3}y(y-1)^{3}$.

Table 9: Error profiles in terms of the pointwise $\ell^{2}$-norm of $\nabla\hat{u}(x,y;\theta)-\nabla u(x,y)$ for example (32). Coefficients | DeepDDM | DNLA (PINNs) | DNLA (Deep Ritz)
---|---|---|---
$(1,1)$ with $\rho=0.5$ |  $(K=15)$ |  $(K=6)$ |  $(K=10)$

Thanks to the variational treatment of the flux jump condition (26), the
erroneous flux prediction in Dirichlet subproblem (see Table 9) is a key
threat to the DeepDDM scheme [38] but not our proposed algorithms, which
demonstrates that our methods are robust with respect to the varying
coefficients. Furthermore, DNLA (PINNs) achieves a smaller approximation error
$\nabla\hat{u}_{1}^{[k]}-\nabla u_{1}^{[k]}$ inside the subdomain $\Omega_{1}$
(see Table 9), and therefore performs better than that of DNLA (deep Ritz)
(see Table 8).

##  5 Conclusion

Motivated by the observation that domain decomposition methods are essentially
continuous schemes, a novel mesh-free method is developed in this work to
realize the classical Dirichlet-Neumann algorithm using neural networks, which
paves the way for effectively solving the elliptic interface problem with
high-contrast coefficients and irregular interfaces. Instead of a direct
substitution of the local problem solver with modern deep learning tools, a
rigorous error analysis is first established to obtain the error bound of
boundary penalty treatment for both the Dirichlet and Neumann subproblems, and
therefore sheds light on the design of training loss functions that can
alleviate the burden of erroneous flux prediction along subdomain interfaces.
In other words, the accuracy of flux transmission (or derivatives of the
network solution at interfaces plays a key role in combining non-overlapping
domain decomposition methods with techniques from the deep learning community,
which has not yet been well addressed by the existing deep learning-based
algorithms and is the focus of this work. Moreover, the tedious tuning of
penalty coefficients can also get insights from our theoretical study. A wide
variety of numerical examples are carried out to validate the effectiveness
and robustness of our proposed methods, achieving promising performance
especially in the presence of erroneous flux data during outer iterations.

We believe that our theoretical and experimental studies can also be
generalized to the Robin-Robin algorithm [47, 45], while substantial
improvements can be made by employing coarse grid correction [42], adaptive
sampling strategy [19], special network architectures [37], and more
comprehensive error analysis [27, 11].

##  Appendix A Detailed Proof of Step 1) in Theorem 3.1

Recall that the function $\hat{u}_{1}\in H^{1}(\Omega_{1})$ is decomposed as a
sum of two local functions, i.e., $\hat{u}_{1}=\hat{u}_{1}^{[k]}+g$, then

$\begin{array}[]{cl}\mathcal{L}_{2}(\hat{u}_{2})\\!\\!\\!\\!\\!&=\displaystyle\frac{1}{2}b_{1}(\hat{u}_{1}^{[k]}+g,\hat{u}_{1}^{[k]}+g)-(f,\hat{u}_{1}^{[k]}+g)_{1}+\frac{\beta_{D}}{2}\Big{(}\lVert\hat{u}_{1}^{[k]}+g\rVert_{L^{2}(\partial\Omega_{1}\cap\partial\Omega)}^{2}\\\
&\displaystyle\ \ \ +\
\lVert\hat{u}_{1}^{[k]}+g-u_{\Gamma}^{[k]}\rVert_{L^{2}(\Gamma)}^{2}\Big{)}\\\
&\displaystyle=\mathcal{L}_{1}(\hat{u}_{1}^{[k]})+b_{1}(\hat{u}_{1}^{[k]},g)-(f,g)_{1}+\frac{1}{2}b_{1}(g,g)+\frac{\beta_{D}}{2}\lVert
g\rVert_{L^{2}(\partial\Omega_{1})}^{2}\\\ &\displaystyle\ \ \ +\
\beta_{D}\left((\hat{u}_{1}^{[k]},g)_{L^{2}(\partial\Omega_{1}\cap\partial\Omega)}+(\hat{u}_{1}^{[k]}-u_{\Gamma}^{[k]},g)_{L^{2}(\Gamma)}\right)\end{array}$

Note that the function $\hat{u}_{1}^{[k]}\in H^{1}(\Omega_{1})$ is required to
satisfy the equation (13) in the sense of distributions, that is, for any
$g\in H^{2}(\Omega_{1})$,

$b_{1}(\hat{u}_{1}^{[k]},g)+\beta_{D}\left((\hat{u}_{1}^{[k]},g)_{L^{2}(\partial\Omega_{1}\cap\partial\Omega)}+(\hat{u}_{1}^{[k]}-u_{\Gamma}^{[k]},g)_{L^{2}(\Gamma)}\right)=(f,g)_{1}$

and therefore we arrive at for any $\hat{u}_{1}\in H^{1}(\Omega_{1})$

$\mathcal{L}_{1}(\hat{u}_{1})=\mathcal{L}_{1}(\hat{u}_{1}^{[k]})+\int_{\Omega_{1}}\left(\frac{c_{1}}{2}|\nabla
g|^{2}+\frac{1}{2}|g|^{2}\right)dx+\frac{\beta_{D}}{2}\int_{\partial\Omega_{1}}|g|^{2}ds\geq\mathcal{L}_{1}(\hat{u}_{1}^{[k]}).$

##  Appendix B Detailed Proof of Step 1) in Theorem 3.3

Recall that the function $\hat{u}_{2}\in H^{1}(\Omega)$ is decomposed as a sum
of two global functions, i.e., $\hat{u}_{2}=\hat{u}_{2}^{[k]}+g$, then

$\begin{array}[]{cl}\mathcal{L}_{2}(\hat{u}_{2})\\!\\!\\!\\!\\!&=\displaystyle\frac{1}{2}b_{2}(\hat{u}_{2}^{[k]}+g,\hat{u}_{2}^{[k]}+g)-(f,\hat{u}_{2}^{[k]}+g)_{2}+b_{1}(\hat{u}_{1}^{[k]},\hat{u}_{2}^{[k]}+g)-(f,\hat{u}_{2}^{[k]}+g)_{1}\\\
&\displaystyle\ \ \ +\
(q,\hat{u}_{2}^{[k]}+g)_{L^{2}(\Gamma)}+\frac{\beta_{N}}{2}\lVert\hat{u}_{2}^{[k]}+g\rVert_{L^{2}(\partial\Omega)}^{2}\\\
&\displaystyle=\mathcal{L}_{2}(\hat{u}_{2}^{[k]})+b_{2}(\hat{u}_{2}^{[k]},g)-(f,g)_{2}+b_{1}(\hat{u}_{1}^{[k]},g)-(f,g)_{1}+(q,g)_{L^{2}(\Gamma)}\\\
&\displaystyle\ \ \ +\
\beta_{N}(\hat{u}_{2}^{[k]},g)_{L^{2}(\partial\Omega)}+\frac{1}{2}b_{2}(g,g)+\frac{\beta_{N}}{2}\lVert
g\rVert_{L^{2}(\partial\Omega)}^{2}\end{array}$

where $g\in H^{1}(\Omega)$ is arbitrary and is defined over the entire domain.
Note that $\hat{u}_{1}^{[k]}\in H^{1}(\Omega_{1})$ and
$\hat{u}_{2}^{[k]}|_{\Omega_{2}}\in H^{1}(\Omega_{2})$ are required to satisfy
the equations (13) and (20) in the sense of distributions, namely,

$b_{1}(\hat{u}_{1}^{[k]},g_{1})-(c_{1}\nabla\hat{u}_{1}^{[k]}\cdot\bm{n}_{1},g_{1})_{L^{2}(\partial\Omega_{1})}=(f,g_{1})_{1}\
\ \ \textnormal{for any}\ g_{1}\in H^{1}(\Omega_{1})$

and

$b_{2}(\hat{u}_{2}^{[k]},g_{2})-(c_{2}\nabla\hat{u}_{2}^{[k]}\cdot\bm{n}_{2},g_{2})_{L^{2}(\partial\Omega_{2})}=(f,g_{2})_{2}\
\ \ \textnormal{for any}\ g_{2}\in H^{1}(\Omega_{2})$

respectively, we then have by the boundary conditions imposed in (21) and (20)
that for any $\hat{u}_{2}\in H^{1}(\Omega)$

$\begin{array}[]{cl}\mathcal{L}_{2}(\hat{u}_{2})\\!\\!\\!\\!\\!&=\displaystyle\mathcal{L}_{2}(\hat{u}_{2}^{[k]})+(c_{1}\nabla\hat{u}_{1}^{[k]}\cdot\bm{n}_{1}+\beta_{N}\hat{u}_{2}^{[k]},g)_{L^{2}(\partial\Omega_{1}\cap\partial\Omega)}+(c_{2}\nabla\hat{u}_{2}^{[k]}\cdot\bm{n}_{2}+\beta_{N}\hat{u}_{2}^{[k]},g)_{L^{2}(\partial\Omega_{2}\cap\partial\Omega)}\\\
&\displaystyle\ \ \ \ +\
(c_{1}\nabla\hat{u}_{1}^{[k]}\cdot\bm{n}_{1}+c_{2}\nabla\hat{u}_{2}^{[k]}\cdot\bm{n}_{2}+q,g)_{L^{2}(\Gamma)}+\frac{1}{2}b_{2}(g,g)+\frac{\beta_{N}}{2}\lVert
g\rVert_{L^{2}(\partial\Omega)}^{2}\\\
&\displaystyle=\mathcal{L}_{2}(\hat{u}_{2}^{[k]})+\int_{\Omega_{2}}\left(\frac{c_{2}}{2}|\nabla
g|^{2}+\frac{1}{2}|g|^{2}\right)dx+\frac{\beta_{N}}{2}\int_{\partial\Omega}|g|^{2}ds\geq\mathcal{L}_{2}(\hat{u}_{2}^{[k]}).\end{array}$

##  Appendix C Detailed Experimental Results

The hyperparameter configuration used for our numerical experiments is briefly
summarized in Table 10, together with the iterative network solutions using
different methods for each example.

Table 10: Hyperparameter configuration used for our experiments. |   Train Datasets ($N_{\Omega_{i}}$, $N_{D_{i}}$, $N_{\Gamma}$)  |   Penalty Coeffs ($\beta_{D}$, $\beta_{N}$)  |   Network (depth, width)
---|---|---|---
$c_{2}=10^{3}$ | DeepDDM | $(20k,5k,5k)$ | 400 | $(6,50)$
DNLA (PINNs) | $(20k,5k,5k)$ | $(800,800k)$ | $(6,50)$
 DNLA (deep Ritz)  | $(20k,5k,5k)$ | $(800,800k)$ | $(6,50)$
$c_{2}=1$ | DeepDDM | $(20k,5k,5k)$ | 400 | $(6,50)$
DNLA (PINNs) | $(20k,5k,5k)$ | $(800,800)$ | $(6,50)$
 DNLA (deep Ritz)  | $(20k,5k,5k)$ | $(800,800)$ | $(6,50)$
Table 11: The iterative solutions $\hat{u}^{[k]}(x,y;\theta)$ for numerical
example (29) with $(c_{1},c_{2})=(1,10^{3})$.

DeepDDM DNLA (PINNs) DNLA (deep Ritz)

Table 12: The iterative solutions $\hat{u}^{[k]}(x,y;\theta)$ for numerical
example (29) with $(c_{1},c_{2})=(1,1)$.

DeepDDM DNLA (PINNs) DNLA (deep Ritz)

Table 13: The iterative solutions $\hat{u}^{[k]}(x,y;\theta)$ for numerical
example (31) with $(c_{1},c_{2})=(1,10^{3})$.

DeepDDM DNLA (PINNs) DNLA (deep Ritz)

Table 14: The iterative solutions $\hat{u}^{[k]}(x,y;\theta)$ for numerical
example (31) with $(c_{1},c_{2})=(1,1)$.

DeepDDM DNLA (PINNs) DNLA (deep Ritz)

Table 15: The iterative solutions $\hat{u}^{[k]}(x,y;\theta)$ for numerical
example (32) with $(c_{1},c_{2})=(1,10^{3})$.

DeepDDM DNLA (PINNs) DNLA (deep Ritz)

Table 16: The iterative solutions $\hat{u}^{[k]}(x,y;\theta)$ for numerical
example (32) with $(c_{1},c_{2})=(1,1)$.

DeepDDM DNLA (PINNs) DNLA (deep Ritz)

## Acknowledgments

The computations were done on the high performance computers of School of
Mathematical Sciences, Tongji University.

## References

  * [1] C. Bajaj, L. McLennan, T. Andeen, and A. Roy, Robust learning of physics informed neural networks, arXiv preprint arXiv:2110.13330, (2021).
  * [2] N. Baker, F. Alexander, T. Bremer, A. Hagberg, Y. Kevrekidis, H. Najm, M. Parashar, A. Patra, J. Sethian, S. Wild, et al., Workshop report on basic research needs for scientific machine learning: Core technologies for artificial intelligence, tech. report, 2019.
  * [3] J. H. Bramble, Interpolation between Sobolev spaces in Lipschitz domains with an application to multigrid theory, Mathematics of Computation, 64 (1995), pp. 1359–1365.
  * [4] S. C. Brenner, L. R. Scott, and L. R. Scott, The mathematical theory of finite element methods, vol. 3, Springer, 2008.
  * [5] L. Chen, H. Wei, and M. Wen, An interface-fitted mesh generator and virtual element methods for elliptic interface problems, Journal of Computational Physics, 334 (2017), pp. 327–348.
  * [6] Z. Chen and J. Zou, Finite element methods and their convergence for elliptic and parabolic interface problems, Numerische Mathematik, 79 (1998), pp. 175–202.
  * [7] A. Daw, J. Bu, S. Wang, P. Perdikaris, and A. Karpatne, Rethinking the importance of sampling in physics-informed neural networks, arXiv preprint arXiv:2207.02338, (2022).
  * [8] T. Dockhorn, A discussion on solving partial differential equations using neural networks, arXiv preprint arXiv:1904.07200, (2019).
  * [9] V. Dolean, P. Jolivet, and F. Nataf, An introduction to domain decomposition methods: algorithms, theory, and parallel implementation, SIAM, 2015.
  * [10] S. Dong and Z. Li, Local extreme learning machines and domain decomposition for solving linear and nonlinear partial differential equations, Computer Methods in Applied Mechanics and Engineering, 387 (2021), p. 114129.
  * [11] C. Duan, Y. Jiao, Y. Lai, X. Lu, Q. Quan, and J. Z. Yang, Analysis of deep Ritz methods for laplace equations with dirichlet boundary conditions, arXiv preprint arXiv:2111.02009, (2021).
  * [12] V. Dwivedi, N. Parashar, and B. Srinivasan, Distributed learning machines for solving forward and inverse problems in partial differential equations, Neurocomputing, 420 (2021), pp. 299–316.
  * [13] L. C. Evans, Partial Differential Equations, vol. 19, American Mathematical Society, 2010.
  * [14] T.-P. Fries and T. Belytschko, The extended/generalized finite element method: an overview of the method and its applications, International Journal for Numerical Methods in Engineering, 84 (2010), pp. 253–304.
  * [15] M. J. Gander and O. Dubois, Optimized Schwarz methods for a diffusion problem with discontinuous coefficient, Numerical Algorithms, 69 (2015), pp. 109–144.
  * [16] D. Gilbarg, N. S. Trudinger, D. Gilbarg, and N. Trudinger, Elliptic partial differential equations of second order, vol. 224, Springer, 1977.
  * [17] I. Goodfellow, Y. Bengio, and A. Courville, Deep learning, MIT Press, 2016.
  * [18] H. Guo and X. Yang, Deep unfitted Nitsche method for elliptic interface problems, arXiv preprint arXiv:2107.05325, (2021).
  * [19] C. He, X. Hu, and L. Mu, A mesh-free method using piecewise deep neural network for elliptic interface problems, Journal of Computational and Applied Mathematics, 412 (2022), pp. 114–358.
  * [20] A. Heinlein, A. Klawonn, M. Lanser, and J. Weber, Combining machine learning and domain decomposition methods for the solution of partial differential equations—A review, GAMM-Mitteilungen, 44 (2021), p. e202100001.
  * [21] W.-F. Hu, T.-S. Lin, and M.-C. Lai, A discontinuity capturing shallow neural network for elliptic interface problems, Journal of Computational Physics, 469 (2022), p. 111576.
  * [22] W.-F. Hu, T.-S. Lin, Y.-H. Tseng, and M.-C. Lai, A hybrid neural-network and finite-difference method for solving poisson equation with jump discontinuities on interfaces, arXiv preprint arXiv:2210.05523, (2022).
  * [23] Z. Hu, A. D. Jagtap, G. E. Karniadakis, and K. Kawaguchi, When do extended physics-informed neural networks (XPINNs) improve generalization?, arXiv preprint arXiv:2109.09444, (2021).
  * [24] Z. Hu, A. D. Jagtap, G. E. Karniadakis, and K. Kawaguchi, Augmented physics-informed neural networks (APINNs): A gating network-based soft domain decomposition methodology, arXiv preprint arXiv:2211.08939, (2022).
  * [25] A. D. Jagtap and G. E. Karniadakis, Extended physics-informed neural networks (XPINNs): A generalized space-time domain decomposition based deep learning framework for nonlinear partial differential equations, Communications in Computational Physics, 28 (2020), pp. 2002–2041.
  * [26] A. D. Jagtap, E. Kharazmi, and G. E. Karniadakis, Conservative physics-informed neural networks on discrete domains for conservation laws: Applications to forward and inverse problems, Computer Methods in Applied Mechanics and Engineering, 365 (2020), pp. 113–028.
  * [27] Y. Jiao, Y. Lai, Y. Lo, Y. Wang, and Y. Yang, Error analysis of deep Ritz methods for elliptic equations, arXiv preprint arXiv:2107.14478, (2021).
  * [28] G. E. Karniadakis, I. G. Kevrekidis, L. Lu, P. Perdikaris, S. Wang, and L. Yang, Physics-informed machine learning, Nature Reviews Physics, 3 (2021), pp. 422–440.
  * [29] D. P. Kingma and J. Ba, Adam: A method for stochastic optimization, arXiv preprint arXiv:1412.6980, (2014).
  * [30] I. E. Lagaris, A. Likas, and D. I. Fotiadis, Artificial neural networks for solving ordinary and partial differential equations, IEEE transactions on neural networks, 9 (1998), pp. 987–1000.
  * [31] I. E. Lagaris, A. C. Likas, and D. G. Papageorgiou, Neural-network methods for boundary value problems with irregular boundaries, IEEE Transactions on Neural Networks, 11 (2000), pp. 1041–1049.
  * [32] M.-C. Lai, C.-C. Chang, W.-S. Lin, W.-F. Hu, and T.-S. Lin, A shallow Ritz method for elliptic problems with singular sources, Journal of Computational Physics, (2022), p. 111547.
  * [33] R. J. LeVeque, Finite difference methods for ordinary and partial differential equations: steady-state and time-dependent problems, SIAM, 2007\.
  * [34] R. J. LeVeque et al., Finite volume methods for hyperbolic problems, vol. 31, Cambridge university press, 2002.
  * [35] R. J. LeVeque and Z. Li, The immersed interface method for elliptic equations with discontinuous coefficients and singular sources, SIAM Journal on Numerical Analysis, 31 (1994), pp. 1019–1044.
  * [36] K. Li, K. Tang, T. Wu, and Q. Liao, D3M: A deep domain decomposition method for partial differential equations, IEEE Access, 8 (2019), pp. 5283–5294.
  * [37] S. Li, Y. Xia, Y. Liu, and Q. Liao, A deep domain decomposition method based on Fourier features, arXiv preprint arXiv:2205.01884, (2022).
  * [38] W. Li, X. Xiang, and Y. Xu, Deep domain decomposition method: Elliptic problems, in Mathematical and Scientific Machine Learning, PMLR, 2020, pp. 269–286.
  * [39] Z. Li and K. Ito, The immersed interface method: Numerical solutions of PDEs involving interfaces and irregular domains, SIAM, 2006.
  * [40] J. L. Lions and E. Magenes, Non-homogeneous boundary value problems and applications: Vol. 1, vol. 181, Springer Science & Business Media, 2012\.
  * [41] D. Mercier, Minimal regularity of the solutions of some transmission problems, Mathematical Methods in the Applied Sciences, 26 (2003), pp. 321–348.
  * [42] V. Mercier, S. Gratton, and P. Boudier, A coarse space acceleration of deep-DDM, arXiv preprint arXiv:2112.03732, (2021).
  * [43] N. Metropolis and S. Ulam, The Monte Carlo method, Journal of the American Statistical Association, 44 (1949), pp. 335–341.
  * [44] J. Müller and M. Zeinhofer, Error estimates for the deep Ritz method with boundary penalty, in Mathematical and Scientific Machine Learning, PMLR, 2022, pp. 215–230.
  * [45] X. Na and X. Xu, Domain decomposition methods for elliptic problems with high contrast coefficients revisited, arXiv preprint arXiv:2212.12216, (2022).
  * [46] A. Paszke, S. Gross, S. Chintala, G. Chanan, E. Yang, Z. DeVito, Z. Lin, A. Desmaison, L. Antiga, and A. Lerer, Automatic differentiation in pytorch, (2017).
  * [47] A. Quarteroni and A. Valli, Domain decomposition methods for partial differential equations, Oxford University Press, 1999.
  * [48] M. Raissi, P. Perdikaris, and G. E. Karniadakis, Physics-informed neural networks: A deep learning framework for solving forward and inverse problems involving nonlinear partial differential equations, Journal of Computational physics, 378 (2019), pp. 686–707.
  * [49] H. Sheng and C. Yang, PFNN-2: A domain decomposed penalty-free neural network method for solving partial differential equations, arXiv preprint arXiv:2205.00593, (2022).
  * [50] K. Shukla, A. D. Jagtap, and G. E. Karniadakis, Parallel physics-informed neural networks via domain decomposition, Journal of Computational Physics, 447 (2021), p. 110683.
  * [51] J. Sirignano and K. Spiliopoulos, DGM: A deep learning algorithm for solving partial differential equations, Journal of Computational Physics, 375 (2018), pp. 1339–1364.
  * [52] Q. Sun, X. Xu, and H. Yi, Domain decomposition learning methods for solving elliptic problems, arXiv preprint arXiv:2207.10358, (2022).
  * [53] A. Taghibakhshi, N. Nytko, T. Zaman, S. MacLachlan, L. Olson, and M. West, Learning interface conditions in domain decomposition solvers, arXiv preprint arXiv:2205.09833, (2022).
  * [54] A. Toselli and O. Widlund, Domain decomposition methods: Algorithms and theory, vol. 34, Springer Science & Business Media, 2004.
  * [55] Y.-H. Tseng, T.-S. Lin, W.-F. Hu, and M.-C. Lai, A cusp-capturing PINN for elliptic interface problems, arXiv preprint arXiv:2210.08424, (2022).
  * [56] S. Wang, Y. Teng, and P. Perdikaris, Understanding and mitigating gradient flow pathologies in physics-informed neural networks, SIAM Journal on Scientific Computing, 43 (2021), pp. A3055–A3081.
  * [57] Z. Wang and Z. Zhang, A mesh-free method for interface problems using the deep learning approach, Journal of Computational Physics, 400 (2020), p. 108963.
  * [58] S. Wu and B. Lu, INN: Interfaced neural networks as an accessible meshless approach for solving interface pde problems, Journal of Computational Physics, 470 (2022), p. 111588.
  * [59] W. Wu, X. Feng, and H. Xu, Improved deep neural networks with domain decomposition in solving partial differential equations, Journal of Scientific Computing, 93 (2022), pp. 1–34.
  * [60] D. Yang, A parallel nonoverlapping Schwarz domain decomposition method for elliptic interface problems, (1997).
  * [61] D. Yang, Finite elements for elliptic problems with wild coefficients, Mathematics and computers in simulation, 54 (2000), pp. 383–395.
  * [62] B. Yu et al., The deep Ritz method: a deep learning-based numerical algorithm for solving variational problems, Communications in Mathematics and Statistics, 6 (2018), pp. 1–12.
  * [63] Y. Zang, G. Bao, X. Ye, and H. Zhou, Weak adversarial networks for high-dimensional partial differential equations, Journal of Computational Physics, 411 (2020), p. 109409.