# A Local Discontinuous Galerkin method for the Benajamin-Ono equation

Mukul Dwivedi
Department of Mathematics, Indian Institute of Technology Jammu, Jagti, NH-44
Bypass Road, Post Office Nagrota, Jammu - 181221, India
<EMAIL_ADDRESS>and  Tanmay Sarkar
Department of Mathematics, Indian Institute of Technology Jammu, Jagti, NH-44
Bypass Road, Post Office Nagrota, Jammu - 181221, India
<EMAIL_ADDRESS>
###### Abstract.

The main purpose of this paper is to design a local discontinuous Galerkin
(LDG) method for the Benjamin-Ono equation. We analyze the stability and error
estimates for the semi-discrete LDG scheme. We prove that the scheme is
$L^{2}$-stable and it converges at a rate $\mathcal{O}(h^{k+1/2})$ for general
nonlinear flux. Furthermore, we develop a fully discrete LDG scheme using the
four-stage fourth order Runge-Kutta method and ensure the devised scheme is
strongly stable in case of linear flux using two-step and three-step stability
approach under an appropriate time step constraint. Numerical examples are
provided to validate the efficiency and accuracy of the method.

###### Key words and phrases:

Local Discontinuous Galerkin method, Benjamin-Ono equation, Hilbert transform.

###### 2020 Mathematics Subject Classification:

Primary: 65M60, 35R09; Secondary: 65M12.

##  1\. Introduction

We consider the following Cauchy problem associated to generalized Benjamin-
Ono (BO) equation:

$\begin{cases}U_{t}+f(U)x-\mathcal{H}U{xx}=0,\qquad&(x,t)\in\mathbb{R}\times(0,T],\\\
U(x,0)=U_{0}(x),\qquad&x\in\mathbb{R},\end{cases}$ (1.1)

where $T>0$ is fixed, $U_{0}$ represents the prescribed initial data, $f$ is
any given function of $U$ and $U:\mathbb{R}\times(0,T]\rightarrow\mathbb{R}$
is the unknown, and $\mathcal{H}$ denotes the Hilbert transform [11, 35],
defined by the principle value integral

$\displaystyle\mathcal{H}U(x):=\text{P.V.}\frac{1}{\pi}\int_{\mathbb{R}}\frac{U(x-y)}{y}~{}dy.$

The Benjamin-Ono equation (1.1) is a nonlinear, non-local partial differential
equation that finds application in various physical phenomena [22, 27]. In
particular, the propagation of weakly nonlinear internal long waves in a fluid
with a thin region of stratification can be represented by the BO equation.
Originating from the modeling of waves in shallow water, it offers insights
into the behavior of these waves, including their propagation and interaction.
Moreover, the BO equation emerges as a fundamental model in diverse fields
such as fluid dynamics, nonlinear optics, and plasma physics, reflecting its
broad relevance in understanding wave dynamics and nonlinear wave
interactions. Furthermore, we mention that it defines a Hamiltonian system,
and with the help of inverse scattering method (see [1]), families of
localized solitary wave solutions, called as soliton solutions [7], can be
obtained under the appropriate assumptions on the initial data. Since the BO
equation is completely integrable, it admits infinitely many conserved
quantities [7].

The investigation into the well-posedness of the Cauchy problem (1.1)
associated to the BO equation has been the subject of extensive research over
the years. Pioneering work in the local well-posedness was conducted by Iório
[21] for the initial data in $H^{s}(\mathbb{R})$, $s>3/2$ and making use of
the conserved quantities, the global well-posedness for data in
$H^{s}(\mathbb{R})$, $s\geq 2$ is demonstrated. Afterwards, the local well-
posedness result is improved by Ponce [29] for $H^{3/2}(\mathbb{R})$ data
along with the global result for any solution in $H^{s}(\mathbb{R})$, $s\geq
3/2$ by using the smoothing properties and energy estimates. Further
improvements in the local well-posedness were given by Koch et al. [25] for
data in $H^{s}(\mathbb{R})$, $s>5/4$ and using the argument given in [25],
Kenig at al. [23] extended the result for $s>9/8$. Tao in [34] also obtained
the global well-posedness in $H^{1}(\mathbb{R})$ by introducing the Gauge
transformation. The idea of Gauge transformation given in [34] was further
improved by Burq et al. [5] to carry out the local well-posedness to
$H^{s}(\mathbb{R})$ for $s>1/4$ and by Kenig et al. [20] to extend the same to
$H^{s}(\mathbb{R})$ for $s\geq 0$. Molinet [28] has also obtained the global
well-posedness for the periodic data in $L^{2}(\mathbb{R})$.

It is well-known that due to the effects of dispersion and nonlinear
convection, finding a reliable method for the BO equation is quite challenging
task. However, various numerical methods have been developed and are employed
in practice for the equation (1.1), although we refer only to the relevant
literature here. Thomée et al. [35] introduced a fully implicit finite
difference scheme and recently Dutta et al. [11] established the convergence
of the fully discrete Crank-Nicolson scheme including discretization of
Hilbert transform and Galtung [16, 17] devised a convergent Galerkin scheme.

The discontinuous Galerkin (DG) approach, a class of finite element methods,
was first introduced by Reed and Hill [30] in the framework of neutron
transport equation. It employs discontinuous piecewise polynomial spaces for
approximating both the solution and test functions. While this method demands
more computational resources due to the discontinuities at the interfaces, it
offers flexibility in selecting fluxes at these interfaces, ensuring high
accuracy and stability. Afterwards, a significant breakthrough emerged in DG
methodology given by Cockburn et al., for details, we refer to [9] and
references therein, where the DG method is applied to nonlinear conservation
laws. In their approach, the space DG discretization is incorporated with the
Runge-Kutta (RK) time discretization and hence this method is often called as
the Runge-Kutta discontinuous Galerkin (RKDG) method. The RKDG method is
highly parallelizable due to the local character of DG schemes and the
explicit time discretization. Nevertheless, equations containing higher-order
derivative terms may not be directly tackled using the DG method. In such
instances, higher-order derivatives might introduce instability and
inconsistency as emphasized in [18]. To address these issues, Bassi and Rebay
[4] made substantial progress by adapting the RKDG method for compressible
Navier-Stokes equations. In their approach, both the solution and its gradient
are treated as independent variables, marking a notable stride in the
evolution of DG methodology.

Cockburn and Shu in [10] introduced a generalization of the Bassi and Rebay
method called the local discontinuous Galerkin (LDG) method, stands as a
cornerstone in addressing intricate higher-order problems. The methodology
revolves around the transformation of equations into first-order systems,
achieved by introducing auxiliary variables to approximate lower derivatives.
The term “local” in LDG underscores the provisional nature of these auxiliary
variables, which can be locally eliminated. Success with LDG method relies on
designing numerical fluxes at the interfaces. These fluxes must be chosen to
ensure stability and allow for solving the extra variables locally. This
method is tailored to handle these equations efficiently and accurately.

The LDG method has been developed to deal with equations that have higher-
order derivative terms. For instance, Yan and Shu [41] devised a LDG method
for KdV type equations, which have third-order spatial derivatives. They
obtained the error estimates with order of convergence $k+1/2$ in the linear
case. Afterwards, Xu and Shu [38, 39] studied nonlinear equations and observed
that the LDG method still provides similar levels of accuracy and order of
convergence. Furthermore, Xu and Shu [40] also expanded the LDG method to
handle equations with fourth and fifth order spatial derivatives. Levy et al.
[26] also worked on adapting the LDG method for equations with compactly
supported traveling wave solutions appearing in nonlinear dispersive
equations.

In more recent times, the LDG method has become popular for dealing with
partial differential equations that involve the non-local operator. Xu and
Hesthaven [36] came up with an LDG method that breaks down the fractional
Laplacian of order $\alpha$ $(1<\alpha<2)$ into second-order derivatives and
fractional integrals of order $2-\alpha$. This method turned out to be very
effective, giving the optimal rates of $k+1$ in the linear case and $k+1/2$ in
the nonlinear setup. Similarly, Aboelenen [2] and Dwivedi et al. [15]
developed a LDG method specifically designed for fractional Schrödinger-type
equations and fractional Korteweg-de Vries equation respectively.

In times there has been advancement in developing a stable fully-discrete LDG
method known as the high-order RK-LDG method. Several research studies have
delved into assessing the stability and accuracy of this method in the context
of various problems, see for reference [31, 32, 33]. Very recently, Hunter et
al. [19] determined the stability using the Fourier method of fully-discrete
implicit-explicit RK Method for the linearized KdV equation endowed with
periodic initial data. Moreover, using the generalized Gauss-Radau projection,
Ai et al. [3] deduced an error estimates for the fully-discrete RK-LDG method
developed for solving the scalar nonlinear conservation laws. However, up to
our knowledge, LDG method has not been developed for the Benjamin-Ono
equation.

In this paper, our approach to design the LDG scheme for the BO equation (1.1)
involves introduction of auxiliary variables to represent (1.1) into the
system of lower order derivatives. Since the equation involves a non-local
operator, the original problem is defined over the entire real line. However,
for numerical approach, we restrict the problem to a sufficiently large
bounded domain. To design the scheme, we employ the usual DG method within
each element to all equations of the system. A critical aspect of this process
is the construction of appropriate numerical fluxes at the interior
interfaces. Boundary numerical fluxes are to be chosen from the provided
boundary conditions. The main ingredients of the paper are enlisted below:

  1. (1)

We design a LDG scheme for the BO equation (1.1) and establish the stability
of the devised semi-discrete scheme with a general nonlinear flux. The
stability analysis has also been extended for fully-discrete four-stage fourth
order RK scheme.

  2. (2)

We also carry out the convergence analysis of the semi-discrete scheme and
obtain a convergence rate of $\mathcal{O}(h^{k+1/2})$ considering any general
nonlinear flux.

  3. (3)

The theoretical convergence rates are validated through the numerical
experiments. We demonstrate the rates obtained by numerical illustrations are
optimal and it preserves the conserved quantity like mass and momentum in
discrete set up.

The rest of the paper is organized as follows. We commence our investigation
by introducing a few preliminary lemmas and semi-discrete LDG scheme in
Section 2. Section 3 outlines the stability analysis of the LDG scheme for
semi-discrete and fully-discrete LDG schemes incorporating the four-stage
fourth order Runge-Kutta scheme for time discretization. We present the
convergence analysis of the proposed scheme and prove the order of convergence
$\mathcal{O}(h^{k+1/2})$ for general nonlinear flux. The efficiency of the
scheme is validated through some numerical examples presented in Section 5.
Concluding remarks and a few remarks about future work are given in Section 6.

