# Adaptive Learning on the Grids for Elliptic Hemivariational Inequalities

Jianguo Huang<EMAIL_ADDRESS>
School of Mathematical Sciences, and MOE-LSC,
Shanghai Jiao Tong University, Shanghai 200240, China
Chunmei Wang<EMAIL_ADDRESS>
Department of Mathematics & Statistics,
Texas Tech University, 1108 Memorial Circle, Lubbock, TX 79409, USA
Haoqin Wang<EMAIL_ADDRESS>
School of Mathematical Sciences, and MOE-LSC,
Shanghai Jiao Tong University, Shanghai 200240, China

###### Abstract

This paper introduces a deep learning method for solving an elliptic
hemivariational inequality (HVI). In this method, an expectation minimization
problem is first formulated based on the variational principle of underlying
HVI, which is solved by stochastic optimization algorithms using three
different training strategies for updating network parameters. The method is
applied to solve two practical problems in contact mechanics, one of which is
a frictional bilateral contact problem and the other of which is a
frictionless normal compliance contact problem. Numerical results show that
the deep learning method is efficient in solving HVIs and the adaptive mesh-
free multigrid algorithm can provide the most accurate solution among the
three learning methods.

Keywords. Deep Learning; Elliptic Hemivariational Inequalities; Contact
Problems; Mesh-Free Methods; Multigrid.

##  1 Introduction

This paper considers an efficient numerical method for solving an elliptic
hemivariational inequality (HVI), which was first introduced by
Panagiotopoulos in the context of engineering in the early 1980s [50]. After
Panagiotopoulos’ pioneering work, a lot of mathematical effort has been
devoted to the study of HVIs (e.g., [51, 49, 48, 44, 22]), and HVIs have
become a powerful mathematical tool in contact mechanics. In practice, the
solution to elliptic HVIs is only available numerically. There have been
various numerical methods to approximate the solution of elliptic HVIs (e.g.,
see [26, 23, 24, 19, 59]). After discretization, we often need to solve a non-
convex and non-smooth optimization problem, which is challenging and technical
to attack. One of the most popular methods for this non-convex optimization is
the iterative convexification approach [58, 45], where a sequence of convex
problems approaching the original non-convex problem is constructed and
solved. Based on this iterative convexification approach, many HVIs describing
contact problems were solved, such as the frictional bilateral contact problem
[3, 25], the frictionless normal compliance contact problem [4, 25], and the
frictionless unilateral contact problem [25]. Except for the iterative
convexification approach, other alternative choices are the proximal bundle
method (e.g., see [42]), the bundle Newton method (e.g., see [2]), and the
primal-dual active-set algorithm [36]. Recently, Feng, Han and Huang [19] used
the double bundle method (cf. [33]) to solve the discrete non-convex and non-
smooth problem arising from discretization of some HVIs. As a whole, the
existing numerical approaches mentioned above are rather involved to implement
in practice.

With the advance of deep learning techniques originated in computer science,
considerable attention in applied mathematics and computational mathematics
has been drawn to apply deep learning to scientific computing, especially for
numerical methods to solve the partial differential equations (PDEs). Neural
network-based numerical methods for solving PDEs date back to the 1990s [39,
13, 38] and achieve significant improvement recently (e.g., see [43, 54, 53,
57, 21, 34]). In those methods, deep neural networks (DNNs) are applied to
parametrize PDE solutions and appropriate parameters are identified by
minimizing an optimization problem formulated from the given PDE. The key to
the success of neural networks-based methods is the universal approximation
property of DNNs (cf. [29, 12, 28, 5, 37, 61, 62, 16, 55, 46, 47, 41, 17,
56]). Though deep learning has made great achievements in solving PDEs,
especially in high-dimensional cases [21], it is also challenging to obtain a
highly accurate solution. Recently, several deep learning frameworks have been
developed along this line. For instance, E and Yu proposed the Deep Ritz
method for high-dimensional variational problems [18]. Karniadakis and his
team designed several neural networks by incorporating physical information
[52, 32, 31]. Gu, Yang and Zhou developed the SelectNet to adaptively choose
training samples in the learning propocess [20]. Liang et al. discussed how to
design data-driven activation functions in [40]. Bao and his team introduced
weak adversarial networks to find the weak solution of PDEs in [63]. Chen et
al. proposed Friedrichs learning strategies to solve various PDEs in a unified
way [9]. Huang, Wang and Yang combined deep learning with traditional
iteration to devise the Int-Deep for solving low-dimensional PDEs with a
finite element accuracy [30]. Dong and Li combined the ideas of domain
decomposition and extreme learning to form a new deep learning method for PDEs
[14].

Though deep learning-based methods have been proposed to handle variational
inequalities in [30], to the best of our knowledge, there is no study on HVIs
in the literature. In this paper, we propose a deep learning-based method to
solve an elliptic HVI based on its equivalent variational form [22] and
compare the numerical performance of three different training strategies for
updating parameters. In our method, the solution space of the HVI is
parameterized via DNNs and an approximation is found by minimizing an
unconstrained expectation minimization problem, which can be solved by
stochastic gradient descent methods or its variants (cf. [6]). In particular,
the unconstrained expectation minimization problem is reformulated based on
the variational principle of the HVI. Therefore, the resulting deep learning
optimization problem has a clear physical meaning. In the meantime, we compare
three algorithms for training networks: the classical stochastic algorithm
with basic neural network parametrization, a blockwise training algorithm for
multi-block neural networks, and an adaptive mesh-free multigrid algorithm. As
applications, we employ the deep learning method to a frictional bilateral
contact problem and a frictionless contact problem with normal compliance. As
we shall see in the numerical experiments, the deep learning method is
efficient in solving HVIs and the adaptive mesh-free multigrid algorithm can
provide the most accurate solution among the three learning methods discussed.
Finally, it deserves to emphasize that our method is easy to be implemented in
programming and suitable for practical applications.

The rest of this paper is organized as follows. In Section 2, an elliptic HVI
and its applications in contact mechanics are introduced briefly. In Section
3, the deep learning method as well as three training algorithms for solving
the HVI are presented in details. In Section 4, two numerical examples are
provided to demonstrate the efficiency of the deep learning methods. Finally,
we summarize our work with a short conclusion in Section 5.

##  2 The elliptic hemivariational inequality with applications

In this section, we first introduce some notations. For a real Banach space
$X$ equipped with a norm $\|\cdot\|_{X}$, denote by $X^{*}$ its dual space
equipped with a norm $\|\cdot\|_{X^{*}}$. The notation
$\langle\cdot,\cdot\rangle$ stands for the dual pairing between $X^{*}$ and
$X$. These notations apply to a Hilbert space $H$ naturally. Denote by
$\mathcal{L}(X,Y)$ the space of all continuous linear operators from one
normed linear space $X$ to another normed linear space $Y$. We also use the
standard notations for Sobolev spaces and norms or seminorms (cf. [1]). For
any locally Lipschitz continuous functional $j$ on a Banach space $X$, denote
by $j^{0}(u;w)$ the generalized (Clarke) directional derivative of $j$ at $u$
in the direction $w$ (cf. [10, 11]); i.e.,

$j^{0}(u;w)=\lim_{v\to u}\sup_{t\downarrow 0}\frac{j(v+tw)-j(v)}{t},\quad u\in
X,w\in X.$

###  2.1 The elliptic hemivariational inequality and its equivalent
minimization problem

In this paper, we consider an elliptic HVI on a spatial domain $\Omega$ in a
finite-dimensional Euclidean space. For simplicity, we denote the boundary or
part of the boundary of the domain as $\Gamma$. Let $H$ be a Hilbert space.
The elliptic HVI can be described as follows.

Problem 1 Find $u\in H$ such that

$\langle Au,v\rangle+\int_{\Gamma}j^{0}(\gamma_{j}u;\gamma_{j}v)\,{\rm
d}s\geq\langle f,v\rangle,\quad v\in H,$ (2.1)

where $j:\Gamma\times\mathbb{R}^{m}\to\mathbb{R}$ is a locally Lipschitz
continuous functional for some positive integer $m$,
$\gamma_{j}\in\mathcal{L}(H,L^{2}(\Gamma;\mathbb{R}^{m}))$. In the study of
Problem 1, we need the following assumptions [22]:

$\displaystyle(A1)\;A:H\to H^{*}\ \text{is Lipschitz continuous and strongly
monotone; i.e., for a constant}\ m_{A}>0,$
$\displaystyle\qquad\qquad\qquad\qquad\qquad\langle
Av_{1}-Av_{2},v_{1}-v_{2}\rangle\geq m_{A}\|v_{1}-v_{2}\|_{H}^{2},\quad
v_{1},v_{2}\in H.$ $\displaystyle(A2)\;j(\cdot,\boldsymbol{z})\ \text{is
measurable on}\ \Gamma,\ \boldsymbol{z}\in\mathbb{R}^{m};\text{there exists}\
\boldsymbol{z}\in L^{2}(\Gamma;\mathbb{R}^{m})\ \text{such that}\
j(\cdot,\boldsymbol{z}(\cdot))\in L^{1}(\Gamma);$
$\displaystyle\qquad\text{there exist non-negative constants}\ c_{0},c_{1},\
\text{and}\ \alpha_{j}\ \text{such that}$
$\displaystyle\qquad\qquad\qquad\qquad\qquad\|\partial
j(\boldsymbol{z})\|_{\mathbb{R}^{m}}\leq
c_{0}+c_{1}\|\boldsymbol{z}\|_{\mathbb{R}^{m}},\quad\boldsymbol{z}\in\mathbb{R}^{m},$
$\displaystyle\qquad\qquad\qquad\qquad\qquad
j^{0}(\boldsymbol{z}_{1};\boldsymbol{z}_{2}-\boldsymbol{z}_{1})+j^{0}(\boldsymbol{z}_{2};\boldsymbol{z}_{1}-\boldsymbol{z}_{2})\leq\alpha_{j}\|\boldsymbol{z}_{1}-\boldsymbol{z}_{2}\|_{\mathbb{R}^{m}}^{2},\quad\boldsymbol{z}_{1},\boldsymbol{z}_{2}\in\mathbb{R}^{m}.$
$\displaystyle(A3)\;f\in H^{*}.$ $\displaystyle(A4)\;\text{Denote by}\
c_{\Gamma}\ \text{an upper bound of the norm of the operator}\ \gamma_{j}.\
\text{There holds}$
$\displaystyle\qquad\qquad\qquad\qquad\qquad\|\gamma_{j}v\|_{L^{2}(\Gamma;\mathbb{R}^{m})}\leq
c_{\Gamma}\|v\|_{H},\quad v\in H.$

We recall an important result on the solution existence and uniqueness for
Problem 1 (cf. [25]).

######  Theorem 2.1.

Assume (A1)-(A4) and $\alpha_{j}c_{\Gamma}^{2}<m_{A}$ hold true. For any given
$f\in H^{*}$, Problem 1 has a unique solution.

Next, assume $A\in\mathcal{L}(H,H^{*})$ is symmetric, i.e.,

$\langle Av_{1},v_{2}\rangle=\langle Av_{2},v_{1}\rangle,\quad v_{1},v_{2}\in
H.$

We turn to consider the following minimization problem:

Problem 2 Find $u\in H$ such that

$u=\operatorname*{arg\,min}_{v\in H}E(v),$ (2.2)

where

$E(v)=\frac{1}{2}\langle Av,v\rangle+\int_{\Gamma}j(\gamma_{j}v)\,{\rm
d}s-\langle f,v\rangle,\quad v\in H.$

As given in [22], we have the following equivalence result.

######  Theorem 2.2.

Assume (A1)-(A4), $\alpha_{j}c_{\Gamma}^{2}<m_{A}$, and
$A\in\mathcal{L}(H,H^{*})$ is symmetric. Then Problem 2 is equivalent to
Problem 1.

Note that the functional $E(v)$ in (2.2) is reformulated from (2.1) based on
the variational principle. Hence, it has clear physical meaning in mechanics.

###  2.2 Some applications in contact mechanics

Let $\Omega\subset\mathbb{R}^{2}$ be the reference configuration of a linear
elastic body. We assume that $\Omega$ is an open, bounded, and connected
domain with Lipschitz continuous boundary $\Gamma=\partial\Omega$. The
boundary is made of three disjoint and measurable parts: $\Gamma_{D}$,
$\Gamma_{T}$, and $\Gamma_{C}$ such that $\text{meas}(\Gamma_{D})>0$ and
$\text{meas}(\Gamma_{C})>0$, Denote by $\cdot$ and $|\cdot|$ the canonical
inner product and the induced norm, respectively. For a vector field
$\boldsymbol{v}\in\mathbb{R}^{2}$, we use
$v_{\nu}=\boldsymbol{v}\cdot\boldsymbol{\nu}$ for its normal component and
$\boldsymbol{v}_{\tau}=\boldsymbol{v}-v_{\nu}\boldsymbol{\nu}$ for its
tangential component, where $\boldsymbol{\nu}$ is the unit outward normal
vector to $\Gamma$. The linearized strain tensor associated with a
displacement field $\boldsymbol{u}:\Omega\to\mathbb{R}^{2}$ is denoted by
$\boldsymbol{\varepsilon}(\boldsymbol{u})$ and the stress field is denoted by
$\boldsymbol{\sigma}:\Omega\to\mathbb{S}^{2}$. In addition, we assume a volume
force of density $\boldsymbol{f}_{0}\in L^{2}(\Omega,\mathbb{R}^{2})$ acting
in $\Omega$. Besides, the body is assumed to be fixed on $\Gamma_{D}$, is
subject to an action of the surface traction of density $\boldsymbol{f}_{2}\in
L^{2}(\Gamma_{T},\mathbb{R}^{2})$ on $\Gamma_{T}$, and is in contact on
$\Gamma_{C}$. We usually drop the spatial variable $\boldsymbol{x}$ for
simplicity when the dependence is clear without any confusion.

