11institutetext: Interdisciplinary Center for Security, Reliability and Trust
University of Luxembourg, Luxembourg
11email<EMAIL_ADDRESS>

# Probabilistic Inferences in Bayesian Networks

Jianguo Ding

###### Abstract

Bayesian network is a complete model for the variables and their
relationships, it can be used to answer probabilistic queries about them. A
Bayesian network can thus be considered a mechanism for automatically applying
Bayes’ theorem to complex problems. In the application of Bayesian networks,
most of the work is related to probabilistic inferences. Any variable updating
in any node of Bayesian networks might result in the evidence propagation
across the Bayesian networks. This paper sums up various inference techniques
in Bayesian networks and provide guidance for the algorithm calculation in
probabilistic inference in Bayesian networks.

##  1 Introduction

Because a Bayesian network is a complete model for the variables and their
relationships, it can be used to answer probabilistic queries about them. For
example, the network can be used to find out updated knowledge of the state of
a subset of variables when other variables (the evidence variables) are
observed. This process of computing the posterior distribution of variables
given evidence is called probabilistic inference. A Bayesian network can thus
be considered a mechanism for automatically applying Bayes’ theorem to complex
problems.

In the application of Bayesian networks, most of the work is related to
probabilistic inferences. Any variable updating in any node of Bayesian
networks might result in the evidence propagation across the Bayesian
networks. How to examine and execute various inferences is the important task
in the application of Bayesian networks.

This chapter will sum up various inference techniques in Bayesian networks and
provide guidance for the algorithm calculation in probabilistic inference in
Bayesian networks. Information systems are of discrete event characteristics,
this chapter mainly concerns the inferences in discrete events of Bayesian
networks.

##  2 The Semantics of Bayesian Networks

The key feature of Bayesian networks is the fact that they provide a method
for decomposing a probability distribution into a set of local distributions.
The independence semantics associated with the network topology specifies how
to combine these local distributions to obtain the complete joint probability
distribution over all the random variables represented by the nodes in the
network. This has three important consequences.

Firstly, naively specifying a joint probability distribution with a table
requires a number of values exponential in the number of variables. For
systems in which interactions among the random variables are sparse, Bayesian
networks drastically reduce the number of required values.

Secondly, efficient inference algorithms are formed in that work by
transmitting information between the local distributions rather than working
with the full joint distribution.

Thirdly, the separation of the qualitative representation of the influences
between variables from the numeric quantification of the strength of the
influences has a significant advantage for knowledge engineering. When
building a Bayesian network model, one can focus first on specifying the
qualitative structure of the domain and then on quantifying the influences.
When the model is built, one is guaranteed to have a complete specification of
the joint probability distribution.

The most common computation performed on Bayesian networks is the
determination of the posterior probability of some random variables, given the
values of other variables in the network. Because of the symmetric nature of
conditional probability, this computation can be used to perform both
diagnosis and prediction. Other common computations are: the computation of
the probability of the conjunction of a set of random variables, the
computation of the most likely combination of values of the random variables
in the network and the computation of the piece of evidence that has or will
have the most influence on a given hypothesis.

A detailed discussion of inference techniques in Bayesian networks can be
found in the book by Pearl [Pearl, 2000].

  * •

Probabilistic semantics. Any complete probabilistic model of a domain must,
either explicitly or implicitly, represent the joint distribution which the
probability of every possible event as defined by the values of all the
variables. There are exponentially many such events, yet Bayesian networks
achieve compactness by factoring the joint distribution into local,
conditional distributions for each variable given its parents. If $x_{i}$
denotes some value of the variable $X_{i}$ and $\pi(x_{i})$ denotes some set
of values for $X_{i}$’s parents $\pi(x_{i})$, then $P(x_{i}|\pi(x_{i}))$
denotes this conditional distribution. For example, $P(x_{4}|x_{2},x_{3})$ is
the probability of wetness given the values of sprinkler and rain. Here
$P(x_{4}|x_{2},x_{3})$ is the brief of $P(x_{4}|\\{x_{2},x_{3}\\})$. The set
parentheses are omitted for the sake of readability. We use the same
expression in this thesis. The global semantics of Bayesian networks specifies
that the full joint distribution is given by the product

$P(x_{1},\ldots,x_{n})=\prod_{i}P(x_{i}|\pi(x_{i}))$ (1)

Equation 1 is also called the chain rule for Bayesian networks.

Figure 1: Causal Influences in A Bayesian Network.

In the example Bayesian network in Figure 1, we have

$P(x_{1},x_{2},x_{3},x_{4},x_{5})=P(x_{1})P(x_{2}|x_{1})P(x_{3}|x_{1})P(x_{4}|x_{2},x_{3})P(x_{5}|x_{4})$
(2)

Provided the number of parents of each node is bounded, it is easy to see that
the number of parameters required grows only linearly with the size of the
network, whereas the joint distribution itself grows exponentially. Further
savings can be achieved using compact parametric representations, such as
noisy-OR models, decision tress, or neural networks, for the conditional
distributions [Pearl, 2000].

There are also entirely equivalent local semantics, which assert that each
variable is independent of its non-descendants in the network given its
parents. For example, the parents of $X_{4}$ in Figure 1 are $X_{2}$ and
$X_{3}$ and they render $X_{4}$ independent of the remaining non-descendant,
$X_{1}$. That is,

$P(x_{4}|x_{1},x_{2},x_{3})=P(x_{4}|x_{2},x_{3})$ (3)

The collection of independence assertions formed in this way suffices to
derive the global assertion in Equation 2, and vice versa. The local semantics
are most useful in constructing Bayesian networks, because selecting as
parents the direct causes of a given variable automatically satisfies the
local conditional independence conditions. The global semantics lead directly
to a variety of algorithms for reasoning.

  * •

Evidential reasoning. From the product specification in Equation 2, one can
express the probability of any desired proposition in terms of the conditional
probabilities specified in the network. For example, the probability that the
sprinkler was on, given that the pavement is slippery, is

$\displaystyle P(X_{3}=on|X_{5}=true)$ (4)
$\displaystyle=\frac{P(X_{3}=on,X_{5}=true)}{P(X_{5}=true)}$
$\displaystyle=\frac{\sum_{x_{1},x_{2},x_{4}}P(x_{1},x_{2},X_{3}=on,x_{4},X_{5}=true)}{\sum_{x_{1},x_{2},x_{3},x_{4}}P(x_{1},x_{2},x_{3},x_{4},X_{5}=true)}$
$\displaystyle=\frac{\sum_{x_{1},x_{2},x_{4}}P(x_{1})P(x_{2}|x_{1})P(X_{3}=on|x_{1})P(x_{4}|x_{2},X_{3}=on)P(X_{5}=true|x_{4})}{\sum_{x_{1},x_{2},x_{3},x_{4}}P(x_{1})P(x_{2}|x_{1})P(x_{3}|x_{1})P(x_{4}|x_{2},x_{3})P(X_{5}=true|x_{4})}$

These expressions can often be simplified in the ways that reflect the
structure of the network itself.

It is easy to show that reasoning in Bayesian networks subsumes the
satisfiability problem in propositional logic and hence reasoning is NP-hard
[Cooper, 1990]. Monte Carlo simulation methods can be used for approximate
inference [Pearl, 1987], given that estimates are gradually improved as the
sampling proceeds. (Unlike join-tree methods, these methods use local message
propagation on the original network structure.) Alternatively, variational
methods [Jordan et al., 1998] provide bounds on the true probability.

  * •

Functional Bayesian networks. The networks discussed so far are capable of
supporting reasoning about evidence and about actions. Additional refinement
is necessary in order to process counterfactual information. For example, the
probability that ”the pavement would not have been slippery had the sprinkler
been OFF, given that the sprinkler is in fact ON and that the pavement is in
fact slippery” cannot be computed from the information provided in Figure 1
and Equation 2. Such counterfactual probabilities require a specification in
the form of functional networks, where each conditional probability
$P(x_{i}|\pi(i))$ is replaced by a functional relationship
$x_{i}=f_{i}(\pi(i),\epsilon_{i})$, where $\epsilon_{i}$ is a stochastic
(unobserved) error term. When the functions $f_{i}$ and the distributions of
$\epsilon_{i}$ are known, all counterfactual statements can be assigned unique
probabilities, using evidence propagation in a structure called a ”twin
network”. When only partial knowledge about the functional form of $f_{i}$ is
available, bounds can be computed on the probabilities of counterfactual
sentences [Balke & Pearl, 1995] [Pearl, 2000].

  * •

Causal discovery. One of the most exciting prospects in recent years has been
the possibility of using Bayesian networks to discover causal structures in
raw statistical data [Pearl & Verma, 1991] [Spirtes et al., 1993] [Pearl,
2000], which is a task previously considered impossible without controlled
experiments. Consider, for example, the following pattern of dependencies
among three events: $A$ and $B$ are dependent, $B$ and $C$ are dependent, yet
$A$ and $C$ are independent. If you ask a person to supply an example of three
such events, the example would invariably portray $A$ and $C$ as two
independent causes and $B$ as their common effect, namely, $A\rightarrow
B\leftarrow C$. Fitting this dependence pattern with a scenario in which $B$
is the cause and $A$ and $C$ are the effects is mathematically feasible but
very unnatural, because it must entail fine tuning of the probabilities
involved; the desired dependence pattern will be destroyed as soon as the
probabilities undergo a slight change.

Such thought experiments tell us that certain patterns of dependency, which
are totally void of temporal information, are conceptually characteristic of
certain causal directionalities and not others. When put together
systematically, such patterns can be used to infer causal structures from raw
data and to guarantee that any alternative structure compatible with the data
must be less stable than the one(s) inferred; namely, slight fluctuations in
parameters will render that structure incompatible with the data.

  * •

Plain beliefs. In mundane decision making, beliefs are revised not by
adjusting numerical probabilities but by tentatively accepting some sentences
as ”true for all practical purposes”. Such sentences, called plain beliefs,
exhibit both logical and probabilistic characters. As in classical logic, they
are propositional and deductively closed; as in probability, they are subject
to retraction and to varying degrees of entrenchment. Bayesian networks can be
adopted to model the dynamics of plain beliefs by replacing ordinary
probabilities with non-standard probabilities, that is, probabilities that are
infinitesimally close to either zero or one [Goldszmidt & Pearl, 1996].

  * •

Models of cognition. Bayesian networks may be viewed as normative cognitive
models of propositional reasoning under uncertainty [Pearl, 2000]. They handle
noise and partial information by using local, distributed algorithm for
inference and learning. Unlike feed forward neural networks, they facilitate
local representations in which nodes correspond to propositions of interest.
Recent experiments [Tenenbaum & Griffiths, 2001] suggest that they capture
accurately the causal inferences made by both children and adults. Moreover,
they capture patterns of reasoning that are not easily handled by any
competing computational model. They appear to have many of the advantages of
both the “symbolic” and the “subsymbolic” approaches to cognitive modelling.

Two major questions arise when we postulate Bayesian networks as potential
models of actual human cognition.

Firstly, does an architecture resembling that of Bayesian networks exist
anywhere in the human brain? No specific work had been done to design neural
plausible models that implement the required functionality, although no
obvious obstacles exist.

Secondly, how could Bayesian networks, which are purely propositional in their
expressive power, handle the kinds of reasoning about individuals, relations,
properties, and universals that pervades human thought? One plausible answer
is that Bayesian networks containing propositions relevant to the current
context are constantly being assembled as needed to form a more permanent
store of knowledge. For example, the network in Figure 1 may be assembled to
help explain why this particular pavement is slippery right now, and to decide
whether this can be prevented. The background store of knowledge includes
general models of pavements, sprinklers, slipping, rain, and so on; these must
be accessed and supplied with instance data to construct the specific Bayesian
network structure. The store of background knowledge must utilize some
representation that combines the expressive power of first-order logical
languages (such as semantic networks) with the ability to handle uncertain
information.

##  3 Reasoning Structures in Bayesian Networks

###  3.1 Basic reasoning structures

####  3.1.1 d-Separation in Bayesian Networks

d-Separation is one important property of Bayesian networks for inference.
Before we define d-separation, we first look at the way that evidence is
transmitted in Bayesian Networks. There are two types of evidence:

  * •

Hard Evidence (instantiation) for a node $A$ is evidence that the state of $A$
is definitely a particular value.

  * •

Soft Evidence for a node $A$ is any evidence that enables us to update the
prior probability values for the states of $A$.

d-Separation (Definition):

Two distinct variables $X$ and $Z$ in a causal network are d-separated if, for
all paths between $X$ and $Z$, there is an intermediate variable $V$ (distinct
from $X$ and $Z$) such that either

  * •

the connection is serial or diverging and $V$ is instantiated or

  * •

the connection is converging, and neither $V$ nor any of $V$’s descendants
have received evidence.

If $X$ and $Z$ are not d-separated, we call them d-connected.

####  3.1.2 Basic structures of Bayesian Networks

Based on the definition of d-seperation, three basic structures in Bayesian
networks are as follows:

  1. 1.

Serial connections

Consider the situation in Figure 2. $X$ has an influence on $Y$, which in turn
has an influence on $Z$. Obviously, evidence on $Z$ will influence the
certainty of $Y$, which then influences the certainty of $Z$. Similarly,
evidence on $Z$ will influence the certainty on $X$ through $Y$. On the other
hand, if the state of $Y$ is known, then the channel is blocked, and $X$ and
$Z$ become independent. We say that $X$ and $Z$ are d-separated given $Y$, and
when the state of a variable is known, we say that it is instantiated (hard
evidence).

We conclude that evidence may be transmitted through a serial connection
unless the state of the variable in the connection is known.

Figure 2: Serial Connection. When $Y$ is Instantiated, it blocks the
communication between $X$ and $Z$.

  2. 2.

Diverging connections

The situation in Figure 3 is called a diverging connection. Influence can pass
between all the children of $X$ unless the state of $X$ is known. We say that
$Y_{1},Y_{2},\ldots,Y_{n}$ are d-separated given $X$.

Evidence may be transmitted through a diverging connection unless it is
instantiated.

Figure 3: Diverging Connection. If $X$ is instantiated, it blocks the
communication between its children.

  3. 3.

Converging connections

Figure 4: Converging Connection. If $Y$ changes certainty, it opens for the
communication between its parents.

A description of the situation in Figure 4 requires a little more care. If
nothing is known about $Y$ except what may be inferred from knowledge of its
parents $X_{1},\ldots,X_{n}$, then the parents are independent: evidence on
one of the possible causes of an event does not tell us anything about other
possible causes. However, if anything is known about the consequences, then
information on one possible cause may tell us something about the other
causes.

This is the explaining away effect illustrated in Figure 1. $X_{4}$ (pavement
is wet) has occurred, and $X_{3}$ (the sprinkler is on) as well as $X_{2}$
(it’s raining) may cause $X_{4}$. If we then get the information that $X_{2}$
has occurred, the certainty of $X_{3}$ will decrease. Likewise, if we get the
information that $X_{2}$ has not occurred, then the certainty of $X_{3}$ will
increase.