##  2\. Semi-discrete LDG scheme

###  2.1. Preliminary results

Hereby we describe a few relevant properties of the Hilbert transform through
the following lemma. It is worthwhile to mention that these properties are
instrumental for subsequent analysis.

######  Lemma 2.1.

(See [24, Chapter 15]) Let $\phi\in L^{2}(\mathbb{R})$ be a sufficiently
smooth function. Then the Hilbert transform $\mathcal{H}$ satisfies the
following properties:

  1. i)

Skew symmetric:

$(\mathcal{H}\phi_{1},\phi_{2})=-(\phi_{1},\mathcal{H}\phi_{2}),\qquad\forall\phi_{1},\phi_{2}\in
L^{2}(\mathbb{R}).$

  2. ii)

Commutes with derivatives:

$\mathcal{H}\phi_{x}=(\mathcal{H}\phi)_{x}.$

  3. iii)

$L^{2}$-isometry property:

$\left\|\mathcal{H}\phi\right\|_{L^{2}(\mathbb{R})}=\left\|\phi\right\|_{L^{2}(\mathbb{R})}.$

  4. iv)

Orthogonality:

$(\mathcal{H}\phi,\phi)=0.$

where $(\cdot,\cdot)$ is the standard $L^{2}$-inner product.

Given that the original problem is defined over the entire real line due to
the involvement of a non-local operator $\mathcal{H}$, but for the numerical
purpose, we restrict it to a bounded domain $\Omega:=[a,b]$, where $a<b$. We
assume that $U$ has compact support within $\Omega$. Hence it becomes
imperative to impose the boundary conditions:

$U(a,t)=0=U(b,t),\quad\text{for all }t<T.$

Moreover, the properties of the Hilbert transform introduced in Lemma 2.1
remain applicable for a bounded domain $\Omega$, provided $\phi$ has a compact
support within $\Omega$.

We partition the domain $\Omega$ into intervals
$I_{i}=(x_{i-\frac{1}{2}},x_{i+\frac{1}{2}})$ with
$a=x_{\frac{1}{2}}<x_{3/2}<\cdots<x_{N+\frac{1}{2}}=b$, where $N$ represents
the number of elements. This partition creates a mesh of elements denoted by
$\mathcal{I}$, with each element having a spatial step size
$h_{i}=x_{i+\frac{1}{2}}-x_{i-\frac{1}{2}}$ and a maximum step size
$h=\max\limits_{1\leq i\leq N}\\{h_{i}\\}$. In conjunction with this mesh, we
define the broken Sobolev spaces as follows:

$H^{1}(\Omega,\mathcal{I}):=\\{v:\Omega\to\mathbb{R}\,|\,v|_{I_{i}}\in
H^{1}(I_{i}),\,i=1,\cdots,N\\};$

and

$L^{2}(\Omega,\mathcal{I}):=\\{v:\Omega\to\mathbb{R}\,|\,v|_{I_{i}}\in
L^{2}(I_{i}),\,i=1,\cdots,N\\}.$

Within this framework, we introduce the notation $v_{i+\frac{1}{2}}$ to
represent the value of $v$ at the nodes $\\{x_{i+\frac{1}{2}}\\}$, and denote
the one-sided limits as

$v^{\pm}_{i+\frac{1}{2}}=v(x^{\pm}_{i+\frac{1}{2}}):=\lim\limits_{x\to
x^{\pm}_{i+\frac{1}{2}}}v(x).$

We define the local inner product and local $L^{2}(I_{i})$ norm as follows:

$(u,v)_{I_{i}}=\int_{I_{i}}uv\,dx,\quad\text{and
}\quad\|u\|_{I_{i}}=(u,u)_{I_{i}}^{\frac{1}{2}}.$

With all these preparation we introduce the auxiliary variables $P$ and $Q$
such that

$P=\mathcal{H}Q,\qquad Q=U_{x}.$

As a consequence, the equation (1.1) can be represented in the following
equivalent form of first order differential system

$\begin{split}U_{t}&=-(f(U)-P)_{x},\\\ P&=\mathcal{H}Q,\\\
Q&=U_{x}.\end{split}$ (2.1)

Prior to introducing the LDG scheme, we assume that the exact solution
$(U,P,Q)$ of the system (2.1) belongs to

$\mathcal{T}_{3}\times
K(\Omega,\mathcal{I}):=H^{1}(0,T;H^{1}(\Omega,\mathcal{I}))\times
L^{2}(0,T;H^{1}(\Omega,\mathcal{I}))\times
L^{2}(0,T;L^{2}(\Omega,\mathcal{I})).$

This implies that the solution $(U,P,Q)$ of (2.1) satisfies the following
system:

$\begin{split}\left(U_{t},v\right)_{I_{i}}&=\left(f(U)-P,v_{x}\right)_{I_{i}}-\left(fv-
Pv\right)|_{x_{i-\frac{1}{2}}^{+}}^{x_{i+\frac{1}{2}}^{-}},\\\
\left(P,w\right)_{I_{i}}&=\left(\mathcal{H}Q,w\right)_{I_{i}},\\\
\left(Q,z\right)_{I_{i}}&=-\left(U,z_{x}\right)_{I_{i}}+\left(Uz\right)|_{x_{i-\frac{1}{2}}^{+}}^{x_{i+\frac{1}{2}}^{-}},\end{split}$
(2.2)

for all $w\in L^{2}(\Omega,\mathcal{I})$, $v,z\in H^{1}(\Omega,\mathcal{I})$,
and for $i=1,\cdots,N$.

We define the finite dimensional discontinuous piecewise polynomial space
$V_{h}^{k}\subset H^{1}(\Omega,\mathcal{I})$ by

$V_{h}^{k}=\\{v\in L^{2}(\Omega):v|_{I_{i}}\in P^{k}(I_{i}),~{}\forall
i=1,2,\cdots,N\\},$ (2.3)

where $P^{k}(I_{i})$ is space of polynomials of degree up to order $k$ $(\geq
1)$ on $I_{i}$.

###  2.2. LDG scheme

Applying DG approach in all the equations of the above system (2.2), we design
the scheme as follows: we seek an approximations

$(u_{h},p_{h},q_{h})\in H^{1}(0,T;V_{h}^{k})\times L^{2}(0,T;V_{h}^{k})\times
L^{2}(0,T;V_{h}^{k})=:\mathcal{T}_{3}\times\mathcal{V}_{h}^{k}$

to $(U,P,Q)$, where $U$ is an exact solution of (1.1) with
$P=\mathcal{H}Q,~{}Q=U_{x},$ such that for all test functions
$(v,w,z)\in\mathcal{T}_{3}\times\mathcal{V}_{h}^{k}$ and $i=1,\cdots,N$, the
following system of equations holds:

$\begin{split}\left((u_{h})_{t},v\right)_{I_{i}}&=\left(\left(f(u_{h})-p_{h}\right),v_{x}\right)_{I_{i}}-\left(\hat{f}_{h}v-\hat{p}_{h}v\right)|_{x_{i-\frac{1}{2}}^{+}}^{x_{i+\frac{1}{2}}^{-}},\\\
\left(p_{h},w\right)_{I_{i}}&=\left(\mathcal{H}q_{h},w\right)_{I_{i}},\\\
\left(q_{h},z\right)_{I_{i}}&=-\left(u_{h},z_{x}\right)_{I_{i}}+\left(\hat{u}_{h}z\right)|_{x_{i-\frac{1}{2}}^{+}}^{x_{i+\frac{1}{2}}^{-}},\\\
\left(u_{h}^{0},v\right)_{I_{i}}&=\left(U_{0},v\right)_{I_{i}}.\end{split}$
(2.4)

To complete the LDG scheme (2.4), it is necessary to define the numerical
fluxes $\hat{u}_{h},$ $\hat{p}_{h}$ and the nonlinear flux $\hat{f}_{h}.$ We
introduce the following notations:

$\\{\\!\\!\\{u\\}\\!\\!\\}=\frac{u^{-}+u^{+}}{2},\qquad\llbracket
u\rrbracket=u^{+}-u^{-}.$

We choose the alternative numerical flux which is given by

$\displaystyle\hat{p}_{h}=p_{h}^{+},\qquad\hat{u}_{h}=u_{h}^{-},$ (2.5)

or alternatively

$\displaystyle\hat{p}_{h}=p_{h}^{-},\qquad\hat{u}_{h}=u_{h}^{+},$

at interface $x_{i+\frac{1}{2}}$, $i=1,2,\cdots,N-1$. We set the boundary flux
as

$\begin{split}(\hat{u}_{h})_{N+\frac{1}{2}}=U(b,t)=0,&\qquad(\hat{u}_{h})_{\frac{1}{2}}=U(a,t)=0,~{}\text{
for all }t<T,\\\
(\hat{p}_{h})_{N+\frac{1}{2}}=(p_{h})^{-}_{N+\frac{1}{2}},&\qquad(\hat{p}_{h})_{\frac{1}{2}}=(p_{h})^{+}_{\frac{1}{2}}.\\\
\end{split}$ (2.6)

For the nonlinear flux $\hat{f}_{h}$, we can use any monotone flux [41]. In
particular, we consider the following Lax-Friedrichs flux

$\hat{f}_{h}=\hat{f}(u_{h}^{-},u_{h}^{+})=\frac{1}{2}\big{(}f(u_{h}^{-})+f(u_{h}^{+})-\delta\llbracket
u_{h}\rrbracket\big{)},\qquad\delta=\max\limits_{u}|f^{\prime}(u)|,$ (2.7)

where the maximum is taken over a range of $u$ in a relevant element.

The proposed LDG scheme (2.4)-(2.6) for the BO equation works in the following
way: given $u_{h}$, we use the third equation of (2.4) to obtain $q_{h}$
locally; more precisely, $q_{h}$ in the cell $I_{i}$ can be computed with the
information of $u_{h}$ in the cells $I_{i-1}$ and $I_{i}$. Afterwards, with
the help of $q_{h}$ in the cell $I_{i}$, one can obtain $p_{h}$ locally in the
cell $I_{i}$. Finally, we update the approximate solution $u_{h}$ in the cell
$I_{i}$ incorporating $p_{h}$ and $u_{h}$ in the cells $I_{i-1}$ and $I_{i}$.
In a similar way, for the choice of alternative fluxes
$\hat{p}_{h}=p_{h}^{-},~{}\hat{u}_{h}=u_{h}^{+}$ the algorithm can be adopted
accordingly.