To discuss the contact problem, we introduce a space

$V=\\{\boldsymbol{v}\in
H^{1}(\Omega;\mathbb{R}^{2})|\boldsymbol{v}=\boldsymbol{0}\ \text{a.e. on}\
\Gamma_{D}\\}$

equipped with the inner product

$(\boldsymbol{u},\boldsymbol{v})_{V}=\int_{\Omega}\boldsymbol{\varepsilon}(\boldsymbol{u})\cdot\boldsymbol{\varepsilon}(\boldsymbol{v})\,{\rm
d}x,\quad\boldsymbol{u},\boldsymbol{v}\in V,$

and the associated norm
$\|\boldsymbol{v}\|_{V}=\sqrt{(\boldsymbol{v},\boldsymbol{v})_{V}}$. Thanks to
$\text{meas}(\Gamma_{D})>0$ and Korn’s inequality [7], we know $V$ is a
Hilbert space with the norm $\|\cdot\|_{V}$. Besides, we need to introduce
another Hilbert space $Q=L^{2}(\Omega;\mathbb{S}^{2})$, which is equipped with
the inner product

$(\boldsymbol{\sigma},\boldsymbol{\tau})_{Q}=\int_{\Omega}{\sigma}_{ij}(\boldsymbol{x}){\tau}_{ij}(\boldsymbol{x})\,{\rm
d}x.$

Here and below, we use the Einstein summation convention, which means the
summation is implied for an index exactly appeared two times in a quantity.

First of all, we introduce an HVI to describe the frictional bilateral contact
problem (cf. [19, 25]). We let

$V_{1}=\\{\boldsymbol{v}\in V|\boldsymbol{v}_{\nu}=0\ \text{on}\
\Gamma_{C}\\},$

and define

$\displaystyle\langle
A\boldsymbol{u},\boldsymbol{v}\rangle=(\boldsymbol{\mathcal{F}}(\boldsymbol{\varepsilon}(\boldsymbol{u})),\boldsymbol{\varepsilon}(\boldsymbol{v}))_{Q};$
$\displaystyle\gamma_{j_{\tau}}:V_{1}\to L^{2}(\Gamma_{C};\mathbb{R}^{2})\
\text{such that}\ \gamma_{j_{\tau}}(\boldsymbol{v})=\boldsymbol{v}_{\tau};$
$\displaystyle\int_{\Gamma}j^{0}(\gamma_{j}\boldsymbol{u};\gamma_{j}\boldsymbol{v})\,{\rm
d}s=\int_{\Gamma_{C}}j_{\tau}^{0}(\gamma_{j_{\tau}}\boldsymbol{u};\gamma_{j_{\tau}}\boldsymbol{v})\,{\rm
d}s;$
$\displaystyle\langle\boldsymbol{f},\boldsymbol{v}\rangle=\int_{\Omega}\boldsymbol{f}_{0}\cdot\boldsymbol{v}\,{\rm
d}x+\int_{\Gamma_{T}}\boldsymbol{f}_{2}\cdot\boldsymbol{v}\,{\rm d}s,\
\boldsymbol{v}\in V;$

where the linear elasticity operator
$\boldsymbol{\mathcal{F}}=(F_{ijkl})_{1\leq i,j,k,l\leq
2}:\Omega\times\mathbb{S}^{2}\to\mathbb{S}^{2}$ is symmetric, bounded, and
satisfies the following property [22]

$\boldsymbol{\mathcal{F}}(\boldsymbol{\sigma})\cdot\boldsymbol{\sigma}\geq
m_{\boldsymbol{\mathcal{F}}}|\boldsymbol{\sigma}|^{2},\quad
m_{\boldsymbol{\mathcal{F}}}>0,\ \boldsymbol{\sigma}\in\mathbb{S}^{2},$ (2.3)

$j_{\tau}:\Gamma_{C}\times\mathbb{R}^{2}\to\mathbb{R}$ is locally Lipschitz on
$\mathbb{R}^{2}$ for a.e. $\boldsymbol{x}\in\Gamma_{C}$ and satisfies the
assumption $(A2)$ with constants $c_{0}$, $c_{1}$ and $\alpha_{j_{\tau}}$(cf.
[19, 25]). Then the frictional bilateral contact problem can be described as a
HVI in the form

$(\boldsymbol{\mathcal{F}}(\boldsymbol{\varepsilon}(\boldsymbol{u})),\boldsymbol{\varepsilon}(\boldsymbol{v}))_{Q}+\int_{\Gamma_{C}}j_{\tau}^{0}(\boldsymbol{u}_{\tau};\boldsymbol{v}_{\tau})\,{\rm
d}s\geq\langle\boldsymbol{f},\boldsymbol{v}\rangle,\quad\boldsymbol{v}\in
V_{1}.$ (2.4)

Let $\boldsymbol{u}\in V_{1}$ and $\lambda_{1}$ be the smallest positive
eigenvalue of the eigenvalue problem

$\int_{\Omega}\boldsymbol{\varepsilon}(\boldsymbol{u})\cdot\boldsymbol{\varepsilon}(\boldsymbol{v})\,{\rm
d}x=\lambda_{1}\int_{\Gamma_{C}}\boldsymbol{u}_{\tau}\cdot\boldsymbol{v}_{\tau}\,{\rm
d}s,\quad\boldsymbol{v}\in V_{1}.$

Then the assumption $(A4)$ is satisfied if
$c_{\Gamma}\geq\sqrt{1/\lambda_{1}}$ holds true.

The assumption $(A1)$ is satisfied with $m_{A}=m_{\boldsymbol{\mathcal{F}}}$.
Following Theorem 2.1 and Theorem 2.2, if
$\alpha_{j_{\tau}}<\lambda_{1}m_{\boldsymbol{\mathcal{F}}}$, (2.4) has a
unique solution and is equivalent to the following optimization problem:

$\boldsymbol{u}=\operatorname*{arg\,min}_{\boldsymbol{v}\in
V_{1}}E(\boldsymbol{v}),$ (2.5)

where

$E(\boldsymbol{v})=\frac{1}{2}(\boldsymbol{\mathcal{F}}(\boldsymbol{\varepsilon}(\boldsymbol{v})),\boldsymbol{\varepsilon}(\boldsymbol{v}))_{Q}+\int_{\Gamma_{C}}j_{\tau}(\boldsymbol{v}_{\tau})\,{\rm
d}s-\langle\boldsymbol{f},\boldsymbol{v}\rangle.$

Next, we introduce an HVI to describe the frictionless normal compliance
contact problem (cf. [19, 25]). Define

$\displaystyle\langle
A{\boldsymbol{u}},{\boldsymbol{v}}\rangle=(\boldsymbol{\mathcal{F}}(\boldsymbol{\varepsilon}(\boldsymbol{u})),\boldsymbol{\varepsilon}(\boldsymbol{v}))_{Q};$
$\displaystyle\gamma_{j_{\nu}}:V\to L^{2}(\Gamma_{C})\ \text{such that}\
\gamma_{j_{\nu}}(\boldsymbol{v})=\boldsymbol{v}_{\nu};$
$\displaystyle\int_{\Gamma}j^{0}(\gamma_{j}\boldsymbol{u};\gamma_{j}\boldsymbol{v})\,{\rm
d}s=\int_{\Gamma_{C}}j_{\nu}^{0}(\gamma_{j_{\nu}}\boldsymbol{u};\gamma_{j_{\nu}}\boldsymbol{v})\,{\rm
d}s;$
$\displaystyle\langle{\boldsymbol{f}},\boldsymbol{v}\rangle=\int_{\Omega}\boldsymbol{f}_{0}\cdot\boldsymbol{v}\,{\rm
d}x+\int_{\Gamma_{T}}\boldsymbol{f}_{2}\cdot\boldsymbol{v}\,{\rm d}s,\
\boldsymbol{v}\in V;$

where the linear elasticity operator
$\boldsymbol{\mathcal{F}}=(F_{ijkl})_{1\leq i,j,k,l\leq
2}:\Omega\times\mathbb{S}^{2}\to\mathbb{S}^{2}$ is symmetric, bounded, and
satisfies (2.3)(cf. [22]). $j_{\nu}:\Gamma_{C}\times\mathbb{R}\to\mathbb{R}$
is locally Lipschitz on $\mathbb{R}$ for a.e. $\boldsymbol{x}\in\Gamma_{C}$
and satisfies the assumption $(A2)$ with constants $c_{0}$, $c_{1}$ and
$\alpha_{j_{\nu}}$ (cf. [19, 25]).

Therefore, the frictionless normal compliance contact problem reads

$(\boldsymbol{\mathcal{F}}(\boldsymbol{\varepsilon}(\boldsymbol{u})),\boldsymbol{\varepsilon}(\boldsymbol{v}))_{Q}+\int_{\Gamma_{C}}j_{\nu}^{0}(\boldsymbol{u}_{\nu};\boldsymbol{v}_{\nu})\,{\rm
d}s\geq\langle\boldsymbol{f},\boldsymbol{v}\rangle,\quad\boldsymbol{v}\in V.$
(2.6)

Let $\boldsymbol{u}\in V$ and $\lambda_{2}$ be the smallest positive
eigenvalue of the eigenvalue problem

$\int_{\Omega}\boldsymbol{\varepsilon}(\boldsymbol{u})\cdot\boldsymbol{\varepsilon}(\boldsymbol{v})\,{\rm
d}x=\lambda_{2}\int_{\Gamma_{C}}\boldsymbol{u}_{\nu}\cdot\boldsymbol{v}_{\nu}\,{\rm
d}s,\quad\boldsymbol{v}\in V.$

Then the assumption $(A4)$ is satisfied if
$c_{\Gamma}\geq\sqrt{1/\lambda_{2}}$ holds true. The assumption $(A1)$ is
satisfied with $m_{A}=m_{\boldsymbol{\mathcal{F}}}$. Following Theorem 2.1 and
Theorem 2.2, if $\alpha_{j_{\nu}}<\lambda_{2}m_{\boldsymbol{\mathcal{F}}}$,
(2.6) has a unique solution and is equivalent to the following optimization
problem:

$\boldsymbol{u}=\operatorname*{arg\,min}_{\boldsymbol{v}\in
V}E(\boldsymbol{v}),$ (2.7)

where

$E(\boldsymbol{v})=\frac{1}{2}(\boldsymbol{\mathcal{F}}(\boldsymbol{\varepsilon}(\boldsymbol{v})),\boldsymbol{\varepsilon}(\boldsymbol{v}))_{Q}+\int_{\Gamma_{C}}j_{\nu}(\boldsymbol{v}_{\nu})\,{\rm
d}s-\langle\boldsymbol{f},\boldsymbol{v}\rangle.$

##  3 Deep learning-based methods for HVIs

###  3.1 The deep learning method

The main idea of deep learning-based HVIs solvers is to treat DNNs as an
efficient parametrization of the solution space of an HVI. The HVI solution is
identified via seeking a DNN ${\phi}(\boldsymbol{x},\boldsymbol{\theta})$ with
input $\boldsymbol{x}$ and parameters $\boldsymbol{\theta}$ that minimizes the
variational minimization problem related to the HVI. From the discussion in
Section 2, we know that Problem 1 is equivalent to Problem 2. This motivates
the following problem.

Problem 3 Find $\boldsymbol{\theta}^{*}$ such that

$\boldsymbol{\theta}^{*}=\operatorname*{arg\,min}_{\boldsymbol{\theta}}E({\phi}(\boldsymbol{x};\boldsymbol{\theta})),$
(3.1)

where

$E({\phi}(\boldsymbol{x};\boldsymbol{\theta}))=\frac{1}{2}\langle
A{\phi}(\boldsymbol{x};\boldsymbol{\theta}),{\phi}(\boldsymbol{x};\boldsymbol{\theta})\rangle+\int_{\Gamma}j(\gamma_{j}{\phi}(\boldsymbol{x};\boldsymbol{\theta}))\,{\rm
d}s-\langle f,{\phi}(\boldsymbol{x};\boldsymbol{\theta})\rangle.$ (3.2)

In contact mechanics, the first and the third terms on the right-hand side of
(3.2) usually can be formulated as the integrals. Then the objective function
in (3.2) can be viewed as a sum of expectations of several random variables,
which can be solved by stochastic gradient descent methods or its variants
(cf. [6]). We refer to Subsection 4.1 for details along this line.

In this paper, we use two neural networks to approximate the solution of the
HVI introduced in Section 2. The first one is the residual neural network
(ResNet) proposed in [27].

Mathematically, the ResNet can be formulated as follow [15]:

$\boldsymbol{h}_{0}=\boldsymbol{V}\boldsymbol{x},\
\boldsymbol{g}_{\ell}=\sigma(\boldsymbol{W}_{\ell}\boldsymbol{h}_{\ell-1}+\boldsymbol{b}_{\ell}),\
\boldsymbol{h}_{\ell}=\boldsymbol{h}_{\ell-1}+\boldsymbol{U}_{\ell}\boldsymbol{g}_{\ell},\
\ell=1,2,\dots,L,\
{\phi}(\boldsymbol{x};\boldsymbol{\theta})=\boldsymbol{a}^{T}\boldsymbol{h}_{L},$

