Least-Squares Neural Network (LSNN) Method For Linear Advection-Reaction
  Equation: Non-constant Jumps
The least-squares ReLU neural network (LSNN) method was introduced and
studied for solving linear advection-reaction equation with discontinuous
solution in \cite{Cai2021linear,cai2023least}. The method is based on an
equivalent least-squares formulation and \cite{cai2023least} employs ReLU
neural network (NN) functions with $\lceil \log_2(d+1)\rceil+1$-layer
representations for approximating solutions. In this paper, we show
theoretically that the method is also capable of accurately approximating
non-constant jumps along discontinuous interfaces that are not necessarily
straight lines. Theoretical results are confirmed through multiple numerical
examples with $d=2,3$ and various non-constant jumps and interface shapes,
showing that the LSNN method with $\lceil \log_2(d+1)\rceil+1$ layers
approximates solutions accurately with degrees of freedom less than that of
mesh-based methods and without the common Gibbs phenomena along discontinuous
interfaces having non-constant jumps.
Zhiqiang Cai, Junpyo Choi, Min Liu
# Least-Squares Neural Network (LSNN) Method

###### Abstract

The least-squares ReLU neural network method (LSNN) was introduced and studied for solving linear advection-reaction equation with discontinuous solution in [4, 5]. The method is based on an equivalent least-squares formulation and employs ReLU neural network (NN) functions with a \(\lceil\log_{2}(d+1)\rceil+1\) layer representation for approximating the solution. In this paper, we show theoretically that the method is also capable of approximating a non-constant jump along the discontinuous interface of the underlying problem that is not necessarily a straight line. Numerical results for test problems with various non-constant jumps and interfaces show that the LSNN method with \(\lceil\log_{2}(d+1)\rceil+1\) layers approximates the solution accurately with DoFs less than that of mesh-based methods and without the common Gibbs phenomena along the discontinuous interface.

Least-Squares Method, ReLU Neural Network, Linear Advection-Reaction Equation, Discontinuous Solution

65N15, 65N99

## 1 Introduction

Let \(\Omega\) be a bounded domain in \(\mathbb{R}^{d}\) (\(d\geq 2\)) with Lipschitz boundary \(\partial\Omega\), and denote the advective velocity field by \(\boldsymbol{\beta}(\mathbf{x})=(\beta_{1},\cdots,\beta_{d})^{T}\in C^{0}( \bar{\Omega})^{d}\). Define the inflow part of the boundary \(\Gamma=\partial\Omega\) by

\[\Gamma_{-}=\{\mathbf{x}\in\Gamma:\,\boldsymbol{\beta}(\mathbf{x})\cdot \boldsymbol{n}(\mathbf{x})<0\} \tag{1}\]

with \(\boldsymbol{n}(\mathbf{x})\) being the unit outward normal vector to \(\Gamma\) at \(\mathbf{x}\in\Gamma\). Consider the linear advection-reaction equation