##  3\. Stability estimate

Hereby we analyze the stability and accuracy of the proposed LDG scheme
(2.4)-(2.6) for (1.1).

###  3.1. Stability of semi-discrete scheme

We prove the following stability lemma for general flux function using an
appropriate compact form:

######  Lemma 3.1.

($L^{2}$-stability) The LDG scheme (2.4)-(2.6) for Benjamin-Ono equation (1.1)
is $L^{2}$-stable. We have the following estimate

$\left\|u_{h}(\cdot,T)\right\|_{L^{2}(\Omega)}\leq
C\left\|u_{h}^{0}\right\|_{L^{2}(\Omega)},$ (3.1)

for any $T>0$ and a constant $C$.

###### Proof.

To carry out the stability estimate of the LDG scheme (2.4)-(2.6), we define
the corresponding compact form

$\displaystyle\mathcal{B}(u,p,q;v,w,z)$
$\displaystyle=\sum\limits_{i=1}^{N}\Big{[}\left(u_{t},v\right)_{I_{i}}-\left(f(u)-p,v_{x}\right)_{I_{i}}+\left(p,w\right)_{I_{i}}-\left(\mathcal{H}q,w\right)_{I_{i}}$
$\displaystyle\quad+\left(q,z\right)_{I_{i}}+\left(u,z_{x}\right)_{I_{i}}+\left((\hat{f}-\hat{p})v\right)|_{x_{i-\frac{1}{2}}^{+}}^{x_{i+\frac{1}{2}}^{-}}-\left(\hat{u}z\right)|_{x_{i-\frac{1}{2}}^{+}}^{x_{i+\frac{1}{2}}^{-}}\Big{]},$
(3.2)

for all $(u,p,q)\in\mathcal{T}_{3}\times K(\Omega,\mathcal{I})$ and
$(v,w,z)\in\mathcal{T}_{3}\times\mathcal{V}_{h}^{k}$. We represent the
numerical fluxes in the following way for subsequent analysis:

$\displaystyle\sum\limits_{i=1}^{N}(\hat{f}v)|_{x_{i-\frac{1}{2}}^{+}}^{x_{i+\frac{1}{2}}^{-}}=-\hat{f}_{\frac{1}{2}}v^{+}_{\frac{1}{2}}+\hat{f}_{N+\frac{1}{2}}v^{-}_{N+\frac{1}{2}}-\sum\limits_{i=1}^{N-1}\hat{f}_{i+\frac{1}{2}}\llbracket
v\rrbracket_{i+\frac{1}{2}},$
$\displaystyle\sum\limits_{i=1}^{N}(\hat{p}v)|_{x_{i-\frac{1}{2}}^{+}}^{x_{i+\frac{1}{2}}^{-}}=-p_{\frac{1}{2}}^{+}v^{+}_{\frac{1}{2}}+p_{N+\frac{1}{2}}^{-}v^{-}_{N+\frac{1}{2}}-\sum\limits_{i=1}^{N-1}p^{+}_{i+\frac{1}{2}}\llbracket
v\rrbracket_{i+\frac{1}{2}},$

and similarly

$\displaystyle\sum\limits_{i=1}^{N}(\hat{u}z)|_{x_{i-\frac{1}{2}}^{+}}^{x_{i+\frac{1}{2}}^{-}}=-\hat{u}_{\frac{1}{2}}z^{+}_{\frac{1}{2}}+\hat{u}_{N+\frac{1}{2}}z^{-}_{N+\frac{1}{2}}-\sum\limits_{i=1}^{N-1}u^{-}_{i+\frac{1}{2}}\llbracket
z\rrbracket_{i+\frac{1}{2}}=-\sum\limits_{i=1}^{N-1}u^{-}_{i+\frac{1}{2}}\llbracket
z\rrbracket_{i+\frac{1}{2}},$ (3.3)

where we have incorporated the boundary numerical fluxes defined in equations
(2.5)-(2.6). Using the above estimates, the compact form (3.1) can be
rewritten as

$\begin{split}\mathcal{B}(u,p,q;v,w,z)&=\sum\limits_{i=1}^{N}\Big{[}\left(u_{t},v\right)_{I_{i}}-\left(\left(f(u)-p\right),v_{x}\right)_{I_{i}}+\left(p,w\right)_{I_{i}}-\left(\mathcal{H}q,w\right)_{I_{i}}\\\
&\qquad+\left(q,z\right)_{I_{i}}+\left(u,z_{x}\right)_{I_{i}}\Big{]}+\mathcal{IF}(u,p;v,z),\end{split}$
(3.4)

where numerical flux at interfaces $\mathcal{IF}$ is given by

$\begin{split}\mathcal{IF}(u,p;v,z)&:=-\left(\hat{f}_{\frac{1}{2}}-p^{+}_{\frac{1}{2}}\right)v^{+}_{\frac{1}{2}}+\left(\hat{f}_{N+\frac{1}{2}}-p^{-}_{N+\frac{1}{2}}\right)v^{-}_{N+\frac{1}{2}}\\\
&\quad\quad-\sum\limits_{i=1}^{N-1}(\hat{f}_{i+\frac{1}{2}}-p^{+}_{i+\frac{1}{2}})\llbracket
v\rrbracket_{i+\frac{1}{2}}+\sum\limits_{i=1}^{N-1}u^{-}_{i+\frac{1}{2}}\llbracket
z\rrbracket_{i+\frac{1}{2}}.\end{split}$ (3.5)

By choosing $(v,w,z)=(u,-q,p)$ in (3.4), it yields

$\displaystyle\mathcal{B}(u,p,q;u,-q,p)$
$\displaystyle=\sum\limits_{i=1}^{N}\Big{[}\left(u_{t},u\right)_{I_{i}}-\left(f(u),u_{x}\right)_{I_{i}}+\left(p,u_{x}\right)_{I_{i}}$
$\displaystyle\quad-\left(p,q\right)_{I_{i}}+\left(\mathcal{H}q,q\right)_{I_{i}}+\left(q,p\right)_{I_{i}}+\left(u,p_{x}\right)_{I_{i}}\Big{]}+\mathcal{IF}(u,p;u,p)$
$\displaystyle=\sum\limits_{i=1}^{N}\Big{[}\left(u_{t},u\right)_{I_{i}}-\left(f(u),u_{x}\right)_{I_{i}}+\left(p,u_{x}\right)_{I_{i}}+\left(u,p_{x}\right)_{I_{i}}\Big{]}+\mathcal{IF}(u,p;u,p).$
(3.6)

Using the integration by parts, we obtain
$(p,u_{x})_{I_{i}}+(u,p_{x})_{I_{i}}=(up)|_{x_{i-\frac{1}{2}}^{+}}^{x_{i+\frac{1}{2}}^{-}},$
and consequently, the compact form reduces to

$\begin{split}\mathcal{B}(u,p,q;u,-q,p)=\sum\limits_{i=1}^{N}\Big{[}\left(u_{t},u\right)_{I_{i}}-\left(f(u),u_{x}\right)_{I_{i}}\Big{]}+\underbrace{\sum\limits_{i=1}^{N}(up)|_{x_{i-\frac{1}{2}}^{+}}^{x_{i+\frac{1}{2}}^{-}}+\mathcal{IF}(u,p;u,p)}_{\mathcal{E}_{2}}.\end{split}$
(3.7)

Furthermore, we observe that

$\sum\limits_{i=1}^{N}(up)|_{x_{i-\frac{1}{2}}^{+}}^{x_{i+\frac{1}{2}}^{-}}=-u^{+}_{\frac{1}{2}}p^{+}_{\frac{1}{2}}+u^{-}_{N+\frac{1}{2}}p^{-}_{N+\frac{1}{2}}-\sum\limits_{i=1}^{N-1}u^{-}_{i+\frac{1}{2}}\llbracket
p\rrbracket_{i+\frac{1}{2}}-\sum\limits_{i=1}^{N-1}p^{+}_{i+\frac{1}{2}}\llbracket
u\rrbracket_{i+\frac{1}{2}},$

and using the above identity, we deduce

$\begin{split}\mathcal{E}_{2}&=\sum\limits_{i=1}^{N}(up)|_{x_{i-\frac{1}{2}}^{+}}^{x_{i+\frac{1}{2}}^{-}}+\mathcal{IF}(u,p;u,p)\\\
&=-u^{+}_{\frac{1}{2}}p^{+}_{\frac{1}{2}}+u^{-}_{N+\frac{1}{2}}p^{-}_{N+\frac{1}{2}}-\sum\limits_{i=1}^{N-1}u^{-}_{i+\frac{1}{2}}\llbracket
p\rrbracket_{i+\frac{1}{2}}-\sum\limits_{i=1}^{N-1}p^{+}_{i+\frac{1}{2}}\llbracket
u\rrbracket_{i+\frac{1}{2}}\\\
&\quad-\left(\hat{f}_{\frac{1}{2}}-p^{+}_{\frac{1}{2}}\right)u^{+}_{\frac{1}{2}}+\left(\hat{f}_{N+\frac{1}{2}}-p^{-}_{N+\frac{1}{2}}\right)u^{-}_{N+\frac{1}{2}}-\sum\limits_{i=1}^{N-1}(\hat{f}_{i+\frac{1}{2}}-p^{+}_{i+\frac{1}{2}})\llbracket
u\rrbracket_{i+\frac{1}{2}}\\\
&\quad+\sum\limits_{i=1}^{N-1}u^{-}_{i+\frac{1}{2}}\llbracket
p\rrbracket_{i+\frac{1}{2}}\\\
&=-\hat{f}_{\frac{1}{2}}u^{+}_{\frac{1}{2}}+\hat{f}_{N+\frac{1}{2}}u^{-}_{N+\frac{1}{2}}-\sum\limits_{i=1}^{N-1}\hat{f}_{i+\frac{1}{2}}\llbracket
u\rrbracket_{i+\frac{1}{2}}.\end{split}$ (3.8)

Let us define $F(u)=\int^{u}f(u)\,du$. Then we have

$\sum\limits_{i=1}^{N}\left(f(u),u_{x}\right)_{I_{i}}=\sum\limits_{i=1}^{N}F(u)|_{u^{+}_{i-\frac{1}{2}}}^{u^{-}_{i+\frac{1}{2}}}=-\sum\limits_{i=1}^{N-1}\llbracket
F(u)\rrbracket_{i+\frac{1}{2}}-F(u)_{\frac{1}{2}}+F(u)_{N+\frac{1}{2}}.$ (3.9)