where $\boldsymbol{V}\in\mathbb{R}^{N_{0}\times d}$,
$\boldsymbol{W}_{\ell}\in\mathbb{R}^{N_{\ell}\times N_{0}}$,
$\boldsymbol{U}_{\ell}\in\mathbb{R}^{N_{0}\times N_{\ell}}$,
$\boldsymbol{b}_{\ell}\in\mathbb{R}^{N_{\ell}}$ for $\ell=1,\dots,L$,
$\boldsymbol{a}\in\mathbb{R}^{{N_{0}}\times m}$. $\sigma(x)$ is a non-linear
activation function. For the purpose of simplicity, we consider
$N_{0}=N_{\ell}=N$ and $\boldsymbol{U}_{\ell}$ is set as the identity matrix.
Here, $L$ is the depth of the ResNet, and $N$ is the width of the network, and
$\boldsymbol{\theta}=\\{\boldsymbol{V},\boldsymbol{a},\boldsymbol{W}_{\ell},\boldsymbol{b}_{\ell}:1\leq\ell\leq
L\\}$ denotes the set of all parameters in $\boldsymbol{\phi}$, which uniquely
determines the neural network.

The other one is a special neural network consisting of an input block and a
few parallel blocks as visualized in Figure 1. Precisely speaking, the whole
network is denoted by $B(\boldsymbol{x},\boldsymbol{\theta})$. The input block
is a ResNet with width $N_{in}$ and depth $L_{in}$ denoted by
$B_{in}(\boldsymbol{x},\boldsymbol{\theta}_{in})$. Following the input block
are $P$ parallel blocks as independent ResNets with width $N_{p}$ and depth
$L_{p}$ denoted by $B_{p}(\boldsymbol{x},\boldsymbol{\theta}_{p})$ for $1\leq
p\leq P$. Then the network $B(\boldsymbol{x},\boldsymbol{\theta})$ can be
formulated as

$B(\boldsymbol{x},\boldsymbol{\theta})=\sum_{p=1}^{P}B_{p}(B_{in}(\boldsymbol{x},\boldsymbol{\theta}_{in}),\boldsymbol{\theta}_{p})$

In this network structure, different ResNets in parallel are trained with
samples at different levels of discretization grids to obtain a PDE solution.
The detailed training algorithm will be introduced later.

Figure 1: The illustration of the network
$B(\boldsymbol{x},\boldsymbol{\theta})$ with $P=5$.

###  3.2 Three strategies for updating network parameters

Now let us introduce three different strategies to update the parameters of
the networks. First of all, we introduce a basic training algorithm to obtain
a numerical solution to the HVI in the form of a ResNet
${\phi}(\boldsymbol{x};\boldsymbol{\theta})$ or the network
${B}(\boldsymbol{x};\boldsymbol{\theta})$. This algorithm is summarized in
Algorithm 1 below.

Algorithm 1 The basic training algorithm.

0: The desired HVI.

0: The parameter set $\boldsymbol{\theta}$ in the solution ResNet
${\phi}(\boldsymbol{x};\boldsymbol{\theta})$ or the solution network
${B}(\boldsymbol{x};\boldsymbol{\theta})$.

Set $Epoch$ as the total iteration number, and assign $N$, $N_{\Gamma}$ as the
sample sizes in the domain $\Omega$, the boundary $\Gamma$, respectively.

Initialize ${\phi}(\boldsymbol{x};\boldsymbol{\theta})$ or
${B}(\boldsymbol{x};\boldsymbol{\theta})$ following the default random
initialization of PyTorch.

for $k=1,\cdots,Epoch$ do

Generate uniformly distributed samples
$\\{\boldsymbol{x}_{i}\\}_{i=1}^{N}\subset\Omega$,
$\\{\boldsymbol{x}_{i}\\}_{i=1}^{N_{\Gamma}}\subset\Gamma$.

Update $\boldsymbol{\theta}$ using the gradient of (3.2) evaluated at the
generated samples.

end for

Next, we introduce a blockwise training algorithm without utilizing multigrid
sampling to obtain the approximate solution to the HVI in the form of the
network ${B}(\boldsymbol{x},\boldsymbol{\theta})$ introduced before. This
algorithm can serve as a comparison to the adaptive multigrid training
algorithm to be introduced later. It consists of two phases: the
initialization phase and the refinement phase. In the initialization phase, we
train the network ${B}(\boldsymbol{x},\boldsymbol{\theta})$ by Algorithm 1 for
$Epoch_{int}$ iterations. In the refinement phase, we train only one block
${B}_{p}(\boldsymbol{x},\boldsymbol{\theta}_{p})$ together with the input
block (i.e., the parameters in other blocks are fixed) in each iteration. We
use $Epoch_{re}$ to denote the number of blockwise refinement steps and
$Epoch_{b}$ to specify how many times a certain block is trained contiguously.
Thus, the total training iterations in the refinement phase is
$Epoch_{re}\times P\times Epoch_{b}$. The blockwise training algorithm is
summarized in Algorithm 2 below.

Algorithm 2 The blockwise training algorithm.

0: The desired HVI.

0: The parameter set $\boldsymbol{\theta}$ in the solution network
${B}(\boldsymbol{x};\boldsymbol{\theta})$.

Set parameters $Epoch_{int}$, $Epoch_{re}$, and $Epoch_{b}$ as iteration
numbers. Assign parameters $N$, $N_{\Gamma}$ for sample sizes in the domain
$\Omega$, the boundary $\Gamma$, respectively. Let $P$ be the number of
blocks.

Initialize ${B}(\boldsymbol{x};\boldsymbol{\theta})$ following the default
random initialization of PyTorch.

for $k=1,\cdots,Epoch_{int}$ do

Generate uniformly distributed samples
$\\{\boldsymbol{x}_{i}\\}_{i=1}^{N}\subset\Omega$,
$\\{\boldsymbol{x}_{i}\\}_{i=1}^{N_{\Gamma}}\subset\Gamma$.

Update $\boldsymbol{\theta}$ using the gradient of (3.2) evaluated at the
generated samples.

end for

for $j=1,\cdots,Epoch_{re}$ do

for $p=1,\cdots,P$ do

for $\ell=1,\cdots,Epoch_{b}$ do

Generate uniformly distributed samples
$\\{\boldsymbol{x}_{i}\\}_{i=1}^{N}\subset\Omega$,
$\\{\boldsymbol{x}_{i}\\}_{i=1}^{N_{\Gamma}}\subset\Gamma$.

Update all the parameters of the $p$-th block and the input block using the
gradient of (3.2) evaluated at the sampled grid points just above.

end for

end for

end for

Since deep learning is challenging to obtain a very accurate solution to PDEs,
we use some technical strategies for training the networks to improve the
accuracy. It is well known that multigrid methods improve solution accuracy
via computation on a hierarchy of discretization meshes. Multigrid computation
traverses from fine grids to coarse grids and finally back to fine grids so as
to efficiently improve the solution error at different levels. The special
order of traverse in multigrid methods is motivated by the fact that
traditional iterative methods tend to reduce high-frequency errors better than
low-frequency errors. In contrast, deep learning optimization reduces low-
frequency errors faster than high-frequency errors as discussed in [60, 8].
Based on this observation, we propose a mesh-free multigrid method for
training the parameters of the networks. A straightforward design of such a
method would move the computation among different levels of grids in the
opposite way of multigrid. However, the numerical performance is not
satisfactory. To overcome this limitation, we adopt the adaptive strategy,
that means, the computation moves from the current grids to the grids
associated with the smallest loss function evaluated at the corresponding
grids. This idea leads to an adaptive multigrid method described below.

Suppose there are $P$ levels of grids generated by uniform discretization with
step sizes $H,2H$,$\cdots$, $2^{P-1}H$, where $H$ is the step size of the
finest grids. Let $Epoch_{int}$, $Epoch_{re}$, and $Epoch_{b}$ be the same
parameters as in the blockwise training algorithm. This algorithm utilizing
multigrid sampling also consists of an initialization phase and a refinement
phase. The initialization phase of the adaptive multigrid training is the same
as the one of the blockwise training algorithm. In the refinement phase, we
first evaluate a loss $E_{p}$ according to (3.2) using $N$ samples from
$\Omega_{p}^{H}$, $N_{\Gamma}$ samples from $\Gamma_{p}^{H}$, where
$\Omega^{H}_{p}$ is the grids in the domain with a step size $2^{P-1}H$,
$\Gamma_{p}^{H}$ is the grids on the boundary $\Gamma$ with a step size
$2^{P-1}H$. If $E_{k}$ is the smallest one among $\\{E_{p}\\}_{p=1}^{P}$, we
then train the $k$-th block and the input block with $N$ samples from
$\Omega_{p}^{H}$, $N_{\Gamma}$ samples from $\Gamma_{p}^{H}$. Thus, the total
training iterations in the refinement phase is $Epoch_{re}\times Epoch_{b}$.
Such a training algorithm chooses the grids corresponding to the smallest loss
as training samples to refine the network parameters and, hence, it is called
an adaptive multigrid training algorithm. See Algorithm 3 below for a detailed
description.

Algorithm 3 The adaptive multigrid training.

0: The desired HVI.

0: The parameter set $\theta$ in the solution network
${B}(\boldsymbol{x};\boldsymbol{\theta})$.

Set parameters $Epoch_{int}$, $Epoch_{b}$, and $Epoch_{re}$ as iteration
numbers, parameters $N$ and $N_{\Gamma}$ as sample sizes in the domain
$\Omega$ and the boundary $\Gamma$, respectively. Let $P$ be the number of
blocks.

Initialize ${B}(\boldsymbol{x};\boldsymbol{\theta})$ following the default
random initialization of PyTorch.

for $k=1,\cdots,Epoch_{int}$ do

Generate uniformly distributed samples
$\\{\boldsymbol{x}_{i}\\}_{i=1}^{N}\subset\Omega$,
$\\{\boldsymbol{x}_{i}\\}_{i=1}^{N_{\Gamma}}\subset\Gamma$.

Update $\boldsymbol{\theta}$ using the gradient of (3.2) evaluated at the
generated samples.

end for

for $j=1,\cdots,Epoch_{re}$ do

for $p=1,\cdots,P$ do

Generate uniformly distributed samples
$\\{\boldsymbol{x}_{i}\\}_{i=1}^{N}\subset\Omega^{H}_{p}$,
$\\{\boldsymbol{x}_{i}\\}_{i=1}^{N_{\Gamma}}\subset\Gamma_{p}^{H}$. Evaluate
the loss $E_{p}$ according to (3.2) using the samples just above.

end for

Let $p=\operatorname*{arg\,min}_{1\leq p\leq P}E_{p}$.

for $\ell=1,\cdots,Epoch_{b}$ do

Generate uniformly distributed samples
$\\{\boldsymbol{x}_{i}\\}_{i=1}^{N}\subset\Omega^{H}_{p}$,
$\\{\boldsymbol{x}_{i}\\}_{i=1}^{N_{\Gamma}}\subset\Gamma_{p}^{H}$.Update all
the parameters of the $p$-th block and the input block using the gradient of
(3.2) evaluated at the sampled grid points just above.

end for

end for

##  4 Numerical experiments

In this section, we shall illustrate the performance of the deep learning
method to a frictional bilateral contact problem and a frictionless normal
compliance contact problem. Numerical comparisons of the deep learning method
and the virtual element method (VEM) are provided. As we will see, the
approximation accuracy in the form of DNNs is nearly the same as that of VEM
in a fine meshsize, but the deep learning method is much easier to implement.
In the meantime, we investigate the numerical performance in terms of
different network structures and different training algorithms. As we will
see, the adaptive mesh-free multigrid algorithm can provide a more accurate
approximation to HVIs than other deep learning methods.

###  4.1 The algorithm description of deep learning methods for contact
problems

From Section 2, we know the frictional bilateral contact problem (2.4) is
equivalent to optimization problem (2.5) and the frictionless normal
compliance contact problem (2.6) is equivalent to (2.7). In order to
numerically solve these problems by deep learning method, we need to
parametrize the solution space of HVIs as introduced in Section 3 firstly. The
solution spaces are $V_{1}$ for the frictional bilateral contact problem (2.5)
and $V$ for the frictionless normal compliance contact problem (2.7). To deal
with the constraints in the admissible spaces $V_{1}$ or $V$, a natural way is
the penalty method that penalizes the loss function with extra terms to
enforce these constraints. However, tuning parameters in the penalty method
may be tedious in practice. Therefore, we construct DNNs satisfying these
constraints automatically as follows:

$\boldsymbol{\phi}(\boldsymbol{x};\boldsymbol{\theta})=\boldsymbol{b}(\boldsymbol{x})*\boldsymbol{\psi}(\boldsymbol{x};\boldsymbol{\theta}),$

where $\boldsymbol{b}(\boldsymbol{x})$ is a known smooth vector-valued
function such that $\boldsymbol{\phi}=\boldsymbol{0}$ on the Dirichlet
boundary $\Gamma_{D}$ and $\phi_{\nu}=0$ on the contact boundary $\Gamma_{C}$
for the frictional bilateral contact problem, or
$\boldsymbol{\phi}=\boldsymbol{0}$ on the Dirichlet boundary $\Gamma_{D}$ for
the frictionless normal compliance contact problem;
$\boldsymbol{\psi}(\boldsymbol{x};\boldsymbol{\theta})$ is an arbitrary DNN in
$\\{\boldsymbol{\psi}(\boldsymbol{x};\boldsymbol{\theta})\\}_{\boldsymbol{\theta}}\approx
H^{1}(\Omega;\mathbb{R}^{2})$; and “$*$” stands for the Hadamard product.
Thus, (3.1) can be formulated as follows