The three preceding cases cover all ways in which evidence may be transmitted
through a variable.

##  4 Classification of Inferences in Bayesian Networks

In Bayesian networks, 4 popular inferences are identified as:

  1. 1.

Forward Inference

Forward inferences is also called predictive inference (from causes to
effects). The inference reasons from new information about causes to new
beliefs about effects, following the directions of the network arcs. For
example, in Figure 2, $X\rightarrow Y\rightarrow Z$ is a forward inference.

  2. 2.

Backward Inference

Backward inferences is also called diagnostic inference (from effects to
causes). The inference reasons from symptoms to cause, Note that this
reasoning occurs in the opposite direction to the network arcs. In Figure 2 ,
$Z\rightarrow Y$ is a backward inference. In Figure 3 , $Y_{i}\rightarrow
X(i\in[1,n])$ is a backward inference.

  3. 3.

Intercausal Inference

Intercausal inferences is also called explaining away (between parallel
variables). The inference reasons about the mutual causes (effects) of a
common effect (cause). For example, in Figure 4, if the $Y$ is instantiated,
$X_{i}$ and $X_{j}(i,j\in[1,n])$ are dependent. The reasoning
$X_{i}\leftrightarrow X_{j}(i,j\in[1,n])$ is an intercausal inference. In
Figure 3, if $X$ is not instantiated, $Y_{i}$ and $Y_{j}(i,j\in[1,n])$ are
dependent. The reasoning $Y_{i}\leftrightarrow Y_{j}(i,j\in[1,n])$ is an
intercausal inference.

  4. 4.

Mixed inference

Mixed inferences is also called combined inference. In complex Bayesian
networks, the reasoning does not fit neatly into one of the types described
above. Some inferences are a combination of several types of reasoning.

###  4.1 Inference in Bayesian Networks

####  4.1.1 inference in basic models

  * •

in Serial Connections

    * –

the forward inference executes with the evidence forward propagation. For
example, in Figure 5, consider the inference $X\rightarrow Y\rightarrow Z$.
111Note: In this chapter, $P(X^{+})$ is the abbreviation of $P(X=true)$,
$P(X^{-})$ is the abbreviation of $P(|X=false)$. For simple expression, we use
$P(Y|X)$ to denote $P(Y=true|X=true)$ by default. But in express $P(Y^{+}|X)$,
$X$ denotes both situations $X=true$ and $X=false$.

Figure 5: Inference in Serial Connection

If Y is instantiated, X and Z are independent, then we have following example:

$P(Z|XY)=P(Z|Y)$;

$P(Z^{+}|Y^{+})=0.95$;

$P(Z^{-}|Y^{+})=0.05$;

$P(Z^{+}|Y^{-})=0.01$;

$P(Z^{-}|Y^{-})=0.99$;

if Y is not instantiated, X and Z are dependent, then

$P(Z^{+}|X^{+}Y)=P(Z^{+}|Y^{+})P(Y^{+}|X^{+})+P(Z^{+}|Y^{-})P(Y^{-}|X^{+})$

$=0.95*0.85+0.01*0.15=0.8075+0.0015=0.809$;

$P(Z^{-}|X^{-}Y)=P(Z^{-}|Y^{+})P(Y^{+}|X^{-})+P(Z^{-}|Y^{-})P(Y^{-}|X^{-})$

$=0.05*0.03+0.99*0.97=0.0015+0.9603=0.9618$.

    * –

the backward inference executes the evidence backward propagation. For
example, in Figure 5, consider the inference $Z\rightarrow Y\rightarrow X$.

      1. 1.

If $Y$ is instantiated ($P(Y^{+})=1$ or $P(Y^{-})=1)$, $X$ and $Z$ are
independent, then

$\displaystyle P(X|YZ)=P(X|Y)=\frac{P(X)P(Y|X)}{P(Y)}$ (5)

$P(X^{+}|Y^{+}Z)=P(X^{+}|Y^{+})=\frac{P(X^{+})P(Y^{+}|X^{+})}{P(Y^{+})}=\frac{09*0.85}{1}=0.765$;

$P(X^{+}|Y^{-}Z)=P(X^{+}|Y^{-})=\frac{P(X^{+})P(Y^{-}|X^{+})}{P(Y^{-})}=\frac{09*0.15}{1}=0.135$.

      2. 2.

If $Y$ is not instantiated, $X$ and $Z$ are dependent (See the dashed lines in
Figure 5). Suppose $P(Z^{+})=1$ then

$P(X^{+}|YZ^{+})=\frac{P(X^{+}YZ^{+})}{P(YZ^{+})}=\frac{P(X^{+}YZ^{+})}{\sum_{X}P(XYZ^{+})}$;

$P(X^{+}YZ^{+})=P(X^{+}Y^{+}Z^{+})+P(X^{+}Y^{-}Z^{+})=0.9*0.85*0.95+0.9*0.15*0.05=0.72675+0.00675=0.7335$;

$\sum_{X}P(XYZ^{+})=P(X^{+}Y^{+}Z^{+})+P(X^{+}Y^{-}Z^{+})+P(X^{-}Y^{+}Z^{+})+P(X^{-}Y^{-}Z^{+})\\\
=0.9*0.85*0.95+0.9*0.15*0.99+0.1*0.03*0.95+0.1*0.97*0.01\\\
=0.72675+0.13365+0.00285+0.00097=0.86422$;

$P(X^{+}|YZ^{+})=\frac{P(X^{+}YZ^{+})}{\sum_{X}P(XYZ^{+})}=\frac{0.7335}{0.86422}=0.8487.$

In serial connections, there is no intercausal inference.

  * •

in Diverging Connections

    * –

the forward inference executes with the evidence forward propagation. For
example, in Figure 6, consider the inference $Y\rightarrow X$ and
$Y\rightarrow Z$, the goals are easy to obtain by nature.

Figure 6: Inference in Diverging Connection

    * –

the backward inference executes with the evidence backward propagation, see
the dashed line in Figure 6, consider the inference $(XZ)\rightarrow Y$, $X$
and $Z$ are instantiated by assumption, suppose $P(X^{+}=1)$, $P(Z^{+}=1)$.
Then,

$\displaystyle
P(Y^{+}|X^{+}Z^{+})=\frac{P(Y^{+}X^{+}Z^{+})}{P(X^{+}Z^{+})}=\frac{P(Y^{+})P(X^{+}|Y^{+})P(Z^{+}|Y^{+})}{P(X^{+}Z^{+})}$
$\displaystyle=\frac{0.98*0.95*0.90}{1}=0.8379$ (6)

    * –

the intercausal inference executes between effects with a common cause. In
Figure 6, if $Y$ is not instantiated, there exists intercausal inference in
diverging connections. Consider the inference $X\rightarrow Z$,

$P(X^{+}|YZ^{+})=\frac{P(X^{+}YZ^{+})}{P(YZ^{+})}=\frac{P(X^{+}Y^{+}Z^{+})+P(X^{+}Y^{-}Z^{+})}{P(Y^{+}Z^{+})+P(Y^{-}Z^{+})}$;

$=\frac{0.98*0.95*0.90+0.02*0.01*0.03}{0.98*0.90+0.02*0.03}=0.94936$.

  * •

in Converging Connections,

    * –

the forward inference executes with the evidence forward propagation. For
example, in Figure 7, consider the inference $(XZ)\rightarrow Y$, $P(Y|XZ)$ is
easy to obtain by the definition of Bayesian Network in by nature.

Figure 7: Inference in Converging Connection

    * –

the backward inference executes with the evidence backward propagation. For
example, in Figure 7, consider the inference $Y\rightarrow(XZ)$.

$P(Y)=\sum_{XZ}P(XYZ)=\sum_{XZ}(P(Y|XZ)P(XZ))$,

$P(XZ|Y)=\frac{P(Y|XZ)P(XZ)}{P(Y)}=\frac{P(Y|XZ)P(X)P(Z)}{\sum_{XZ}(P(Y|XZ)P(XZ))}$.

Finally,

$P(X|Y)=\sum_{Z}P(XZ|Y)$,

$P(Z|Y)=\sum_{X}P(XZ|Y)$.

    * –

the intercausal inference executes between causes with a common effect, and
the intermediate node is instantiated, then $P(Y^{+})=1$ or $P(Y^{-})=1$. In
Figure 7, consider the inference $X\rightarrow Z$, suppose $P(Y^{+})=1$,

$P(Z^{+}|X^{+}Y^{+})=\frac{P(Z^{+}X^{+}Y^{+})}{P(X^{+}Y^{+})}=\frac{P(Z^{+}X^{+}Y^{+})}{\sum_{Z}P(X^{+}Y^{+}Z)}$;

$P(Z^{+}X^{+}Y^{+})=P(X^{+})P(Z^{+})P(Y^{+}|X^{+}Z^{+})$;

$\sum_{Z}P(X^{+}YZ)=P(X^{+}Y^{+}Z^{+})+P(X^{+}Y^{+}Z^{-})$;

$P(Z^{+}|X^{+}Y^{+})=\frac{P(Z^{+}X^{+}Y^{+})}{\sum_{Z}P(X^{+}Y^{+}Z)}=\frac{P(X^{+})P(Z^{+})P(Y^{+}|X^{+}Z^{+})}{P(X^{+}Y^{+}Z^{+})+P(X^{+}Y^{+}Z^{-})}.$

####  4.1.2 inference in complex model

For complex models in Bayesian networks, there are single-connected networks,
multiple-connected, or event looped networks. It is possible to use some
methods, such as Triangulated Graphs, Clustering and Join Trees [Bertele &
Brioschi, 1972] [Finn & Thomas, 2007] [Golumbic, 1980], etc., to simplify them
into a polytree. Once a polytree is obtained, the inference can be executed by
the following approaches.

Polytrees have at most one path between any pair of nodes; hence they are also
referred to as singly-connected networks.

Suppose $X$ is the query node, and there is some set of evident nodes
$E,X\notin E$. The posterior probability (belief) is denoted as
$\mathbb{B}(X)=P(X|E)$, see Figure 8.

Figure 8: Evidence Propagation in Polytree

$E$ can be splitted into 2 parts: $E^{+}$ and $E^{-}$. $E^{-}$ is the part
consisting of assignments to variables in the subtree rooted at $X$, $E^{+}$
is the rest of it.

$\pi_{X}(E^{+})=P(X|E^{+})$

$\lambda_{X}(E^{-})=P(E^{-}|X)$

$\mathbb{B}(X)=P(X|E)=P(X|E^{+}E^{-})=\frac{P(E^{-}|XE^{+})P(X|E^{+})}{P(E^{-}|E^{+})}=\frac{P(E^{-}|X)P(X|E^{+})}{P(E^{-}|E^{+})}=\alpha\pi_{X}(E^{+})\lambda_{X}(E^{-})$
(7)

$\alpha$ is a constant independent of $X$.

where