Incorporating (3.8) and (3.9) in equation (3.7), we have

$\displaystyle\mathcal{B}(u,p,q;u,-q,p)=\left(u_{t},u\right)_{L^{2}(\Omega)}$
$\displaystyle+\sum\limits_{i=1}^{N-1}\llbracket
F(u)\rrbracket_{i+\frac{1}{2}}+F(u)_{\frac{1}{2}}-F(u)_{N+\frac{1}{2}}$
$\displaystyle-\hat{f}_{\frac{1}{2}}u^{+}_{\frac{1}{2}}+\hat{f}_{N+\frac{1}{2}}u^{-}_{N+\frac{1}{2}}-\sum\limits_{i=1}^{N-1}\hat{f}_{i+\frac{1}{2}}\llbracket
u\rrbracket_{i+\frac{1}{2}}.$ (3.10)

Since $(u_{h},p_{h},q_{h})$ satisfies the LDG scheme (2.4)-(2.6), we have

$\mathcal{B}(u_{h},p_{h},q_{h};v,w,z)=0,$

for any $(v,w,z)\in V_{h}^{k}$. As a result, we have

$\displaystyle\left((u_{h})_{t},u_{h}\right)_{L^{2}(\Omega)}$
$\displaystyle+\sum\limits_{i=1}^{N-1}\llbracket
F(u_{h})\rrbracket_{i+\frac{1}{2}}+F(u_{h})_{\frac{1}{2}}-F(u_{h})_{N+\frac{1}{2}}$
$\displaystyle-(\hat{f}_{h})_{\frac{1}{2}}(u_{h})^{+}_{\frac{1}{2}}+(\hat{f}_{h})_{N+\frac{1}{2}}(u_{h})^{-}_{N+\frac{1}{2}}-\sum\limits_{i=1}^{N-1}(\hat{f}_{h})_{i+\frac{1}{2}}\llbracket
u_{h}\rrbracket_{i+\frac{1}{2}}=0.$

Since the numerical flux $\hat{f}_{h}=\hat{f}(u_{h}^{-},u_{h}^{+})$ is
monotone, it is non-decreasing in its first argument and non-increasing in its
second argument. As a consequence, we have

$\llbracket
F(u_{h})\rrbracket_{i+\frac{1}{2}}-\hat{f}_{i+\frac{1}{2}}\llbracket
u_{h}\rrbracket_{i+\frac{1}{2}}>0,$

for all $i=1,2,\cdots,N-1$. Dropping the positive term from left hand side of
equation (3.1) and using the boundary conditions, we end up with

$\begin{split}\frac{1}{2}\frac{d}{dt}\left\|u_{h}\right\|^{2}_{L^{2}(\Omega)}\leq
0.\end{split}$ (3.11)

Applying the Gronwall’s inequality, we obtain

$\left\|u_{h}(\cdot,T)\right\|_{L^{2}(\Omega)}\leq
C\left\|u_{h}^{0}\right\|_{L^{2}(\Omega)}.$

Hence the result follows. ∎

###  3.2. Stability of fully-discrete scheme

Hereby we concentrate on the fully-discrete scheme of (1.1). In particular, we
derive stability estimates for an explicit fourth-order Runge-Kutta time
discretization assuming that the LDG scheme remains stable in spatial
discretization. It is worth noting that proving stability under higher-order
time discretizations is quite challenging. Currently, we have focused only on
the linear case, specifically choosing $f(U)=0$ for simplicity. We follow the
similar approach presented in [33] and it was also analyzed for the fractional
KdV equation [15].

We define the associated LDG operators $D^{h}_{\pm}$ and $D^{h}_{\mathcal{H}}$
for approximating first order derivative and nonlocal terms in the integral
form

$\begin{split}(D^{h}_{\pm}v,w)&=-(v,w_{x})-\sum_{i=1}^{N-1}v^{\pm}_{i+\frac{1}{2}}\llbracket
w\rrbracket_{i+\frac{1}{2}},\qquad\forall v,w\in V_{h}^{k},\\\
(D^{h}_{\mathcal{H}}v,w)&=\left(\mathcal{H}v,w\right),\qquad\forall v,w\in
V_{h}^{k}.\end{split}$ (3.12)

We are assuming either $v^{+}_{\frac{1}{2}}=0=w^{-}_{N+\frac{1}{2}}$ or
$v^{-}_{N+\frac{1}{2}}=0=w^{+}_{\frac{1}{2}}$ in the above definitions.
Moreover, the LDG operators have the following properties.

######  Proposition 3.2 (Antisymmetry [15]).

The LDG operators $D^{h}_{\pm}$ and $D^{h}_{\mathcal{H}}$ are antisymmetric.
In particular, we have

$(D^{h}_{+}e,g)=-(e,D^{h}_{-}g),\qquad\forall e,g\in V_{h}^{k},$

and we denote ${D^{h}_{+}}^{T}=-D^{h}_{-}$.

Using the LDG operators (3.12), the semi-discrete scheme (2.4)-(2.6) can be
represented as

$\left(\frac{d}{dt}u_{h},v_{h}\right)=(L_{h}u_{h},v_{h}),\qquad\forall
v_{h}\in V_{h}^{k},$

where $L_{h}=-D^{h}_{+}D^{h}_{\mathcal{H}}D^{h}_{-}$. Then the semi-discrete
scheme (2.4)-(2.6) corresponds to an autonomous ODE system

$\frac{d}{dt}u_{h}=L_{h}u_{h}.$ (3.13)

With this, we consider the standard four stage explicit fourth order RK-LDG
scheme which can be presented as follows:

$\displaystyle u_{h}^{n+1}=P_{4}(\Delta tL_{h})u^{n}_{h},$ (3.14)

where the operator $P_{4}$ is defined by

$\displaystyle P_{4}(\Delta tL_{h})=I+\Delta tL_{h}+\frac{1}{2}(\Delta
tL_{h})^{2}+\frac{1}{6}(\Delta tL_{h})^{3}+\frac{1}{24}(\Delta tL_{h})^{4}.$

We say the operator $L_{h}$ is semi-negative if
$(v_{h},(L_{h}+L_{h}^{T})v_{h}))\leq 0$, $\forall v_{h}\in V_{h}^{k}$, and we
denote it by $L_{h}+L_{h}^{T}\leq 0$. We define $\mathcal{L}_{h}:=\Delta
tL_{h}$ and a bilinear form on $V_{h}^{k}$ as

$[u_{h},v_{h}]:=-\big{(}u_{h},(L_{h}+L_{h}^{T})v_{h}\big{)}.$

######  Definition 3.3 (Strongly stable and monotonically stable).

Let $u_{h}^{n}$ be the approximate solution obtained by the fully-discrete
scheme (3.14) and $u^{0}_{h}$ is discretized initial data (cf. (2.4)). Then
the scheme (3.14) is strongly stable if there is an integer $n_{0}$ such that

$\left\|u_{h}^{n}\right\|_{L^{2}(\Omega)}\leq\left\|u_{h}^{0}\right\|_{L^{2}(\Omega)},\qquad\forall
n\geq n_{0}.$ (3.15)

Moreover, the scheme (3.14) is monotonically stable if we have

$\left\|u_{h}^{n+1}\right\|_{L^{2}(\Omega)}\leq\left\|u_{h}^{n}\right\|_{L^{2}(\Omega)},\qquad\forall
n\geq 0,$ (3.16)

provided CFL number is sufficiently small. Clearly, monotonically stable
implies the strongly stable.

We state the following result related to energy equality. This equality
provides the information of the solution after one time step.

######  Lemma 3.4 (Energy equality [33]).

Let $u_{h}^{n}$ be the solution of fully-discrete scheme (3.14). Then there
holds

$\left\|u_{h}^{n+1}\right\|^{2}_{L^{2}(\Omega)}-\left\|u_{h}^{n}\right\|^{2}_{L^{2}(\Omega)}=\mathcal{Q}(u_{h}^{n}),\qquad\forall
n\geq 1,$

where

$\mathcal{Q}(u_{h}^{n})=\frac{1}{576}\left\|\mathcal{L}_{h}^{4}u_{h}^{n}\right\|^{2}_{L^{2}(\Omega)}-\frac{1}{72}\left\|\mathcal{L}_{h}^{3}u_{h}^{n}\right\|^{2}_{L^{2}(\Omega)}+\Delta
t\sum\limits_{i,j=0}^{3}\alpha_{ij}[\mathcal{L}_{h}^{i}u_{h}^{n},\mathcal{L}_{h}^{j}u_{h}^{n}],$

and

$A=\left(\alpha_{ij}\right)_{i,j=0}^{3}=-\begin{pmatrix}1&1/2&1/6&1/24\\\
1/2&1/3&1/8&1/24\\\ 1/6&1/8&1/24&1/48\\\ 1/24&1/24&1/48&1/144\end{pmatrix}.$

The next lemma establishes the negativity of the quadratic form. Consequently,
it further provides conditions for strong stability.

######  Lemma 3.5.

(See [33, Lemma 2.4]) Let $L_{h}$ be a semi-negative operator and

$\mathcal{Q}_{1}(u)=\zeta\left\|\mathcal{L}_{h}^{3}(u)\right\|_{L^{2}(\Omega)}+\Delta
t\sum\limits_{i,j=0}^{m}\bar{\alpha}_{ij}[\mathcal{L}_{h}^{i}u,\mathcal{L}_{h}^{j}u],$
(3.17)

where $\bar{\alpha}_{ij}=\bar{\alpha}_{ji}$, $m\geq 2$ and
$\mathcal{L}_{h}=\Delta tL_{h}$. If $\zeta<0$ and
$\bar{A}=\left(\bar{\alpha}\right)_{i,j=0}^{2}$ is negative definite, then
there exists a constant $c_{0}>0$, such that $\mathcal{Q}_{1}(u)\leq 0$,
provided $\left\|\mathcal{L}\right\|\leq c_{0}$.

Since the semi-discrete scheme (3.13) honours the spatial stability from Lemma
3.1, that is

$\left(\frac{d}{dt}u_{h},u_{h}\right)=(L_{h}u_{h},u_{h})\leq 0,$

the LDG operator $L_{h}$ is semi-negative [32], i.e. $L_{h}$ satisfies:
$L_{h}+L_{h}^{T}\leq 0$.