$\boldsymbol{\theta}^{*}=\operatorname*{arg\,min}_{\boldsymbol{\theta}}E(\boldsymbol{b}(\boldsymbol{x})*\boldsymbol{\psi}(\boldsymbol{x};\boldsymbol{\theta})),$
(4.1)

and
$\boldsymbol{u}^{DL}(\boldsymbol{x};\boldsymbol{\theta}^{*})=\boldsymbol{b}(\boldsymbol{x})*\boldsymbol{\psi}(\boldsymbol{x};\boldsymbol{\theta}^{*})$
is the approximate solution to the target HVI.

Next, we discuss the specific form $E(\boldsymbol{b}*\boldsymbol{\psi})$ in
(4.1) for the frictional bilateral contact problem (2.5) and the frictionless
normal compliance contact problem (2.7), respectively.

For the frictional bilateral contact problem (2.5),

$\displaystyle\quad
E(\boldsymbol{b}(\boldsymbol{x})*\boldsymbol{\psi}(\boldsymbol{x};\boldsymbol{\theta}))$
$\displaystyle=\frac{1}{2}(\boldsymbol{\mathcal{F}}(\boldsymbol{\varepsilon}(\boldsymbol{b}(\boldsymbol{x})*\boldsymbol{\psi}(\boldsymbol{x};\boldsymbol{\theta})),\boldsymbol{\varepsilon}(\boldsymbol{b}(\boldsymbol{x})*\boldsymbol{\psi}(\boldsymbol{x};\boldsymbol{\theta}))_{Q}+\int_{\Gamma_{C}}j_{\tau}((\boldsymbol{b}*\boldsymbol{\psi})_{\tau}(\boldsymbol{x};\boldsymbol{\theta}))\,{\rm
d}s-\langle\boldsymbol{f},\boldsymbol{b}(\boldsymbol{x})*\boldsymbol{\psi}(\boldsymbol{x};\boldsymbol{\theta})\rangle$
$\displaystyle=|\Omega|\mathbb{E}_{\boldsymbol{\xi}_{1}}\left[\frac{1}{2}\boldsymbol{\mathcal{F}}(\boldsymbol{\varepsilon}(\boldsymbol{b}(\boldsymbol{\xi_{1}})*\boldsymbol{\psi}(\boldsymbol{\xi}_{1};\boldsymbol{\theta})))_{ij}\boldsymbol{\varepsilon}(\boldsymbol{b}(\boldsymbol{\xi_{1}})*\boldsymbol{\psi}(\boldsymbol{\xi}_{1};\boldsymbol{\theta}))_{ij}-\boldsymbol{f}_{0}(\boldsymbol{\xi}_{1})\cdot\boldsymbol{b}(\boldsymbol{\xi_{1}})*\boldsymbol{\psi}(\boldsymbol{\xi}_{1};\boldsymbol{\theta})\right]$
$\displaystyle\quad-|\Gamma_{T}|\mathbb{E}_{\boldsymbol{\xi}_{2}}\big{[}\boldsymbol{f}_{2}(\boldsymbol{\xi}_{2})\cdot\boldsymbol{b}(\boldsymbol{\xi_{2}})*\boldsymbol{\psi}(\boldsymbol{\xi}_{2};\boldsymbol{\theta})\big{]}+|\Gamma_{C}|\mathbb{E}_{\boldsymbol{\xi}_{3}}\big{[}j_{\tau}((\boldsymbol{b}*\boldsymbol{\psi})_{\tau}(\boldsymbol{\xi}_{3};\boldsymbol{\theta}))\big{]},$
(4.2)

with $\boldsymbol{\xi}_{1}$, $\boldsymbol{\xi}_{2}$, and
$\boldsymbol{\xi}_{3}$ being random vectors following the uniform distribution
over $\Omega$, $\Gamma_{T}$, and $\Gamma_{C}$, respectively. In practice, the
minimization problem (4.1) is solved by the stochastic gradient descent method
[6] or its variants (e.g. Adam[35]) by randomly sampling the integral domains
in the loss function. In each iteration of the optimization algorithm, a
stochastic loss function defined below is minimized instead of the original
loss function in (4.1):

$\displaystyle\quad\hat{E}(\boldsymbol{b}(\boldsymbol{x})*\boldsymbol{\psi}(\boldsymbol{x};\boldsymbol{\theta}))$
$\displaystyle=\frac{|\Omega|}{N}\sum_{l=1}^{N}\left[\frac{1}{2}\boldsymbol{\mathcal{F}}(\boldsymbol{\varepsilon}(\boldsymbol{b}(\boldsymbol{\xi}_{l})*\boldsymbol{\psi}(\boldsymbol{\xi}_{l};\boldsymbol{\theta})))_{ij}\boldsymbol{\varepsilon}(\boldsymbol{b}(\boldsymbol{\xi}_{l})*\boldsymbol{\psi}(\boldsymbol{\xi}_{l};\boldsymbol{\theta}))_{ij}-\boldsymbol{f}_{0}(\boldsymbol{\xi}_{l})\cdot\boldsymbol{b}(\boldsymbol{\xi}_{l})*\boldsymbol{\psi}(\boldsymbol{\xi}_{l};\boldsymbol{\theta})\right]$
$\displaystyle-\frac{|\Gamma_{T}|}{N_{T}}\sum_{l=1}^{N_{T}}\big{[}\boldsymbol{f}_{2}(\boldsymbol{\eta}_{l})\cdot\boldsymbol{b}(\boldsymbol{\eta}_{l})*\boldsymbol{\psi}(\boldsymbol{\eta}_{l};\boldsymbol{\theta})\big{]}+\frac{|\Gamma_{C}|}{N_{T}}\sum_{l=1}^{N_{C}}\big{[}j_{\tau}((\boldsymbol{b}*\boldsymbol{\psi})_{\tau}(\boldsymbol{\zeta}_{l};\boldsymbol{\theta}))\big{]},$
(4.3)

where $\\{\boldsymbol{\xi}_{i}\\}_{i=1}^{N}$,
$\\{\boldsymbol{\eta}_{i}\\}_{i=1}^{N_{T}}$, and
$\\{\boldsymbol{\zeta}_{i}\\}_{i=1}^{N_{C}}$ are independent random vectors
following the uniform distribution over $\Omega$, $\Gamma_{T}$, and
$\Gamma_{C}$, respectively.

Similarly, for the frictionless normal compliance contact problem (2.7),

$\displaystyle\quad
E(\boldsymbol{b}(\boldsymbol{x})*\boldsymbol{\psi}(\boldsymbol{x};\boldsymbol{\theta}))$
$\displaystyle=\frac{1}{2}(\boldsymbol{\mathcal{F}}(\boldsymbol{\varepsilon}(\boldsymbol{b}(\boldsymbol{x})*\boldsymbol{\psi}(\boldsymbol{x};\boldsymbol{\theta})),\boldsymbol{\varepsilon}(\boldsymbol{b}(\boldsymbol{x})*\boldsymbol{\psi}(\boldsymbol{x};\boldsymbol{\theta}))_{Q}+\int_{\Gamma_{C}}j_{\nu}((b*\psi)_{\nu}(\boldsymbol{x};\boldsymbol{\theta}))\,{\rm
d}s-\langle\boldsymbol{f},\boldsymbol{b}(\boldsymbol{x})*\boldsymbol{\psi}(\boldsymbol{x};\boldsymbol{\theta}))\rangle$
$\displaystyle=|\Omega|\mathbb{E}_{\boldsymbol{\xi}_{1}}\left[\frac{1}{2}\boldsymbol{\mathcal{F}}(\boldsymbol{\varepsilon}(\boldsymbol{b}(\boldsymbol{\xi}_{1})*\boldsymbol{\psi}(\boldsymbol{\xi}_{1};\boldsymbol{\theta})))_{ij}\boldsymbol{\varepsilon}(\boldsymbol{b}(\boldsymbol{\xi}_{1})*\boldsymbol{\psi}(\boldsymbol{\xi}_{1};\boldsymbol{\theta}))_{ij}-\boldsymbol{f}_{0}(\boldsymbol{\xi}_{1})\cdot\boldsymbol{b}({\boldsymbol{\xi}_{1}})*\boldsymbol{\psi}(\boldsymbol{\xi}_{1};\boldsymbol{\theta})\right]$
$\displaystyle\quad-|\Gamma_{T}|\mathbb{E}_{\boldsymbol{\xi}_{2}}\big{[}\boldsymbol{f}_{2}(\boldsymbol{\xi}_{2})\cdot\boldsymbol{b}(\boldsymbol{\xi}_{2})*\boldsymbol{\psi}(\boldsymbol{\xi}_{2};\boldsymbol{\theta})\big{]}+|\Gamma_{C}|\mathbb{E}_{\boldsymbol{\xi}_{3}}\big{[}j_{\nu}((b*\psi)_{\nu}(\boldsymbol{\xi}_{3};\boldsymbol{\theta}))\big{]},$
(4.4)

with $\boldsymbol{\xi}_{1}$, $\boldsymbol{\xi}_{2}$, and
$\boldsymbol{\xi}_{3}$ being random vectors following the uniform distribution
over $\Omega$, $\Gamma_{T}$, and $\Gamma_{C}$, respectively. In the
implementation of stochastic optimization algorithm, the loss function
evaluated in each iteration is defined as

$\displaystyle\quad\hat{E}(\boldsymbol{b}(\boldsymbol{x})*\boldsymbol{\psi}(\boldsymbol{x};\boldsymbol{\theta}))$
$\displaystyle=\frac{|\Omega|}{N}\sum_{l=1}^{N}\left[\frac{1}{2}\boldsymbol{\mathcal{F}}(\boldsymbol{\varepsilon}(\boldsymbol{b}(\boldsymbol{\xi}_{l})*\boldsymbol{\psi}(\boldsymbol{\xi}_{l};\boldsymbol{\theta})))_{ij}\boldsymbol{\varepsilon}(\boldsymbol{b}(\boldsymbol{\xi}_{l})*\boldsymbol{\psi}(\boldsymbol{\xi}_{l};\boldsymbol{\theta}))_{ij}-\boldsymbol{f}_{0}(\boldsymbol{\xi}_{l})\cdot\boldsymbol{b}(\boldsymbol{\xi}_{l})*\boldsymbol{\psi}(\boldsymbol{\xi}_{l};\boldsymbol{\theta})\right]$
$\displaystyle-\frac{|\Gamma_{T}|}{N_{T}}\sum_{l=1}^{N_{T}}\big{[}\boldsymbol{f}_{2}(\boldsymbol{\eta}_{l})\cdot\boldsymbol{b}(\boldsymbol{\eta}_{l})*\boldsymbol{\psi}(\boldsymbol{\eta}_{l};\boldsymbol{\theta})\big{]}+\frac{|\Gamma_{C}|}{N_{C}}\sum_{l=1}^{N_{C}}\big{[}j_{\nu}(({b}*{\psi})_{\nu}(\boldsymbol{\zeta}_{l};\boldsymbol{\theta}))\big{]},$
(4.5)

where $\\{\boldsymbol{\xi}_{i}\\}_{i=1}^{N}$,
$\\{\boldsymbol{\eta}_{i}\\}_{i=1}^{N_{T}}$, and
$\\{\boldsymbol{\zeta}_{i}\\}_{i=1}^{N_{C}}$ are independent random vectors
following the uniform distribution over $\Omega$, $\Gamma_{T}$, and
$\Gamma_{C}$, respectively.

###  4.2 Numerical results

In our numerical experiments, we apply the ResNet
$\boldsymbol{\phi}(\boldsymbol{x};\boldsymbol{\theta})$ or the special block
ResNet $\boldsymbol{B}(\boldsymbol{x};\boldsymbol{\theta})$ introduced in
Section 3 as the solution ansatz to HVIs. In the ResNet, we set the depth of
$\boldsymbol{\phi}$ as $L=8$ and the width as $N=50$. In the block ResNet
$\boldsymbol{B}$, we take a ResNet with depth $L=4$ and width $N=50$ as the
input block and take a ResNet with depth $L=4$ and width $N=10$ as each part
in the parallel block. Therefore, the total number of parameters in
$\boldsymbol{\phi}(\boldsymbol{x};\boldsymbol{\theta})$ is about $20,400$ and
the total number of $\boldsymbol{B}(\boldsymbol{x};\boldsymbol{\theta})$ is
about $12,400$. Note that if we choose $\sigma(x)=\text{ReLU}(x)$, the
derivative of the corresponding DNN is a constant almost everywhere, which
leads to DNNs not being able to capture the feature of the PDE solution. In
the literature, the activation function $\text{ReLU}^{\alpha}$ or $\tanh$ are
used for different problem, where $\text{ReLU}(x)=\max\\{x,0\\}$ and $\alpha$
is a positive integer, to overcome this difficulty. All neural networks are
trained by Adam optimizer [35] with a default learning rate $\eta=0.001$ and
exponential decay rates $\beta_{1}=0.9$ and $\beta_{2}=0.999$. The batch size
in the domain is 1024 and the number of training dates on each boundary is 256
for all problems. All numerical experiments are implemented in Python 3.7
using Pytorch 1.3 in an Nvidia GEFORCE RTX 2080 Ti GPU card.

Before reporting numerical results, let us summarize notations used in this
section. Denote by $\boldsymbol{u}^{DL}=(u^{DL}_{1},u^{DL}_{2})^{T}$ as the
approximate solution estimated by the deep learning method. Suppose
$\mathcal{T}_{h}$ is a uniform triangulation of $\Omega$ into $K$ and
$h=$diam$(K)$. Since the true solution is unavailable, for the bilateral
contact problem, we use $\boldsymbol{u}^{ref}$ as the reference solution
evaluated by VEM with $h=2^{-7}$ (see [19]). Similarly, for the frictionless
contact with normal compliance problem, we take $\boldsymbol{u}^{ref}$ as the
reference solution evaluated by VEM with $h=2^{-9}$ (see [19]). Denote the
relative difference between the deep learning solution and the reference
solution as

$\mathcal{E}_{r}=\frac{\|\boldsymbol{u}^{DL}-\boldsymbol{u}^{ref}\|_{E}}{\|\boldsymbol{u}^{ref}\|_{E}},$

where the energy norm $\|\cdot\|_{E}$ is defined by

$\|\boldsymbol{v}\|_{E}=\frac{1}{\sqrt{2}}\left(\boldsymbol{\mathcal{F}}(\boldsymbol{\varepsilon}(\boldsymbol{v})),\boldsymbol{\varepsilon}(\boldsymbol{v})\right)_{Q}^{1/2}.$

######  Example 4.1.

Consider the bilateral contact problem (2.4). The domain
$\Omega=(0,4)\times(0,4)$ is the cross section of a three-dimensional linearly
elastic body and the plane stress condition is imposed. The body is clamped on
$\Gamma_{D}=\\{4\\}\times(0,4)$ and vertical tractions act on
$\Gamma_{T}=(\\{0\\}\times(0,4))\cup((0,4)\times\\{4\\})$. The frictional
contact happens on the boundary $\Gamma_{C}=(0,4)\times\\{0\\}$. The linear
elasticity tensor $\boldsymbol{\mathcal{F}}$ is

$\mathcal{F}_{ij}=\frac{E\kappa}{(1-\kappa^{2})}(\varepsilon_{11}+\varepsilon_{22})\delta_{ij}+\frac{E}{1+\kappa}\varepsilon_{ij},\quad
1\leq i,j,\leq 2,$

where $E$ is the Young modulus, $\kappa$ is the Poisson ratio of the material,
and $\delta_{ij}$ is the Kronecker symbol. In numerical simulations, we use
the following data

$\displaystyle E=2000\ \text{daN}/\text{mm}^{2},\quad\kappa=0.4,$
$\displaystyle\boldsymbol{f}_{0}(x,y)=(0,0)^{T}\ \text{daN}/\text{mm}^{2}\
\text{in}\ \Omega,$
$\displaystyle\boldsymbol{f}_{2}(x,y)=\begin{cases}(200(5-y),-200)^{T}\
\text{daN}/\text{mm}^{2}\quad\text{on}\ \\{0\\}\times(0,4)\\\ (0,0)^{T}\
\text{daN}/\text{mm}^{2}\quad\text{on}\ (0,4)\times\\{4\\},\end{cases}$
$\displaystyle
j_{\tau}(\boldsymbol{z})=\int_{0}^{\|\boldsymbol{z}\|}450e^{-2000t}+450\,{\rm
d}t,\quad\boldsymbol{z}=(x,y).$

Since that $\text{ReLU}^{\alpha}$ activation function is still not appropriate
in this problem because of the numerical overflow in the training process. We
choose the tanh activation function in this example. Figure 2 shows the
numerical solution in the form of ResNet after $50,000$ epochs. The first row
of Figure 2 displays each component of $\boldsymbol{u}^{DL}$ in the domain
whereas the second row shows the approximation on the contact boundary.

Figure 2: The numerical solution in the domain (upper) and on the contact
boundary (bottom) of the bilateral contact problem.

For the purpose of quantifying the accuracy of deep learning methods and
comparing the efficiency of different networks and algorithms, we evaluate the
relative error between the numerical solution and a reference solution
obtained by VEM. For Algorithm 1, we train the network with $50,000$ epochs
either using the ResNet $\boldsymbol{\phi}$ or the block ResNet
$\boldsymbol{B}$. For Algorithm 2, we set $Epoch_{int}=9000$,
$Epoch_{b}=1000$, and $Epoch_{re}=9$. In order to make fair comparisons, we
stop Algorithm 2 after $50,000$ total epochs to update network parameters. For
Algorithm 3, we set the discretization step size $H=1/50$ and use
$Epoch_{int}=9000$, $Epoch_{re}=41$, and $Epoch_{b}=1000$.

From Table 1, we see that the relative errors are reduced to $5\%$ nearly by
all training algorithms. Besides, the number of parameters of ResNet
$\boldsymbol{\phi}$ is almost twice of the one of block ResNet
$\boldsymbol{B}$, but the solution in the form of block ResNet
$\boldsymbol{B}$ is more accurate than the one of ResNet $\boldsymbol{\phi}$.
Moreover, the adaptive mesh-free multigrid algorithm can provide a more
accurate solution than others after the same number of epochs: $41.4\%$ higher
accuracy than the basic algorithm and $35.5\%$ higher accuracy than the
blockwise algorithm.

Table 1: The relative error of different algorithms for the bilateral contact problem. Algorithm 1: the basic algorithm in the literature. Algorithm 2: blockwise training. Algorithm 3: adaptive mesh-free multigrid training.  Algorithm |  Algorithm 1 (ResNet $\boldsymbol{\phi}$) |  Algorithm 1 (Block ResNet $\boldsymbol{B}$) |  Algorithm 2 |  Algorithm 3
---|---|---|---|---
$\mathcal{E}_{r}$ | 0.0481 | 0.0412 | 0.0437 | 0.0282

######  Example 4.2.

Consider the frictionless normal compliance contact problem (2.6). The domain
$\Omega=(0,1)\times(0,1)$ is the cross section of a three-dimensional linearly
elastic body and the plane stress condition is imposed. The body is clamped on
$\Gamma_{D}=(\\{0\\}\times(0,1))\cup(\\{1\\}\times(0,1))$ and the vertical
traction acts on $\Gamma_{T}=\ (0,1)\times\\{1\\}$. The frictional contact
happens on the boundary $\Gamma_{C}=(0,1)\times\\{0\\}$. The linear elasticity
tensor $\boldsymbol{\mathcal{F}}$ is

$\mathcal{F}_{ij}=\frac{E\kappa}{(1+\kappa)(1-2\kappa)}(\varepsilon_{11}+\varepsilon_{22})\delta_{ij}+\frac{E}{1+\kappa}\varepsilon_{ij},\quad
1\leq i,j\leq 2,$

where $E$ is the Young modulus, $\kappa$ is the Poisson ratio of the material,
and $\delta_{ij}$ is the Kronecker symbol. In numerical simulations, we use
the following data

$\displaystyle E=70\ \text{GPa},\quad\kappa=0.3,$
$\displaystyle\boldsymbol{f}_{0}(x,y)=(0,0)^{T}\ \text{GPa}\quad\text{in}\
\Omega$ $\displaystyle\boldsymbol{f}_{2}(x,y)=(0,-52)^{T}\
\text{GPa}\quad\text{on}\ \Gamma_{T}$ $\displaystyle
j_{\nu}(u_{\nu})=\begin{cases}0,\quad u_{\nu}\leq 0,\\\
50u_{\nu}^{2}+0.1u_{\nu},\quad u_{\nu}\in(0,0.1],\\\
20.1u_{\nu}-50u_{\nu}^{2}-1,\quad u_{\nu}\in(0.1,0.15),\\\
200u_{\nu}^{2}-54.9u_{\nu}+4.625,\quad u_{\nu}\geq 0.15.\end{cases}$

The activation function $\text{ReLU}^{2}=\max\\{x^{2},0\\}$ is chosen in this
example based on numerical experience in this example. Figure 3 shows the
numerical solution in the form of ResNet with activation function
$\sigma(x)=\text{ReLU}^{2}$ after $50,000$ epochs. The first row of Figure 3
displays each component of $\boldsymbol{u}^{DL}$ in the domain whereas the
second row shows the numerical solution on the contact boundary.

Figure 3: The numerical solution in the domain (upper) and on the contact
boundary (bottom) of the frictionless normal compliance contact problem.

Similarly, we evaluate the relative error between the numerical solution and
the reference solution to compare the numerical performance of different
networks and different training algorithms. For Algorithm 1, we train the
network with $50,000$ epochs either using the ResNet $\boldsymbol{\phi}$ or
the block ResNet $\boldsymbol{B}$. For Algorithm 2, we set $Epoch_{int}=9000$,
$Epoch_{b}=1000$, and $Epoch_{re}=9$. In order to make a fair comparison, we
stop Algorithm 2 after $50,000$ total epochs. For Algorithm 3, we set the
discretization step size $H=1/200$ and use $Epoch_{int}=9000$,
$Epoch_{re}=41$, and $Epoch_{b}=1000$.

From Table 2, we see that relative errors are reduced to $7\%$ by all deep
learning methods. However, only the adaptive mesh-free multigrid algorithm can
provide an approximate solution with a relative error less than $5\%$ after
$50,000$ epochs. In particular, the adaptive mesh-free multigrid algorithm
improves the accuracy by $38.4\%$ compared to the basic algorithm and by
$22.0\%$ compared to the blockwise algorithm. At last, recall that the number
of parameters of the ResNet $\boldsymbol{\phi}$ is twice the one of the block
ResNet $\boldsymbol{B}$, but the approximate solution in the form of the block
ResNet $\boldsymbol{B}$ is more accurate than the one of the ResNet
$\boldsymbol{\phi}$.

Table 2: The relative error of different algorithms for the frictionless normal compliance contact problem. Algorithm 1: the basic algorithm in the literature. Algorithm 2: blockwise training. Algorithm 3: adaptive mesh-free multigrid training.  Algorithm |  Algorithm 1 (ResNet $\boldsymbol{\phi}$) |  Algorithm 1 (Block ResNet $\boldsymbol{B}$) |  Algorithm 2 |  Algorithm 3
---|---|---|---|---
$\mathcal{E}_{r}$ | 0.0706 | 0.0556 | 0.0558 | 0.0435

##  5 Conclusion

This paper focuses on developing a deep learning method for solving HVIs based
on their variational problems. First of all, the solution space is
parameterized via DNNs, and then the HVI is reformulated as an expectation
minimization problem, which are therefore worked out by stochastic gradient
descent method or its variants (e.g. Adam) combined with three different
training strategies for updating network parameters. As applications to
contact mechanics, a frictional bilateral contact problem and a frictionless
normal compliance contact problem are carried out in details. Numerical
results show that the deep learning method is efficient in solving HVIs and
the adaptive mesh-free multigrid algorithm can provide the best accuracy
solution among the three learning methods discussed. It deserves to mention
that the proposed method is easy to realize in programming.

## Acknowledgments

J. H. was partially supported by the National Key Research and Development
Project (2020YFA0709800) and NSFC (Grant No. 12071289). C. W. was partially
supported by National Science Foundation Award DMS-1849483. We would like to
thank Dr. Fang Feng for offering the comparison data of numerical solutions
for contact problems by the virtual element method used in numerical
experiments.

## References

  * [1] R.A. Adams.  Sobolev Spaces.  Academic Press, New York-London, 1975.
  * [2] C.C. Baniotopoulos, J. Haslinger, Z. Morávková.  Mathematical modeling of delamination and nonmonotone friction problems by hemivariational inequalities.  Appl. Math., 50:1–25, 2005.
  * [3] M. Barboteu, K. Bartosz, P. Kalita.  An analytical and numerical approach to a bilateral contact problem with nonmonotone friction.  Int. J. Appl. Math. Comput. Sci., 23:263–276, 2013.
  * [4] M. Barboteu, K. Bartosz, P. Kalita, A. Ramadan.  Analysis of a contact problem with normal compliance, finite penetration and nonmonotone slip dependent friction.  Commun. Contemp. Math., 16:1350016, 29, 2014.
  * [5] A. R. Barron.  Universal approximation bounds for superpositions of a sigmoidal function.  IEEE Transactions on Information Theory, 39:930–945, 1993.
  * [6] L. Bottou, F. E. Curtis, J. Nocedal.  Optimization methods for large-scale machine learning.  SIAM Rev., 60:223–311, 2018.
  * [7] S. C. Brenner.  Korn’s inequalities for piecewise $H^{1}$ vector fields.  Math. Comp., 73:1067–1087, 2004.
  * [8] Y. Cao, Z. Fang, Y. Wu, D.-X. Zhou, Q. Gu.  Towards Understanding the Spectral Bias of Deep Learning.  arXiv e-prints, page arXiv:1912.01198, 2019.
  * [9] F. Chen, J. Huang, C. Wang, H. Yang.  Friedrichs Learning: Weak Solutions of Partial Differential Equations via Deep Learning.  arXiv e-prints, page arXiv:2012.08023, 2020.
  * [10] F.H. Clarke.  Generalized gradients and applications.  Trans. Amer. Math. Soc., 205:247–262, 1975.
  * [11] F.H. Clarke.  Optimization and nonsmooth analysis.  John Wiley & Sons, Inc., New York, 1983.
  * [12] G. Cybenko.  Approximation by superpositions of a sigmoidal function.  Math. Control Signals Systems, 2:303–314, 1989.
  * [13] M. W. M. G. Dissanayake, N. Phan-Thien.  Neural-network-based approximations for solving partial differential equations.  Comm. Numer. Methods Engrg., 10:195–201, 1994.
  * [14] S. Dong, Z. Li.  Local extreme learning machines and domain decomposition for solving linear and nonlinear partial differential equations.  arXiv e-prints, page arXiv:2012.02895, 2020.
  * [15] W. E, C. Ma, Q. Wang.  A Priori Estimates of the Population Risk for Residual Networks.  arXiv e-prints, page arXiv:1903.02154, 2019.
  * [16] W. E, C. Ma, L. Wu.  A Priori Estimates of the Population Risk for Two-layer Neural Networks.  arXiv e-prints, page arXiv:1810.06397, 2018.
  * [17] W. E, S. Wojtowytsch.  Representation formulas and pointwise properties for barron functions.  arXiv e-prints, page arXiv:2006.05982, 2020.
  * [18] W. E, B. Yu.  The deep ritz method: a deep learning-based numerical algorithm for solving variational problems.  Commun. Math. Stat., 6:1–12, 2018.
  * [19] F. Feng, W. Han, J. Huang.  Virtual element method for an elliptic hemivariational inequality with applications to contact mechanics.  J. Sci. Comput., 81:2388–2412, 2019.
  * [20] Y. Gu, H. Yang, C. Zhou.  Selectnet: Self-paced learning for high-dimensional partial differential equations.  arXiv e-prints, page arXiv:2001.04860, 2020.
  * [21] J. Han, A. Jentzen, W. E.  Solving high-dimensional partial differential equations using deep learning.  Proc. Natl. Acad. Sci. USA, 115:8505–8510, 2018.
  * [22] W. Han.  Minimization principles for elliptic hemivariational inequalities.  Nonlinear Anal. Real World Appl., 54:103114, 13, 2020.
  * [23] W. Han, S. Migórski, M. Sofonea.  A class of variational-hemivariational inequalities with applications to frictional contact problems.  SIAM J. Math. Anal., 46:3891–3912, 2014.
  * [24] W. Han, M. Sofonea.  Numerical analysis of hemivariational inequalities in contact mechanics.  Acta Numer., 28:175–286, 2019.
  * [25] W. Han, M. Sofonea, M. Barboteu.  Numerical analysis of elliptic hemivariational inequalities.  SIAM J. Numer. Anal., 55:640–663, 2017.
  * [26] J. Haslinger, M. Miettinen, P. D. Panagiotopoulos.  Finite element method for hemivariational inequalities.  Kluwer Academic Publishers, Dordrecht, 1999.
  * [27] K. He, X. Zhang, S. Ren, J. Sun.  Deep residual learning for image recognition.  In 2016 IEEE Conference on Computer Vision and Pattern Recognition (CVPR), pages 770–778, 2016.
  * [28] K. Hornik.  Approximation capabilities of multilayer feedforward networks.  Neural Networks, 4:251 – 257, 1991.
  * [29] K. Hornik, M. Stinchcombe, H. White.  Multilayer feedforward networks are universal approximators.  Neural Networks, 2:359 – 366, 1989.
  * [30] J. Huang, H. Wang, H. Yang.  Int-deep: A deep learning initialized iterative method for nonlinear problems.  J. Comput. Phys., 419:109675, 24, 2020.
  * [31] A. D. Jagtap, G. E. Karniadakis.  Extended physics-informed neural networks (xpinns): A generalized space-time domain decomposition based deep learning framework for nonlinear partial differential equations.  Commun. Comput. Phys., 28:2002–2041, 2020.
  * [32] A. D. Jagtap, K. Kawaguchi, G. E. Karniadakis.  Adaptive activation functions accelerate convergence in deep and physics-informed neural networks.  J. Comput. Phys., 404:109136, 23, 2020.
  * [33] K. Joki, A.M. Bagirov, N. Karmitsa, M.M. Mäkelä, S. Taheri.  Double bundle method for finding Clarke stationary points in nonsmooth DC programming.  SIAM J. Optim., 28:1892–1919, 2018.
  * [34] Y. Khoo, J. Lu, L. Ying.  Solving for high-dimensional committor functions using artificial neural networks.  Res. Math. Sci., 6:1–13, 2019.
  * [35] D. P. Kingma, J. Ba.  Adam: a Method for Stochastic Optimization.  arXiv e-prints, page arXiv:1412.6980, 2014.
  * [36] V.A. Kovtunenko.  A hemivariational inequality in crack problems.  Optimization, 60:1071–1089, 2011.
  * [37] V. Kůrková.  Kolmogorov’s theorem and multilayer neural networks.  Neural Networks, 5:501–506, 1992.
  * [38] I.E. Lagaris, A. Likas, D. I. Fotiadis.  Artificial neural networks for solving ordinary and partial differential equations.  IEEE Trans. Neural Networks, 9:987–1000, 1998.
  * [39] H. Lee, I. S. Kang.  Neural algorithm for solving differential equations.  J. Comput. Phys., 91:110–131, 1990.
  * [40] S. Liang, L. Lyu, C. Wang, H. Yang.  Reproducing activation function for deep learning.  arXiv e-prints, page arXiv:2101.04844, 2021.
  * [41] J. Lu, Z. Shen, H. Yang, S. Zhang.  Deep Network Approximation for Smooth Functions.  arXiv e-prints, page arXiv:2001.03040, 2020.
  * [42] M.M. Mäkelä.  Survey of bundle methods for nonsmooth optimization.  Optim. Methods Softw., 17:1–29, 2002.
  * [43] K.S. McFall, J. R. Mahan.  Artificial neural network method for solution of boundary value problems with exact satisfaction of arbitrary boundary conditions.  IEEE Trans. Neural Networks, 20:1221–1233, 2009.
  * [44] S. Migórski, A. Ochal, M. Sofonea.  Nonlinear inclusions and hemivariational inequalities.  Springer, New York, 2013.
  * [45] E.S. Mistakidis and P.D. Panagiotopoulos.  Numerical treatment of problems involving nonmonotone boundary or stress-strain laws.  Computers & Structures, 64:553 – 565, 1997.
  * [46] H. Montanelli, Q. Du.  New error bounds for deep relu networks using sparse grids.  SIAM J. Math. Data Sci., 1:78–92, 2019.
  * [47] H. Montanelli, H. Yang.  Error bounds for deep Relu networks using the Kolmogorov-Arnold superposition theorem.  arXiv e-prints, page arXiv:1906.11945, 2019.
  * [48] D. Motreanu, P. D. Panagiotopoulos.  Minimax theorems and qualitative properties of the solutions of hemivariational inequalities.  Kluwer Academic Publishers, Dordrecht, 1999.
  * [49] Z. Naniewicz, P. D. Panagiotopoulos.  Mathematical theory of hemivariational inequalities and applications.  Marcel Dekker, Inc., New York, 1995.
  * [50] P. D. Panagiotopoulos.  Nonconvex energy functions. hemivariational inequalities and substationarity principles.  Acta Mech., 48:111–130, 1983.
  * [51] P. D. Panagiotopoulos.  Hemivariational inequalities. Applications in Mechanics and Engineering.  Springer, Berlin, 1993.
  * [52] M. Raissi, P. Perdikaris, G.E. Karniadakis.  Physics-informed neural networks: a deep learning framework for solving forward and inverse problems involving nonlinear partial differential equations.  J. Comput. Phys., 378:686 – 707, 2019.
  * [53] K. Rudd, S. Ferrari.  A constrained integration (CINT) approach to solving partial differential equations using artificial neural networks.  Neurocomputing, 155:277 – 285, 2015.
  * [54] R. Shekari Beidokhti, A. Malek.  Solving initial-boundary value problems for systems of partial differential equations using neural networks and optimization techniques.  J. Franklin Inst., 346:898 – 913, 2009.
  * [55] Z. Shen, H. Yang, S. Zhang.  Deep Network Approximation Characterized by Number of Neurons.  arXiv e-prints, page arXiv:1906.05497, 2019.
  * [56] J.W. Siegel, J. Xu.  Approximation rates for neural networks with general activation functions.  Neural Networks, 128:313 – 321, 2020.
  * [57] J. Sirignano, K. Spiliopoulos.  DGM: a deep learning algorithm for solving partial differential equations.  J. Comput. Phys., 375:1339 – 1364, 2018.
  * [58] M.A. Tzaferopoulos, E.S. Mistakidis, C.D. Bisbos, P.D. Panagiotopoulos.  Comparison of two methods for the solution of a class of nonconvex energy problems using convex minimization algorithms.  Comput. & Structures, 57:959–971, 1995.
  * [59] F. Wang, H. Qi.  A discontinuous Galerkin method for an elliptic hemivariational inequality for semipermeable media.  Appl. Math. Lett., 109:106572, 8, 2020.
  * [60] Z.-Q.J. Xu, Y. Zhang, T. Luo, Y Xiao, Z Ma.  Frequency principle: Fourier analysis sheds light on deep neural networks.  Commun. Comput. Phys., 28:1746–1767, 2020\.
  * [61] D. Yarotsky.  Error bounds for approximations with deep relu networks.  Neural Networks, 94:103–114, 2017.
  * [62] D. Yarotsky.  Optimal approximation of continuous functions by very deep relu networks.  In 31st Annual Conference on Learning Theory, volume 75, pages 1–11. 2018.
  * [63] Y. Zang, G. Bao, X. Ye, H. Zhou.  Weak adversarial networks for high-dimensional partial differential equations.  J. Comput. Phys., 411:109409, 14, 2020.


# Comparative Analysis of Retrieval Systems
in the Real World

Dmytro Mozolevskyi Waseem AlShikh Writer, Inc.
<EMAIL_ADDRESS>

##  1 Introduction

In the rapidly evolving landscape of information retrieval and natural
language processing, the integration of advanced language models with search
and retrieval systems has become a cornerstone for enhancing the quality and
efficiency of data processing. This report presents a comprehensive analysis
of various state-of-the-art methods that combine cutting-edge language models
with sophisticated retrieval techniques. Our objective is to evaluate and
compare these methods based on their performance in two critical aspects:
accuracy, as measured by the RobustQA average score (Han et al., 2023), and
efficiency, determined by the average response time.

We have explored a diverse range of methods including Azure Cognitive Search
Retriever integrated with GPT-4,111https://github.com/Azure-Samples/azure-
search-openai-demo Pinecone’s Canopy
framework,222https://www.pinecone.io/blog/canopy-rag-framework/ various
implementations of Langchain with Pinecone and different language models
(OpenAI,
Cohere),333https://python.langchain.com/docs/integrations/vectorstores/pinecone,444https://python.langchain.com/docs/integrations/retrievers/cohere
LlamaIndex with Weaviate Vector Store’s hybrid
search,555https://docs.llamaindex.ai/en/stable/examples/vector_stores/WeaviateIndexDemo-
Hybrid.html Google’s RAG implementation on Cloud VertexAI-
Search,666https://python.langchain.com/docs/templates/rag-google-cloud-
vertexai-search Amazon SageMaker’s
RAG,777https://aws.amazon.com/blogs/machine-learning/question-answering-using-
retrieval-augmented-generation-with-foundation-models-in-amazon-sagemaker-
jumpstart/ and a novel approach combining a graph search algorithm with a
language model and retrieval awareness (Writer Retrieval).

The impetus for this analysis stems from the increasing demand for robust and
responsive question-answering systems in various domains, from customer
service to academic research. As the complexity of queries and the volume of
information grows, it becomes imperative to not only retrieve relevant
information quickly but also to ensure the precision and adaptability of the
responses. The RobustQA metric (Han et al., 2023), a pivotal element in our
analysis, offers a nuanced view of how well these systems perform under
diverse paraphrasing of questions, reflecting real-world querying scenarios.

This report aims to provide insights into the strengths and weaknesses of each
method, drawing comparisons that would help in understanding which
combinations of technologies and approaches are best suited for specific
applications. By delving into performance data, we aim to shed light on the
current state of AI-driven search and retrieval systems and to pave the way
for informed decisions in the deployment and development of these
technologies.

##  2 Background

In the realm of natural language processing (NLP) and question-answering (QA)
systems, the need for robust and reliable evaluation metrics has always been
paramount. The RobustQA (Han et al., 2023) metric is a significant step
towards addressing this need.

RobustQA is an innovative framework designed to evaluate the robustness of QA
systems. It focuses on assessing how well these systems handle diverse
paraphrasings of questions. This approach is crucial because real-world
queries often come in various forms and styles, and a truly robust QA system
must handle this variability effectively.

The framework employs a novel methodology that goes beyond traditional
evaluation metrics, which often rely on a limited set of question variations.
RobustQA introduces a broader and more challenging set of paraphrased
questions, thereby providing a more comprehensive and realistic assessment of
a QA system’s performance.

By incorporating RobustQA into our evaluation of different language model
integrations with search and retrieval systems, we aim to gain deeper insights
into their real-world effectiveness. This approach allows us to understand not
only the accuracy of these systems in responding to queries but also their
adaptability and resilience to varied linguistic expressions. Various
statistics about the data in RobustQA are shown in Table 1.

Domain Label # Test Questions # Documents # Passages Data Source NQ Wikipedia
[NQ] 3,610 - 21,015,324 NQ RobustQA Web-search [SE] 31,760 13,791,373
13,791,592 SearchQA Biomedical [BI] 1,956 15,559,026 37,406,880 BioASQ Finance
[FI] 3,669 57,638 105,777 FiQA Lifestyle [LI] 2,214 119,461 241,780 LoTTE
Recreation [RE] 2,096 166,975 315,203 LoTTE Technology [TE] 2,115 638,509
1,252,402 LoTTE Science [SC] 1,426 1,694,164 3,063,916 LoTTE Writing [WR]
2,696 199,994 347,322 LoTTE

Table 1: Table and caption recreated from (Han et al., 2023, Table 2). Data
summary: Natural Questions (Kwiatkowski et al., 2019, NQ) (top) v.s. RobustQA
(bottom). # Documents for NQ is missing because we directly use the passage
split provided by (Karpukhin et al., 2020). Passages consists of 100
continuous tokens at most from the original documents.

##  3 Experiments and Results

In our empirical evaluation, we rigorously assessed the efficacy of eight
different retrieval system configurations by scrutinizing two primary
dimensions: the RobustQA metric score and the response time latency. Both
metrics are salient indicators of the systems’ proficiency in handling real-
world query scenarios with swiftness and precision. The technical
configurations tested ranged between different combinations of search
architectures, integrated language models, and indexing techniques.

RobustQA Score Estimation: The RobustQA score computation leverages a
precision-oriented methodology that quantifies the systems’ accuracy in
retrieving the correct response against a corpus of paraphrased queries. This
measure is calculated through an analytic comparison, where each system-
generated response is benchmarked against a gold-standard dataset that
encapsulates a broad spectrum of question paraphrases, reflecting realistic
variance in linguistic expression. A high-resolution scoring mechanism, it
accommodates shuffled question sequences to mitigate position bias and ensure
comprehensiveness in evaluation.

Response Time Measurement: This metric is quantified by recording the interval
from the initiation of the query request to the point when a complete response
set is rendered to the end user. Precision timers embedded within the system’s
operational pipeline log these intervals, capturing the raw computational
throughput and data retrieval cadence. This temporal metric critically
reflects the operational efficiency and scalability of the system, with
adjustments made for network latency effects to ensure the metric purity.

Moving to the experimental setup, each retrieval system configuration was
rigorously evaluated in a controlled test environment to ensure repeatability
and consistency. The subsequent report meticulously tabulates these findings:

Pipeline |  RobustQA Avg. score |  Avg. response time (secs)
---|---|---
Azure Cognitive Search Retriever + GPT4 + Ada |  72.36 |  >1.0s
Canopy (Pinecone) |  59.61 |  >1.0s
Langchain + Pinecone + OpenAI |  61.42 |  <0.8s
Langchain + Pinecone + Cohere |  69.02 |  <0.6s
LlamaIndex + Weaviate Vector Store - Hybrid Search |  75.89 |  <1.0s
RAG Google Cloud VertexAI-Search + Bison |  51.08 |  >0.8s
RAG Amazon SageMaker |  32.74 |  <2.0s
Graph search algorithm + LLM + Retrieval awareness (Writer Retrieval) |  86.31 |  <0.6s
Table 2: Comparison of different pipelines on RobustQA in terms of both
performance and response time. Best scores are in bold.

### Azure Cognitive Search Retriever + GPT4 + Ada888https://github.com/Azure-
Samples/azure-search-openai-demo

This combination shows a strong performance in terms of accuracy, indicating
good integration between Azure’s search capabilities and GPT-4’s Ada model.
The response time is over one second, suggesting a balance between
thoroughness and speed.

### Canopy Pinecone)999https://www.pinecone.io/blog/canopy-rag-framework/

The lower score here could indicate limitations in either the Canopy
framework’s retrieval ability or its integration with Pinecone. The response
time is similar to the Azure Cognitive Search setup.

### Langchain + Pinecone +
OpenAI101010https://python.langchain.com/docs/integrations/vectorstores/pinecone

This method shows a modest improvement in accuracy over Canopy, with a notably
faster response time. The Langchain framework might be more efficient in
querying or processing results.

### Langchain + Pinecone +
Cohere111111https://python.langchain.com/docs/integrations/retrievers/cohere

This setup shows a significant jump in accuracy and a further reduction in
response time. The Cohere model’s integration with Langchain and Pinecone
seems particularly effective.

### LlamaIndex + Weaviate Vecor Score - Hybrid
Search121212https://docs.llamaindex.ai/en/stable/examples/vector_stores/WeaviateIndexDemo-
Hybrid.html

This method is second best in terms of accuracy, suggesting that LlamaIndex’s
hybrid search approach effectively leverages Weaviate’s vector store
capabilities. The sub-one-second response time also indicates efficiency.

### RAG Google Cloud VertexAI-Search +
Bison131313https://python.langchain.com/docs/templates/rag-google-cloud-
vertexai-search

This combination shows the lowest accuracy, which might be due to how RAG and
Bison interact or how they leverage Google Cloud’s VertexAI-Search. The
response time is moderate.

### RAG Amazon SageMaker141414https://aws.amazon.com/blogs/machine-
learning/question-answering-using-retrieval-augmented-generation-with-
foundation-models-in-amazon-sagemaker-jumpstart/

This setup has the lowest score and the longest response time, indicating
possible inefficiencies or mismatches in the integration of RAG with Amazon
SageMaker.

### Graph search algorithm + LLM + Retrieval awareness (Writer Retrieval)

This method outperforms all others in accuracy while maintaining a fast
response time, highlighting the effectiveness of combining a graph search
algorithm with a language model and retrieval awareness.

Azure+GPT4+AdaCanopy
(Pinecone)Langchain+Pinecone+OpenAILangchain+Pinecone+CohereLlamaIndex+WeaviateRAG-
Google-VertexAI+BisonRAG Amazon SageMakerWriter
Retrieval$0$$20$$40$$60$$80$RobustQA Avg.ScoreRobustQA Avg. ScoreRuntime
(secs)$0$$0.5$$1$$1.5$$2$Runtime (secs)

Figure 1: Visual representation of different pipelines in terms of both
RobustQA Average score and Runtime.

##  4 Conclusion

Figure 1 visually represents the comparative analysis of different search and
retrieval methods integrated with language models. It juxtaposes the RobustQA
average score (in blue bars) and average response time (in red line with
markers) for each method. This visualization aids in understanding how each
method performs in terms of accuracy and efficiency. The Graph search
algorithm with LLM and Retrieval awareness (Writer Retrieval) stands out as
the most effective method, balancing high accuracy with quick response times.
LlamaIndex with Weaviate Vector Store is also notable for its high accuracy.
On the other hand, RAG implementations (both Google Cloud and Amazon
SageMaker) lag in performance. This analysis suggests a trend where
specialized retrieval-aware methods combined with efficient language models
lead to better performance in both accuracy and response time.

## References

  * Han et al. (2023) Rujun Han, Peng Qi, Yuhao Zhang, Lan Liu, Juliette Burger, William Yang Wang, Zhiheng Huang, Bing Xiang, and Dan Roth. 2023.  RobustQA: Benchmarking the robustness of domain adaptation for open-domain question answering.  In _Findings of the Association for Computational Linguistics: ACL 2023_ , pages 4294–4311, Toronto, Canada. Association for Computational Linguistics.
  * Karpukhin et al. (2020) Vladimir Karpukhin, Barlas Oguz, Sewon Min, Patrick Lewis, Ledell Wu, Sergey Edunov, Danqi Chen, and Wen-tau Yih. 2020.  Dense passage retrieval for open-domain question answering.  In _Proceedings of the 2020 Conference on Empirical Methods in Natural Language Processing (EMNLP)_ , pages 6769–6781, Online. Association for Computational Linguistics.
  * Kwiatkowski et al. (2019) Tom Kwiatkowski, Jennimaria Palomaki, Olivia Redfield, Michael Collins, Ankur Parikh, Chris Alberti, Danielle Epstein, Illia Polosukhin, Jacob Devlin, Kenton Lee, Kristina Toutanova, Llion Jones, Matthew Kelcey, Ming-Wei Chang, Andrew M. Dai, Jakob Uszkoreit, Quoc Le, and Slav Petrov. 2019.  Natural Questions: A Benchmark for Question Answering Research.  _Transactions of the Association for Computational Linguistics_ , 7:453–466.


Solutions of the (3 + 1) dimensional

Charney - Obukhov equation for the ocean,

Part I: case of “separation” of variables

V.A. Goloveshkin, A.G. Kudryavtsev, N.N. Myagkov

Institute of Applied Mechanics,

Russian Academy of Sciences, Moscow 125040, Russia

###### Abstract

Exact solutions describing Rossby waves and vortices in ocean propagating
along the zonal direction at a constant velocity are considered for the (3 +
1) -dimensional nonlinear Charney - Obukhov equation. In the first part of our
work, we give examples of solutions of the Charney-Obukhov equation that
satisfy the nonlinear boundary conditions of the ocean, for a special case
that corresponds to the ”separation” of variables.

In recent papers [1], [2], [3], exact solutions of the Charney–Obukhov
equation for the ocean were obtained in the form of synoptic-scale waves and
eddies against the background of a zonal flow depending on the vertical
coordinate $z$. The zonal flow in these solutions is itself a solution to the
Charney–Obukhov equation, while waves and vortices are not. These solutions
are of physical interest, since they demonstrate the inseparable connection
between the emerging motion of the ocean and the zonal flow. It should be
noted that for solutions describing drift in the zonal direction with a
constant velocity $V$, Charney–Obukhov equation can be reduced to solving a
linear Helmholtz-type equation. This reduction is well known in hydrodynamics
and was used to obtain solutions in [1], [2], [3]. To describe the processes
in the ocean, the Charney–Obukhov equations are supplemented with nonlinear
boundary conditions at the bottom and on the surface of the ocean. Taking into
account the connection of the Charney-Obukhov equation with a linear equation,
the Charney–Obukhov equation itself without taking into account the boundary
conditions has a large set of solutions and a superposition of solutions. An
interesting property of the considered model is the possibility of fulfilling
nonlinear boundary conditions for a wide class of solutions of the Charney-
Obukhov equation. In the first part of our work, we give examples of solutions
of the Charney-Obukhov equation that satisfy the nonlinear boundary conditions
of the ocean, for a special case that corresponds to the ”separation” of
variables.

It is well known that a large number of large scale circulations of the
Earth’s atmosphere and ocean is described by nonlinear quasi-geostrophic
potential vorticity equation (1), which is also known as the Charney - Obukhov
equation. We consider the (3+1)-dimensional nonlinear Charney–Obukhov equation
in the form

${\frac{\partial^{3}}{\partial{x}^{2}\partial
t}}p\left(t,x,y,z\right)+{\frac{\partial^{3}}{\partial{y}^{2}\partial
t}}p\left(t,x,y,z\right)+{\frac{\partial^{3}}{\partial{z}^{2}\partial
t}}p\left(t,x,y,z\right)\\\ +\left({\frac{\partial}{\partial
x}}p\left(t,x,y,z\right)\right)\left({\frac{\partial^{3}}{\partial
y\partial{x}^{2}}}p\left(t,x,y,z\right)+{\frac{\partial^{3}}{\partial{y}^{3}}}p\left(t,x,y,z\right)+{\frac{\partial^{3}}{\partial{z}^{2}\partial
y}}p\left(t,x,y,z\right)\right)\\\ -\left({\frac{\partial}{\partial
y}}p\left(t,x,y,z\right)\right)\left({\frac{\partial^{3}}{\partial{x}^{3}}}p\left(t,x,y,z\right)+{\frac{\partial^{3}}{\partial{y}^{2}\partial
x}}p\left(t,x,y,z\right)+{\frac{\partial^{3}}{\partial{z}^{2}\partial
x}}p\left(t,x,y,z\right)\right)\\\ +\beta\,{\frac{\partial}{\partial
x}}p\left(t,x,y,z\right)=0$ (1)

Where $p$ is the dimensionless geostrophic stream function, bearing the sense
of relative pressure perturbation; $\beta$ is the dimensionless meridional
(northern) gradient of Coriolis parameter. As usual, we assume that the
x-coordinate is east, the y-coordinate is north, and the z-coordinate is up.

We consider for equation (1) the boundary conditions for an ocean with a flat
bottom and a rigid lid. Such boundary conditions are often used to model waves
and vortices in the ocean, a detailed discussion of these boundary conditions
is given, for example, in the book [4]. The boundary condition operator has
the form

${\frac{\partial^{2}}{\partial z\partial
t}}p\left(t,x,y,z\right)-\left({\frac{\partial}{\partial
y}}p\left(t,x,y,z\right)\right){\frac{\partial^{2}}{\partial z\partial
x}}p\left(t,x,y,z\right)\\\ +\left({\frac{\partial}{\partial
x}}p\left(t,x,y,z\right)\right){\frac{\partial^{2}}{\partial z\partial
y}}p\left(t,x,y,z\right)$ (2)

We require that the boundary condition operator be equal to zero at $z=0$ and
$z=H$, where H is the ocean depth ( $0\leq z\leq H$).

Solution 1.

The solution that satisfies the boundary conditions is

$F\left(x-Vt,y\right)+\left(\cos\left({\frac{\pi\,nz}{H}}\right)P-{\frac{\beta\,{H}^{2}}{{\pi}^{2}{n}^{2}}}-V\right)y$
(3)

where $P$ \- arbitrary constant, $n=1,2,3,...$, the function $F$ satisfies the
equation

${\frac{\partial^{2}}{\partial{s}^{2}}}F\left(s,y\right)+{\frac{\partial^{2}}{\partial{y}^{2}}}F\left(s,y\right)+\frac{{\pi}^{2}{n}^{2}}{{H}^{2}}F\left(s,y\right)=0$
(4)

The boundary condition operator on this solution has the form

$-\left({\frac{\partial}{\partial
x}}F\left(x-Vt,y\right)\right)\sin\left({\frac{\pi\,nz}{H}}\right)\pi\,nP{H}^{-1}$
(5)

Solution 2.

The solution that satisfies the boundary conditions is

$\sin\left({\it
k_{z}}\,z\right)F\left(x-Vt,y\right)+\left({\frac{\beta\,\left(\sin\left(Kz\right)M+\cos\left(Kz\right)-1\right)}{{K}^{2}}}-V\right)y$
(6)

where ${K}^{2}={{\it k_{z}}}^{2}+{{\it K_{r}}}^{2}$, ${\it k_{z}},{\it K_{r}}$
\- arbitrary constants,

$M={\frac{{\it k_{z}}\,\cos\left({\it k_{z}}\,H\right)\cos\left(KH\right)-{\it
k_{z}}\,\cos\left({\it k_{z}}\,H\right)+K\sin\left({\it
k_{z}}\,H\right)\sin\left(KH\right)}{K\sin\left({\it
k_{z}}\,H\right)\cos\left(KH\right)-{\it k_{z}}\,\cos\left({\it
k_{z}}\,H\right)\sin\left(KH\right)}}$ (7)

The function $F$ satisfies the equation

${\frac{\partial^{2}}{\partial{s}^{2}}}F\left(s,y\right)+{\frac{\partial^{2}}{\partial{y}^{2}}}F\left(s,y\right)+{{\it
K_{r}}}^{2}F\left(s,y\right)=0$ (8)

The boundary condition operator on this solution has the form

$\left({\frac{\partial}{\partial
x}}F\left(x-Vt,y\right)\right)\beta{K}^{-2}\,\left(M\left(K\sin\left({\it
k_{z}}\,z\right)\cos\left(Kz\right)-{\it k_{z}}\cos\left({\it
k_{z}}\,z\right)\sin\left(Kz\right)\right)\right.\\\ \left.-{\it
k_{z}}\cos\left({\it k_{z}}\,z\right)\cos\left(Kz\right)+{\it
k_{z}}\cos\left({\it k_{z}}\,z\right)-K\sin\left({\it
k_{z}}\,z\right)\sin\left(Kz\right)\right)$ (9)

Solution 3.

The solution that satisfies the boundary conditions is

$\cos\left({\it
k_{z}}\,z\right)F\left(x-Vt,y\right)+\left({\frac{\beta\,\left(\cos\left(Kz\right)M-1\right)}{{K}^{2}}}-V\right)y$
(10)

where ${K}^{2}={{\it k_{z}}}^{2}+{{\it K_{r}}}^{2}$, ${\it k_{z}},{\it K_{r}}$
\- arbitrary constants,

$M={\frac{{\it k_{z}}\,\sin\left({\it k_{z}}\,H\right)}{{\it
k_{z}}\sin\left({\it k_{z}}\,H\right)\cos\left(KH\right)-K\cos\left({\it
k_{z}}\,H\right)\sin\left(KH\right)}}$ (11)

The function $F$ satisfies the equation (8).

The boundary condition operator on this solution has the form

$\left({\frac{\partial}{\partial
x}}F\left(x-Vt,y\right)\right)\beta{K}^{-2}\left(M\left({\it
k_{z}}\,\sin\left({\it k_{z}}\,z\right)\cos\left(Kz\right)-K\cos\left({\it
k_{z}}\,z\right)\sin\left(Kz\right)\right)\right.\\\ \left.-{\it
k_{z}}\sin\left({\it k_{z}}\,z\right)\right)$ (12)

Solution 4.

The solution that satisfies the boundary conditions is

$\sinh\left({\it
k_{z}}\,z\right)F\left(x-Vt,y\right)+\left({\frac{\beta\,\left(\sin\left(Kz\right)M+\cos\left(Kz\right)-1\right)}{{K}^{2}}}-V\right)y$
(13)

where ${K}^{2}=-{{\it k_{z}}}^{2}+{{\it K_{r}}}^{2}$, ${\it k_{z}},{\it
K_{r}}$ \- arbitrary constants, ${\it k_{z}}<{\it K_{r}}$,

$M={\frac{\cosh\left({\it k_{z}}\,H\right){\it
k_{z}}\,\cos\left(KH\right)-\cosh\left({\it k_{z}}\,H\right){\it
k_{z}}+\sinh\left({\it k_{z}}\,H\right)K\sin\left(KH\right)}{-\cosh\left({\it
k_{z}}\,H\right){\it k_{z}}\,\sin\left(KH\right)+K\sinh\left({\it
k_{z}}\,H\right)\cos\left(KH\right)}}$ (14)

The function $F$ satisfies the equation (8).

The boundary condition operator on this solution has the form

$\left({\frac{\partial}{\partial
x}}F\left(x-Vt,y\right)\right)\beta{K}^{-2}\,\left(-\cosh\left({\it
k_{z}}\,z\right){\it k_{z}}\,\sin\left(Kz\right)M\right.\\\
\left.+\cosh\left({\it k_{z}}\,z\right){\it k_{z}}-\cosh\left({\it
k_{z}}\,z\right){\it k_{z}}\,\cos\left(Kz\right)\right.\\\
\left.+\sinh\left({\it k_{z}}\,z\right)K\cos\left(Kz\right)M-\sinh\left({\it
k_{z}}\,z\right)K\sin\left(Kz\right)\right)$ (15)

Solution 5.

The solution that satisfies the boundary conditions is

$\sinh\left({\it
k_{z}}\,z\right)F\left(x-Vt,y\right)+\left({\frac{\beta\,\left(\sinh\left(Kz\right)M-\cosh\left(Kz\right)+1\right)}{{K}^{2}}}-V\right)y$
(16)

where ${K}^{2}={{\it k_{z}}}^{2}-{{\it K_{r}}}^{2}$, ${\it k_{z}},{\it K_{r}}$
\- arbitrary constants, ${\it k_{z}}>{\it K_{r}}$,

$M={\frac{-\cosh\left({\it k_{z}}\,H\right){\it
k_{z}}\,\cosh\left(KH\right)+\cosh\left({\it k_{z}}\,H\right){\it
k_{z}}+\sinh\left({\it k_{z}}\,H\right)K\sinh\left(KH\right)}{-\cosh\left({\it
k_{z}}\,H\right){\it k_{z}}\,\sinh\left(KH\right)+K\sinh\left({\it
k_{z}}\,H\right)\cosh\left(KH\right)}}$ (17)

The function $F$ satisfies the equation (8).

The boundary condition operator on this solution has the form

$\left({\frac{\partial}{\partial
x}}F\left(x-Vt,y\right)\right)\beta{K}^{-2}\,\left(-\cosh\left({\it
k_{z}}\,z\right){\it k_{z}}\,\sinh\left(Kz\right)M\right.\\\
\left.-\cosh\left({\it k_{z}}\,z\right){\it k_{z}}+\cosh\left({\it
k_{z}}\,z\right){\it k_{z}}\,\cosh\left(Kz\right)\right.\\\
\left.+\sinh\left({\it k_{z}}\,z\right)K\cosh\left(Kz\right)M-\sinh\left({\it
k_{z}}\,z\right)K\sinh\left(Kz\right)\right)$ (18)

Solution 6.

The solution that satisfies the boundary conditions is

$\sinh\left({\it
k_{z}}\,z\right)F\left(x-Vt,y\right)+\left(-1/2\,\beta\,{z}^{2}+1/2\,\beta\,HMz-V\right)y$
(19)

where ${\it k_{z}}$ \- arbitrary constant,

$M={\frac{2\,\sinh\left({\it k_{z}}\,H\right)-\cosh\left({\it
k_{z}}\,H\right){\it k_{z}}\,H}{\sinh\left({\it
k_{z}}\,H\right)-\cosh\left({\it k_{z}}\,H\right){\it k_{z}}\,H}}$ (20)

The function $F$ satisfies the equation

${\frac{\partial^{2}}{\partial{s}^{2}}}F\left(s,y\right)+{\frac{\partial^{2}}{\partial{y}^{2}}}F\left(s,y\right)+{{\it
k_{z}}}^{2}F\left(s,y\right)=0$ (21)

The boundary condition operator on this solution has the form

$-1/2\,\left({\frac{\partial}{\partial
x}}F\left(x-Vt,y\right)\right)\beta\,\left(-\cosh\left({\it
k_{z}}\,z\right){\it k_{z}}\,{z}^{2}+\cosh\left({\it k_{z}}\,z\right){\it
k_{z}}\,HMz\right.\\\ \left.+2\,\sinh\left({\it
k_{z}}\,z\right)z-\sinh\left({\it k_{z}}\,z\right)HM\right)$ (22)

Solution 7.

The solution that satisfies the boundary conditions is

$\cosh\left({\it
k_{z}}\,z\right)F\left(x-Vt,y\right)+\left({\frac{\beta\,\left(\cos\left(Kz\right)M-1\right)}{{K}^{2}}}-V\right)y$
(23)

where ${K}^{2}=-{{\it k_{z}}}^{2}+{{\it K_{r}}}^{2}$, ${\it k_{z}},{\it
K_{r}}$ \- arbitrary constants, ${\it k_{z}}<{\it K_{r}}$,

$M={\frac{\sinh\left({\it k_{z}}\,H\right){\it k_{z}}}{\sinh\left({\it
k_{z}}\,H\right){\it k_{z}}\,\cos\left(KH\right)+K\cosh\left({\it
k_{z}}\,H\right)\sin\left(KH\right)}}$ (24)

The function $F$ satisfies the equation (8).

The boundary condition operator on this solution has the form

$-\left({\frac{\partial}{\partial
x}}F\left(x-Vt,y\right)\right)\beta{K}^{-2}\,\left(\sinh\left({\it
k_{z}}\,z\right){\it k_{z}}\,\cos\left(Kz\right)M-\sinh\left({\it
k_{z}}\,z\right){\it k_{z}}\right.\\\ \left.+\cosh\left({\it
k_{z}}\,z\right)K\sin\left(Kz\right)M\right)$ (25)

Solution 8.

The solution that satisfies the boundary conditions is

$\cosh\left({\it
k_{z}}\,z\right)F\left(x-Vt,y\right)+\left({\frac{\beta\,\left(\cosh\left(Kz\right)M+1\right)}{{K}^{2}}}-V\right)y$
(26)

where ${K}^{2}={{\it k_{z}}}^{2}-{{\it K_{r}}}^{2}$, ${\it k_{z}},{\it K_{r}}$
\- arbitrary constants, ${\it k_{z}}>{\it K_{r}}$,

$M={\frac{\sinh\left({\it k_{z}}\,H\right){\it k_{z}}}{-\sinh\left({\it
k_{z}}\,H\right){\it k_{z}}\,\cosh\left(KH\right)+K\cosh\left({\it
k_{z}}\,H\right)\sinh\left(KH\right)}}$ (27)

The function $F$ satisfies the equation (8).

The boundary condition operator on this solution has the form

$\left({\frac{\partial}{\partial
x}}F\left(x-Vt,y\right)\right)\beta{K}^{-2}\,\left(-\sinh\left({\it
k_{z}}\,z\right){\it k_{z}}\,\cosh\left(Kz\right)M-\sinh\left({\it
k_{z}}\,z\right){\it k_{z}}\right.\\\ \left.+\cosh\left({\it
k_{z}}\,z\right)K\sinh\left(Kz\right)M\right)$ (28)

Solution 9.

The solution that satisfies the boundary conditions is

$\cosh\left({\it
k_{z}}\,z\right)F\left(x-Vt,y\right)+\left(-1/2\,\beta\,{z}^{2}+1/2\,\beta\,{H}^{2}M-V\right)y$
(29)

where ${\it k_{z}}$ \- arbitrary constant,

$M=1-2\,{\frac{\cosh\left({\it k_{z}}\,H\right)}{\sinh\left({\it
k_{z}}\,H\right)H{\it k_{z}}}}$ (30)

The function $F$ satisfies the equation (21).

The boundary condition operator on this solution has the form

$-1/2\,\left({\frac{\partial}{\partial
x}}F\left(x-Vt,y\right)\right)\beta\,\left(-\sinh\left({\it
k_{z}}\,z\right){\it k_{z}}\,{z}^{2}+\sinh\left({\it k_{z}}\,z\right){\it
k_{z}}\,{H}^{2}M\right.\\\ \left.+2\,\cosh\left({\it k_{z}}\,z\right)z\right)$
(31)

All presented solutions have the form

$f(z)F\left(x-Vt,y\right)+U(z;V)y$ (32)

where $U(z;V)$ is the velocity of the steady zonal background flow which
depends on the velocity $V$. Thus, the first term of the formula (32) can be
represented as the product of the horizontal coordinate function and the
vertical coordinate function, this is a case of ”separation” of variables.

## References

  * [1] A.G.Kudryavtsev, N.N.Myagkov. “On exact solutions of the Charney Obukhov equation for the ocean,” Physics Letters A 446 (2022) 128282 https://doi.org/10.1016/j.physleta.2022.128282
  * [2] A.G. Kudryavtsev, N.N. Myagkov. “New exact spatially localized solutions of the (3 + 1) -dimensional Charney Obukhov equation for the ocean,” Physics of Fluids 34, 126604 (2022); https://doi.org/10.1063/5.0129694
  * [3] A.G. Kudryavtsev, N.N. Myagkov. “On the superposition of solutions of the (3+1) dimensional Charney–Obukhov equation for the ocean,” Physics of Fluids 35, 051701 (2023); https://doi.org/10.1063/5.0150230
  * [4] J. Pedlosky, Geophysical Fluid Dynamics, Springer-Verlag, New York, 1987.