$\lambda_{X}(E^{-})=\\{\begin{array}[]{cc}1&if\ evidence\ is\ X=x_{i}\\\ 0&if\
evidence\ is\ for\ another\ x_{j}\\\ \end{array}$ (8)

$\pi_{X}(E^{+})=\sum_{u_{1},...,u_{m}}P(X|u_{1},...,u_{m})\prod_{i}\pi_{X}(u_{i})$
(9)

  1. 1.

Forward inference in Polytree

Node $X$ sends $\pi$ messages to its children.

$\pi_{X}(U)=\\{\begin{array}[]{cc}1&if\ x_{i}\in X\ is\ entered\\\ 0&if\
evidentce\ is\ for\ another\ value\ x_{j}\\\
\sum_{u_{1},...u_{m}}P(X|u_{1},...u_{m})\prod_{i}\pi_{X}(u_{i})&otherwise\end{array}$
(10)

  2. 2.

Backward inference in Polytree Node $X$ sends new $\lambda$ messages to its
parents.

$\lambda_{X}(Y)=\prod_{y_{j}\in Y}[\sum_{j}P(y_{j}|X)\lambda_{X}(y_{j})]$ (11)

###  4.2 Related Algorithms for Probabilistic Inference

Various types of inference algorithms exist for Bayesian networks [Lauritzen &
Spiegelhalter, 1988] [Pearl, 1988] [Pearl, 2000] [Neal, 1993]. Each class
offers different properties and works better on different classes of problems,
but it is very unlikely that a single algorithm can solve all possible problem
instances effectively. Every resolution is always based on a particular
requirement. It is true that almost all computational problems and
probabilistic inference using general Bayesian networks have been shown to be
NP-hard by Cooper [Cooper, 1990].

In the early 1980’s, Pearl published an efficient message propagation
inference algorithm for polytrees [Kim & Pearl, 1983] [Peal, 1986]. The
algorithm is exact, and has polynomial complexity in the number of nodes, but
works only for singly connected networks. Pearl also presented an exact
inference algorithm for multiple connected networks called loop cutset
conditioning algorithm [Peal, 1986]. The loop cutset conditioning algorithm
changes the connectivity of a network and renders it singly connected by
instantiating a selected subset of nodes referred to as a loop cutset. The
resulting single connected network is solved by the polytree algorithm, and
then the results of each instantiation are weighted by their prior
probabilities. The complexity of this algorithm results from the number of
different instantiations that must be considered. This implies that the
complexity grows exponentially with the size of the loop cutest being
$O(d^{c})$, where $d$ is the number of values that the random variables can
take, and $c$ is the size of the loop cutset. It is thus important to minimize
the size of the loop cutset for a multiple connected network. Unfortunately,
the loop cutset minimization problem is NP-hard. A straightforward application
of Pearl’s algorithm to an acyclic digraph comprising one or more loops
invariably leads to insuperable problems [Koch & Westphall, 2001] [Neal,
1993].

Another popular exact Bayesian network inference algorithm is Lauritzen and
Spiegelhalter’s clique-tree propagation algorithm [Lauritzen & Spiegelhalter,
1988]. It is also called a ”clustering” algorithm. It first transforms a
multiple connected network into a clique tree by clustering the triangulated
moral graph of the underlying undirected graph and then performs message
propagation over the clique tree. The clique propagation algorithm works
efficiently for sparse networks, but still can be extremely slow for dense
networks. Its complexity is exponential in the size of the largest clique of
the transformed undirected graph.

In general, the existent exact Bayesian network inference algorithms share the
property of run time exponentiality in the size of the largest clique of the
triangulated moral graph, which is also called the induced width of the graph
[Lauritzen & Spiegelhalter, 1988].

##  5 Conclusion

This chapter summarizes the popular inferences methods in Bayesian networks.
The results demonstrates that the evidence can propagated across the Bayesian
networks by any links, whatever it is forward or backward or intercausal
style. The belief updating of Bayesian networks can be obtained by various
available inference techniques. Theoretically, exact inferences in Bayesian
networks is feasible and manageable. However, the computing and inference is
NP-hard. That means, in applications, in complex huge Bayesian networks, the
computing and inferences should be dealt with strategically and make them
tractable. Simplifying the Bayesian networks in structures, pruning unrelated
nodes, merging computing, and approximate approaches might be helpful in the
inferences of large scale Bayeisan networks.

## References

  * [Balke & Pearl, 1995] A. Balke and J. Pearl. Counterfactuals and policy analysis in structural models. Proceedings of the Eleventh Conference on Uncertainty in Artificial Intelligence, pages 11-18, 1995. Morgan Kaufmann.
  * [Bertele & Brioschi, 1972] Bertele, U. and Brioschi, F. (1972). Nonserial Dynamic Programming. Academic Press, London, ISBN-13: 978-0120934508.
  * [Cooper, 1990] G. Cooper. Computational complexity of probabilistic inference using Bayesian belief networks. Artificial Intelligence, 42:393-405, 1990.
  * [Finn & Thomas, 2007 ] Finn V. Jensen and Thomas D. Nielsen (2007). Bayesian Networks and Decision Graphs. Springer, ISBN-13:978-0-387-68281-5.
  * [Golumbic, 1980] Golumbic, M. C. (1980). Algorithmic Graph Theory and Perfect Graphs. Academic Press, London,ISBN-13: 978-0122892608.
  * [Goldszmidt & Pearl, 1996] M. Goldszmidt and J. Pearl. Qualitative Probabilities for Default Reasoning, Belief Revision, and Causal Modeling. Artificial Intelligence, 84(1-2): 57-112, July 1996.
  * [Jordan et al., 1998] M. I. Jordan, Z. Ghahramani, T. S. Jaakkola and L. K. Saul. An Introduction to Variational Methods for Graphical Models. M. I. Jordan (Ed.), Learning in Graphical Models. Kluwer, Dordrecht, The Netherlands, 1998.
  * [Kim & Pearl, 1983] Jin H. Kim and Judea Pearl. A computational model for combined causal and diagnostic reasoning in inference systems. In Proceedings of the Eighth International Joint Conference on Artificial Intelligence (IJCAI-83), pages 190-193, 1983. Morgan Kaufmann.
  * [ Koch & Westphall, 2001] F. L. Koch, and C. B. Westphall. Decentralized Network Management Using Distributed Artificial Intelligence. Journal of Network and systems management, Vol. 9, No. 4, December 2001.
  * [Lauritzen & Spiegelhalter, 1988] S. L. Lauritzen and D. J. Spiegelhalter. Local Computations with Probabilities on Graphical Structures and Their Application to Expert Systems. Journal of the Royal Statistical Society, Series B 50:157-224, 1988.
  * [Neal, 1993] R. M. Neal, Probabilistic Inference Using Markov Chain Monte Carlo methods, Tech. Rep. CRG-TR93-1, University of Toronto, Department of Computer Science, 1993.
  * [Peal, 1986] J. Pearl. A constraint-propagation approach to probabilistic reasoning, Uncertainty in Artificial Intelligence. North-Holland, Amsterdam, pages 357-369, 1986.
  * [Pearl, 1987] J. Pearl. Evidential Reasoning Using Stochastic Simulation of Causal Models. Artificial Intelligence, 32:247-257, 1987.
  * [Pearl, 1988] J. Pearl. Probabilistic Reasoning in Intelligent Systems: Networks of Plausible Inference. Morgan Kaufmann, San Mateo, CA, 1988.
  * [Pearl, 2000] J. Pearl. Causality: Models, Reasoning, and Inference. Cambridge, England: Cambridge University Press. New York, NY, 2000.
  * [Pearl & Verma, 1991] J. Pearl and T. Verma. A theory of inferred causation. J. A. Allen, R. Fikes and E. Sandewall (Eds.), Principles of Knowledge Representation and Reasoning. Proceedings of the Second International Conference, pages 441-452. Morgan Kaufmann, San Mateo, CA, 1991.
  * [Spirtes et al., 1993] P. Spirtes, C. Glymour and R. Scheines. Causation, Prediction, and Search. Springer-Verlag, New York, 1993.
  * [Tenenbaum & Griffiths, 2001] J. B. Tenenbaum and T. L. Griffiths. Structure learning in human causal induction. Advances in Neural Information Processing Systems, volume 13, Denver, Colorado, 2001. MIT Press.


Stability analysis and error estimates of local discontinuous Galerkin method
for convection-diffusion equations on overlapping mesh with non-periodic
boundary conditions111Supported by NSF grant DMS-1818467 and KMITL Research
Fund, Research Seed Grant for New Lecturer.

Nattaporn Chuenjarern222Department of Mathematics, Faculty of Science, King
Mongkut’s Institute of Technology Ladkrabang, Bangkok, Thailand 10520. E-mail:
<EMAIL_ADDRESS>Kanognudge Wuttanachamsri333Department of
Mathematics, Faculty of Science, King Mongkut’s Institute of Technology
Ladkrabang, Bangkok, Thailand 10520. E-mail<EMAIL_ADDRESS>Yang
Yang444Department of Mathematical Sciences, Michigan Technological University,
Houghton, MI 49931. E-mail<EMAIL_ADDRESS>
Abstract

A new local discontinuous Galerkin (LDG) method for convection-diffusion
equations on overlapping meshes with periodic boundary conditions was
introduced in [14]. With the new method, the primary variable $u$ and the
auxiliary variable $p=u_{x}$ are solved on different meshes. In this paper, we
will extend the idea to convection-diffusion equations with non-periodic
boundary conditions, i.e. Neumann and Dirichlet boundary conditions. The main
difference is to adjust the boundary cells. Moreover, we study the stability
and suboptimal error estimates. Finally, numerical experiments are given to
verify the theoretical findings.

Keywords: Local discontinuous Galerkin method, Stability, Error analysis,
Overlapping meshes

##  1 Introduction

In this paper, we apply the local discontinuous Galerkin (LDG) method on
overlapping meshes provided in [14] for the convection-diffusion equations

$u_{t}+f(u)_{x}=(a^{2}(u)u_{x})_{x},\>\>\>\>x\in[0,1],\>\>t>0,$ (1.1)

as well as its two dimensional version. We assume that $a(u)\geq 0$.

In 1973, Reed and Hill first introduced the discontinuous Galerkin (DG) method
in the framework of neutron linear transportation [23]. This method gained
even greater popularity for good stability, high-order accuracy, and
flexibility on h-p adaptivity and complex geometry. Subsequently, Cockburn et.
al. proposed in a series of papers [6, 7, 8, 9] the Runge-Kutta discontinuous
Galerkin (RKDG) methods for hyperbolic conservation laws. Later, in [10],
Cockburn and Shu introduced the LDG method for convection-diffusion equations
motivated by successful solving compressible Navier-Stokes equations in [1].

As in traditional LDG method, we introduce an auxiliary variable $p=A(u)_{x}$
with $A(u)=\int^{u}a(s)\ ds$ to represent the derivative of the primary
variable $u$, and rewrite (1.1) into the following system of first order
equations

$\displaystyle\begin{cases}u_{t}+f(u)_{x}=(a(u)p)_{x},\\\
p=A(u)_{x}.\end{cases}$ (1.2)

Then we can solve $u$ and $p$ on the same mesh by using the DG method. The LDG
method shares all the nice features of the DG methods for hyperbolic
equations, and it becomes one of the most popular numerical methods for
solving convection-diffusion equations. However, due to the discontinuity
nature of the numerical approximations, it may not be easy to construct and
analyze the scheme for some specials convection-diffusion equations. For
example, the convection terms of chemotaxis model [19, 22] and miscible
displacements in porous media [11, 12] are products of one of the primary
variable and the derivative of another one. Therefore, the upwind flux for the
convection term may not be easy to obtain. One of the alternatives is to use
other methods, such as mixed finite element method, to obtain continuous
approximations of the derivatives, see e.g. [18]. A more general idea is to
use the Lax-Friedrichs flux, see e.g. [16, 20, 28] for the error estimates for
miscible displacements and chemotaxis models. The main technique is to use the
diffusion term to control the convection term [24, 25, 26]. Moreover, to make
the numerical solutions to be physically relevant, we have to add a
sufficiently large penalty which depends on the numerical approximations of
the derivatives of the primary variables [17, 20, 2]. Another possible way is
to construct flux-free schemes, such as the central discontinuous Galerkin
(CDG) method [21] and the staggered discontinuous Galerkin (SDG) method [4].
However, the CDG scheme doubles the computational cost as we have to solve
each equation in (1.2) on both the primary and dual meshes twice and it is not
easy to apply limiters in SDG method because it requires partial continuity of
the numerical approximations.

Recently, one of the authors in this paper introduced a new LDG method on
overlapping meshes [14] by solving $u$ and $p$ on primitive and dual meshes,
respectively, hence $p$ is continuous across the interfaces on the primitive
mesh. The scheme is proved to be stable under the $L^{2}$-norm and can be used
to construct third-order maximum-principle-preserving schemes [13]. However,
in some special cases, it may not enjoy the optimal convergence rates. The
suboptimal convergence rate can be observed numerically if all the following
three conditions are satisfied: (1) Odd order polynomials are used in the
finite element space, (2) The dual mesh generated by connecting the midpoints
of the primitive mesh, (3) No penalty is added to the numerical scheme. If one
of the conditions is violated, the convergence rate will turn out to be
optimal. Later, in [3], we used Fourier analysis to explicitly write out the
error between the numerical and exact solutions and verify the optimal
convergence rate for linear parabolic equations with periodic boundary
conditions in one space dimension. Moreover, we also found out some
superconvergence points that may depend on the perturbation constant in the
construction of the dual mesh.

Both works given above are for problems with periodic boundary conditions. To
implement the scheme, we need to combine the two boundary cells at the
boundaries into one and find a polynomial approximation on the new cell. It is
impossible to do that for general Dirichlet and Neumann boundary conditions,
which are more realistic in practice, see e.g. [11, 12, 20]. In this paper, we
will discuss the stability and error estimates of the new LDG methods for
problems with Neumann and Dirichlet boundary conditions. The difficulty for
the Neumann and Dirichlet boundary conditions is how to deal with the boundary
cells of the dual mesh since two boundary cells cannot be combined. One
possible way is to leave two boundary cells after generating the dual mesh,
and introduce suitable numerical fluxes are the boundaries. For simplicity of
presentation, we only demonstrate the proof for nonlinear parabolic equations

$\displaystyle\begin{cases}u_{t}=(a(u)p)_{x},\\\ p=A(u)_{x},\end{cases}$ (1.3)

where $\displaystyle A(u)=\int^{u}a(t)\>dt$. The extension to general
nonlinear convection-diffusion equations can be obtained following [14], hence
we only demonstrate the results without proof.

The rest of the paper is organized as follows: we first discuss the LDG scheme
on overlapping mesh in Section 2. In Section 3, we demonstrate the stability
analysis of the scheme for the Neumann and Dirichlet boundary conditions. The
error estimate will be provided in Section 4. The extension to problems in two
space dimensions will be discussed in Section 5. In Section 6, the numerical
experiments will be given to demonstrate the accuracy of the scheme on non-
periodic boundary conditions. We will end in Section 7 with concluding
remarks.

##  2 Preliminary

In this section, we proceed to demonstrate the new LDG method for solving the
one-dimensional diffusion equation (1.3) on overlapping meshes with three
different boundary conditions, i.e. periodic, Neumann and Dirichlet boundary
conditions.

$I_{i}$$I_{i+1}$$P_{i+1/2}$$x_{i-1/2}$$x_{i+1/2}$$x_{i+3/2}$$\tilde{x}_{i}$$\tilde{x}_{i+1}$
Figure 1: Overlapping meshes

###  2.1 Overlapping meshes

The new LDG method solves the variables $u$ and $p$ on two different meshes as
shown in Figure 1. First, we define the primitive mesh on which the primary
variable $u$ is solved. We give a partition of the computational domain
$\Omega=[0,1]$ as
$0=x_{\frac{1}{2}}<x_{\frac{3}{2}}<\cdots<x_{N+\frac{1}{2}}=1$, and denote the
$i-$th cell as

$I_{i}=[x_{i-\frac{1}{2}},x_{i+\frac{1}{2}}],\>\>\>\>i=1,...,N.$

Moreover, we denote

$\Delta
x_{i}=x_{i+\frac{1}{2}}-x_{i-\frac{1}{2}},\>\>\>\>x_{i}=\frac{x_{i+\frac{1}{2}}+x_{i-\frac{1}{2}}}{2}$

as the cell length and the cell center of $I_{i}$, respectively, and define
$\displaystyle\Delta x=\max_{i}\Delta x_{i}$.

We now demonstrate how to create the dual mesh, namely the P-mesh in [14], for
solving the auxiliary variable $p$ for problems with periodic boundary
conditions. We choose a point $\tilde{x}_{i}$ given as

$\tilde{x}_{i}=x_{i}+\frac{\Delta
x_{i}}{2}\xi_{0,i},\>\>\>\>\xi_{0,i}\in[-1,1],\>\>\>i=1,...,N.$ (2.1)

It is easy to see that
$\tilde{x}_{i}\in[x_{i-\frac{1}{2}},x_{i+\frac{1}{2}}]$. Define

$P_{i-\frac{1}{2}}=[\tilde{x}_{i-1},\tilde{x}_{i}],\>\>\>\>i=1,...,N,$

as the ($i-\frac{1}{2})$-th cell of the dual mesh where we denote
$\tilde{x}_{0}=\tilde{x}_{N}-1$. Moreover,

$\Delta\tilde{x}_{i-\frac{1}{2}}=\tilde{x}_{i}-\tilde{x}_{i-1},\>\>\>\>\tilde{x}_{i-\frac{1}{2}}=\frac{\tilde{x}_{i}+\tilde{x}_{i-1}}{2}$

stand for the cell length and the cell center of $P_{i-\frac{1}{2}}$,
respectively. The P-mesh will consist of all $P$ cells. For periodic boundary
conditions, we can also define
$P_{\frac{1}{2}}=[0,\tilde{x}_{1}]\cup[\tilde{x}_{N},1]$ as stated in [14],
see Figure 4.

Next, we demonstrate how to deal with the boundary cells on the P-mesh for
problems with Neumann and Dirichlet boundary conditions. We leave the left and
the right boundary cells and define $P_{\frac{1}{2}}=[0,\tilde{x}_{1}]$ and
$P_{N+\frac{1}{2}}=[\tilde{x}_{N},1]$, see Figure 4. This mesh is called the
L-mesh. The other way is to combine the boundary cells with their neighbour as
$P_{\frac{3}{2}}=[0,\tilde{x}_{2}]$ and
$P_{N-\frac{1}{2}}=[\tilde{x}_{N-1},1]$, see Figure 4. This mesh is called
C-mesh.

If we take $\xi_{0,i}$ as a constant independent of $i$, it is easy to check
that

$\min\\{\Delta x_{i-1},\Delta
x_{i}\\}\leq\Delta\tilde{x}_{i-\frac{1}{2}}\leq\max\\{\Delta x_{i-1},\Delta
x_{i}\\},$

and hence we have $\max_{i}\Delta\tilde{x}_{i-\frac{1}{2}}\leq\Delta x.$

$I_{i}$$I_{i+1}$$I_{1}$$I_{2}$$I_{N-1}$$I_{N}$$P_{i+1/2}$$P_{3/2}$$P_{1/2}$$P_{N-1/2}$$P_{1/2}$$x_{1/2}$$x_{3/2}$$x_{5/2}$$x_{i-1/2}$$x_{i+1/2}$$x_{i+3/2}$$x_{N-3/2}$$x_{N-1/2}$$x_{N+1/2}$$\tilde{x}_{0}$$\tilde{x}_{1}$$\tilde{x}_{2}$$\tilde{x}_{i}$$\tilde{x}_{i+1}$$\tilde{x}_{N-1}$$\tilde{x}_{N}$$\tilde{x}_{0}$

Figure 2: Overlapping meshes for periodic boundary conditions.

$I_{i}$$I_{i+1}$$I_{1}$$I_{2}$$I_{N-1}$$I_{N}$$P_{i+1/2}$$P_{3/2}$$P_{1/2}$$P_{N-1/2}$$P_{N+1/2}$$x_{1/2}$$x_{3/2}$$x_{5/2}$$x_{i-1/2}$$x_{i+1/2}$$x_{i+3/2}$$x_{N-3/2}$$x_{N-1/2}$$x_{N+1/2}$$\tilde{x}_{0}$$\tilde{x}_{1}$$\tilde{x}_{2}$$\tilde{x}_{i}$$\tilde{x}_{i+1}$$\tilde{x}_{N-1}$$\tilde{x}_{N}$$\tilde{x}_{N+1}$

Figure 3: Overlapping mesh 1, L-mesh.

$I_{i}$$I_{i+1}$$I_{1}$$I_{2}$$I_{N-1}$$I_{N}$$P_{i+1/2}$$P_{3/2}$$P_{N-1/2}$$x_{1/2}$$x_{3/2}$$x_{5/2}$$x_{i-1/2}$$x_{i+1/2}$$x_{i+3/2}$$x_{N-3/2}$$x_{N-1/2}$$x_{N+1/2}$$\tilde{x}_{1}$$\tilde{x}_{2}$$\tilde{x}_{i}$$\tilde{x}_{i+1}$$\tilde{x}_{N-1}$$\tilde{x}_{N}$

Figure 4: Overlapping mesh2, C-mesh.

###  2.2 Norms

In this section, we define some norms that will be used throughout the paper.

For any interval $I$, we define $||u||_{I}$ and $||u||_{\infty,I}$ to be the
standard $L^{2}$\- and $L^{\infty}$\- norms of $u$ on $I$, respectively. For
any natural number $\ell$, we consider the norm of Sobolev space
$H^{\ell}(I)$, defined by

$||u||_{\ell,I}=\left\\{\sum_{0\leq\beta\leq\ell}\left\|\frac{\partial^{\beta}u}{\partial
x^{\beta}}\right\|^{2}_{I}\right\\}^{\frac{1}{2}}.$

For convenience, if $I$ is the whole computational domain, then the
corresponding subscript will be omitted.

Moreover, for any $u\in C(I_{i})$, we define

$||u||_{\Gamma_{i}}=|u^{-}_{i+\frac{1}{2}}|+|u^{+}_{i-\frac{1}{2}}|.$

Similarly, for any $u\in C(P_{i-\frac{1}{2}})$, we define

$||u||_{\Gamma_{i-\frac{1}{2}}}=|u^{-}_{i-1}|+|u^{+}_{i}|.$

###  2.3 Numerical schemes with Neumann and Dirichlet boundary conditions

In this section, we consider the LDG method for the following nonlinear
diffusion equation

$\begin{cases}u_{t}=(a(u)p)_{x},\>\>\>x\in[0,1],\>\>t>0,\\\
p=A(u)_{x},\>\>\>x\in[0,1],\\\ \end{cases}$ (2.2)

subject to the following boundary conditions

  1. Neumann boundary condition:

$u_{x}(0,t)=u_{x}(1,t)=0.$ (2.3)

  2. Dirichlet boundary condition:

$u(0,t)=u(1,t)=0.$ (2.4)

We define the finite element spaces to be

$\displaystyle V_{h}$ $\displaystyle=\\{u_{h}:u_{h}|_{I_{i}}\in
P^{k}(I_{i}),i=1,\ldots,N\\},$ $\displaystyle P^{L}_{h}$
$\displaystyle=\\{p_{h}:p_{h}|_{P_{i-\frac{1}{2}}}\in
P^{k}(P_{i-\frac{1}{2}}),i=1,\ldots,N+1\\},$ $\displaystyle P^{C}_{h}$
$\displaystyle=\\{p_{h}:p_{h}|_{P_{i-\frac{1}{2}}}\in
P^{k}(P_{i-\frac{1}{2}}),i=2,\ldots,N\\},$

where $P^{L}_{h}$ and $P^{C}_{h}$ are finite element spaces for the L-mesh and
C-mesh, respectively, and $P^{k}$ is the space of polynomials of degree up to
$k$.

Multiplying (2.2) with the test functions and using the integration by parts,
we obtain

$\displaystyle\int_{I_{i}}(u_{h})_{t}vdx$
$\displaystyle=-\int_{I_{i}}a(u_{h})p_{h}v_{x}dx+\hat{a}_{i+\frac{1}{2}}\widehat{p_{h}}_{i+\frac{1}{2}}v^{-}_{i+\frac{1}{2}}-\hat{a}_{i-\frac{1}{2}}\widehat{p_{h}}_{i-\frac{1}{2}}v^{+}_{i-\frac{1}{2}},$
(2.5) $\displaystyle\int_{P_{i-\frac{1}{2}}}p_{h}wdx$
$\displaystyle=-\int_{P_{i-\frac{1}{2}}}A(u_{h})w_{x}dx+A(u_{h}(\tilde{x}_{i}))w^{-}_{i}-A(u_{h}(\tilde{x}_{i-1}))w^{+}_{i-1}.$
(2.6)

Then the LDG method on overlapping meshes is defined as follows:

  * •

L-mesh: find $(u_{h},p_{h})\in V_{h}\times P^{L}_{h}$, such that for any test
functions $(v,w)\in V_{h}\times P^{L}_{h}$ we have (2.5) and (2.6),

  * •

C-mesh: find $(u_{h},p_{h})\in V_{h}\times P^{C}_{h}$, such that for any test
functions $(v,w)\in V_{h}\times P^{C}_{h}$ we have (2.5) and (2.6),

where $v^{-}_{i+\frac{1}{2}}=v^{-}(x_{i+\frac{1}{2}})$ and
$w^{-}_{i}=w^{-}(\tilde{x}_{i})$. Likewise for $v^{+}_{i-\frac{1}{2}}$ and
$w^{+}_{i-1}$.

We denote the jump of the function $s$ across the cell interface
$x=x_{i-\frac{1}{2}}$ as
$[s]_{i-\frac{1}{2}}=s^{+}_{i-\frac{1}{2}}-s^{-}_{i-\frac{1}{2}}$. Similarly,
$[w]_{i}=w^{+}_{i}-w^{-}_{i}$ denotes the jump of the function $w$ across the
cell interface $x=\tilde{x}_{i}$ on the P-mesh. Due to the boundary conditions
we also define the numerical fluxes of $u$ and $p$ on L-meshes as

  1. Neumann boundary condition: $u^{-}_{\frac{1}{2}}=u^{+}_{\frac{1}{2}}$, $u^{+}_{N+\frac{1}{2}}=u^{-}_{N+\frac{1}{2}}$, $p^{-}_{1}=p^{+}_{1}$ and $p^{+}_{N}=p^{-}_{N}$.

  2. Dirichlet boundary condition: $u^{-}_{\frac{1}{2}}=0$, $u^{+}_{N+\frac{1}{2}}=0$, $p^{-}_{1}=0$ and $p^{+}_{N}=0$.

The definition of the numerical fluxes at the boundary for C-mesh is similar.
The numerical flux $\hat{a}$ at the point $x_{i-\frac{1}{2}}$ is defined as

$\hat{a}_{i-\frac{1}{2}}=\frac{[A(u_{h})]_{i-\frac{1}{2}}}{[u_{h}]_{i-\frac{1}{2}}}.$
(2.7)

Also, we choose the numerical flux $\widehat{p_{h}}_{i-\frac{1}{2}}$ as the
value of $p_{h}$ evaluated at $x=x_{i-\frac{1}{2}}$ with the penalty term

$\widehat{p_{h}}_{i-\frac{1}{2}}=p_{h}(x_{i-\frac{1}{2}})+\frac{\alpha_{i-\frac{1}{2}}}{\Delta\tilde{x}_{i-\frac{1}{2}}}[u_{h}]_{i-\frac{1}{2}}.$
(2.8)

Notice that $p_{h}$ is continuous at the interfaces of the primitive cells and
hence $p_{h}(x_{i-\frac{1}{2}})$ is well defined.

Finally, we define

$\displaystyle H_{u}(u_{h},p_{h},v)$
$\displaystyle=-\sum_{i=1}^{N}\int_{I_{i}}a(u_{h})p_{h}v_{x}dx+\sum_{i=1}^{N}\left(\hat{a}_{i+\frac{1}{2}}\hat{p}_{i+\frac{1}{2}}v^{-}_{i+\frac{1}{2}}-\hat{a}_{i-\frac{1}{2}}\hat{p}_{i-\frac{1}{2}}v^{+}_{i-\frac{1}{2}}\right),$
(2.9) $\displaystyle H_{p}(u_{h},w)$
$\displaystyle=-\sum_{i=1}^{N+1}\int_{P_{i-\frac{1}{2}}}A(u_{h})w_{x}dx+\sum_{i=1}^{N+1}\left(A(u_{h}(\tilde{x}_{i}))w^{-}_{i}-A(u_{h}(\tilde{x}_{i-1}))w^{+}_{i-1}\right),$
(2.10)

for L-mesh. For C-mesh, the two summations in (2.10) are from $i=2$ to $N$
Then the LDG scheme can be rewritten as

$\displaystyle\int_{\Omega}(u_{h})_{t}vdx$
$\displaystyle=H_{u}(u_{h},p_{h},v),$ (2.11)
$\displaystyle\int_{\Omega}p_{h}wdx$ $\displaystyle=H_{p}(u_{h},w).$ (2.12)

##  3 Stability analysis

In this section, we demonstrate the stability of the new LDG method on
overlapping meshes with non-periodic boundary conditions.

###  3.1 Neumann boundary condition

In this subsection, the stability of the new LDG method on overlapping meshes
with the Neumann boundary conditions will be demonstrated.

######  Lemma 3.1.

Suppose $H_{u}$ and $H_{p}$ are defined in (2.9) and (2.10), respectively, and
the dual mesh is given as either the L-mesh or the C-mesh, then we have

$\displaystyle
H_{u}(u_{h},p_{h},v)+H_{p}(u_{h},w)=-\sum_{i=2}^{N}\frac{[A(u_{h})]_{i-\frac{1}{2}}}{[u_{h}]_{i-\frac{1}{2}}}\frac{\alpha_{i-\frac{1}{2}}}{\Delta\tilde{x}_{i-\frac{1}{2}}}[u_{h}]^{2}_{i-\frac{1}{2}}.$
(3.1)

###### Proof.

We only prove for L-mesh. The case for C-mesh is basically the same. Taking
$w=p_{h}$ in (2.10), using integration by parts and applying the Neumann
boundary conditions, we obtain

$\displaystyle H_{p}(u_{h},p_{h})$
$\displaystyle=-\sum_{i=1}^{N+1}\int_{P_{i-\frac{1}{2}}}A(u_{h})(p_{h})_{x}dx+\sum_{i=1}^{N+1}\left(A(u_{h}(\tilde{x}_{i}))(p_{h}^{-})_{i}-A(u_{h}(\tilde{x}_{i-1}))(p_{h}^{+})_{i-1}\right)$
$\displaystyle=-\int_{P_{\frac{1}{2}}}A(u_{h})(p_{h})_{x}dx+A(u_{h}(\tilde{x}_{1}))(p_{h}^{-})_{1}-A(u_{h}(\tilde{x}_{0}))(p_{h}^{+})_{0}$
$\displaystyle\>\>\>\>\>-\sum_{i=2}^{N}\int^{x_{i-\frac{1}{2}}}_{\tilde{x}_{i-1}}A(u_{h})(p_{h})_{x}dx-\sum_{i=2}^{N}\int_{x_{i-\frac{1}{2}}}^{\tilde{x}_{i}}A(u_{h})(p_{h})_{x}dx$
$\displaystyle\>\>\>\>\>+\sum_{i=2}^{N}\left(A(u_{h}(\tilde{x}_{i}))(p_{h}^{-})_{i}-A(u_{h}(\tilde{x}_{i-1}))(p_{h}^{+})_{i-1}\right)$
$\displaystyle\>\>\>\>\>-\int_{P_{N+\frac{1}{2}}}A(u_{h})(p_{h})_{x}dx+A(u_{h}(\tilde{x}_{N+1}))(p_{h}^{-})_{N+1}-A(u_{h}(\tilde{x}_{N}))(p_{h}^{+})_{N}$
$\displaystyle=\int_{P_{\frac{1}{2}}}a(u_{h})(u_{h})_{x}p_{h}dx$
$\displaystyle\>\>\>\>\>+\sum_{i=2}^{N}\int^{x_{i-\frac{1}{2}}}_{\tilde{x}_{i-1}}a(u_{h})(u_{h})_{x}p_{h}dx+\sum_{i=2}^{N}\int_{x_{i-\frac{1}{2}}}^{\tilde{x}_{i}}a(u_{h})(u_{h})_{x}p_{h}dx+\sum_{i=2}^{N}[A(u_{h})]_{i-\frac{1}{2}}p_{h}(x_{i-\frac{1}{2}})$
$\displaystyle\>\>\>\>\>+\int_{P_{N+\frac{1}{2}}}a(u_{h})(u_{h})_{x}p_{h}dx$
$\displaystyle=\sum_{i=1}^{N}\int_{I_{i}}a(u_{h})(u_{h})_{x}p_{h}dx+\sum_{i=2}^{N}[A(u_{h})]_{i-\frac{1}{2}}p_{h}(x_{i-\frac{1}{2}}).$
(3.2)

Taking $v=u_{h}$ in (2.9), we obtain

$\displaystyle H_{u}(u_{h},p_{h},v)$
$\displaystyle=-\sum_{i=1}^{N}\int_{I_{i}}p_{h}a(u_{h})(u_{h})_{x}dx-\sum_{i=2}^{N}\hat{a}_{i-\frac{1}{2}}\hat{p}_{i-\frac{1}{2}}[u_{h}]_{i-\frac{1}{2}}$
$\displaystyle\>\>\>\>\>-\hat{a}_{\frac{1}{2}}\hat{p}_{\frac{1}{2}}(u_{h})^{+}_{\frac{1}{2}}+\hat{a}_{N+\frac{1}{2}}\hat{p}_{N+\frac{1}{2}}(u_{h})^{-}_{N+\frac{1}{2}}$
$\displaystyle=-\sum_{i=1}^{N}\int_{I_{i}}p_{h}a(u_{h})(u_{h})_{x}dx-\sum_{i=2}^{N}\frac{[A(u_{h})]_{i-\frac{1}{2}}}{[u_{h}]_{i-\frac{1}{2}}}\left(p_{h}(x_{i-\frac{1}{2}})+\frac{\alpha_{i-\frac{1}{2}}}{\Delta\tilde{x}_{i-\frac{1}{2}}}\right)[u_{h}]_{i-\frac{1}{2}}.$
(3.3)

where the last step, we obtain it by applying the Neumann boundary conditions.
Summing (3.2) and (3.3), we have (3.1) which further leads to the $L^{2}$
stability of the LDG method on overlapping meshes with Neumann boundary
conditions. ∎

###  3.2 Dirichlet boundary conditions

In this subsection, we will describe the stability of the new LDG method on
overlapping meshes for problems with Dirichlet boundary conditions.

With a minor change of the previous proof, we can obtain the following lemma.

######  Lemma 3.2.

Suppose $H_{u}$ and $H_{p}$ defined in (2.9) and (2.10), respectively, then we
have

$\displaystyle H_{u}(u_{h},p_{h},v)+H_{p}(u_{h},w)$
$\displaystyle=-\sum_{i=1}^{N+1}\frac{[A(u_{h})]_{i-\frac{1}{2}}}{[u_{h}]_{i-\frac{1}{2}}}\frac{\alpha_{i-\frac{1}{2}}[u_{h}]^{2}_{i-\frac{1}{2}}}{\Delta\tilde{x}_{i-\frac{1}{2}}}.$
(3.4)

The proof is very similar to that of Lemma 3.1, so we omit it and only
demonstrate the result as the following theorem.

######  Theorem 3.3.

The LDG method introduced (2.5) and (2.6) with the boundary conditions (2.3)
and (2.4) are stable and

$\frac{1}{2}\frac{d}{dt}||u_{h}||^{2}+||p_{h}||^{2}\leq 0.$

##  4 Error estimates

In this section, we demonstrate the error estimates. We will consider linear
equations. Moreover, for simplicity, we will discuss the problem with the
Neumann boundary condition on L-meshes only. For the other cases, we can apply
the same procedure.

First, we use $e$ to denote the error between the exact and numerical
solutions i.e. $e_{u}=u-u_{h}$ and $e_{p}=p-p_{h}$, then we can get the error
equations from (2.5) and (2.6) as

$\displaystyle\int_{I_{i}}(e_{u})_{t}vdx$
$\displaystyle=-\int_{I_{i}}e_{p}v_{x}dx+\widehat{e_{p}}_{i+\frac{1}{2}}v^{-}_{i+\frac{1}{2}}-\widehat{e_{p}}_{i-\frac{1}{2}}v^{+}_{i-\frac{1}{2}},$
(4.1) $\displaystyle\int_{P_{i-\frac{1}{2}}}e_{p}wdx$
$\displaystyle=-\int_{P_{i-\frac{1}{2}}}e_{u}w_{x}dx+e_{u}(\tilde{x}_{i})w^{-}_{i}-e_{u}(\tilde{x}_{i-1})w^{+}_{i-1}.$
(4.2)

Next, we introduce some basic properties of the finite element space that will
be used.

######  Lemma 4.1.

Assuming $u\in V_{h}$, there exists constant $C>0$ independent of $\Delta x$
and $u$ such that for $\beta\geq 1$

$||\partial^{\beta}_{x}u||_{I_{i}}\leq C\Delta
x_{i}^{-\beta}||u||_{I_{i}},\>\>\>\>||u||_{\Gamma_{i}}\leq C\Delta
x_{i}^{-1/2}||u||_{I_{i}}.$

Similarly, for any $u\in P^{L}_{h}$, there exists constant $C>0$ independent
of $\Delta\tilde{x}$ and $u$ such that for $\beta\geq 1$

$||\partial^{\beta}_{x}u||_{P_{i-\frac{1}{2}}}\leq
C\Delta\tilde{x}_{i-\frac{1}{2}}^{-\beta}||u||_{P_{i-\frac{1}{2}}},\>\>\>\>||u||_{\Gamma_{i-\frac{1}{2}}}\leq
C\Delta\tilde{x}_{i-\frac{1}{2}}^{-1/2}||u||_{P_{i-\frac{1}{2}}}.$

We also introduce the standard $L^{2}$ projections $P^{1}_{k}$ into $V_{h}$
and $P^{2}_{k}$ into $P^{L}_{h}$ by:

$\int_{I_{i}}P^{1}_{k}uvdx=\int_{I_{i}}uvdx,\>\>\forall v\in P^{k}(I_{i}),$

and

$\int_{P_{i-\frac{1}{2}}}P^{2}_{k}uvdx=\int_{P_{i-\frac{1}{2}}}uvdx,\>\>\forall
v\in P^{k}({P_{i-\frac{1}{2}}}),$

respectively. By the scaling argument, we obtain the following lemma [5]

######  Lemma 4.2.

Suppose the function $u(x)\in C^{k+1}(I_{i})$, then there exists positive
constant $C$ independent of $\Delta x$ and $u$, such that

$||u-P^{1}_{k}u||_{I_{i}}+\Delta x_{i}||(u-P^{1}_{k}u)_{x}||_{I_{i}}+\Delta
x^{1/2}_{i}||u-P^{1}_{k}u||_{\infty,{I_{i}}}\leq C\Delta
x^{k+1}_{i}||u||_{k+1,I_{i}}.$

Moreover, if $u(x)\in C^{k+1}(P_{i-\frac{1}{2}})$, then there exists positive
constant $C$ independent of $\Delta\tilde{x}$ and $u$, such that

$||u-P^{2}_{k}u||_{P_{i-\frac{1}{2}}}+\Delta\tilde{x}_{i-\frac{1}{2}}||(u-P^{2}_{k}u)_{x}||_{P_{i-\frac{1}{2}}}+\Delta\tilde{x}^{1/2}_{i-\frac{1}{2}}||u-P^{2}_{k}u||_{\infty,P_{i-\frac{1}{2}}}\leq
C\Delta\tilde{x}^{k+1}_{i-\frac{1}{2}}||u||_{k+1,P_{i-\frac{1}{2}}},$

where $||u||_{k+1,I}$ is the standard $H^{k+1}$-norm over the interval $I$.

As the general treatment of the finite element method, we split the errors
into two terms as

$e_{u}=\eta_{u}-\xi_{u},\>\>\>\>e_{p}=\eta_{p}-\xi_{p},$

where

$\eta_{u}=u-P^{1}_{k}u,\>\>\>\xi_{u}=u_{h}-P^{1}_{k}u,\>\>\>\eta_{p}=p-P^{2}_{k}p,\>\>\>\xi_{p}=p_{h}-P^{2}_{k}p.$

With the above notations, the equations (4.1) and (4.2) can be rewritten as

$\displaystyle\int_{I_{i}}(\xi_{u})_{t}vdx$
$\displaystyle=\int_{I_{i}}e_{p}v_{x}dx-\widehat{e_{p}}_{i+\frac{1}{2}}v^{+}_{i+\frac{1}{2}}-\widehat{e_{p}}_{i-\frac{1}{2}}v^{+}_{i-\frac{1}{2}}$
(4.3) $\displaystyle\int_{P_{i-\frac{1}{2}}}\xi_{p}wdx$
$\displaystyle=\int_{P_{i-\frac{1}{2}}}e_{u}w_{x}dx-
e_{u}(\tilde{x}_{i})w^{-}_{i}+e_{u}(\tilde{x}_{i-1})w^{+}_{i-1}$ (4.4)

Now, we can state the main theorem.

######  Theorem 4.3.

Suppose the exact solution $u\in C^{k+2}(\Omega)$ and the finite element space
is made up of piecewise polynomials of degree $k$. Moreover, the numerical
solutions satisfy (2.5) and (2.6). Then the error between the numerical and
exact solutions satisfies

$||e_{u}||+\int_{0}^{T}||e_{p}||dt\leq C\Delta x^{k},$

where $C$ is independent of $\Delta x$.

###### Proof.

Sum up (4.3) and (4.4) with $v=\xi_{u}$ and $w=\xi_{p}$, and then sum up over
$i$ to obtain

$\displaystyle\frac{1}{2}\frac{d}{dt}||\xi_{u}||^{2}+||\xi_{p}||^{2}$
$\displaystyle=\sum_{i=1}^{N}\int_{I_{i}}(\eta_{p}-\xi_{p})(\xi_{u})_{x}dx+\sum_{i=2}^{N}\left(\eta_{p_{i-\frac{1}{2}}}-\xi_{p_{i-\frac{1}{2}}}+\alpha_{i-\frac{1}{2}}\frac{[\eta_{u}-\xi_{u}]_{i-\frac{1}{2}}}{\Delta\tilde{x}_{i-\frac{1}{2}}}\right)[\xi_{u}]_{i-\frac{1}{2}}$
$\displaystyle\>\>\>\>\>-\sum_{i=1}^{N}\int_{I_{i}}(\eta_{u}-\xi_{u})(\xi_{p})_{x}dx-\sum_{i=2}^{N}[\eta_{u}-\xi_{u}]_{i-\frac{1}{2}}\xi_{p}(x_{i-\frac{1}{2}})$
$\displaystyle=\sum_{i=1}^{N}\int_{I_{i}}\eta_{p}(\xi_{u})_{x}dx+\sum_{i=2}^{N}\left(\eta_{p_{i-\frac{1}{2}}}+\alpha_{i-\frac{1}{2}}\frac{[\eta_{u}]_{i-\frac{1}{2}}}{\Delta\tilde{x}_{i-\frac{1}{2}}}\right)[\xi_{u}]_{i-\frac{1}{2}}$
$\displaystyle\>\>\>\>\>-\sum_{i=1}^{N}\int_{I_{i}}\eta_{u}(\xi_{p})_{x}dx-\sum_{i=2}^{N}[\eta_{u}]_{i-\frac{1}{2}}\xi_{p}(x_{i-\frac{1}{2}})+H_{u}^{N}(\xi_{u},\xi_{p},\xi_{u})+H_{p}^{N}(\xi_{u},\xi_{p})$
$\displaystyle=R_{1}+R_{2}+R_{3},$ (4.5)

where

$\displaystyle R_{1}$
$\displaystyle=\sum_{i=1}^{N}\int_{I_{i}}\eta_{p}(\xi_{u})_{x}dx-\sum_{i=1}^{N}\int_{I_{i}}\eta_{u}(\xi_{p})_{x}dx,$
$\displaystyle R_{2}$
$\displaystyle=\sum_{i=2}^{N}\left(\alpha_{i-\frac{1}{2}}\frac{[\eta_{u}]_{i-\frac{1}{2}}}{\Delta\tilde{x}_{i-\frac{1}{2}}}\right)[\xi_{u}]_{i-\frac{1}{2}}+H_{u}^{N}(\xi_{u},\xi_{p},\xi_{u})+H_{p}^{N}(\xi_{u},\xi_{p}),$
$\displaystyle R_{3}$
$\displaystyle=\sum_{i=2}^{N}\eta_{p_{i-\frac{1}{2}}}[\xi_{u}]_{i-\frac{1}{2}}-\sum_{i=2}^{N}[\eta_{u}]_{i-\frac{1}{2}}\xi_{p}(x_{i-\frac{1}{2}}).$

Now we estimate $R_{i}$ where $i=1,2,3$ term by term.

$\displaystyle R_{1}$
$\displaystyle\leq\sum_{i=1}^{N}\left(||\eta_{p}||_{I_{i}}||(\xi_{u})_{x}||_{I_{i}}+||(\eta_{u})_{x}||_{I_{i}}||\xi_{p}||_{I_{i}}\right)$
$\displaystyle\leq\sum_{i=1}^{N}\left(||\eta_{p}||_{P_{i-\frac{1}{2}}\cup
P_{i+\frac{1}{2}}}||\xi_{u}||_{I_{i}}+||(\eta_{u})_{x}||_{I_{i}}||\xi_{p}||_{I_{i}}\right)$
$\displaystyle\leq C\Delta
x^{k}\sum_{i=1}^{N}\left[\left(||p||_{k+1,P_{i-\frac{1}{2}}}+||p||_{k+1,P_{i+\frac{1}{2}}}\right)||\xi_{u}||_{I_{i}}+||u||_{k+1,I_{i}}||\xi_{p}||_{I_{i}}\right]$
$\displaystyle\leq C\Delta x^{k}\left(||\xi_{u}||+||\xi_{p}||\right),$ (4.6)

where we applied the Cauchy-Schwarz inequality to the first step. In the
second step, we used Lemmas 4.1 and 4.2. Also, the Cauchy-Schwarz inequality
was used again in the last step. Applying Lemma 3.1, we obtain the estimate of
$R_{2}$

$\displaystyle R_{2}$
$\displaystyle\leq\sum_{i=2}^{N}\frac{\alpha_{i-\frac{1}{2}}}{\Delta\tilde{x}_{i-\frac{1}{2}}}\left([\eta_{u}]_{i-\frac{1}{2}}[\xi_{u}]_{i-\frac{1}{2}}-[\xi_{u}]^{2}_{i-\frac{1}{2}}\right)$
$\displaystyle\leq C\sum_{i=2}^{N}\frac{\alpha_{i-\frac{1}{2}}}{\Delta
x}[\eta_{u}]^{2}_{i-\frac{1}{2}}$ $\displaystyle\leq
C\sum_{i=2}^{N}\alpha_{i-\frac{1}{2}}\Delta
x^{2k}\left(||u||^{2}_{I_{i-1}}+||u||^{2}_{I_{i}}\right)$ $\displaystyle\leq
C\Delta x^{2k},$ (4.7)

where Lemma 4.2 was applied in step 2. While steps 2 and 4 follow from direct
computation. Finally, we estimate $R_{3}$.

$\displaystyle R_{3}$
$\displaystyle\leq\sum_{i=2}^{N}||\eta_{p}||_{\infty,P_{i-\frac{1}{2}}}\left(||\xi_{u}||_{\Gamma_{i-1}}+||\xi_{u}||_{\Gamma_{i}}\right)+\sum_{i=2}^{N}\left(||\eta_{u}||_{\Gamma_{i-1}}+||\eta_{u}||_{\Gamma_{i}}\right)||\xi_{p}||_{\infty,I_{i}}$
$\displaystyle\leq C\Delta
x^{k}\sum_{i=2}^{N}\left[||p||_{k+1,P_{i-\frac{1}{2}}}\left(||\xi_{u}||_{i-1}+||\xi_{u}||_{i}\right)+\left(||u||_{i-1}+||u||_{i}\right)||\xi_{p}||_{k+1,I_{i}}\right]$
$\displaystyle\leq C\Delta
x^{k}\left(||p||_{k+1}||\xi_{u}||+||u||_{k+1}||\xi_{p}||\right),$ (4.8)

where step 1 is straightforward, step 2 follows from Lemmas 4.1 and 4.2 and
the last step we applied the Cauchy-Schwarz inequality. Substitute (4.6)-(4.8)
into (4.5) to obtain

$\frac{1}{2}\frac{d}{dt}||\xi_{u}||^{2}+||\xi_{p}||^{2}\leq C\Delta
x^{2k}+C\Delta x^{k}(||\xi_{u}||+||\xi_{p}||),$

which further yields

$\frac{1}{2}\frac{d}{dt}||\xi_{u}||^{2}+||\xi_{p}||^{2}\leq C\Delta
x^{2k}+||\xi_{u}||^{2}.$

Finally, we can apply the Gronwall’s inequality and complete the proof. ∎

##  5 Numerical scheme for two-dimensional case

In this section, we will construct the scheme in two-dimensional case for the
following problem over the domain $\Omega=[0,1]\times[0,1]$

$\begin{cases}u_{t}=(a(u)p)_{x}+(b(u)q)_{x},\\\ p=A(u)_{x},\\\
q=B(u)_{y}.\end{cases}$ (5.1)

where $\displaystyle A(u)=\int^{u}a(t)\>dt$ and $\displaystyle
B(u)=\int^{u}b(t)\>dt.$ We consider the following boundary conditions

  1. Neumann boundary condition:

$u_{x}(0,y,t)=u_{x}(1,y,t)=u_{y}(x,0,t)=u_{y}(x,1,t)=0.$ (5.2)

  2. Dirichlet boundary condition:

$u(0,y,t)=u(1,y,t)=u(x,0,t)=u(x,1,t)=0.$ (5.3)

First, we give a rectangular decomposition of $\Omega$ which is the primitive
mesh for the primary variable $u$. Let
$0=x_{\frac{1}{2}}<x_{\frac{3}{2}}<\cdots<x_{N_{x}+\frac{1}{2}}=1$, and
$0=y_{\frac{1}{2}}<y_{\frac{3}{2}}<\cdots<y_{N_{y}+\frac{1}{2}}=1$ be grid
points in $x$ and $y$ directions, respectively, and denote the $i,j-th$ cell
as

$I_{ij}=I_{i}\times
J_{j}=[x_{i-\frac{1}{2}},x_{i+\frac{1}{2}}]\times[y_{j-\frac{1}{2}},y_{j+\frac{1}{2}}]$

for all $i=1,\ldots,N_{x}$ and $j=1,\ldots,N_{y}$. Also, we denote

$\Delta
x_{i}=x_{i+\frac{1}{2}}-x_{i-\frac{1}{2}},\>\>\>\>x_{i}=\frac{x_{i+\frac{1}{2}}+x_{i-\frac{1}{2}}}{2},\>\>\>\>\Delta
y_{i}=y_{j+\frac{1}{2}}-y_{j-\frac{1}{2}},\>\>\>\>y_{j}=\frac{y_{j+\frac{1}{2}}+y_{j-\frac{1}{2}}}{2},$

and

$\Delta x=\max_{i}\Delta x_{i},\>\>\>\>\Delta y=\max_{j}\Delta
y_{j},\>\>\>\>h=\max\\{\Delta x,\Delta y\\}.$

Next, we define the P-mesh and Q-mesh for solving the auxiliary variable $p$
and $q$. we choose $\tilde{x}_{i}$ given as

$\tilde{x}_{i}=x_{i}+\frac{\Delta x_{i}}{2}\xi_{0},\>\>\>\>\xi_{0}\in[-1,1]$
(5.4)

to define the P-mesh as

$P_{i-\frac{1}{2},j}=[\tilde{x}_{i-1},\tilde{x}_{i}]\times[y_{j-\frac{1}{2}},y_{j+\frac{1}{2}}],\>\>\>i=1,...,N_{x},$

with $\tilde{x}_{0}=\tilde{x}_{N_{x}-1}$. Similarly, we pick a point
$\tilde{y}_{j}$ given as

$\tilde{y}_{j}=y_{j}+\frac{\Delta
y_{j}}{2}\eta_{0},\>\>\>\>\eta_{0}\in[-1,1],$ (5.5)

and define Q-mesh as

$Q_{i,j-\frac{1}{2}}=[x_{i-\frac{1}{2}},x_{i+\frac{1}{2}}]\times[\tilde{y}_{j-1},\tilde{y}_{j}],\>\>\>j=1,...,N_{y},$

with $\tilde{y}_{0}=\tilde{y}_{N_{y}-1}$. Similar to the problem in one space
dimension, we need to deal with the boundary cells for problems with the
Neumann and Dirichlet boundary conditions. We can leave the left and right
boundary cells as $P_{\frac{1}{2},j}=[0,\tilde{x}_{1}]\times J_{j}$,
$P_{N_{x}+\frac{1}{2},j}=[\tilde{x}_{N_{x}},1]\times J_{j}$,
$Q_{i,\frac{1}{2}}=I_{i}\times[0,\tilde{y}_{1}]$ and
$Q_{i,N_{y}+\frac{1}{2}}=I_{i}\times[\tilde{y}_{N_{y}},1]$ . This mesh is
called the L-mesh. The other way is to combine the boundary cells with their
neighbour as $P_{\frac{3}{2},j}=[0,\tilde{x}_{2}]\times
J_{j},P_{N_{x}-\frac{1}{2}}=[\tilde{x}_{N_{x}-1},1]\times
J_{j},Q_{i,\frac{3}{2}}=I_{i}\times[0,\tilde{x}_{2}]$ and
$Q_{N_{y}-\frac{1}{2}}=I_{i}\times[\tilde{y}_{N_{y}-1},1]$. This mesh is
called C-mesh.

We define the finite element spaces for the L-mesh, $P^{L}_{h}$ and
$Q^{L}_{h}$, and for the C-mesh, $P^{C}_{h}$ and $Q^{C}_{h}$, to be

$\displaystyle V_{h}$ $\displaystyle=\\{u_{h}:u_{h}|_{I_{i}j}\in
P^{k}(I_{i}j),i=1,\ldots,N_{x},j=1,\ldots,N_{y}\\},$ $\displaystyle P^{L}_{h}$
$\displaystyle=\\{p_{h}:p_{h}|_{P_{i-\frac{1}{2}},j}\in
P^{k}(P_{i-\frac{1}{2},j}),i=1,\ldots,N_{x}+1,j=1,\ldots,N_{y}\\},$
$\displaystyle P^{C}_{h}$
$\displaystyle=\\{p_{h}:p_{h}|_{P_{i-\frac{1}{2}},j}\in
P^{k}(P_{i-\frac{1}{2},j}),i=2,\ldots,N_{x},j=1,\ldots,N_{y}\\},$
$\displaystyle Q^{L}_{h}$
$\displaystyle=\\{p_{h}:p_{h}|_{Q_{i,j-\frac{1}{2}}}\in
P^{k}(Q_{i,j-\frac{1}{2}}),i=1,\ldots,N_{x},j=2,\ldots,N_{y}+1\\},$
$\displaystyle Q^{C}_{h}$
$\displaystyle=\\{p_{h}:p_{h}|_{Q_{i,j-\frac{1}{2}}}\in
P^{k}(Q_{i,j-\frac{1}{2}}),i=1,\ldots,N_{x},j=2,\ldots,N_{y}\\},$

where $P^{k}$ is the space of polynomials of degree up to $k$.

Now, we can introduce the LDG method on overlapping mesh for (5.1).
Multiplying the test functions and using the integration by parts, we obtain

$\displaystyle\int_{I_{ij}}(u_{h})_{t}vdxdy$
$\displaystyle=-\int_{I_{i}}a(u_{h})p_{h}v_{x}dxdy+\int_{J_{j}}\hat{a}_{i+\frac{1}{2},j}\widehat{p_{h}}_{i+\frac{1}{2},j}v^{-}_{i+\frac{1}{2},j}dy-\int_{J_{j}}\hat{a}_{i-\frac{1}{2},j}\widehat{p_{h}}_{i-\frac{1}{2},j}v^{+}_{i-\frac{1}{2},j}dy,$
$\displaystyle\>\>\>\>\>-\int_{J_{i}}b(u_{h})p_{h}v_{y}dxdy+\int_{I_{i}}\hat{a}_{i,j+\frac{1}{2}}\widehat{p_{h}}_{i,j+\frac{1}{2}}v^{-}_{i,j+\frac{1}{2}}dx-\int_{I_{i}}\hat{a}_{i,j-\frac{1}{2}}\widehat{p_{h}}_{i,j-\frac{1}{2}}v^{+}_{i,j-\frac{1}{2}}dx,$
(5.6) $\displaystyle\int_{P_{i-\frac{1}{2},j}}p_{h}wdxdy$
$\displaystyle=-\int_{P_{i-\frac{1}{2},j}}A(u_{h})w_{x}dxdy+\int_{J_{j}}A(u_{h}(\tilde{x}_{i}))w^{-}_{i}dy-\int_{J_{j}}A(u_{h}(\tilde{x}_{i-1}))w^{+}_{i-1}dy,$
(5.7) $\displaystyle\int_{Q_{i,j-\frac{1}{2}}}q_{h}zdxdy$
$\displaystyle=-\int_{Q_{i,j-\frac{1}{2}}}B(u_{h})z_{y}dxdy+\int_{I_{i}}B(u_{h}(\tilde{y}_{j}))z^{-}_{j}dx-\int_{I_{i}}B(u_{h}(\tilde{y}_{j-1}))z^{+}_{j-1}dx.$
(5.8)

Then the LDG method on overlapping meshes for (5.1) is defined as follows:

  * •

L-mesh: find $(u_{h},p_{h},q_{h})\in V_{h}\times P^{L}_{h}\times Q^{L}_{h}$,
such that for any test functions $(v,w,z)\in V_{h}\times P^{L}_{h}\times
Q^{L}_{h}$ we have (5.6) - (5.8),

  * •

C-mesh: find $(u_{h},p_{h},q_{h})\in V_{h}\times P^{C}_{h}\times Q^{C}_{h}$,
such that for any test functions $(v,w,z)\in V_{h}\times P^{C}_{h}\times
Q^{C}_{h}$ we have (5.6) - (5.8).

We denoted $\displaystyle
u^{+}_{i-\frac{1}{2},j},u^{-}_{i+\frac{1}{2},j},u^{+}_{i,j-\frac{1}{2}},u^{+}_{i,j+\frac{1}{2}}$
as the traces of $u\in V_{h}$ on the four edges of $I_{ij}$, respectively.
Likewise for the traces of along the vertical edges of $P_{i-\frac{1}{2},j}$
and the horizontal edges of $Q_{i,j-\frac{1}{2}}$. Moreover, we use
$[u]=u^{+}-u^{-}$ and $\displaystyle\\{u\\}=\frac{1}{2}(u^{+}+u^{-})$ as the
jump and average of $u$ at the cells interfaces, receptively.

The numerical flux $\hat{a}$ along the edge $x_{i-\frac{1}{2}}$ is defined as

$\hat{a}_{i-\frac{1}{2},j}=\frac{[A(u_{h})]_{i-\frac{1}{2},j}}{[u_{h}]_{i-\frac{1}{2}},j}.$
(5.9)

Similarly, the numerical flux $\hat{b}$ along the edge $y_{j-\frac{1}{2}}$ is
defined as

$\hat{b}_{j-\frac{1}{2}}=\frac{[B(u_{h})]_{i,j-\frac{1}{2}}}{[u_{h}]_{i,j-\frac{1}{2}}}.$
(5.10)

Also, we choose the numerical flux

$\widehat{p_{h}}_{i-\frac{1}{2},j}=p_{h}(x_{i-\frac{1}{2}},y)+\frac{\alpha_{i-\frac{1}{2},j}}{\Delta\tilde{x}_{i-\frac{1}{2},j}}[u_{h}]_{i-\frac{1}{2},j},$
(5.11)

and

$\widehat{q_{h}}_{i,j-\frac{1}{2}}=q_{h}(x,y_{j-\frac{1}{2}})+\frac{\alpha_{i,j-\frac{1}{2}}}{\Delta\tilde{y}_{i,j-\frac{1}{2},j}}[u_{h}]_{i,j-\frac{1}{2}},$
(5.12)

where $[s]_{i-\frac{1}{2},j}=s^{+}_{i-\frac{1}{2},j}-s^{-}_{i-\frac{1}{2},j}$
stands for the jump of the function $s$ across the cell boundary
$\\{x_{i-\frac{1}{2}}\\}\times J_{j}$. Similarly for $[s]_{i,j-\frac{1}{2}}$.

To obtain the stability analysis and error estimates, we can follow the same
analyses for the problem in one-dimensional space. Therefore, we will omit the
proof and only state the results in the following two.

######  Theorem 5.1.

The LDG method introduced (5.6), (5.7) and (5.8) with the boundary conditions
(5.2) and (5.3) are stable and

$\frac{1}{2}\frac{d}{dt}||u_{h}||^{2}+||p_{h}||^{2}+||q_{h}||^{2}\leq 0.$

######  Theorem 5.2.

Suppose the exact solution for linear parabolic equation (5.1) with
$a(u)=b(u)=1$ satisfies $u\in C^{k+1}(\Omega)$ and the finite element space is
made up of piecewise polynomial of degree k in each cell. The numerical
solutions satisfy (5.6) - (5.8). Then the error between the numerical and
exact solutions satisfies

$||u-u_{h}||+\int_{0}^{T}||p-p_{h}||+||q-q_{h}||dt\leq Ch^{k},$

where $C$ is independent of $h$.

##  6 Numerical Experiments

In this section, several numerical experiments will be given to demonstrate
the stability and accuracy of the new LDG method on overlapping mesh for non-
periodic boundary conditions.

######  Example 6.1.

To compare the CFL number for two different meshes, we consider third-order
SSP Runge-Kutta time discretization [15] to solve the following heat equation
in one space dimension

$u_{t}=u_{xx},\>\>\>\>x\in[0,2\pi],$ (6.1)

subject to Neumann and Dirichlet boundary conditions

$\displaystyle u_{x}(0,t)$ $\displaystyle=u_{x}(2\pi,t)=0,$ $\displaystyle
u(0,t)$ $\displaystyle=u(2\pi,t)=0,$

respectively.

We consider a uniform mesh and use linear polynomial $P^{1}$ with the dual
mesh is generated by using the midpoint of the primitive mesh to compare the
CFL number to the two different meshes. The results in Table 1 shows that the
CFL number for C-mesh is larger than that for L-mesh. This is mainly because
the small cell effect does not works for C-mesh as we combine the small cells
near the boundary with its neighbor.

CFL | N | L-mesh | C-mesh
---|---|---|---
Neumann | Dirichlet | Neumann | Dirichlet
$L^{2}$ norm | order |  $L^{2}$ norm | order |  $L^{2}$ norm | order |  $L^{2}$ norm | order
| 10 | 1.65E-01 | - | 9.22E+00 | - | 4.68E-02 | - | 6.47E-02 | -
| 20 | 1.37E+01 | - | 1.45E+07 | - | 1.83E-02 | 1.36 | 2.36E-02 | 1.46
0.25 | 40 | 2.06E+12 | - | 2.27E+33 | - | 8.48E-03 | 1.11 | 9.70E-03 | 1.28
| 80 | 1.17E+59 | - | 1.31E+143 | - | 4.23E-03 | 1.00 | 4.52E-03 | 1.10
| 160 | INF | - | Nan | - | 2.12E-03 | 0.99 | 2.19E-03 | 1.05
Table 1: Example 6.1: The results of CFL testing.

For all the following numerical experiments, we use piecewise polynomials of
degree $k=1,2.$ Moreover, we consider third-order SSP Runge-Kutta time
discretization [15] with $\Delta t=0.01\Delta x^{2}$ to reduce the time error.

######  Example 6.2.

We solve the following heat equation in one space dimension

$\begin{cases}u_{t}=u_{xx},\>\>\>\>x\in[0,2\pi],\\\
u(x,0)=\cos(x),\end{cases}$ (6.2)

with Neumann boundary condition

$u_{x}(0,t)=u_{x}(2\pi,t)=0$ (6.3)

Clearly, the exact solution is

$u(x,t)=e^{-t}\cos(x).$

We consider a uniform mesh and take $\xi_{0}=0$ in (2.1), i.e. the dual mesh
is generated by using the midpoint of the primitive mesh. We take the final
time $T=0.5$ and compute the $L^{2}$-norm of the error between the numerical
and exact solutions. In Table 2, the results show that we can only obtain
suboptimal accuracy if $k$ is an odd number with the penalty parameter
$\alpha=0$. One way to obtain the optimal accuracy is to choose $\alpha\neq
0$.

k | N | L-mesh | C-mesh
---|---|---|---
no penalty | $\alpha=1.0$ | no penalty | $\alpha=1.0$
$L^{2}$ norm | order |  $L^{2}$ norm | order |  $L^{2}$ norm | order |  $L^{2}$ norm | order
| 10 | 9.51E-02 | - | 2.12E-02 | - | 3.91E-02 | - | 2.78E-02 | -
| 20 | 4.66E-02 | 1.03 | 4.61E-03 | 2.20 | 1.8E-02 | 1.08 | 6.98E-03 | 1.99
1 | 40 | 2.30E-02 | 1.02 | 1.08E-03 | 2.09 | 8.66E-03 | 1.10 | 1.64E-03 | 2.10
| 80 | 1.14E-02 | 1.01 | 2.63E-04 | 2.04 | 4.24E-03 | 1.03 | 3.91E-04 | 2.06
| 160 | 5.67E-03 | 1.00 | 6.49E-05 | 2.02 | 2.13E-03 | 0.99 | 9.52E-05 | 2.04
| 10 | 1.29E-03 | - | 9.37E-04 | - | 2.26E-03 | - | 1.87E-03 | -
| 20 | 1.60E-04 | 3.01 | 1.14E-04 | 3.04 | 3.56E-04 | 2.67 | 1.76E-04 | 3.41
2 | 40 | 1.99E-05 | 3.00 | 1.41E-05 | 3.01 | 5.63E-05 | 2.66 | 2.00E-05 | 3.14
| 80 | 2.49E-06 | 3.00 | 1.76E-06 | 3.00 | 9.26E-06 | 2.60 | 2.46E-06 | 3.03
| 160 | 3.12E-07 | 3.00 | 2.20E-07 | 3.00 | 1.57E-06 | 2.56 | 3.06E-07 | 3.00
Table 2: Example 6.2: midpoint with Neumann boundary condition.

In addition, we also take $\alpha=0$ and $\xi_{0}=0.1$ which is closed to 0
and $\xi_{0}=\sqrt{3}/3$ which is away from 0. The results in Table 3
demonstrated that we may recover the optimal convergence rates if we take
$\xi_{0}\neq 0$. However, we may observe only the $\displaystyle
k+\frac{1}{2}$ order of accuracy when we consider the C-mesh treatment without
the penalty term.

k | N | L-mesh | C-mesh
---|---|---|---
$\xi_{0}=0.1$ | $\xi_{0}=\sqrt{3}/3$ | $\xi_{0}=0.1$ | $\xi_{0}=\sqrt{3}/3$
$L^{2}$ norm | order |  $L^{2}$ norm | order |  $L^{2}$ norm | order |  $L^{2}$ norm | order
| 10 | 4.18E-02 | - | 1.87E-02 | - | 4.09E-02 | - | 3.08E-02 | -
| 20 | 9.24E-03 | 2.18 | 4.05E-03 | 2.09 | 1.71E-02 | 1.26 | 6.99E-03 | 2.14
1 | 40 | 2.25E-03 | 2.04 | 1.05E-03 | 2.07 | 7.31E-03 | 1.22 | 1.62E-03 | 2.11
| 80 | 5.65E-04 | 1.99 | 2.55E-04 | 2.04 | 2.21E-03 | 1.72 | 3.87E-04 | 2.06
| 160 | 1.42E-04 | 1.99 | 6.28E-05 | 2.02 | 5.66E-04 | 1.98 | 9.45E-05 | 2.03
| 10 | 1.51E-03 | - | 1.29E-03 | - | 2.29E-03 | - | 2.76E-03 | -
| 20 | 1.79E-04 | 3.07 | 1.55E-04 | 3.05 | 6.59E-04 | 2.68 | 4.60E-04 | 2.59
2 | 40 | 2.22E-05 | 3.02 | 1.93E-05 | 3.01 | 5.65E-05 | 2.67 | 7.73E-05 | 2.57
| 80 | 2.76E-06 | 3.00 | 2.41E-06 | 3.00 | 9.26E-06 | 2.61 | 1.33E-06 | 2.54
| 160 | 3.45E-07 | 3.00 | 3.01E-07 | 3.00 | 1.57E-06 | 2.56 | 2.30E-06 | 2.52
Table 3: Example 6.2: $\xi_{0}=0.1$ and $\xi_{0}=\sqrt{3}/3$ with Neumann
boundary condition.

######  Example 6.3.

We solve

$\begin{cases}u_{t}=u_{xx},\>\>\>\>x\in[0,2\pi],\\\
u(x,0)=\sin(x),\end{cases}$ (6.4)

subject to Dirichlet boundary condition

$u(0,t)=u(2\pi,t)=0.$ (6.5)

Clearly, the exact solution is

$u(x,t)=e^{-t}\sin(x).$

We also consider a uniform mesh and take $\xi_{0}=0$ in (2.1). Then, we choose
$\alpha=0$, and take $\xi_{0}=0.1$ which is closed to 0 and
$\xi_{0}=\sqrt{3}/3$ which is away from 0. According to Table 4 and Table 5,
we can observe the same outcomes discussed for the Neumann boundary condition
in Example 6.2. However, we can obtain the optimal accuracy when we consider
this case without the penalty term.

k | N | L-mesh | C-mesh
---|---|---|---
no penalty | $\alpha=1.0$ | no penalty | $\alpha=1.0$
$L^{2}$ norm | order |  $L^{2}$ norm | order |  $L^{2}$ norm | order |  $L^{2}$ norm | order
| 10 | 7.19E-02 | - | 1.82E-02 | - | 5.04E-02 | - | 3.16E-02 | -
| 20 | 3.54E-02 | 1.02 | 4.26E-03 | 2.10 | 2.21E-02 | 1.29 | 6.45E-03 | 2.30
1 | 40 | 1.76E-02 | 1.00 | 1.04E-03 | 2.04 | 9.51E-03 | 1.22 | 1.50E-03 | 2.10
| 80 | 8.81E-03 | 1.00 | 2.57E-04 | 2.01 | 4.50E-03 | 1.08 | 3.70E-04 | 2.02
| 160 | 4.40E-03 | 1.00 | 6.42E-05 | 2.00 | 2.18E-03 | 1.04 | 9.23E-05 | 2.00
| 10 | 1.32E-03 | - | 9.75E-04 | - | 1.96E-03 | - | 1.59E-03 | -
| 20 | 1.63E-04 | 3.01 | 1.16E-04 | 3.08 | 2.41E-04 | 3.03 | 2.03E-04 | 2.97
2 | 40 | 2.02E-05 | 3.01 | 1.42E-05 | 3.03 | 2.99E-05 | 3.01 | 2.30E-05 | 3.14
| 80 | 2.51E-06 | 3.01 | 1.76E-06 | 3.01 | 3.73E-06 | 3.00 | 2.68E-06 | 3.10
| 160 | 3.13E-07 | 3.00 | 2.20E-07 | 3.00 | 4.66E-07 | 3.00 | 3.20E-07 | 3.06
Table 4: Example 6.3 midpoint with Dirichlet boundary condition. k | N | L-mesh | C-mesh
---|---|---|---
$\xi_{0}=0.1$ | $\xi_{0}=\sqrt{3}/3$ | $\xi_{0}=0.1$ | $\xi_{0}=\sqrt{3}/3$
$L^{2}$ norm | order |  $L^{2}$ norm | order |  $L^{2}$ norm | order |  $L^{2}$ norm | order
| 10 | 3.85E-02 | - | 1.58E-02 | - | 5.34E-02 | - | 3.91E-02 | -
| 20 | 9.56E-03 | 2.01 | 3.95E-03 | 2.00 | 2.21E-02 | 1.27 | 1.26E-02 | 1.63
1 | 40 | 2.33E-03 | 2.04 | 9.88E-04 | 2.00 | 8.97E-03 | 1.30 | 4.02E-03 | 1.65
| 80 | 5.77E-04 | 2.01 | 2.47E-04 | 2.00 | 2.73E-03 | 1.71 | 1.37E-03 | 1.55
| 160 | 1.44E-04 | 2.00 | 6.18E-05 | 2.00 | 7.52E-04 | 1.86 | 4.74E-04 | 1.53
| 10 | 1.61E-03 | - | 1.36E-03 | - | 1.99E-03 | - | 2.30E-03 | -
| 20 | 1.87E-04 | 3.11 | 1.61E-04 | 3.08 | 2.46E-04 | 3.02 | 2.61E-04 | 3.14
2 | 40 | 2.27E-05 | 3.04 | 1.96E-05 | 3.03 | 3.05E-05 | 3.01 | 3.07E-05 | 3.09
| 80 | 2.80E-06 | 3.02 | 2.43E-06 | 3.01 | 3.81E-06 | 3.00 | 3.71E-06 | 3.04
| 160 | 3.47E-07 | 3.01 | 3.02E-07 | 3.01 | 4.76E-07 | 3.00 | 4.57E-07 | 3.02
Table 5: Example 6.3 $\xi_{0}=0.1$ , and $\xi_{0}=\sqrt{3}/3$ with Dirichlet
boundary condition.

######  Example 6.4.

We solve the following heat equation in two space dimension

$\begin{cases}u_{t}=u_{xx}+u_{yy},\>\>\>\>(x,y)\in[0,2\pi]\times[0,2\pi],\\\
u(x,y,0)=\cos(x)\cos(y),\end{cases}$ (6.6)

subject to Neumann boundary condition

$u_{x}(0,y,t)=u_{x}(2\pi,y,t)=0\>\>\text{and}\>\>u_{y}(x,0,t)=u_{y}(x,2\pi,t)=0.$
(6.7)

Clearly, the exact solution is

$u(x,y,t)=e^{-2t}\cos(x)\cos(y).$

We also consider a uniform mesh, and take $\xi_{0}=0$ and $\eta_{0}=0$ to
examine the dual meshes which were generated by the midpoint in both
directions. We choose the final time to be $T=0.1$. Table 6 shows suboptimal
accuracy if $k$ is an odd number with the penalty parameter $\alpha=0$. Also,
one way to obtain the optimal accuracy is choosing $\alpha\neq 0$. However,
when the C-mesh was used we can only obtain the $\displaystyle k+\frac{1}{2}$
order of accuracy.

k | N | L-mesh | C-mesh
---|---|---|---
no penalty | $\alpha=1.0$ | no penalty | $\alpha=1.0$
$L^{2}$ norm | order |  $L^{2}$ norm | order |  $L^{2}$ norm | order |  $L^{2}$ norm | order
| 16 | 3.72E-01 | - | 1.34E-01 | - | 5.51E-01 | - | 5.52E-01 | -
1 | 64 | 1.36E-01 | 1.45 | 3.98E-02 | 1.78 | 1.92E-01 | 1.52 | 1.96E-01 | 1.50
| 256 | 5.51E-02 | 1.30 | 1.12E-02 | 1.83 | 6.99E-02 | 1.45 | 6.14E-02 | 1.67
| 1024 | 2.54E-02 | 1.11 | 2.81E-03 | 2.00 | 2.17E-02 | 1.69 | 2.00E-02 | 1.62
| 16 | 2.06E-01 | - | 2.26E-01 | - | 7.85E-01 | - | 6.23E-01 | -
| 64 | 3.40E-02 | 2.59 | 4.51E-02 | 2.34 | 1.44E-01 | 2.39 | 1.12E-01 | 2.47
2 | 256 | 5.20E-03 | 2.71 | 7.91E-03 | 2.51 | 2.59E-02 | 2.48 | 2.02E-02 | 2.48
| 1024 | 7.51E-04 | 2.79 | 1.14E-03 | 2.79 | 4.63E-03 | 2.48 | 3.65E-03 | 2.47
| 4096 | 9.52E-05 | 2.97 | 1.58E-04 | 2.85 | 8.21E-04 | 2.49 | 6.54E-04 | 2.48
| 16384 | 1.19E-05 | 2.99 | 1.98E-05 | 2.99 | 1.47E-04 | 2.49 | 1.18E-04 | 2.48
Table 6: Example 6.4: midpoint with Neumann boundary condition.

Moreover, we take $\alpha=0$, and consider the combination of $\xi_{0}=0,0.5$
and $\eta_{0}=0,0.5$ The results in Table 7 demonstrated that we may recover
optimal convergence rates we if take $\xi_{0}\neq 0$ and $\eta_{0}\neq 0$.
However, when we use the C-mesh, we can only obtain $\displaystyle
k+\frac{1}{2}$ order of accuracy as shown in Table 8.

k | N | L-mesh
---|---|---
$\xi_{0}=0$, $\eta_{0}=0.5$ |  $\xi_{0}=0.5$, $\eta_{0}=0$ |  $\xi_{0}=0.5$, $\eta_{0}=0.5$
$L^{2}$ norm | order |  $L^{2}$ norm | order |  $L^{2}$ norm | order
| 16 | 3.48E-01 | - | 3.47E-01 | - | 3.22E-01 | -
1 | 64 | 1.15E-01 | 1.59 | 1.16E-01 | 1.57 | 9.04E-02 | 1.83
| 256 | 4.22E-02 | 1.44 | 4.27E-02 | 1.45 | 2.28E-02 | 1.99
| 1024 | 1.84E-02 | 1.19 | 1.85E-02 | 1.20 | 5.55E-03 | 2.04
| 16 | 1.76E-01 | - | 1.78E-01 | - | 1.69E-01 | -
| 64 | 3.76E-02 | 2.66 | 2.78E-02 | 2.68 | 2.92E-02 | 2.53
2 | 256 | 3.74E-03 | 2.88 | 4.06E-03 | 2.78 | 4.63E-03 | 2.66
| 1024 | 5.03E-04 | 2.90 | 5.78E-04 | 2.81 | 5.93E-04 | 2.96
| 4096 | 6.35E-05 | 2.98 | 7.35E-05 | 2.98 | 7.59E-05 | 2.97
| 16384 | 7.99E-06 | 2.99 | 9.29E-06 | 2.99 | 9.52E-06 | 2.99
Table 7: Example 6.4: combination of $\xi_{0}=0,0.5$ and $\eta_{0}=0,0.5$ with Neumann boundary condition on L-mesh. k | N | C-mesh
---|---|---
$\xi_{0}=0$, $\eta_{0}=0.5$ |  $\xi_{0}=0.5$, $\eta_{0}=0$ |  $\xi_{0}=0.5$, $\eta_{0}=0.5$
$L^{2}$ norm | order |  $L^{2}$ norm | order |  $L^{2}$ norm | order
| 16 | 5.54E-01 | - | 5.56E-01 | - | 5.57E-01 | -
1 | 64 | 1.94E-01 | 1.51 | 1.93E-01 | 1.53 | 1.84E-02 | 1.60
| 256 | 6.13E-02 | 1.66 | 5.70E-02 | 1.76 | 5.97E-02 | 1.62
| 1024 | 1.99E-02 | 1.62 | 1.83E-02 | 1.64 | 1.89E-03 | 1.66
| 16 | 5.48E-01 | - | 4.85E-01 | - | 3.53E-01 | -
| 64 | 8.88E-02 | 2.62 | 7.85E-02 | 2.63 | 5.42E-02 | 2.70
2 | 256 | 1.49E-02 | 2.57 | 1.25E-02 | 2.65 | 8.57E-02 | 2.66
| 1024 | 2.59E-03 | 2.53 | 2.06E-03 | 2.60 | 1.46E-03 | 2.56
| 4096 | 4.48E-04 | 2.53 | 3.45E-04 | 2.58 | 2.58E-04 | 2.50
| 16384 | 8.02E-05 | 2.48 | 6.10E-05 | 2.50 | 4.57E-05 | 2.50
Table 8: Example 6.4: combination of $\xi_{0}=0,0.5$ and $\eta_{0}=0,0.5$ with
Neumann boundary condition on C-mesh.

######  Example 6.5.

We solve

$\begin{cases}u_{t}=u_{xx}+u_{yy},\>\>\>\>(x,y)\in[0,2\pi]\times[0,2\pi],\\\
u(x,y,0)=\sin(x)\sin(y),\end{cases}$ (6.8)

subject to Dirichlet boundary condition

$u(0,y,t)=u(2\pi,y,t)=0\>\>\text{and}\>\>u(x,0,t)=u(x,2\pi,t)=0$ (6.9)

Clearly, the exact solution is

$u(x,y,t)=e^{-2t}\sin(x)\sin(y).$

We also explore the problem on a uniform mesh and take $\xi_{0}=0$ and
$\eta_{0}=0$ to examine the dual meshes that were generated by the midpoint in
both directions. We take $\alpha=0$ and consider the combination of
$\xi_{0}=0,0.5$ and $\eta_{0}=0,0.5$. As stated in Table 9, Table 10 and Table
11, we can observe the same results discussed for problems with Neumann
boundary condition in Example 6.4.

k | N | L-mesh | C-mesh
---|---|---|---
no penalty | $\alpha=1.0$ | no penalty | $\alpha=1.0$
$L^{2}$ norm | order |  $L^{2}$ norm | order |  $L^{2}$ norm | order |  $L^{2}$ norm | order
| 16 | 3.55E-01 | - | 8.83E-02 | - | 6.62E-01 | - | 6.63E-01 | -
1 | 64 | 1.08E-01 | 1.73 | 2.18E-02 | 2.01 | 18.2E-01 | 1.86 | 1.84E-01 | 1.85
| 256 | 4.33E-02 | 1.30 | 5.38E-03 | 2.01 | 4.93E-02 | 1.86 | 5.20E-02 | 1.82
| 1024 | 2.05E-02 | 1.08 | 1.33E-03 | 2.01 | 1.57E-02 | 1.64 | 1.45E-02 | 1.85
| 16 | 5.12E-02 | - | 1.69E-01 | - | 6.17E-01 | - | 5.29E-01 | -
| 64 | 9.81E-03 | 2.39 | 3.24E-02 | 2.39 | 1.27E-01 | 2.28 | 1.10E-01 | 2.27
2 | 256 | 1.58E-03 | 2.63 | 5.37E-03 | 2.59 | 2.22E-02 | 2.51 | 1.97E-02 | 2.48
| 1024 | 2.25E-04 | 2.81 | 8.54E-04 | 2.65 | 3.82E-03 | 2.54 | 3.45E-03 | 2.51
| 4096 | 3.00E-05 | 2.90 | 1.07E-04 | 2.99 | 6.65E-04 | 2.52 | 6.07E-04 | 2.51
| 16384 | 3.75E-06 | 3.00 | 1.34E-05 | 2.99 | 1.17E-04 | 2.51 | 1.07E-04 | 2.50
Table 9: Example 6.5: midpoint with Dirichlet boundary condition. k | N | L-mesh
---|---|---
$\xi_{0}=0$, $\eta_{0}=0.5$ |  $\xi_{0}=0.5$, $\eta_{0}=0$ |  $\xi_{0}=0.5$, $\eta_{0}=0.5$
$L^{2}$ norm | order |  $L^{2}$ norm | order |  $L^{2}$ norm | order
| 16 | 3.51E-01 | - | 3.54E-01 | - | 3.47E-01 | -
1 | 64 | 9.62E-01 | 1.86 | 1.07E-01 | 1.72 | 8.04E-02 | 2.05
| 256 | 3.40E-02 | 1.50 | 4.27E-02 | 1.33 | 2.09E-02 | 2.01
| 1024 | 1.49E-02 | 1.18 | 1.69E-02 | 1.33 | 5.23E-03 | 2.00
| 16 | 5.20E-02 | - | 5.20E-02 | - | 5.30E-02 | -
| 64 | 9.66E-03 | 2.42 | 9.66E-03 | 2.43 | 9.54E-03 | 2.47
2 | 256 | 1.54E-03 | 2.64 | 1.54E-03 | 2.65 | 1.50E-03 | 2.67
| 1024 | 2.12E-04 | 2.85 | 2.12E-04 | 2.86 | 1.99E-04 | 2.91
| 4096 | 2.79E-05 | 2.93 | 2.79E-05 | 2.93 | 2.56E-05 | 2.96
| 16384 | 3.50E-06 | 2.99 | 3.49E-06 | 3.00 | 3.20E-06 | 3.00
Table 10: Example 6.5: combination of $\xi_{0}=0,0.5$ and $\eta_{0}=0,0.5$ with Dirichlet boundary condition on L-mesh. k | N | C-mesh
---|---|---
$\xi_{0}=0$, $\eta_{0}=0.5$ |  $\xi_{0}=0.5$, $\eta_{0}=0$ |  $\xi_{0}=0.5$, $\eta_{0}=0.5$
$L^{2}$ norm | order |  $L^{2}$ norm | order |  $L^{2}$ norm | order
| 16 | 6.62E-01 | - | 6.62E-01 | - | 6.62E-01 | -
1 | 64 | 1.81E-01 | 1.86 | 1.81E-01 | 1.87 | 1.80E-02 | 1.88
| 256 | 4.91E-02 | 1.88 | 4.81E-02 | 1.91 | 4.77E-02 | 1.92
| 1024 | 1.59E-02 | 1.62 | 1.32E-02 | 1.87 | 1.28E-03 | 1.91
| 16 | 6.48E-01 | - | 5.81E-01 | - | 6.14E-01 | -
| 64 | 1.31E-01 | 2.30 | 1.12E-01 | 2.37 | 1.17E-01 | 2.39
2 | 256 | 2.34E-02 | 2.50 | 1.85E-02 | 2.60 | 1.96E-02 | 2.58
| 1024 | 4.07E-03 | 2.53 | 3.03E-03 | 2.60 | 3.24E-03 | 2.60
| 4096 | 7.10E-04 | 2.52 | 5.12E-04 | 2.57 | 5.48E-04 | 2.56
| 16384 | 1.25E-05 | 2.51 | 8.83E-05 | 2.54 | 9.45E-05 | 2.54
Table 11: Example 6.5: combination of $\xi_{0}=0,0.5$ and $\eta_{0}=0,0.5$
with Dirichlet boundary condition on C-mesh.

##  7 Conclusion

In this paper, we demonstrate the algorithms of the LDG method on overlapping
mesh for Neumann and Dirichlet boundary conditions. The scheme for each
boundary condition is stable with adjusted boundary cells but the order of
accuracy may not be optimal. We observed that C-mesh could yield larger time
step than L-mesh. However, in some cases, C-mesh cannot yield optimal
convergence rates. With a positive penalty parameter, both meshes can result
in optimal convergence rates.

Acknowledgements This work is supported by NSF grant DMS-1818467 KMITL
Research Fund, Research Seed Grant for New Lecturer.

## References

  * [1] F. Bassi, 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), 267-279.
  * [2] N. Chuenjarern, Z. Xu and Y. Yang, High-order bound-preserving discontinuous Galerkin methods for compressible miscible displacements in porous media on triangular meshes, Journal of Computational Physics, 378 (2019), 110-128.
  * [3] N. Chuenjarern, Y. Yang, Fourier analysis of local discontinuous Galerkin method for linear parabolic equations on overlapping meshes, Journal of Scitific Computing, 81 (2019), 671-688.
  * [4] E. Chung and C.S. Lee, A staggered discontinuous Galerkin method for convection-diffusion equations, Journal of Numerical Mathematics, 20 (2012), 1-31.
  * [5] P.G. Ciarlet, Finite element method for elliptic problems, North-Holland, Amsterda, 1978.
  * [6] B. Cockburn, S. Hou, C. W. Shu, The Runge-Kutta local projection discontinuous Galerkin finite element method for conservation laws. IV: The multidimensional case, Mathematics Computation, 54 (1990), 545-581.
  * [7] B. Cockburn, S. Y. Lin, C. W. Shu, TVB Runge-Kutta local projection discontinuous Galerkin finite element method for conservation laws. III: One-dimensional systems, Journal of Computational Physics, 84 (1989), 90-113.
  * [8] B. Cockburn, C. W. Shu, TVB Runge-Kutta local projection discontinuous Galerkin finite element method for conservation laws. II: General framework, Mathematics Computation, 52 (1989), 411-435.
  * [9] B. Cockburn, C. W. Shu, The Runge-Kutta discontinuous Galerkin method for conservation laws. V: Multidimensional systems, Journal of Computational Physics, 141 (1998), 199-224.
  * [10] B. Cockburn, C.-W. Shu, The local discontinuous Galerkin method for time-dependent convection-diffusion systems, SIAM Journal on Numerical Analysis, 35 (1998), 2440-2463.
  * [11] J. Douglas, Jr., R.E. Ewing and M.F. Wheeler, _A time-discretization procedure for a mixed finite element approximation of miscible displacement in porous media_ , R.A.I.R.O. Analyse numérique, 17 (1983), 249-256.
  * [12] J. Douglas, Jr., R.E. Ewing and M.F. Wheeler, _The approximation of the pressure by a mixed method in the simulation of miscible displacement_ , R.A.I.R.O. Analyse numérique, 17 (1983), 17-33.
  * [13] J. Du and Y. Yang, Maximum-principle-preserving third-order local discontinuous Galerkin methods on overlapping meshes, Journal of computational physics, 377 (2019), 117-141.
  * [14] J. Du, Y. Yang and E. Chung, Stability analysis and error estimates of local discontinuous Galerkin method for convection-diffusion equations on overlapping meshes, BIT Numerical Mathematics, 59 (2019), 853-876.
  * [15] S. Gottlieb, C.-W. Shu and E. Tadmor, Strong stability-preserving high-order time discretization method, SIAM Review, 43(2001), 89-112.
  * [16] H. Guo, F. Yu and Y. Yang, Local discontinuous Galerkin method for incompressible miscible displacement problem in porous media, Journal of Scientific Computing, 71 (2017), 615-633.
  * [17] H. Guo and Y. Yang, Bound-preserving discontinuous Galerkin method for compressible miscible displacement problem in porous media, SIAM Journal on Scientific Computing, 39 (2017), A1969-A1990.
  * [18] H. Guo, Q. Zhang and Y. Yang, A combined mixed finite element method and local discontinuous Galerkin method for miscible displacement problem in porous media, Science China Mathematics, 57 (2014), 2301-2320.
  * [19] E. F. Keller and L. A. Segal, Initiation on slime mold aggregation viewed as instability, Journal of Theoretical of Biology, 26 (1970), 399-415
  * [20] X. Li, C.-W. Shu and Y. Yang, Local discontinuous Galerkin method for the Keller-Segel chemotaxis model, Journal of Scientific Computing, 73 (2017), 943-967.
  * [21] Y. Liu, C.-W. Shu, E. Tadmor and M. Zhang, Central local discontinuous Galerkin method on overlapping cells for diffusion equations, ESAIM: Mathematical Modeling and Numerical Analysis (M2AN), 45 (2011), 1009-1032.
  * [22] C. Patlak, Random walk with persistence and external bias, The bulletin of mathematical biophysics, 15 (1953), 311-338.
  * [23] W.H. Reed and T. R. Hill, Triangular mesh method for the neutron transport equation, Los Alamos Scientific Laboratory Report LA-UR-73-479, Los Alamos, NM, 1973.
  * [24] H. Wang, C.-W. Shu, and Q. Zhang, Stability and error estimates of local discontinuous Galerkin methods with implicit-explicit time-marching for advection-diffusion problems, SIAM Journal on Numerical Analysis, 53 (2015), 206-227.
  * [25] H. Wang, C.-W. Shu and Q. Zhang, Stability analysis and error estimates of local discontinuous Galerkin methods with implicit-explicit time-marching for nonlinear convection-diffusion problems, Applied Mathematics and Computation, 272 (2016), 237-258.
  * [26] H. Wang, S. Wang, Q. Zhang and C.-W. Shu, Local discontinuous Galerkin methods with implicit-explicit time marching for multi-dimensional convectiondiffusion problems, ESAIM:M2AN, 50 (2016), 1083-1105.
  * [27] T. Xiong, J.-M. Qiu and Z. Xu, High order maximum-principle-preserving discontinuous Galerkin method for convection-diffusion equations, SIAM Journal on Scientific Computing, 37 (2015), A583-A608.
  * [28] F. Yu, H. Guo, N. Chuenjarern and Y. Yang, Conservative local discontinuous Galerkin method for compressible miscible displacements in porous media, Journal of Scientific Computing, 73 (2017), 1249-1275.
  * [29] M. Zhang and C.-W. Shu, An analysis of three different formulations of the discontinuous Galerkin method for diffusion equations, Mathematical Models and Methods in Applied Sciences, 13 (2003), 395-413.