The stability analysis of the classical fourth-order RK time discretization of
equation (3.13) has been a topic of considerable interest. Recent
investigations by Sun and Shu [33] have shed light on this matter, where they
introduced a counterexample demonstrating that this method may lack strong
stability for semi-negative operators. However, it is essential to highlight
that the semi-negative operator featured in the counterexample provided by Sun
and Shu [33] is not a DG operator. Furthermore, an additional counterexample
was presented in [37, Example 1 in Section 6], illustrating the absence of
monotonicity stability for semi-negative DG operators within this method.
Moreover, in [33], it was demonstrated that applying the classical fourth-
order Runge-Kutta scheme for two and three consecutive time steps, which can
be viewed as an eight-stage and twelve-stage respectively, fourth-order Runge-
Kutta method, is strongly stable for semi-negative operators.

######  Theorem 3.6 (Two step strong stability [33]).

Let $L_{h}+L_{h}^{T}\leq 0$. Then the four-stage fourth order RK-LDG scheme
(3.14) is strongly stable in two steps and the following estimate holds

$\left\|u_{h}^{n+2}\right\|_{L^{2}(\Omega)}\leq\left\|u_{h}^{n}\right\|_{L^{2}(\Omega)},$

provided $\Delta t\left\|L_{h}\right\|\leq c_{0}$, where $c_{0}$ is a constant
and $\left\|\cdot\right\|$ is an operator norm.

######  Theorem 3.7 (Three step strong stability [15]).

Let $L_{h}+L_{h}^{T}\leq 0$. Then the four-stage fourth order RK-LDG scheme
(3.14) is strongly stable in three steps. Therefore, we have

$\left\|u_{h}^{n+3}\right\|_{L^{2}(\Omega)}\leq\left\|u_{h}^{n}\right\|_{L^{2}(\Omega)},\qquad
n\geq 0,$

provided $\Delta t\left\|L_{h}\right\|\leq c_{0}$, where $c_{0}$ is a
constant.

Stability in two and three steps of the fully discrete LDG scheme (3.14)
combined together proves the following theorem.

######  Theorem 3.8.

Let $L_{h}+L_{h}^{T}\leq 0$. Then the four-stage fourth order RK-LDG scheme
(3.14) is strongly stable. That is, we have

$\left\|u_{h}^{n}\right\|_{L^{2}(\Omega)}\leq\left\|u_{h}^{0}\right\|_{L^{2}(\Omega)},\qquad\forall
n\geq 2,$

provided $\Delta t\left\|L_{h}\right\|\leq c_{0}$, where $c_{0}$ is a
constant.

##  4\. Convergence analysis

To proceed with the error estimates, we define the projection operators into
$V_{h}^{k}$ as follows. For any $g\in L^{2}(\Omega)$, which is sufficiently
smooth, we define

$\begin{split}\int_{I_{i}}\big{(}\mathcal{P}^{-}g(x)-g(x)\big{)}y(x)\,dx&=0\quad\forall~{}y\in
P^{k-1}(I_{i}),\quad i=1,2,\cdots,N,~{}\text{ and
}~{}(\mathcal{P}^{-}u)^{-}_{i+1/2}=u(x^{-}_{i+1/2}),\\\
\int_{I_{i}}\big{(}\mathcal{P}g(x)-g(x)\big{)}y(x)\,dx&=0\quad\forall~{}y\in
P^{k}(I_{i}),\quad i=1,2,\cdots,N.\end{split}$ (4.1)

Here $\mathcal{P}^{-}$ is a special projection, and $\mathcal{P}$ is the
standard $L^{2}$ projection. Let $U$ be an exact solution of (1.1) and $u_{h}$
be an approximate solution obtained by the LDG scheme (2.4)-(2.6), we denote

$\mathcal{P}^{-}_{h}u=\mathcal{P}^{-}U-u_{h},\quad\mathcal{P}_{h}q=\mathcal{P}Q-q_{h},$

and

$\mathcal{P}^{-}_{e}U=\mathcal{P}^{-}U-U,\quad\mathcal{P}_{e}Q=\mathcal{P}Q-Q.$

In the process of formulating the error estimate for the equation (1.1), we
introduce few lemmas concerning the measure between physical flux $f$ and
numerical flux $\hat{f}_{h}$.

######  Lemma 4.1 (See Lemma 3.1 in [42]).

Let $\xi\in L^{2}(\Omega)$ be any piecewise smooth function. On each
interfaces of elements and on the boundary points we define

$\displaystyle\beta(\hat{f};\xi):=\beta(\hat{f};\xi^{-},\xi^{+})=\begin{cases}\llbracket\xi\rrbracket^{-1}(f(\\{\\!\\!\\{\xi\\}\\!\\!\\})-\hat{f}(\xi)),\qquad&if~{}\llbracket\xi\rrbracket\neq
0,\\\
\frac{1}{2}|f^{\prime}(\\{\\!\\!\\{\xi\\}\\!\\!\\})|,\qquad&if~{}\llbracket\xi\rrbracket=0,\end{cases}$

where $\hat{f}(\xi)=\hat{f}(\xi^{-},\xi^{+})$ is a monotone numerical flux
consistent with the given physical flux $f$. Then $\beta(\hat{f};\xi)$ is
bounded and nonnegative for any $\xi^{-},\xi^{+}\in\mathbb{R}$. Moreover, we
have

$\displaystyle\frac{1}{2}|f^{\prime}(\\{\\!\\!\\{\xi\\}\\!\\!\\})|$
$\displaystyle\leq\beta(\hat{f};\xi)+C_{\ast}|\llbracket\xi\rrbracket|,$
$\displaystyle-\frac{1}{8}f^{\prime\prime}(\\{\\!\\!\\{\xi\\}\\!\\!\\})|\llbracket\xi\rrbracket$
$\displaystyle\leq\beta(\hat{f};\xi)+C_{\ast}|\llbracket\xi\rrbracket|^{2}.$

Borrowing the idea from [39], we find the estimate for nonlinear part $f(u)$,
by defining

$\displaystyle\sum\limits_{i=1}^{N}\mathcal{F}_{i}(f;U,u_{h},v)$
$\displaystyle=\sum\limits_{i=1}^{N}\int_{I_{i}}\big{(}f(U)-f(u_{h})\big{)}v_{x}\,dx+\sum\limits_{i=1}^{N}\Big{(}\big{(}f(U)-f(\\{\\!\\!\\{u_{h}\\}\\!\\!\\})\big{)}\llbracket
v\rrbracket\Big{)}_{i+\frac{1}{2}}$
$\displaystyle\quad+\sum\limits_{i=1}^{N}\Big{(}(f(\\{\\!\\!\\{u_{h}\\}\\!\\!\\})-\hat{f})\llbracket
v\rrbracket\Big{)}_{i+\frac{1}{2}}.$ (4.2)

######  Lemma 4.2 (See Corollary 3.6 in [39]).

Let the operator $\mathcal{F}_{i}$, defined by (4.2). Then we have the
following estimate:

$\displaystyle\sum\limits_{i=1}^{N}\mathcal{F}_{i}(f;U,u_{h},v)$
$\displaystyle\leq-\frac{1}{4}\beta(\hat{f};u_{h})\sum\limits_{i=1}^{N}\llbracket
v\rrbracket^{2}_{i+\frac{1}{2}}+\big{(}C+C_{\ast}(\left\|v\right\|_{L^{\infty}(\Omega)}+h^{-1}\left\|U-u_{h}\right\|_{L^{\infty}(\Omega)}^{2})\big{)}\left\|v\right\|_{L^{2}(\Omega)}^{2}$
$\displaystyle\quad+\big{(}C+C_{\ast}h^{-1}\left\|U-u_{h}\right\|_{L^{\infty}(\Omega)}^{2}\big{)}h^{2k+1}.$
(4.3)

We deal with the nonlinear flux $f(u)$ by making an $a~{}priori$ assumption
[39]. Let $h$ be small enough and $k\geq 1$, then there holds

$\left\|U-u_{h}\right\|_{L^{2}(\Omega)}\leq h.$ (4.4)

The above assumption is unnecessary for linear flux $f(u)=u$. We define the
bilinear operator $\mathcal{B}_{0}$ by the following

$\displaystyle\mathcal{B}_{0}(u,p,q;v,w,z)$
$\displaystyle=\sum\limits_{i=1}^{N}\Big{[}\left(u_{t},v\right)_{I_{i}}+\left(p,v_{x}\right)_{I_{i}}+\left(p,w\right)_{I_{i}}-\left(\mathcal{H}q,w\right)_{I_{i}}$
$\displaystyle\quad+\left(q,z\right)_{I_{i}}+\left(u,z_{x}\right)_{I_{i}}\Big{]}+\mathcal{IF}_{0}(u,p;v,z),$
(4.5)

where the term $\mathcal{IF}_{0}$ is given by

$\begin{split}\mathcal{IF}_{0}(u,p;v,z):=p^{+}_{\frac{1}{2}}v^{+}_{\frac{1}{2}}-p^{-}_{N+\frac{1}{2}}v^{-}_{N+\frac{1}{2}}+\sum\limits_{i=1}^{N-1}p^{+}_{i+\frac{1}{2}}\llbracket
v\rrbracket_{i+\frac{1}{2}}+\sum\limits_{i=1}^{N-1}u^{-}_{i+\frac{1}{2}}\llbracket
z\rrbracket_{i+\frac{1}{2}}.\end{split}$ (4.6)

Note that $\mathcal{B}_{0}=\mathcal{B}$ and $\mathcal{IF}_{0}=\mathcal{IF}$,
if we take $f=0$ in the definition of $\mathcal{B}$, that is,
$\mathcal{B}_{0}$ is the linear part of $\mathcal{B}$. Incorporating (3.1) in
(4), we have

$\begin{split}\mathcal{B}_{0}(u,p,q;u,-q,p)&=\left(u_{t},u\right)_{L^{2}(\Omega)}.\end{split}$
(4.7)

######  Theorem 4.3.

Let $U\in C^{k+1}(\Omega)$ be a sufficiently smooth exact solution of (1.1).
Assume that the nonlinear flux $f\in C^{3}(\Omega)$. Let $u_{h}$ be an
approximate solution obtained by the LDG scheme (2.4)-(2.6) with auxiliary
variables $p_{h}$ and $q_{h}$. Let $V_{h}^{k}$ be a space of piecewise
polynomials of degree $k\geq 1$ defined by (2.3). Then for small enough $h$,
there holds the following error estimate

$\left\|U-u_{h}\right\|_{L^{2}(\Omega)}\leq Ch^{k+1/2},$ (4.8)