\[\left\{\begin{array}{rcl}u_{\boldsymbol{\beta}}+\gamma\,u&=&f,&\mbox{in}\; \;\Omega,\\ u&=&g,&\mbox{on}\;\;\Gamma_{-},\end{array}\right. \tag{2}\]

where \(u_{\boldsymbol{\beta}}\) denotes the directional derivative of \(u\) along \(\boldsymbol{\beta}\). Assume that \(\gamma\in C^{0}(\bar{\Omega})\), \(f\in L^{2}(\Omega)\), and \(g\in L^{2}(\Gamma_{-})\) are given scalar-valued functions.

A major challenge in numerical simulation is that the solution of (2) is discontinuous along an interface because of a discontinuous inflow boundary condition, where the discontinuous interface can be the streamline from the inflow boundary. Traditional mesh-based numerical methods often exhibit oscillations near a discontinuity (called the Gibbs phenomena) and may not be extended to nonlinear hyperbolic conservation laws.

The least-squares ReLU neural network (LSNN) method for solving (2) with discontinuous solution was introduced and studied in [4, 5]. The method is based on an equivalent least-squares formulation studied in ([2, 7]) and employs ReLU neural network (NN) functions with a \(\lceil\log_{2}(d+1)\rceil+1\) layer representation for approximating the solution. The LSNN method is capable of automatically approximating the discontinuous solution since the free hyperplanes of ReLU NN functions adapt to the solution (see [3, 4, 5]). Compared to various adaptive mesh refinement (AMR) algorithms that locate the discontinuous interface through local mesh refinement (see, e.g., [6, 8, 10]), the LSNN method is much more effective in terms of the number of the degrees of freedom.
Approximation property of ReLU NN functions to a step function was recently studied in [4, 5]. In particular, we showed theoretically that two or \(\lceil\log_{2}(d+1)\rceil+1\) layer ReLU NN functions are necessary and sufficient to approximate a step function with any given accuracy \(\varepsilon>0\) when the discontinuous interface is a hyperplane or general hyper-surface, respectively. This approximation property was used to establish _a priori_ error estimates of the LSNN method.

The jump of the discontinuous solution of (2) is generally non-constant when the reaction coefficient \(\gamma\) is non-zero. The main purpose of this paper is to establish _a priori_ error estimates (see Theorem 3) for the LSNN method without making the assumption that the jump is constant. To this end, we decompose the solution as the sum of discontinuous and continuous parts (see (3)), so that the discontinuous part of the solution can be described as a cylindrical surface on one subdomain and zero otherwise. Then we construct a continuous piecewise linear (CPWL) function with sharp transition layers along the discontinuous interface to approximate the discontinuous part accurately by \(\mathcal{O}(J_{2}\sqrt{\varepsilon_{2}})+\mathcal{O}(J_{3}\sqrt{\varepsilon_ {3}})\) error, where \(\varepsilon_{2}\) and \(\varepsilon_{3}\) are, respectively, the width of the transition layers and a given positive number bounding the differences of the function values and directional derivatives between the shifted surface and piecewise plane functions constructed to approximate the surface on the corresponding subdomain (see Lemmas 5 and 5 and Theorem 3). From [1, 5, 11, 12], we know that the CPWL function is a ReLU NN function \(\mathbb{R}^{d}\to\mathbb{R}\) with a \(\lceil\log_{2}(d+1)\rceil+1\) layer representation, from which it follows that the discontinuous part of the solution can be approximated by this class of functions for any prescribed accuracy.

The rest of the paper is organized as follows. In Section 2, we briefly review and discuss properties of ReLU NN functions and the LSNN method in [5]. Then theoretical convergence analysis is conducted in Section 3, showing that discretization error of the method mainly depends on the continuous part of the decomposition of the solution. Finally, to demonstrate the effectiveness of the method, we provide numerical results for test problems with various non-constant jumps in Section 4.

## 2 ReLU NN Functions and LSNN Method

This section briefly reviews properties of ReLU neural network (NN) functions and the least-squares ReLU neural network (LSNN) method in [5]. For a given positive integer \(n\), denote the collection of all ReLU NN functions from \(\mathbb{R}^{d}\) to \(\mathbb{R}\) that have a representation with depth \(L\) and total number of hidden neurons \(n\) by \(\mathcal{M}(d,1,L,n)\), and the collection of all ReLU NN functions from \(\mathbb{R}^{d}\) to \(\mathbb{R}\) with a \(L\) layer representation by \(\mathcal{M}(d,1,L)\). Then

\[\mathcal{M}(d,1,L)=\bigcup_{n\in\mathbb{N}}\mathcal{M}(d,1,L,n). \tag{4}\]

**Proposition 2** **(**see [1, 5])**: _The collection of all continuous piecewise linear (CPWL) functions from \(\mathbb{R}^{d}\) to \(\mathbb{R}\) is equal to \(\mathcal{M}(d,1,\lceil\log_{2}(d+1)\rceil+1)\), i.e., the collection of all ReLU NN functions from \(\mathbb{R}^{d}\) to \(\mathbb{R}\) that have a representation with depth \(\lceil\log_{2}(d+1)\rceil+1\)._

**Proposition 3**: _\(\mathcal{M}(d,1,L,n)\subset\mathcal{M}(d,1,L,n+1)\)._

Define the least-squares (LS) functional

\[\mathcal{L}(v;\mathbf{f})=\|v_{\boldsymbol{\beta}}+\gamma\,v-f\|_{0,\Omega}^{ 2}+\|v-g\|_{-\boldsymbol{\beta}}^{2}, \tag{5}\]

where \(\mathbf{f}=(f,g)\) and \(\|\cdot\|_{-\boldsymbol{\beta}}\) is given by

\[\|v\|_{-\boldsymbol{\beta}}=\langle v,v\rangle_{-\boldsymbol{\beta}}^{1/2}= \left(\int_{\Gamma_{-}}|\boldsymbol{\beta}\boldsymbol{\cdot}\boldsymbol{n}|\,v ^{2}\,ds\right)^{1/2}.\]

The LS formulation of problem (2) is to seek \(u\in V_{\boldsymbol{\beta}}\) such that

\[\mathcal{L}(u;\mathbf{f})=\min_{v\in V_{\boldsymbol{\beta}}}\mathcal{L}(v; \mathbf{f}), \tag{6}\]
where \(V_{\boldsymbol{\beta}}=\{v\in L^{2}(\Omega):v_{\boldsymbol{\beta}}\in L^{2}(\Omega)\}\) is a Hilbert space that is equipped with the norm

\[\left|\!\left|\!\left|v\right|\!\right|\!\right|_{\boldsymbol{\beta}}=\left( \left|\!\left|v\right|\!\right|_{0,\Omega}^{2}+\left|\!\left|v_{\boldsymbol{ \beta}}\right|\!\right|\!\right|_{0,\Omega}^{2}\right)^{1/2}.\]

Then the corresponding LS and discrete LS approximations are, respectively, to find \(u_{{}_{N}}\in\mathcal{M}(d,1,L,n)\) such that

\[\mathcal{L}\big{(}u_{{}_{N}};\mathbf{f}\big{)}=\min_{v\in\mathcal{M}(d,1,L,n)} \mathcal{L}\big{(}v;\mathbf{f}\big{)}, \tag{4}\]

and to find \(u_{{}_{T}}^{N}\in\mathcal{M}(d,1,L,n)\) such that

\[\mathcal{L}_{{}_{T}}\big{(}u_{{}_{T}}^{N};\mathbf{f}\big{)}=\min_{v\in \mathcal{M}(d,1,L,n)}\mathcal{L}_{{}_{T}}\big{(}v;\mathbf{f}\big{)}, \tag{5}\]

where \(\mathcal{L}_{{}_{T}}\big{(}v;\mathbf{f}\big{)}\) is the discrete LS functional (see [4]).

## 3 Error Estimate

In this section, we establish error estimate of the LSNN method for the linear advection-reaction equation with a non-constant jump along a discontinuous interface. For simplicity, we restrict ourselves in two dimensions.

To this end, assume the advection velocity field \(\boldsymbol{\beta}\) is piecewise constant. That is, there exists a partition of the domain \(\Omega\) such that \(\boldsymbol{\beta}\) has the same direction but possibly different magnitude at each interior point of subdomain. Without loss of generality, assume that there are only two sub-domains: \(\Omega=\Upsilon_{1}\cup\Upsilon_{2}\) and that the inflow boundary data \(g(\mathbf{x})\) is discontinuous at only one point \(\mathbf{x}_{0}\in\Gamma_{-}\) with \(g(\mathbf{x}_{0}^{-})=\alpha_{1}\) and \(g(\mathbf{x}_{0}^{+})=\alpha_{2}\) from different sides. (Figure 1(a) depicts \(\Upsilon_{1}\) and \(\Upsilon_{2}\) as the left-upper and the right-lower triangles, respectively.)

Let \(I\) be the streamline emanating from \(\mathbf{x}_{0}\), then the discontinuous interface \(I\) divides the domain \(\Omega\) into two sub-domains: \(\Omega=\Omega_{1}\cup\Omega_{2}\), where \(\Omega_{1}\) and \(\Omega_{2}\) are the left-lower and the right-upper sub-domains separated by the discontinuous interface \(I\), respectively (see Figure 1(a)). The corresponding solution \(u\) of (2) is discontinuous across the interface \(I\) and is piecewise smooth with respect to the partition \(\{\Omega_{1},\Omega_{2}\}\). For a given \(\varepsilon_{2}>0\), take a \(\varepsilon_{2}\) neighborhood around the interface \(I\) in the direction of \(\boldsymbol{\beta}\) as in Figure 1(b).

Next, we estimate error in the sub-domain \(\Upsilon_{i}\), say, \(\Upsilon_{2}\). To further simplify error estimate, assume that

\[\Upsilon_{i}=(0,1)\times(0,1),\quad\boldsymbol{\beta}(\mathbf{x})=(0,v_{2}( \mathbf{x}))^{T},\quad\text{and}\quad\mathbf{x}_{0}=(x_{0},0)\text{ for }x_{0}\in(0,1).\]

These assumptions imply that restriction of the interface \(I\) in \(\Upsilon_{i}\) is a vertical line segment

\[I=\{(x_{0},y)\in\Upsilon_{i}:y\in(0,1)\},\]

and that \(\Upsilon_{i}\) is partitioned into two sub-domains

\[\Omega_{1i}=\{\mathbf{x}=(x,y)\in\Upsilon_{i}:\,x<x_{0}\}\quad\text{and}\quad \Omega_{2i}=\{\mathbf{x}=(x,y)\in\Upsilon_{i}:\,x>x_{0}\}\]

(see Figure 1(c)).

In \(\Upsilon_{i}=\Omega_{1i}\cup\Omega_{2i}\cup I\), let \(u_{1}\) and \(u_{2}\) be the solutions of (2) defined only on \(\Omega_{1i}\) with constant inflow boundary conditions \(g=\alpha_{1}\) and \(\alpha_{2}\) on \(\{(x,0):\,x\in[0,x_{0}]\}\), respectively. (When \(\Upsilon_{i}\) is different from \(\Upsilon_{2}\), the discontinuous point is not \(\mathbf{x}_{0}\) but an interior point of the domain \(\Omega\), and the values of the solution \(u\) at that discontinuous point from different sides are taken as constant inflow boundary conditions.) Set \(a(\mathbf{x})=u_{1}(\mathbf{x})-u_{2}(\mathbf{x})\) and let \(\chi(\mathbf{x})\) be a piecewise discontinuous function defined by

\[\chi(\mathbf{x})=\left\{\begin{array}{ll}a(\mathbf{x}),&\mathbf{x}\in\Omega _{1i},\\ 0,&\mathbf{x}\in\Omega_{2i},\end{array}\right. \tag{6}\]
then the solution \(u\) of (2) has the following decomposition (see [5])

\[u(\mathbf{x})=\hat{u}(\mathbf{x})+\,\chi(\mathbf{x})\,\mbox{ in }\Upsilon_{i}. \tag{3}\]

Here, \(\hat{u}(\mathbf{x})=u(\mathbf{x})-\,\chi(\mathbf{x})\) is clearly piecewise smooth; moreover, it is also continuous in \(\Upsilon_{i}\) since \(\hat{u}\big{|}_{I}=u_{2}\big{|}_{I}\) from both sides. Then we have the following error estimate and postpone its proof to section 5.

**Theorem 3.1**: _For any given \(\varepsilon_{2}>\) and \(\varepsilon_{3}>0\), on \(\Upsilon_{i}\), there exists a CPWL function \(p_{i}(\mathbf{x})\) such that_

\[\left|\!\left|\!\left|\chi-p_{i}\right|\!\right|\!\right|_{\boldsymbol{\beta} }\leq D_{1}\sqrt{\varepsilon_{2}}+D_{2}\sqrt{\varepsilon_{3}}. \tag{4}\]

The proof of the theorem is postponed to Section 5.

**Remark 3.2**: _We now construct a CPWL function \(p(\mathbf{x})\) on \(\Omega\) given by_

\[p(\mathbf{x})=p_{i}(\mathbf{x}),\ \mathbf{x}\in\Upsilon_{i},\]

_so that \(p_{i}(\mathbf{x})=p_{i+1}(\mathbf{x})\) on the boundary of \(\Upsilon_{i}\) and \(\Upsilon_{i+1}\). Using the triangle inequality, Theorem 3.1 can be extended to the case that \(\boldsymbol{\beta}\) is piecewise constant to establish the error estimate on the whole domain \(\Omega\)._

**Theorem 3.3**: _Let \(u\) and \(u_{{}_{N}}\) be the solutions of problems (3) and (4), respectively. If the depth of ReLU NNs in (4) is at least \(\lceil\log_{2}(d+1)\rceil+1\), then for a sufficiently large integer \(n\)

Figure 1: Domain decomposition for the case that \(\boldsymbol{\beta}\) is piecewise constant
_there exists an integer \(\hat{n}\leq n\) such that_

\[\left|\!\left|\!\left|u-u_{{}_{N}}\right|\!\right|\!\right|\!\right|_{\boldsymbol{ \beta}}\leq C\,\left(\sqrt{\varepsilon_{2}}+\sqrt{\varepsilon_{3}}+\inf_{v \in\mathcal{M}(d,n-\hat{n})}\left|\!\left|\!\left|\hat{u}-v\right|\!\right|\! \right|_{\boldsymbol{\beta}}\right), \tag{3.4}\]

_where \(\mathcal{M}(d,n-\hat{n})=\mathcal{M}(d,1,\lceil\log_{2}(d+1)\rceil+1,n-\hat{n})\)._

_Proof._ The proof is similar to the one of Theorem 4.4 in [5]. 

**Lemma 3.4**: _Let \(u\), \(u_{{}_{N}}\), and \(u_{{}_{T}}^{N}\) be the solutions of problems (2.3), (2.4), and (2.5), respectively. Then there exist positive constants \(C_{1}\) and \(C_{2}\) such that_

\[\left|\!\left|\!\left|u-u_{{}_{T}}^{\omega}\right|\!\right|\! \right|\!\right|_{\boldsymbol{\beta}}\leq C_{1}\,\left(\left|(\mathcal{L}-\mathcal{L}_{{}_{T }})(u_{{}_{N}}-u_{{}_{T}}^{N},\mathbf{0})\right|+\left|(\mathcal{L}-\mathcal{ L}_{{}_{T}})(u-u_{{}_{N}},\mathbf{0})\right|\right)^{1/2}\] \[+C_{2}\,\left(\sqrt{\varepsilon_{2}}+\sqrt{\varepsilon_{3}}+\inf_ {v\in\mathcal{M}(d,n-\hat{n})}\left|\!\left|\!\left|\hat{u}-v\right|\!\right|\! \right|_{\boldsymbol{\beta}}\right). \tag{3.5}\]

_Proof._ The proof is similar to the one of Lemma 4.7 in [5]. 

## 4 Numerical Experiments

In this section, we present numerical results for two dimensional test problems with constant, piecewise constant, or variable advection velocity fields. The discrete LS functional was minimized by the ADAM optimization algorithm [9] on a uniform mesh with mesh size \(h=10^{-2}\). The directional derivative \(v_{\boldsymbol{\beta}}\) was approximated by the backward finite difference quotient multiplied by \(|\boldsymbol{\beta}|\)

\[v_{\boldsymbol{\beta}}(\mathbf{x}_{{}_{K}})\approx|\boldsymbol{\beta}|\frac{v( \mathbf{x}_{{}_{K}})-v\big{(}\mathbf{x}_{{}_{K}}-\rho\bar{\boldsymbol{\beta} }(\mathbf{x}_{{}_{K}})\big{)}}{\rho}, \tag{4.1}\]

where \(\bar{\boldsymbol{\beta}}=\frac{\boldsymbol{\beta}}{|\boldsymbol{\beta}|}\) and \(\rho=h/4\) (except for the last test problem, which used \(\rho=h/15\)). The LSNN method was implemented with an adaptive learning rate that started with \(0.004\) and was reduced by half for every \(50000\) iterations (except for the fourth test problem, which reduced for every \(100000\) iterations). For each experiment, to avoid local minima, \(10\) ReLU NN functions were trained for \(5000\) iterations each, and then the experiment began with one of the pretrained network functions that gave the minimum loss.

Tables 1 to 5 report the numerical errors in the relative \(L^{2}\), \(V_{\boldsymbol{\beta}}\), and the LS functional with parameters being the total number of weights and biases. Since the input dimension \(d=2\) (see [5]), we employed ReLU NN functions with a \(2\)-\(n_{1}\)-\(n_{2}\)-\(1\) representation or structure, which means the representation has two hidden layers with \(n_{1}\), \(n_{2}\) neurons, respectively. The \(l^{\text{th}}\) layer breaking lines defined by the set

\[\{\mathbf{x}\in\Omega:\boldsymbol{\omega}^{(l)}(N^{(l-1)}\circ\cdots\circ N^{( 2)}\circ N^{(1)}(\mathbf{x}))-\mathbf{b}^{(l)}\text{ has a zero component}\}\text{ with }N^{(0)}=I\]

are depicted in Figures 2 to 6 to follow the behavior of ReLU NN function approximation along the discontinous interface.

All of the test problems are defined on the domain \(\Omega=(0,1)^{2}\) with \(\gamma=1\) (\(f=1\) for the first three test problems and \(f=0\) for the remaining test problems).
respectively. The exact solution of this test problem is

\[u(x,y)=\left\{\begin{array}{rl}1,&(x,y)\in\Omega_{1}=\{(x,y)\in\Omega:x<1/2\}, \\ 1+e^{-y},&(x,y)\in\Omega\setminus\Omega_{1}.\end{array}\right. \tag{10}\]

The LSNN method with a random initialization and 50000 iterations was implemented with 2-20-20-1 ReLU NN functions. The numerical results are presented in Figure 2 and Table 1. The traces (Figure 2(b)) of the exact and numerical solutions on the plane \(y=0.5\) show no difference or oscillation. The exact solution (Figure 2(c)), which has a non-constant jump along the vertical interface (Figure 2(a)) is accurately approximated by a 3 layer ReLU NN function (Figure 2(d) and Table 1). Note that the solution of this test problem takes the same form as \(\chi\) in section 3, which was approximated by a CPWL function constructed by partitioning the domain into rectangles stacking on top of each other. It appears from Figure 2(e) that the 3 layer ReLU NN function approximation has a similar partition and the second layer breaking lines were generated for approximating the jump along the discontinuous interface by sharp transition layers and nonconstat part of the solution, which is consistent with our theoretical analysis on the convergence of the method.

### Problem with a piecewise smooth solution

This example is a modification of subsection 4.1 by changing the inflow boundary condition to

\[\Gamma_{-}=\{(x,0):x\in(0,1)\}\]

\[\text{and}\ \ g(x,y)=\left\{\begin{array}{rl}0,&(x,y)\in\Gamma_{-}^{1}\equiv\{(x,0):x\in(0,1/2)\},\\ 2,&(x,y)\in\Gamma_{-}^{2}=\Gamma_{-}\setminus\Gamma_{-}^{1},\end{array}\right.\]

respectively. The exact solution of this test problem is

\[u(x,y)=\left\{\begin{array}{rl}1-e^{-y},&(x,y)\in\Omega_{1},\\ 1+e^{-y},&(x,y)\in\Omega_{2}.\end{array}\right. \tag{11}\]

The LSNN method with a random initialization and 50000 iterations was implemented with 2-20-20-1 ReLU NN functions. The numerical results are presented in Figure 3 and Table 2. Unlike the previous test problem, the exact solution (Figure 3(c)) consists of two non-constant smooth parts. The LSNN method is capable of approximating the solution accurately without oscillation (Figures 3(b) to 3(d) and Table 2). The 3 layer ReLU NN function approximation has a partition (Figure 3(e)) similar to the one in subsection 4.1 with the second layer breaking lines on both sides for approximating the two non-constant smooth parts of the solution.
respectively. The exact solution of this test problem is

\[u(x,y)=\left\{\begin{array}{rl}1-\sin(2\pi x)e^{-y},&(x,y)\in\Omega_{1},\\ 1+(3/2-x)e^{-y},&(x,y)\in\Omega_{2}.\end{array}\right. \tag{10}\]

The LSNN method with a random initialization and 100000 iterations was implemented with 2-40-40-1 ReLU NN functions. The numerical results are presented in Figure 4 and Table 3. Since the solution on the inflow boundary consists of two non-constant smooth curves, we increased the number of hidden neurons to obtain a more accurate solution. Figures 4(c) and 4(d) and Table 3 show that the approximation is accurate pointwise and in average. The traces (Figure 4(b)) on \(y=0.5\) exhibit no oscillation and note a few corners on the curve, verifying that the ReLU

Figure 2: Approximation results of the problem in subsection 4.1
NN function approximation is a CPWL function. The partition generated by the breaking lines (Figure 4(e)) of the approximation shows how the exact solution was approximated.

### Problem with a piecewise constant advection velocity field

Let \(\bar{\Omega}=\bar{\Upsilon}_{1}\cup\bar{\Upsilon}_{2}\) and

\[\Upsilon_{1}=\{(x,y)\in\Omega:y\geq x\}\text{ and }\Upsilon_{2}=\Omega\setminus \Upsilon_{1}.\]

The advective velocity field is a piecewise constant field given by

\[\mathbf{\beta}(x,y)=\left\{\begin{array}{cl}(-1,\sqrt{2}-1)^{T},&(x,y)\in\Upsilon _{1},\\ (1-\sqrt{2},1)^{T},&(x,y)\in\Upsilon_{2}.\end{array}\right. \tag{4.6}\]

Figure 3: Approximation results of the problem in subsection 4.2
\begin{table}
\begin{tabular}{|l|c|c|c|c|} \hline Network structure & \(\frac{\|u-u_{T}^{N}\|_{0}}{\|u\|_{0}}\) & \(\frac{\|u-u_{T}^{N}\|_{\theta}}{\|u\|_{\theta}}\) & \(\frac{\mathcal{L}^{1/2}(u_{T}^{N},\mathbf{f})}{\mathcal{L}^{1/2}(u_{T}^{N}, \mathbf{0})}\) & Parameters \\ \hline
2-20-20-1 & 0.078036 & 0.013157 & 0.010386 & 501 \\ \hline \end{tabular}
\end{table}
Table 2: Relative errors of the problem in subsection 4.2

Figure 4: Approximation results of the problem in subsection 4.3

\begin{table}
\begin{tabular}{|l|c|c|c|c|} \hline Network structure & \(\frac{\|u-u_{T}^{N}\|_{0}}{\|u\|_{0}}\) & \(\frac{\|u-u_{T}^{N}\|_{\theta}}{\|u\|_{\theta}}\) & \(\frac{\mathcal{L}^{1/2}(u_{T}^{N},\mathbf{f})}{\mathcal{L}^{1/2}(u_{T}^{N}, \mathbf{0})}\) & Parameters \\ \hline
2-40-40-1 & 0.041491 & 0.016480 & 0.012733 & 1801 \\ \hline \end{tabular}
\end{table}
Table 3: Relative errors of the problem in subsection 4.3
The inflow boundary and the inflow boundary condition are given by

\[\Gamma_{-} =\{(1,y):y\in(0,1)\}\cup\{(x,0):x\in(0,1)\}\] \[\text{and}\ \ g(x,y) =\left\{\begin{array}{rl}\exp(x/(\sqrt{2}-1))x,&(x,y)\in\Gamma_{- }^{1}\equiv\{(x,0):x\in(0,1)\},\\ \exp(1/(\sqrt{2}-1))(11+(\sqrt{2}-1)y),&(x,y)\in\Gamma_{-}^{2}=\Gamma_{-} \setminus\Gamma_{-}^{1},\end{array}\right.\]

respectively. Let

\[\widehat{\Upsilon}_{11}=\{(x,y)\in\Upsilon_{1}:y<(1-\sqrt{2})x+1\},\ \widehat{\Upsilon}_{12}=\Upsilon_{1}\setminus\widehat{\Upsilon}_{11},\]

\[\widehat{\Upsilon}_{21}=\{(x,y)\in\Upsilon_{2}:y<\frac{1}{1-\sqrt{2}}(x-1)\},\ \text{and}\ \widehat{\Upsilon}_{22}=\Upsilon_{2}\setminus\widehat{\Upsilon}_{21}.\]

The exact solution of this test problem is

\[u(x,y)=\left\{\begin{array}{rl}\exp(\sqrt{2}x+y)(y+(\sqrt{2}-1)x),&(x,y)\in \widehat{\Upsilon}_{11},\\ \exp(\sqrt{2}x+y)(y+(\sqrt{2}-1)x+10),&(x,y)\in\widehat{\Upsilon}_{12},\\ \exp(x/(\sqrt{2}-1))(x+(\sqrt{2}-1)y),&(x,y)\in\widehat{\Upsilon}_{21},\\ \exp(x/(\sqrt{2}-1))(x+(\sqrt{2}-1)y+10),&(x,y)\in\widehat{\Upsilon}_{22}. \end{array}\right. \tag{4.7}\]

The LSNN method with a random initialization and 300000 iterations was implemented with 2-40-40-1 ReLU NN functions. The numerical results are presented in Figure 5 and Table 4. The exact solution (Figure 5(c)) consists of four non-constant smooth parts and has a non-constant jump along two connected line segments (Figure 5(a)). The traces (Figure 5(b)) of the exact and numerical solutions, the 3 layer ReLU NN function approximation (Figure 5(d)), and Table 4 indicate that the approximation is accurate pointwise and in average. Most of the second layer breaking lines (Figure 5(e)) are along the discontinuous interface, which correspond to sharp transition layers of the approximation for approximating the jump.
The LSNN method with a random initialization and 300000 iterations was implemented with 2-60-60-1 ReLU NN functions. The numerical results are presented in Figure 6 and Table 5. We increased the number of hidden neurons and \(\rho\) was set to \(h/15\) in the finite difference quotient in (4.1) because of the jump along the curved interface (Figure 6(a)). Although theoretical analysis on the convergence of the method in the case of a smooth interface was not conducted, Figures 6(b) to 6(d) and Table 5 show that the LSNN method is still capable of approximating the discontinuous solution with the curved interface accurately without oscillation. Finally, again, most of the second layer breaking lines are along the interface (Figure 6(e)) to approximate the discontinuous jump.

Figure 5: Approximation results of the problem in subsection 4.4
## 5 Proof of Theorem 3.1

In the section, we provide the proof of Theorem 3.1 by constructing a CPWL function to approximate \(\chi(\mathbf{x})\) in (3.1). Note that \(a(\mathbf{x})\) is generally a cylindrical surface and that the jump of \(\chi(\mathbf{x})\) is non-constant with

\[\chi(x,0)=a(x,0)=\alpha_{1}-\alpha_{2}\quad\forall\ x\in[0,x_{0}]\]

(see Figure 7(a)).

For a given \(\varepsilon_{1}>0\), take \(\widehat{\Upsilon}=(0,1)\times(y_{0},y_{0}+\varepsilon_{1})\) (see Figure 7(b)). Without loss of generality, let \(\alpha_{1}=1\), \(\alpha_{2}=0\), and \(\widehat{\Upsilon}=(0,1)\times(0,\varepsilon_{1})\).

Figure 6: Approximation results of the problem in subsection 4.5
Hence,

\[\chi(\mathbf{x})=\chi_{0}(\mathbf{x})+\chi_{1}(\mathbf{x})\text{ on }\widehat{\Upsilon}, \tag{5.1}\]

where \(\chi_{0}(\mathbf{x})\) is a step function and \(\chi_{1}(\mathbf{x})\) vanishes on the inflow boundary given by

\[\chi_{0}(\mathbf{x})=\left\{\begin{array}{ll}1,&\mathbf{x}\in\Omega_{1i} \cap\widehat{\Upsilon},\\ 0,&\mathbf{x}\in\Omega_{2i}\cap\widehat{\Upsilon},\end{array}\right.\quad \text{and}\quad\chi_{1}(\mathbf{x})=\left\{\begin{array}{ll}a(\mathbf{x})-1,&\mathbf{x}\in\Omega_{1i}\cap\widehat{\Upsilon},\\ 0,&\mathbf{x}\in\Omega_{2i}\cap\widehat{\Upsilon}.\end{array}\right.\]

**Lemma 5.1**: _Let_

\[b(\mathbf{x})=\left\{\begin{array}{ll}\boldsymbol{b}\cdot(\mathbf{x}-(x_{0},0)),&\mathbf{x}\in\Omega_{1i}\cap\widehat{\Upsilon},\\ 0,&\mathbf{x}\in\Omega_{2i}\cap\widehat{\Upsilon},\end{array}\right.\]

_where \(\boldsymbol{b}=(0,d)^{T}\) is a constant vector, and \(p_{1}(\mathbf{x})\) be a two layer neural network function on \(\widehat{\Upsilon}\), given by_

\[p_{1}(\mathbf{x})=-c\,\sigma(\mathbf{w}_{1}\cdot\mathbf{x}+x_{0})+c\,\sigma (\mathbf{w}_{2}\cdot\mathbf{x}+x_{0})\]

\begin{table}
\begin{tabular}{|l|c|c|c|c|} \hline Network structure & \(\frac{\|u-u_{T}^{N}\|_{0}}{\|u\|_{0}}\) & \(\frac{\|u-u_{T}^{N}\|_{\boldsymbol{\beta}}}{\|u\|_{\boldsymbol{\beta}}}\) & \(\frac{\mathcal{L}^{1/2}(u_{T}^{N},\mathbf{f})}{\mathcal{L}^{1/2}(u_{T}^{N}, \mathbf{0})}\) & Parameters \\ \hline
2-60-60-1 & 0.046528 & 0.049423 & 0.019995 & 3901 \\ \hline \end{tabular}
\end{table}
Table 5: Relative errors of the problem in subsection 4.5

Figure 7: Illustration of the convergence analysis on one subdomain \(\Upsilon_{i}\)
with weights and coefficient_

\[\mathbf{w}_{1}=\begin{pmatrix}-1\\ -\varepsilon_{2}\end{pmatrix},\quad\mathbf{w}_{2}=\begin{pmatrix}-1\\ \varepsilon_{2}\end{pmatrix},\quad c=\frac{d}{2\varepsilon_{2}}\]

_(see Figure 7(c)). Then we have on \(\widehat{\Upsilon}\),_

\[\left|\!\left|\!\left|b-p_{1}\right|\!\right|\!\right|_{\boldsymbol{\beta}}= \left(\left|\!\left|b-p_{1}\right|\!\right|\!\right|_{0,\widehat{\Upsilon}}^{2 }+\left|\!\left|b_{\boldsymbol{\beta}}-p_{1\boldsymbol{\beta}}\right|\!\right| \!\right|_{0,\widehat{\Upsilon}}^{2}\right)^{1/2}\leq\sqrt{\frac{\varepsilon_{ 1}^{3}}{24}+\frac{B^{2}}{4}}\left|d\right|\sqrt{\varepsilon_{1}\varepsilon_{2}}, \tag{5.2}\]

_where we assume \(\left|v_{2}(\mathbf{x})\right|\leq B\) (recall \(\boldsymbol{\beta}(\mathbf{x})=(0,v_{2}(\mathbf{x}))\))._

Proof.: Set

\[\widehat{\Upsilon}_{\varepsilon_{2}}\equiv\widehat{\Upsilon}_{1,\varepsilon_ {2}}\cup\widehat{\Upsilon}_{2,\varepsilon_{2}}\equiv\{\mathbf{x}\in\Omega_{1 i}\cap\widehat{\Upsilon}:\,\mathbf{w}_{1}\cdot\mathbf{x}+x_{0}<0\}\cup\{ \mathbf{x}\in\Omega_{2i}\cap\widehat{\Upsilon}:\,\mathbf{w}_{2}\cdot\mathbf{x }+x_{0}>0\}.\]

Then we have

\[b(\mathbf{x})-p_{1}(\mathbf{x})=\left\{\begin{array}{ll}-c(\mathbf{w}_{1} \cdot\mathbf{x}+x_{0})&\mathbf{x}\in\widehat{\Upsilon}_{1,\varepsilon_{2}}, \\ -c(\mathbf{w}_{2}\cdot\mathbf{x}+x_{0})&\mathbf{x}\in\widehat{\Upsilon}_{2, \varepsilon_{2}},\\ 0,&\mathbf{x}\in\widehat{\Upsilon}\setminus\widehat{\Upsilon}_{\varepsilon_ {2}}.\end{array}\right.\]

Calculating the double integral,

\[\left|\!\left|b-p_{1}\right|\!\right|_{0,\widehat{\Upsilon}}^{2}=\left|\! \left|b-p_{1}\right|\!\right|_{0,\widehat{\Upsilon}_{1,\varepsilon_{2}}}^{2} +\left|\!\left|b-p_{1}\right|\!\right|_{0,\widehat{\Upsilon}_{2,\varepsilon_ {2}}}^{2}=\frac{d^{2}}{24}\varepsilon_{1}^{4}\varepsilon_{2}, \tag{5.3}\]

and using the directional derivative,

\[\left|\!\left|b_{\boldsymbol{\beta}}-p_{1\boldsymbol{\beta}}\right|\!\right|_ {0,\widehat{\Upsilon}}^{2}=\int_{\widehat{\Upsilon}_{1,\varepsilon_{2}}}(c \mathbf{w}_{1}\cdot\boldsymbol{\beta})^{2}\,d\mathbf{x}+\int_{\widehat{ \Upsilon}_{2,\varepsilon_{2}}}(c\mathbf{w}_{2}\cdot\boldsymbol{\beta})^{2} \,d\mathbf{x}\leq\frac{(dB)^{2}}{4}\varepsilon_{1}\varepsilon_{2}. \tag{5.4}\]

Now (5.2) follows from (5.3) and (5.4). 

**Lemma 5.2**.: _Let \(\widehat{\Upsilon}\), \(I\), \(b(\mathbf{x})\), \(p_{1}(\mathbf{x})\), and \(\boldsymbol{\beta}(\mathbf{x})\) be as in Lemma 5.1 with \(d=\chi_{1}(0,\varepsilon_{1})/\varepsilon_{1}\), and let \(p_{0}(\mathbf{x})\) be a two layer neural network function on \(\widehat{\Upsilon}\) given by_

\[p_{0}(\mathbf{x})=\frac{1}{2\varepsilon_{2}}\left(\sigma(x-x_{0}+\varepsilon_ {2})-\sigma(x-x_{0}-\varepsilon_{2})\right). \tag{5.5}\]

_Then we have on \(\widehat{\Upsilon}\),_

\[\left|\!\left|\!\left|\!\left|\chi-(p_{0}+p_{1})\right|\!\right|\!\right|\! \right|_{\boldsymbol{\beta}}\leq C_{1}\sqrt{\varepsilon_{1}\varepsilon_{2}}+ C_{2}\sqrt{\varepsilon_{1}}, \tag{5.6}\]

_where \(C_{2}\) is given by the square root of_

\[\sup\{(u_{1}(\mathbf{x})-u_{2}(\mathbf{x})-1-b(\mathbf{x}))^{2}: \mathbf{x}\in\Omega_{1i}\cap\widehat{\Upsilon}\}x_{0}\\ +\sup\{2\gamma(\mathbf{x})^{2}(u_{2}(\mathbf{x})-u_{1}(\mathbf{x }))^{2}+2(dB)^{2}:\mathbf{x}\in\Omega_{1i}\cap\widehat{\Upsilon}\}x_{0} \tag{5.7}\]

Proof.: From the triangle inequality,

\[\left|\!\left|\!\left|\chi-(p_{0}+p_{1})\right|\!\right|\!\right|_{\boldsymbol{ \beta}}=\left|\!\left|\!\left|\chi_{0}+\chi_{1}-(p_{0}+p_{1})\right|\!\right|\! \right|_{\boldsymbol{\beta}}\leq\left|\!\left|\!\left|\chi_{0}-p_{0}\right|\! \right|\!\right|_{\boldsymbol{\beta}}+\left|\!\left|\!\left|\chi_{1}-p_{1} \right|\!\right|\!\right|_{\boldsymbol{\beta}}. \tag{5.8}\]

Since \(\left|\!\left|\chi_{0\boldsymbol{\beta}}-p_{0\boldsymbol{\beta}}\right|\!\right|=0\), calculating the double integral,

\[\left|\!\left|\!\left|\chi_{0}-p_{0}\right|\!\right|\!\right|_{ \boldsymbol{\beta}} =\left(\left|\!\left|\chi_{0}-p_{0}\right|\!\right|\!\right|_{0, \widehat{\Upsilon}}^{2}+\left|\!\left|\chi_{0\boldsymbol{\beta}}-p_{0 \boldsymbol{\beta}}\right|\!\right|_{0,\widehat{\Upsilon}}^{2}\right)^{1/2}\] \[=\left|\!\left|\chi_{0}-p_{0}\right|\!\right|_{0,\widehat{\Upsilon}}\] \[=\frac{1}{\sqrt{6}}\sqrt{\varepsilon_{1}\varepsilon_{2}}. \tag{5.9}\]
Next, again, by the triangle inequality,

\[\left|\!\left|\!\left|\chi_{1}-p_{1}\right|\!\right|\!\right|_{\boldsymbol{\beta}} \leq\left|\!\left|\!\left|\chi_{1}-b\right|\!\right|\!\right|_{\boldsymbol{ \beta}}+\left|\!\left|\!\left|b-p_{1}\right|\!\right|\!\right|_{\boldsymbol{ \beta}}.\]

By Lemma 5.1

\[\left|\!\left|\!\left|\!\left|b-p_{1}\right|\!\right|\!\right|\!\right|_{ \boldsymbol{\beta}}\leq\sqrt{\frac{\varepsilon_{1}^{3}}{24}+\frac{B^{2}}{4}} \left|d\right|\sqrt{\varepsilon_{1}\varepsilon_{2}}.\]

To bound \(\left|\!\left|\!\left|\!\left|\chi_{1}-b\right|\!\right|\!\right|_{\boldsymbol{ \beta}}\), recall the definition of the graph norm,

\[\left|\!\left|\!\left|\chi_{1}-b\right|\!\right|\!\right|_{\boldsymbol{\beta} }=\left(\left|\!\left|\chi_{1}-b\right|\!\right|\!\right|_{0,\widehat{\Upsilon }}^{2}+\left|\!\left|\chi_{1\boldsymbol{\beta}}-b_{\boldsymbol{\beta}}\right| \!\right|_{0,\widehat{\Upsilon}}^{2}\right)^{1/2}.\]

First we have

\[\left|\!\left|\chi_{1}-b\right|\!\right|_{0,\widehat{\Upsilon}}^ {2} =\left|\!\left|\chi_{1}-b\right|\!\right|_{0,\Omega_{1i}\cap \widehat{\Upsilon}}^{2}\] \[\leq\sup\{(\chi_{1}(\mathbf{x})-b(\mathbf{x}))^{2}:\mathbf{x}\in \Omega_{1i}\cap\widehat{\Upsilon}\}\varepsilon_{1}x_{0}\] \[=\sup\{(u_{1}(\mathbf{x})-u_{2}(\mathbf{x})-1-b(\mathbf{x}))^{2}: \mathbf{x}\in\Omega_{1i}\cap\widehat{\Upsilon}\}\varepsilon_{1}x_{0}.\]

Next, observing \(b_{\boldsymbol{\beta}}=dv_{2}\) and from (2),

\[\chi_{1\boldsymbol{\beta}}=(u_{1}-u_{2}-1)_{\boldsymbol{\beta}}=(u_{1}-u_{2}) _{\boldsymbol{\beta}}=\gamma(u_{2}-u_{1}),\]

we have

\[\left|\!\left|\chi_{1\boldsymbol{\beta}}-b_{\boldsymbol{\beta}} \right|\!\right|_{0,\widehat{\Upsilon}}^{2} =\left|\!\left|\chi_{1\boldsymbol{\beta}}-b_{\boldsymbol{\beta}} \right|\!\right|_{0,\Omega_{1i}\cap\widehat{\Upsilon}}^{2}\] \[=\left|\!\left|\gamma(u_{2}-u_{1})-dv_{2}\right|\!\right|_{0, \Omega_{1i}\cap\widehat{\Upsilon}}^{2}\] \[\leq\left\|\!\left|\sqrt{2\gamma^{2}(u_{2}-u_{1})^{2}+2(dv_{2})^ {2}}\right\|\!\right|_{0,\Omega_{1i}\cap\widehat{\Upsilon}}^{2}\] \[\leq\sup\{2\gamma(\mathbf{x})^{2}(u_{2}(\mathbf{x})-u_{1}(\mathbf{ x}))^{2}+2(dB)^{2}:\mathbf{x}\in\Omega_{1i}\cap\widehat{\Upsilon}\} \varepsilon_{1}x_{0}.\]

Combining the above inequalities, we obtain (10). 

Given \(\varepsilon_{3}>0\), choose \(\varepsilon_{1}=1/m\) such that

\[\sup\{(u_{1}(\mathbf{x})-u_{2}(\mathbf{x})-\chi(0,j/m)-b_{j}( \mathbf{x}))^{2}:\mathbf{x}\in\Omega_{1i}\cap\widehat{\Upsilon}_{j}\}, \tag{12}\] \[\sup\{2\gamma(\mathbf{x})^{2}(u_{2}(\mathbf{x})-u_{1}(\mathbf{x}) )^{2}+2(d_{j}B)^{2}:\mathbf{x}\in\Omega_{1i}\cap\widehat{\Upsilon}_{j}\}< \varepsilon_{3},\]

where \(\widehat{\Upsilon}_{j}=(0,1)\times(j/m,(j+1)/m)\) for \(j=0,\ldots,m-1\), and \(b_{j}(\mathbf{x})\), \(d_{j}\) as in Lemma 5.2. Then define on each \(\widehat{\Upsilon}_{j}\), \(p_{0j}(\mathbf{x}),p_{1j}(\mathbf{x})\) as in Lemma 5.2, and construct a CPWL function \(p_{i}(\mathbf{x})\) on \(\Upsilon_{i}\) given by

\[p_{i}(\mathbf{x})=p_{0j}(\mathbf{x})+p_{1j}(\mathbf{x}),\;\mathbf{x}\in \widehat{\Upsilon}_{j}.\]
Proof of Theorem 3.1.: By Lemma 5.2 and the given condition,

\[\begin{split}\left\|\chi-p_{i}\right\|_{\boldsymbol{\beta}}& =\left(\sum_{j=1}^{m}\left\|\chi-(p_{0j}+p_{1j})\right\|_{ \boldsymbol{\beta}}^{2}\right)^{1/2}\\ &\leq\left(\sum_{j=1}^{m}\left(D_{1j}\sqrt{\varepsilon_{1} \varepsilon_{2}}+D_{2j}\sqrt{\varepsilon_{1}\varepsilon_{3}}\right)^{2} \right)^{1/2}\\ &\leq\left(m\max_{1\leq j\leq m}\left(D_{1j}\sqrt{\varepsilon_{1} \varepsilon_{2}}+D_{2j}\sqrt{\varepsilon_{1}\varepsilon_{3}}\right)^{2} \right)^{1/2}\\ &=\sqrt{m}\max_{1\leq j\leq m}\left(D_{1j}\sqrt{\varepsilon_{1} \varepsilon_{2}}+D_{2j}\sqrt{\varepsilon_{1}\varepsilon_{3}}\right)\\ &=\max_{1\leq j\leq m}\left(D_{1j}\sqrt{\varepsilon_{2}}+D_{2j} \sqrt{\varepsilon_{3}}\right),\end{split} \tag{5.11}\]

where for the first identity, each norm on the right-hand side is taken over \(\widehat{\Upsilon}_{j}\). Now \(D_{1}=D_{1k}\) and \(D_{2}=D_{2k}\) for some \(1\leq k\leq m\). 