where $C$ is a constant depending on the fixed time $T>0$, $k$ and the bounds
on the derivatives $|f^{(m)}|$, $m=1,2,3$.

###### Proof.

We begin by deriving an error equation. Since $U$ is an exact solution of
(1.1), we define

$\displaystyle P=Q,\qquad Q=U_{x}.$

Then $U$, $P$, and $Q$ satisfy the equation (3.1). Hence for any $(v,w,z)\in
V_{h}^{k}$, we have

$\displaystyle\mathcal{B}(U,P,Q;v,w,z)=\mathcal{B}(u_{h},p_{h},q_{h};v,w,z)=0,$

where $(u_{h},p_{h},q_{h},r_{h})$ is an approximate solution obtained by the
scheme (3.1). By incorporating the bilinear operator $\mathcal{B}_{0}$, we get

$\displaystyle 0$
$\displaystyle=\mathcal{B}(U,P,Q;v,w,z)-\mathcal{B}(u_{h},p_{h},q_{h};v,w,z)$
$\displaystyle=\mathcal{B}_{0}(U,P,Q;v,w,z)-\mathcal{B}_{0}(u_{h},p_{h},q_{h};v,w,z)-\sum\limits_{i=1}^{N}\mathcal{F}_{i}(f;U,u_{h},v)$
$\displaystyle=\mathcal{B}_{0}(U-u_{h},P-p_{h},Q-q_{h},R-r_{h};v,w,z)-\sum\limits_{i=1}^{N}\mathcal{F}_{i}(f;U,u_{h},v).$

Taking into account the projection operators $\mathcal{P}^{-}$ and
$\mathcal{P}$ defined in (4.1), we choose
$(v,w,z)=(\mathcal{P}^{-}_{h}u,-\mathcal{P}_{h}q,\mathcal{P}_{h}p)$. Since
$U-u_{h}=\mathcal{Q}_{h}u-\mathcal{Q}_{e}U$, where
$\mathcal{Q}_{h}=\mathcal{P}_{h}$ or $\mathcal{P}_{h}^{-}$ and
$\mathcal{Q}_{e}=\mathcal{P}_{e}$ or $\mathcal{P}_{e}^{-}$, we have

$\begin{split}\mathcal{B}_{0}(\mathcal{P}^{-}_{h}u,\mathcal{P}_{h}p,\mathcal{P}_{h}q;\mathcal{P}^{-}_{h}u,-\mathcal{P}_{h}q,\mathcal{P}_{h}p)&=\mathcal{B}_{0}(\mathcal{P}^{-}_{e}U,\mathcal{P}_{e}P,\mathcal{P}_{e}Q;\mathcal{P}^{-}_{h}u,-\mathcal{P}_{h}q,\mathcal{P}_{h}p)\\\
&\quad+\sum\limits_{i=1}^{N}\mathcal{F}_{i}(f;U,u_{h},\mathcal{P}^{-}_{h}u).\end{split}$
(4.9)

We estimate the first term on the right-hand side of (4.9). From equation (4),
we have

$\begin{split}\mathcal{B}_{0}(\mathcal{P}^{-}_{e}U,\mathcal{P}_{e}P,\mathcal{P}_{e}Q;\mathcal{P}^{-}_{h}u,-\mathcal{P}_{h}q,\mathcal{P}_{h}p)&=\sum\limits_{i=1}^{N}\Big{[}\left((\mathcal{P}^{-}_{e}U)_{t},\mathcal{P}^{-}_{h}u\right)_{I_{i}}+\left(\mathcal{P}_{e}P,(\mathcal{P}^{-}_{h}u)_{x}\right)_{I_{i}}\\\
&\quad-\left(\mathcal{P}_{e}P,\mathcal{P}_{h}q\right)_{I_{i}}+\left(\mathcal{H}\mathcal{P}_{e}Q,\mathcal{P}_{h}q\right)_{I_{i}}\\\
&\quad+\left(\mathcal{P}_{e}Q,\mathcal{P}_{h}p\right)_{I_{i}}+\left(\mathcal{P}_{e}U,(\mathcal{P}_{h}p)_{x}\right)_{I_{i}}\Big{]}\\\
&\quad+\mathcal{IF}_{0}(\mathcal{P}^{-}_{e}U,\mathcal{P}_{e}P;\mathcal{P}^{-}_{h}u,\mathcal{P}_{h}p).\end{split}$
(4.10)

Since we have

$\displaystyle(\mathcal{P}_{h}p)_{x}\in
P^{k-1}(I_{i}),~{}(\mathcal{P}^{-}_{h}u)_{x}\in
P^{k-1}(I_{i}),\quad\text{and}\quad\mathcal{P}_{h}p,~{}\mathcal{P}_{h}q\in
P^{k}(I_{i}),$

and consequently, from the projection properties defined in (4.1) imply

$\displaystyle\left(\mathcal{P}_{e}P,(\mathcal{P}^{-}_{h}u)_{x}\right)_{I_{i}}=0,$
$\displaystyle\quad\left(\mathcal{P}_{e}Q,\mathcal{P}_{h}p\right)_{I_{i}}=0,\quad\left(\mathcal{P}_{e}P,\mathcal{P}_{h}q\right)_{I_{i}}=0,$
$\displaystyle\left(\mathcal{P}^{-}_{e}U,(\mathcal{P}_{h}p)_{x}\right)_{I_{i}}=0,$
$\displaystyle\quad\left(\mathcal{H}\mathcal{P}_{e}Q,\mathcal{P}_{h}q\right)_{I_{i}}=0,$

for all $i=1,2,\cdots,N$, and $(\mathcal{P}^{-}_{e}U)^{-}_{i+\frac{1}{2}}=0$
for all $i=1,2,\cdots,N-1$. Also note that, from the approximation theory on
the point values associated with projection operators [8, Section 3.2], we
have

$\displaystyle\left(\mathcal{P}_{e}P\right)^{+}_{i+\frac{1}{2}}\leq
Ch^{k+1},\quad(\mathcal{P}^{-}_{e}U)^{-}_{i-\frac{1}{2}}\leq Ch^{k+1},$

for all $i=1,2,\cdots,N-1$, and

$\displaystyle(\mathcal{P}^{-}_{h}u)^{-}_{N+\frac{1}{2}}=0,\quad(\mathcal{P}^{-}_{h}u)^{+}_{\frac{1}{2}}=0,$

by using boundary conditions. Combining these error bounds, we have

$\displaystyle\mathcal{IF}_{0}(\mathcal{P}^{-}_{e}U,\mathcal{P}_{e}P;\mathcal{P}^{-}_{h}u,\mathcal{P}_{h}p)=$
$\displaystyle(\mathcal{P}_{e}P)^{+}_{\frac{1}{2}}(\mathcal{P}^{-}_{h}u)^{+}_{\frac{1}{2}}-(\mathcal{P}_{e}P)^{-}_{N+\frac{1}{2}}(\mathcal{P}^{-}_{h}u)^{-}_{N+\frac{1}{2}}$
$\displaystyle+\sum\limits_{i=1}^{N-1}(\mathcal{P}_{e}P)^{+}_{i+\frac{1}{2}}\llbracket\mathcal{P}^{-}_{h}u\rrbracket_{i+\frac{1}{2}}+\sum\limits_{i=1}^{N-1}(\mathcal{P}^{-}_{e}U)^{-}_{i+\frac{1}{2}}\llbracket\mathcal{P}_{h}p\rrbracket_{i+\frac{1}{2}}$
$\displaystyle\leq$
$\displaystyle\sum\limits_{i=1}^{N-1}\Big{(}C(\epsilon)\big{(}(\mathcal{P}_{e}P)^{+}_{i+\frac{1}{2}}\big{)}^{2}+\varepsilon\llbracket\mathcal{P}^{-}_{h}u\rrbracket_{i+\frac{1}{2}}^{2}\Big{)}$
$\displaystyle\leq$ $\displaystyle
C(\Omega)h^{2k+1}+\varepsilon\sum\limits_{i=1}^{N-1}\llbracket\mathcal{P}^{-}_{h}u\rrbracket_{i+\frac{1}{2}}^{2}.$

Using the above estimates in (4.10) can be written as

$\begin{split}\mathcal{B}_{0}(\mathcal{P}^{-}_{e}U,\mathcal{P}^{+}_{e}P,\mathcal{P}_{e}Q;\mathcal{P}^{-}_{h}u,-\mathcal{P}_{h}q,\mathcal{P}_{h}p)&\leq\left((\mathcal{P}^{-}_{e}U)_{t},\mathcal{P}^{-}_{h}u\right)_{L^{2}(\Omega)}+C(\Omega)h^{2k+1}+\varepsilon\sum\limits_{i=1}^{N-1}\llbracket\mathcal{P}^{-}_{h}u\rrbracket_{i+\frac{1}{2}}^{2}.\end{split}$
(4.11)

Incorporating estimates form (4.7) and (4.11) in (4.9), we have

$\displaystyle\left((\mathcal{P}^{-}_{h}u)_{t},\mathcal{P}^{-}_{h}u\right)_{L^{2}(\Omega)}\leq$
$\displaystyle\left((\mathcal{P}^{-}_{e}U)_{t},\mathcal{P}^{-}_{h}u\right)_{L^{2}(\Omega)}+C(\Omega)h^{2k+1}+\varepsilon\sum\limits_{i=1}^{N-1}\llbracket\mathcal{P}^{-}_{h}u\rrbracket_{i+\frac{1}{2}}^{2}$
$\displaystyle\quad+\sum\limits_{i=1}^{N}\mathcal{F}_{i}(f;U,u_{h},\mathcal{P}^{-}_{h}u).$
(4.12)

Hence estimate (4.2) and equation (4) imply

$\begin{split}&\left((\mathcal{P}^{-}_{h}u)_{t},\mathcal{P}^{-}_{h}u\right)_{L^{2}(\Omega)}+\frac{1}{4}\beta(\hat{f}_{h};\mathcal{P}^{-}_{h}u)\sum\limits_{i=1}^{N}\llbracket\mathcal{P}^{-}_{h}u\rrbracket^{2}_{i+\frac{1}{2}}\\\
&\qquad\leq\left((\mathcal{P}^{-}_{e}U)_{t},\mathcal{P}^{-}_{h}u\right)_{L^{2}(\Omega)}+C(\Omega)h^{2k+1}+\Big{(}C+C_{\ast}h^{-1}\left\|U-u_{h}\right\|_{L^{\infty}(\Omega)}^{2}\Big{)}h^{2k+1}\\\
&\qquad\qquad+\Big{(}C+C_{\ast}\big{(}\left\|\mathcal{P}^{-}_{h}u\right\|_{L^{\infty}(\Omega)}+h^{-1}\left\|U-u_{h}\right\|_{L^{\infty}(\Omega)}^{2}\big{)}\Big{)}\left\|\mathcal{P}^{-}_{h}u\right\|_{L^{2}(\Omega)}^{2}+\varepsilon\sum\limits_{i=1}^{N-1}\llbracket\mathcal{P}^{-}_{h}u\rrbracket_{i+\frac{1}{2}}^{2}.\end{split}$
(4.13)

Utilizing the inverse property $\left\|u\right\|_{L^{\infty}(\Omega)}\leq
h^{-1/2}\left\|u\right\|_{L^{2}(\Omega)}$ and a priori assumption (4.4), we
obtain the estimate

$h^{-1}\left\|U-u_{h}\right\|_{\infty}^{2}h^{2k+1}\leq
h^{-2}\left\|U-u_{h}\right\|^{2}_{L^{2}(\Omega)}h^{2k+1}\leq h^{2k+1}.$

Using the above estimate and the positivity of $\beta$ from Lemma 4.1,
equation (4.13) implies

$\displaystyle\frac{1}{2}\frac{d}{dt}\left\|\mathcal{P}^{-}_{h}u\right\|^{2}_{L^{2}(\Omega)}\leq\left((\mathcal{P}^{-}_{e}U)_{t},\mathcal{P}^{-}_{h}u\right)_{L^{2}(\Omega)}+C_{\Omega}h^{2k+1}+C\left\|\mathcal{P}^{-}_{h}u\right\|^{2}_{L^{2}(\Omega)}.$

Using the standard approximation theory associated to the projection [8,
Section 3.2], we have
$\left\|\mathcal{P}^{-}_{h}u(\cdot,0)\right\|^{2}_{L^{2}(\Omega)}=0$. Finally,
with the help of Gronwall’s inequality, we obtain the estimate

$\left\|U-u_{h}\right\|_{L^{2}(\Omega)}\leq Ch^{k+1/2}.$

This completes the proof. ∎

##  5\. Numerical Experiments

In this section, our aim is to validate the proposed LDG scheme (2.4) for the
BO equation (1.1). We have achieved the stability and error analysis results
for the semi-discrete scheme of the form

$\frac{d}{dt}u_{h}=L_{h}u_{h}.$

Our focus lies on verifying the performance of the scheme using a low storage
explicit Runge-Kutta (LSERK) of fourth order [6, 18] time discretization
method of the form

$\displaystyle r^{(0)}=u_{h}^{n},$ $\displaystyle\text{for }j=1:5$
$\displaystyle\qquad\begin{cases}k^{j}=a_{j}k^{j-1}+\Delta tL_{h}r^{(j-1)},\\\
r^{(j)}=r^{(j-1)}+b_{j}k^{j},\end{cases}$ $\displaystyle u_{h}^{n+1}=r^{(5)},$

where the weighted coefficients $a_{j}$, $b_{j}$ and $c_{j}$ of the LSERK
method given in [18]. The above iteration in comparison with the classical
fourth order method (3.14) is considerably more efficient and accurate than
(3.14) since it has the disadvantage that it requires four extra storage
arrays. We compute the Hilbert transform of the approximation $q_{h}$ of
$Q=U_{x}$ in the scheme (2.4) by using the fast Fourier transform. Hence we
need to use the periodic solution of the Benjamin-Ono equation. We are
providing the steps of algorithm to compute the approximate solution using LDG
scheme (2.4), for more details one may refer to [18]. Order of convergence
[12] for errors is defined for each intermediate step between element numbers
$N_{1}$ and $N_{2}$ as

$\text{order}=\frac{\ln(E(N_{1}))-\ln(E(N_{2}))}{\ln(N_{2})-\ln(N_{1})},$

where error $E$ can be seen as a function of number of elements $N$.

Algorithm 1 Numerical Procedure for Benjamin-Ono Equation

1: Set up problem parameters and initial conditions

2: Compute time step size $dt$ based on CFL condition

3: Determine total number of time steps $Nsteps$ based on Final time and $dt$

4: for $tstep=1$ to $Nsteps$ do

5: for $INTRK=1$ to $5$ do

6: Compute local time $timelocal$ for current stage

7: Evaluate right-hand side $L_{h}u_{h}$ of the differential equation using
Algorithm 2

8: Update solution $u$ using Runge-Kutta stages

9: end for

10: Increment time by $dt$

11: end for

12: return Final solution $u$ and final time $time$

Algorithm 2 Algorithm for Computing RHS of Benjamin-Ono Equation

0: Solution $u$, left boundary value $xL$, right boundary value $xR$, current
time $timelocal$

0: Right-hand side of the equation $rhsu$

1: Compute coefficients and variables based on initial conditions and boundary
values

2: Define field differences at interfaces, including boundary conditions

3: Compute flux for $u$

4: Compute local variable $q$ and its differences using $u$

5: Compute local variable $p$, using Hilbert transform of $q$

6: Compute flux for $p$ ($fluxp$)

7: Compute squared differences of $u$, $du^{2}$ at interfaces, and impose
boundary conditions

8: Compute flux for the convective term ($fluxf$)

9: Compute derivative of the flux ($dfdr$)

10: Compute right-hand side of the semi-discrete PDE ($rhsu$)

11: return $rhsu$

The BO equation (1.1) possesses an infinite number of conserved quantities
[16]. Hereby we consider the first two specific quantities known as $mass$ and
$momentum$ and with normalization, these quantities can be expressed as
follows:

$\displaystyle\mathbf{C}_{1}^{h}$
$\displaystyle:=\frac{\int_{\Omega}u_{h}\,dx}{\int_{\Omega}U_{0}\,dx},\qquad\mathbf{C}_{2}^{h}:=\frac{\|u_{h}\|_{L^{2}(\Omega)}}{\left\|U_{0}\right\|_{L^{2}(\Omega)}}.$

Our aim is to preserve these quantities in the discrete set up.

To verify the proposed LDG scheme (2.4)-(2.6), we compare it with the exact
periodic solution of the Benjamin-Ono equation (1.1) with
$f(U)=\frac{1}{2}U^{2}$. For instance, we consider the solution of (1.1) from
Thomée et. al. [35] (also see [13, 14])

$U(x,t)=\frac{2c\delta}{1-\sqrt{1-\delta^{2}}\cos(c\delta(x-ct))},\qquad\delta=\frac{\pi}{cL},$
(5.1)

where we choose $L=15$, $c=0.25$. We compare the above exact solution at the
final time $T=10$ and $T=20$. This solution characterizes periodic solitary
waves, exhibiting periodicity and amplitude determined by the parameters $L$
and $c$ respectively.

In our numerical implementation, we initially compare the approximate solution
obtained from the LDG scheme (2.4)-(2.6) with the exact solution at time
$t=10$, and further at the final time $t=20$. Figure 5.1 depicts the
comparison of the solution at various times. This validation confirms that the
numerical scheme converges to the exact solution. The Table 5.1 presents the
$L^{2}$-errors of the proposed scheme under polynomial degree up to three. It
is observed that the errors are converging to zero with relatively coarser
grids and optimal orders of convergence are obtained at time $T=10$ for each
polynomial degree up to three and the discrete conserved quantities
$\mathbf{C}_{1}^{h}$ and $\mathbf{C}_{2}^{h}$ are also preserved. Similar
analysis has been carried out at time $T=20$ and details are presented in
Table 5.2 and optimal rates are obtained there as well.

We have also approached to compute the Hilbert transform of the approximation
by using Gaussian quadrature. However, we do not achieve the optimal rates or
the theoretical rates obtained by convergence analysis. Similar observations
are also described in [16, Section 5]. One possible interpretation is that the
principal value integral is not being computed accurately near the
singularity, and it is most likely caused by our approximation of the full
line problem by restricting to a bounded interval while still employing the
full line Hilbert transform.

Figure 5.1. Approximate solution computed by the LDG scheme at $T=10$ and $T=20$ with $N=160$, $k=3$ and initial condition $U(x,0)$ of (1.1). k | k=1 | k=2 | k=3
---|---|---|---
N | error | order | error | order | error | order | $\mathbf{C_{1}^{h}}$ | $\mathbf{C_{2}^{h}}$
40 | 2.91e-01 |  | 3.65e-02 |  | 8.69e-03 |  | 1.01 | 1.01
|  | 2.01 |  | 2.76 |  | 4.37 |  |
80 | 7.15e-02 |  | 5.37e-03 |  | 4.20e-04 |  | 1.00 | 1.00
|  | 2.02 |  | 2.93 |  | 4.07 |  |
160 | 1.77e-02 |  | 7.029e-04 |  | 2.49e-05 |  | 1.00 | 1.00
|  | 2.00 |  | 2.98 |  | 4.00 |  |
320 | 4.40e-03 |  | 8.88e-05 |  | 1.56e-06 |  | 1.00 | 1.00
Table 5.1. $L^{2}$-error and order of convergence for the LDG scheme at time $T=10$ taking $N$ elements and polynomial degree $k$. k | k=1 | k=2 | k=3
---|---|---|---
N | error | order | error | order | error | order | $\mathbf{C_{1}^{h}}$ | $\mathbf{C_{2}^{h}}$
40 | 6.02e-01 |  | 6.015e-01 |  | 5.98e-02 |  | 0.99 | 0.98
|  | 2.05 |  | 2.65 |  | 4.46 |  |
80 | 1.45e-01 |  | 9.54e-02 |  | 2.718e-03 |  | 1.00 | 1.00
|  | 2.02 |  | 2.90 |  | 4.09 |  |
160 | 3.56e-02 |  | 1.27e-02 |  | 1.59e-04 |  | 1.00 | 1.00
|  | 2.01 |  | 2.97 |  | 4.02 |  |
320 | 8.83e-03 |  | 1.62e-03 |  | 9.80e-06 |  | 1.00 | 1.00
Table 5.2. $L^{2}$-error and order of convergence for the LDG scheme at time
$T=20$ taking $N$ elements and polynomial degree $k$.

##  6\. Concluding remarks

We have designed a local discontinuous Galerkin method for the Benjamin-Ono
equation with general nonlinear flux. We have shown that the semi-discrete
scheme is stable and convergent. The stability analysis is also carried out
for four-stage fourth order time marching RK method considering linear flux
function. The stability analysis for nonlinear flux and convergence analysis
for fully-discrete scheme will be addressed in the future work. In the
numerical experiments, it is observed that optimal rates are obtained for
various degrees of polynomials which demonstrates the efficiency and accuracy
of the proposed scheme.

## References

  * [1] M. J. Ablowitz and A. S. Fokas.  The inverse scattering transform for the Benjamin-Ono equation-a pivot to multidimensional problems.  Studies in Applied Mathematics, 68 (1983), no. 1, 1-10.
  * [2] T. Aboelenen.  A high-order nodal discontinuous Galerkin method for nonlinear fractional Schrödinger type equations.  Communications in Nonlinear Science and Numerical Simulation, 54 (2018), 428-452.
  * [3] J. Ai, Y. Xu, C.-W. Shu and Q. Zhang.  $L^{2}$ error estimate to smooth solutions of high order Runge-Kutta discontinuous Galerkin method for scalar nonlinear conservation laws with and without sonic points.  SIAM Journal on Numerical Analysis, 60 (2022), no. 4, 1741-1773.
  * [4] F. Bassi and S. Rebay.  A high-order accurate discontinuous finite element method for the numerical solution of the compressible Navier–Stokes equations.  Journal of Computational Physics, 131 (1997), no. 2, 267-279.
  * [5] N. Burq and F. Planchon.  On well-posedness for the Benjamin–Ono equation.  Mathematische Annalen, 340 (2008), no. 3, 497-542.
  * [6] H. M. Carpenter and A. C. Kennedy.  Fourth-order 2N-storage Runge-Kutta schemes.  NASA TM 109112, NASA Langley Research Center, 1994.
  * [7] K. M. Case.  Benjamin-Ono related equations and their solutions.  Proceedings of the National Academy of Sciences, 76 (1979), no. 1, 1-3.
  * [8] P. G. Ciarlet.  The finite element method for elliptic problems.  Society for Industrial and Applied Mathematics, 2002.
  * [9] B. Cockburn, G. E. Karniadakis and C.-W. Shu.  Discontinuous Galerkin methods: theory, computation and applications.  Springer Science & Business Media, 2012.
  * [10] B. Cockburn and C.-W. Shu.  The local discontinuous Galerkin method for time-dependent convection-diffusion systems.  SIAM Journal on Numerical Analysis, 35 (1998), no. 6, 2440-2463.
  * [11] R. Dutta, H. Holden, U. Koley and N. H. Risebro.  Convergence of finite difference schemes for the Benjamin-Ono equation.  Numerische Mathematik, 134 (2016), no. 2, 249-274.
  * [12] M. Dwivedi and T. Sarkar.  Convergence of a conservative Crank-Nicolson finite difference scheme for the KdV equation with smooth and non-smooth initial data.  arXiv preprint arXiv:2312.14454, (2023).
  * [13] M. Dwivedi and T. Sarkar.  Stability and convergence analysis of a Crank-Nicolson Galerkin scheme for the fractional Korteweg-de Vries equation.  SMAI Journal of Computational Mathematics, to appear.
  * [14] M. Dwivedi and T. Sarkar.  Fully discrete finite difference schemes for the fractional Korteweg-de Vries equation.  arXiv preprint arXiv:2403.08275, (2024).
  * [15] M. Dwivedi and T. Sarkar.  Local discontinuous Galerkin method for fractional Korteweg-de Vries equation.  arXiv preprint arXiv:2404.18069, (2024).
  * [16] S. T. Galtung.  Convergent Crank–Nicolson Galerkin scheme for the Benjamin-Ono equation.  Discrete and Continuous Dynamical Systems, 38 (2018), no. 3, 1243-1268.
  * [17] S. T. Galtung.  Convergence rates of a fully discrete Galerkin scheme for the Benjamin-Ono equation.  XVI International Conference on Hyperbolic Problems: Theory, Numerics, Applications, Springer, (2016), 589-601.
  * [18] J. S. Hesthaven and T. Warburton.  Nodal discontinuous Galerkin methods: algorithms, analysis, and applications.  Springer Science and Business Media, 2007.
  * [19] J. Hunter, Z. Sun and Y. Xing.  Stability and time-Step constraints of implicit-explicit Runge-Kutta methods for the linearized Korteweg-de Vries equation.  Communications on Applied Mathematics and Computation, (2023), 1-30.
  * [20] A. Ionescu and C. E. Kenig.  Global well-posedness of the Benjamin-Ono equation in low-regularity spaces.  Journal of the American Mathematical Society, 20 (2007), no. 3, 753-798.
  * [21] R. J. Iório.  On the Cauchy problem for the Benjamin-Ono equation.  Communications in Partial Differential Equations, 11 (1986), no. 10, 1031-1081.
  * [22] Y. Ishimori.  Solitons in a one-dimensional Lennard-Jones lattice.  Progress of Theoretical Physics, 68 (1982), no. 2, 402-410.
  * [23] C. E. Kenig and K. D. Koenig.  On the local well-posedness of the Benjamin-Ono and modified Benjamin-Ono equations.  Mathematical Research Letters, 10 (2003), no. 6, 879-895.
  * [24] F. W. King.  Hilbert transforms: Volume 2.  Cambridge University Press, 2009.
  * [25] H. Koch and N. Tzvetkov.  On the local well-posedness of the Benjamin-Ono equation in $H^{s}(\mathbb{R})$.  International Mathematics Research Notices, (2003), no. 26, 1449-1464.
  * [26] D. Levy, C.-W. Shu and J. Yan.  Local discontinuous Galerkin methods for nonlinear dispersive equations.  Journal of Computational Physics, 196 (2004), no. 2, 751-772.
  * [27] Y. Matsuno and D. J. Kaup.  Initial value problem of the linearized Benjamin-Ono equation and its applications.  Journal of Mathematical Physics, 38 (1997), no. 10, 5198-5224.
  * [28] L. Molinet.  Global well-posedness in $L^{2}$ for the periodic Benjamin-Ono equation.  American Journal of Mathematics, 130 (2008), no. 3, 635-683.
  * [29] G. Ponce.  On the global well-posedness of the Benjamin-Ono equation.  Differential and Integral Equations, 4 (1991), 527-542.
  * [30] W. H. Reed and T. R. Hill.  Triangular mesh methods for the neutron transport equation.  Technical Report LA-UR-73-479, Los Alamos Scientific Lab, 1973.
  * [31] Z. Sun and C.-W. Shu.  Error analysis of Runge-Kutta discontinuous Galerkin methods for linear time-dependent partial differential equations.  arXiv preprint arXiv:2001.00971, (2020).
  * [32] Z. Sun and C.-W. Shu.  Strong stability of explicit Runge-Kutta time discretizations.  SIAM Journal on Numerical Analysis, 57 (2019), no. 3, 1158-1182.
  * [33] Z. Sun and C.-W. Shu.  Stability of the fourth order Runge-Kutta method for time-dependent partial differential equations.  Annals of Mathematical Sciences and Applications, 2 (2017), no. 2, 255-284.
  * [34] T. Tao.  Global well-posedness of the Benjamin-Ono equation in $H^{1}(\mathbb{R})$.  Journal of Hyperbolic Differential Equations, 1 (2004), no. 1, 27-49.
  * [35] V. Thomée and A. S. V. Murthy.  A numerical method for the Benjamin–Ono equation.  BIT Numerical Mathematics, 38 (1998), 597-611.
  * [36] Q. Xu and J. S. Hesthaven.  Discontinuous Galerkin method for fractional convection-diffusion equations.  SIAM Journal on Numerical Analysis, 52 (2014), no. 1, 405-423.
  * [37] Y. Xu, Q. Zhang, C.-W. Shu and H. wang.  The $L^{2}$-norm stability analysis of Runge-Kutta discontinuous Galerkin methods for linear hyperbolic equations.  SIAM Journal on Numerical Analysis, 57 (2019), no. 4, 1574-1601.
  * [38] Y. Xu and C.-W. Shu.  Local discontinuous Galerkin methods for nonlinear Schrödinger equations.  Journal of Computational Physics, 205 (2005), no. 1, 72-97.
  * [39] Y. Xu and C.-W. Shu.  Error estimates of the semi-discrete local discontinuous Galerkin method for nonlinear convection-diffusion and KdV equations.  Computer Methods in Applied Mechanics and Engineering, 196 (2007), no. 37-40, 3805–3822.
  * [40] Y. Xu and C.-W. Shu.  Local discontinuous Galerkin methods for high-order time-dependent partial differential equations.  Communications in Computational Physics, 7 (2010), no. 1, 1-46.
  * [41] J. Yan and C.-W. Shu.  A local discontinuous Galerkin method for KdV type equations.  SIAM Journal on Numerical Analysis, 40 (2002), no. 2, 769-791.
  * [42] Q. Zhang and C.-W. Shu.  Error estimates to smooth solutions of Runge-Kutta discontinuous Galerkin methods for scalar conservation laws.  SIAM Journal on Numerical Analysis, 42 (2004), no. 2, 641-666.


# Meta Label Propagation for Low-resource Semi-supervised Learning on Graphs

Anonymous Author(s)
Affiliation

###### Abstract

Inspired by the extensive success of deep learning, graph neural networks
(GNNs) have been proposed to learn expressive node representations and
demonstrated promising performance on different graph learning tasks. However,
existing endeavors predominately focus on the conventional semi-supervised
setting where relatively abundant gold-labeled nodes are provided. Such an
assumption is often impractical due to the fact that data labeling is
unbearably laborious and requires intensive domain-knowledge, especially when
considering the heterogeneity of graph-structured data. While under the low-
resource semi-supervised setting, most of the existing GNNs built upon the
message-passing scheme have two major limitations, which can largely undermine
the model performance in practice: (i) the notorious oversmoothing issue when
increasing the message-passing range; and (ii) the over-reliance on the
homophily assumption about the input graph. In this paper, we propose a new
architecture that is able to learn simple yet effective GNN models with
limited labeled data. Specifically, we first propagate the labels of training
data along the graph to generate abundant pseudo labels, then progressively
correct the mislabeled nodes during propagation and learn the prediction model
within a unified meta-learning algorithm. Extensive experiments demonstrate
that our approach offers easy and substantial performance gains compared to
existing techniques on different benchmark datasets.

##  1 Introduction

##  2 Related Work

##  3 Our Approach: Meta Label Propagation

##  4 Experiments

##  5 Conclusion

## Checklist

##  Appendix A Appendix