\documentclass[reqno]{amsart}
\usepackage{hyperref}
\usepackage{graphicx}
\AtBeginDocument{{\noindent\small
2021 UNC Greensboro PDE Conference,\newline
\emph{Electronic Journal of Differential Equations},
Conference 26 (2022), pp. 97--113.\newline
ISSN: 1072-6691. URL: http://ejde.math.txstate.edu or http://ejde.math.unt.edu}
\thanks{\copyright 2022 This work is licensed under a CC BY 4.0 license.}
\vspace{8mm}}
\begin{document} \setcounter{page}{97}
\title[\hfilneg EJDE-2018/conf/26\hfil Semi-Lagrangian forward methods for some PDEs]
{Semi-Lagrangian forward methods for some time-dependent nonlinear partial \\ differential equations}
\author[D. Guo \hfil EJDE-2022/conf/26\hfilneg]
{Daniel X. Guo}
\address{Daniel Guo \newline
Department of Mathematics and Statistics,
University of North Carolina at Wilmington, NC 28403, USA}
\email{guod@uncw.edu}
\thanks{Published August 25, 2022.}
\subjclass[2020]{35Q35, 65M12, 76M20}
\keywords{Semi-Lagrangian forward methods; trajectory; Runge-Kutta method;
\hfill\break\indent time-dependent nonlinear partial differential equations}
\begin{abstract}
In this article, we study one-step Semi-Lagrangian forward method for
computing the numerical solutions of time-dependent nonlinear partial differential
equations with initial and boundary conditions in one space dimension.
Comparing with classic Semi-Lagrangian method, this method is more straight forward
to analyze and implement.
This method is based on Lagrangian trajectory from the departure points (regular nodes)
to the arrival points by Runge-Kutta methods. The arrival points are traced forward
from the departure points along the trajectory of the path. Most likely the arrival
points are not on the regular grid nodes. However, it is convenient to approximate the
high order derivative terms in spatial dimension on regular nodes.
The convergence and stability are studied for the explicit methods.
The numerical examples show that those methods work very efficient for the
time-dependent nonlinear partial differential equations.
\end{abstract}
\maketitle
\numberwithin{equation}{section}
\newtheorem{theorem}{Theorem}[section]
\newtheorem{lemma}[theorem]{Lemma}
\newtheorem{definition}[theorem]{Definition}
\newtheorem{remark}[theorem]{Remark}
\newtheorem{example}[theorem]{Example}
\allowdisplaybreaks
\section{Introduction}
\label{sec:introduction}
Differential equations in mathematical models are very important part,
and depending on how many variables are involved,
the formulations of ordinary differential equations (ODEs), and of
partial differential equations (PDEs).
For some equation, we can use the traditional techniques to obtain an explicit or
implicit solutions. However, most realistic mathematical models cannot be solved
in this way; instead, they must be dealt with by computational methods that
provide approximate solutions in some sense, \cite{BF, GEAR, GH}.
A lot of numerical methods have been investigated to produce the approximations of
the desired solutions. Semi-Lagrangian methods have been introduced at the beginning
of the eighties \cite{ROB} and since then have been extensively used in the numerical
simulation of models for weather forecast, oceanography, and for more general problems
in fluid dynamics \cite{BM, LKM, NSS, SCT}. The basic idea of the semi-Lagrangian method
is to construct the numerical solution over a set of grid points by advancing the
characteristics, or the trajectory path.
In the article \cite{DG4}, a Semi-Lagrangian backward method (traditional method) was
investigated for its stability and convergence. However there are two disadvantages.
One is to compute the departure points by an implicit method. In each iteration at
the approximated departure point, a same order interpolation method comparing to the
iteration method is needed on no-regular grids. The computation cost is increasing
according to the order. The another one is to approximate the terms with spatial
derivatives at the irregular grid points. The approximations are not based on the regular
nodes since the departure points most likely is not on the regular node.
Recently, more applications of semi-Lagrangian method are reported.
Typically they are for the Navier-Stokes equations \cite{BFR}, shallow-water equations
\cite{CAR, BB, MWL}, Advection-diffusion equations \cite{EKSV, BS}, and Burgers
equation \cite{BKK}. Some even tried to use semi-Lagrangian forward method \cite{CRS, MWL}.
The advantages of semi-Lagrangian forward method are the overall order of approximation,
less computations from interpolation, and standard discretization of terms involving
spatial derivatives.
In this article, we study the stability and convergence of semi-Lagrangian forward method.
We consider the time-dependent nonlinear differential equation
$$
u_t + u u_x = f(t,x,u, u_x,...), \quad t > 0, \; x \in \Omega\subset R^n,
$$
with initial-value and boundary conditions. Where $n \ge 1$ could be any positive integer.
On the Lagrangian trajectory path, those equations could be considered as ODEs.
To obtain a numerical approximation of the solution, the arrival and departure points
at $t_{n+1}$ and $t = t_n$ are needed. One could be specified on regular nodes,
but the other one needs to be computed. For example, if the arrival points are on the grid,
then the departure points must be calculated \cite{DG4}.
In this article, the departure points are fixed on regular grids, but the arrival points
must be computed and most likely they are not on the regular grid.
Also Runge-Kutta methods could be applied to different types of partial differential
equations \cite{RBJC, CFN, MS, JGV}. In the paper \cite{DG3}, the truncation error
of the first-order Semi-Lagrangian method was studied.
The complete numerical analysis including stability and convergence is presented in
this paper for the first-order and the second-order Semi-Lagrangian methods.
The numerical analysis and numerical example showed that it is possible to construct
overall order-two/three one-step difference method for the time-dependent partial
differential equations.
This article is organized as follows:
Section 2 shows the specified partial differential equations with given initial-value
and boundary conditions to be studied; In Section 3, the arrival point is computed;
An one-step difference method is constructed and analyzed in Section 4;
In Section 5, the Semi-Lagrangian methods are investigated;
Numerical examples are presented in Section 6;
The final section is the summary and the consideration of future work.
\section{Time-dependent partial differential equations}
In this article, we will focus on the partial differential equation
\[
\frac{du(t,x)}{dt} = f(t,x,u), \quad t > 0, \; x \in \Omega,
\]
with initial-value condition
\[
u(0,x) = u_0(x), \quad x \in \Omega,
\]
where $\Omega$ is a domain in one, two, or three dimensions.
The boundary condition is provided for $u(t,x)$ on $\partial \Omega$ if it is needed.
The total derivative reads
\[
\frac{d}{dt} = \frac{\partial }{\partial t}+u \cdot \frac{\partial }{\partial x},
\]
where $u$ is the velocity field and $x$ is the spatial variable.
In this article, the dimension one is considered for demonstration.
It is more complicated for higher dimensions and it will be reported later.
For examples, the Kortweg-de-Vries (KdV) equation
\[
u_t + u u_x + \delta^2 u_{xxx} = 0,
\]
with boundary and initial conditions. Where $\delta$ is a given constant.
Also the Navier-Stokes equation reads
\begin{gather*}
u_t + (u \cdot \nabla ) u = f + \nu \Delta u - \nabla p, \\
\operatorname{div} u = 0,
\end{gather*}
with boundary and the initial conditions. Even the reformulated Shallow-Water
equations \cite{DGJD1, DGJD2} have the similar form.
To focus on the proposed method, we will study the general time-dependent differential
equation
\begin{equation}
\begin{gathered}
\frac{d u}{dt} = \frac{\partial u}{\partial t}+u \frac{\partial u}{\partial x}
= g(t,x,u),\quad 0 < t < T,\; a \le x \le b \\
u(0,x) = u_0(x), \quad a \le x \le b
\end{gathered} \label{eq:tpde}
\end{equation}
with boundary conditions about $x$ and $g(t,x,u) = f(t, x) + {{L}} u$,
where ${L}$ is a linear derivative operator in spatial $x$, for example
${{L}} u=c u_{xx}$, or $c u_{xxx}$, where $c$ is given constant.
The one-step forward explicit methods for the problem
\eqref{eq:tpde} are proposed. Also we study the stability and convergence of these methods
for problem \eqref{eq:tpde}.
\begin{theorem} \label{thm1}
If $g(t, x, u)$ is a continuous function of $t$ and $x$, and satisfies the Lipschitz
condition in $u$ in the region $0\le t \le T$, $a \le x \le b$, $-\infty < u < \infty$,
and the first order derivative of $u_0(x)$ is continuous on $[a,b]$, then there
exists at least a differentiable solution $u(t,x)$ of the initial-value boundary problem
\eqref{eq:tpde}.
\end{theorem}
This result can be found in many books on partial differential equations,
such as \cite[page 282]{RT}.
\section{Computing the arrival points}
It is well-known that the exact solution of the problem \eqref{eq:tpde}
for most cases is hard to find. The approximations to $u$ will be computed
first at the given points, called grid points, in the domain $[0,T] \times [a,b]$.
At other points in the domain, the approximate solutions are calculated from the
interpolation of the approximate solutions on the grid points.
Now assume that the grid points are equally distributed throughout the domain
$[0,T] \times [a,b]$. Let $N$ and $M$ be two positive integer numbers.
The grid points are formed by calculating
\begin{gather*}
t_k = k\tau, \quad \text{for } k = 0, 1, 2, \dots, N; \\
x_i = a + ih, \quad \text{ for } i = 0, 1, 2, \dots, M;
\end{gather*}
with the step sizes $\tau = T/N$ and $h=(b-a)/M$.
Assume that all solutions $u$ on $(t_k, x^k_i)$ for $i = 0, 1, \dots, M$ are known,
and we need to find the solutions $u$ on $(t_{k+1}, x^{k+1}_i)$ for
$i = 0, 1, 2, \dots, M$. From the view of Lagrangian trajectory, the particles
on arrival points at $t = t_{k+1}$ shall come from some points, called departure
points at $t = t_k$ along the trajectory of path. To compute the arrival
point $x^{k+1}_A$ at $t = t_{k+1}$ for the departure point $x_D^{k}$,
we consider the equation
\begin{equation}
\frac{dx}{dt} = u. \label{eq:path}
\end{equation}
\begin{figure}[htb]
\centering
\includegraphics[width=0.5\textwidth]{fig1} % grid1.eps
\caption{Departure and arrival points on grid.} \label{fig1}
\end{figure}
Note that the arrival point $x^{k+1}_A$ may not be any of $x_i$, $i=0, 1, \dots, M$
at $t = t_{k+1}$ as shown in Figure \ref{fig1}.
Then the arrival point $x^{k+1}_A$ is calculated by integrating the equation along
the trajectory from the departure point $x_D^{k}$ as
\[
x^{k+1}_A = x^{k}_D + \int_{t_k}^{t_{k+1}} u \, dt, \quad
x^k_D \in \{x_i: i = 0, 1, \dots, M\}.
\]
It is obvious that Euler's method will provide the first order approximation as
\[
x^{k+1}_A \approx x^{k}_D + \tau u(t_k, x^k_D).
\]
This is an explicit equation about $x^{k+1}_A$. It is possible to obtain a higher
order approximation by combining the computation of solution $u(t_{k+1}, x^{k+1}_A)$
at the time $t=t_{k+1}$, like, the modified Euler's method.
\section{One-step forward difference method}
For the problem \eqref{eq:tpde}, along the trajectory of the path, we propose a general
one-step forward difference method in the form
\begin{equation}
\begin{gathered}
w^0_D = u_0(x^0_D), \quad x^0_D \in \{x_i: i = 0, 1, \dots, M\},\\
x^{k+1}_A = x^{k}_D + \tau w^k_D, \\
w^{k+1}_A = w^k_D + \tau \phi(t_k,x^k_D,w^k_D,\tau,h), \\
w^{k+1}_D \quad \text{Calculated from linear interpolation of } w^{k+1}_A
\end{gathered}
\label{eq:tdiff}
\end{equation}
for $k = 0, 1, \dots, N-1$,
where $\phi(t,x,u,\tau,h)$ is a continuous function of $t$, $x$, $\tau$, and
$h$ and satisfies the Lipschitz condition in $u$. $x_D^{k}$ is the departure point
on the regular grids; $x_A^{k+1}$ is the arrival point and most likely is not any
of the regular grids like showed in the Figure \ref{fig1}, like $x_{A1}^{k+1}$ or
$x_{A2}^{k+1}$; $w^k_D$ and $w^{k+1}_A$ are approximations of $ u(t_k,x^k_D)$ and
$u(t_{k+1},x^{k+1}_A)$ respectively for $k = 0, 1, \dots, N-1$.
Note that $w^{k+1}_D$ is calculated from linear interpolation of $w^{k+1}_A$.
If we consider a higher order difference method, a higher order interpolation is needed.
Assuming that $\tau = O(h)$ (or $ \tau = \alpha h$, $\alpha > 0$ a constant),
and let $u(t,x)$ be the unique solution of the initial-value problem \eqref{eq:tpde}.
We define the local truncation error for each time step as follows.
\begin{definition} \rm
The local truncation error of the difference method \eqref{eq:tdiff} to the
problem \eqref{eq:tpde} is defined as
\[
\varepsilon_{k+1}(\tau)
= \frac{u^{k+1}_A - [u^k_D+\tau \phi(t_k,x^k_D, u^k_D, \tau, h)]}{\tau}
= \frac{u^{k+1}_A - u^k_D}{\tau} - \phi(t_k,x^k_D, u^k_D, \tau, h),
\]
where $u_D^{k} = u(t_{k}, x_D^{k})$ and $u_A^{k+1} = u(t_{k+1}, x_A^{k+1})$
on the departure and arrival points respectively.
\end{definition}
We need two lemmas whose proofs can be found in any numerical analysis book,
for example \cite{BF}.
\begin{lemma} \label{lem1}
For all $t \ge -1$ and any positive $m$, we have $0 \le (1+t)^m \le e^{mt}$.
\end{lemma}
\begin{lemma} \label{lem2}
If $s$ and $t$ are positive real numbers, $\{a_k\}_{k=0}^n$ is a sequence satisfying
$a_0 \ge -t/s$, and
\[
a_{k+1} \le (1+s) a_k+t, \quad \text{for each } k = 0, 1, 2, \dots, n-1,
\]
then
\[
a_{k+1} \le e^{(k+1)s}\big(a_0+\frac{t}{s}\big) - \frac{t}{s}.
\]
\end{lemma}
\begin{figure}[htb]
\centering
\includegraphics[width=0.5\textwidth]{fig2} % grid2.eps
\caption{Different arrival points on grid.} \label{fig2}
\end{figure}
\begin{theorem} \label{thm2}
The initial-value boundary problem \eqref{eq:tpde} is approximated by a one-step
difference method \eqref{eq:tdiff}. Suppose that there exists a number $\tau_0 > 0$
such that $\phi(t, x, w, \tau, h)$ is continuous and satisfies a Lipschitz condition
in all variables with Lipschitz constant $L$ on the set
\[
D = \{ (t,x,w, \tau,h) |\, 0 \le t \le T, a \le x \le b, -\infty < w < \infty, 0
\le \tau \le \tau_0, 0 \le h \le \alpha \tau_0\}.
\]
Then
\begin{itemize}
\item[(i)] the difference method \eqref{eq:tdiff} is stable;
\item[(ii)] the difference method \eqref{eq:tdiff} is convergent if and only
if it is consistent, which is equivalent to
\[
\phi(t,x,w,0,0) = g(t,x,w), \quad \text{for all } 0 \le t \le T, a \le x \le b,
-\infty < w < \infty;
\]
\item[(iii)] if a function $\varepsilon$ exists and, for each
$k = 0, 1, \dots, N$, the local truncation error $\varepsilon_{k}(\tau)$ satisfies
$| \varepsilon_{k}(\tau)| \le \varepsilon (\tau)$ whenever
$0 \le \tau \le \tau_0$, then there exists a constant $C$ such that
\begin{align*}
\max_{D}|u(t_k,x^k_D) - w^k_D|
&\le e^{T L}\max_{D}|u_0(x^0_D)-w^0_D|+\frac{\varepsilon(\tau)+C h}{L}(e^{T L}-1) \\
&\le \frac{\varepsilon(\tau)+C h}{L}(e^{T L}-1), \quad \text{for } k = 1, 2, \dots, N.
\end{align*}
\end{itemize}
\end{theorem}
\begin{proof}
Let $u$ and $v$ be two solutions of the initial-value problem (1) with the initial values
$u_0$ and $v_0$ respectively. There exist three constants $M$, $K_1$, and $K_2$ such that
\begin{equation} \label{e4}
|u|, |v| \leq M, \quad
|\frac{\partial u}{\partial x}|, |\frac{\partial u}{\partial t}|,
|\frac{\partial v}{\partial x}|, |\frac{\partial v}{\partial t}| \leq K_1,
\quad |\frac{\partial^2 u}{\partial x^2}|, |\frac{\partial^2 v}{\partial x^2}|\le K_2.
\end{equation}
and
\begin{equation} \label{e5}
\tau K_1 \leq 1/2.
\end{equation}
(i) Suppose that $\{u^{k}_{D}\}_{k=0}^{N}$ and $\{v^{k}_{D}\}_{k=0}^{N}$ satisfy
the difference equation \eqref{eq:tdiff} with the initial-values $u_0$ and $v_0$ respectively.
Let
\[
{E}^k_D = u^k_D-v^k_D \quad \text{for the point } (t_k,x^k_D)
\]
and
\[
{E}^k = \max \{|{ E}^k_D|, \text{ for all regular grid points (departure points) at }
t = t_k\}.
\]
Without lose generality, let the arrival points be as in the
Figure \ref{fig2}.
If $x^{k+1}_{D2}$ is located in between $x^{k+1}_{A_{u2}}$ and $x^{k+1}_{A_{v2}}$,
the procedure is very similar. Then
\begin{equation} \label{e6}
u^{k+1}_{D2} = \omega_u u^{k+1}_{A_{u2}}+(1-\omega_u)u^{k+1}_{A_{u1}}, \quad
v^{k+1}_{D2} = \omega_v v^{k+1}_{A_{v2}}+(1-\omega_v)v^{k+1}_{A_{v1}},
\end{equation}
where
\[
\omega_u = \frac{x^{k+1}_{D2}-x^{k+1}_{A_{u1}}}{x^{k+1}_{A_{u2}}-x^{k+1}_{A_{u1}}}, \quad
\omega_v = \frac{x^{k+1}_{D2}-x^{k+1}_{A_{v1}}}{x^{k+1}_{A_{v2}}-x^{k+1}_{A_{v1}}},
\]
and $0 \le \omega_u, \omega_v \le 1$.
Then from \eqref{e6},
\begin{equation} \label{e7}
\begin{aligned}
&|u^{k+1}_{D2} - v^{k+1}_{D2}| \\
& = |\omega_u u^{k+1}_{A_{u2}}+(1-\omega_u)u^{k+1}_{A_{u1}}
- \omega_v v^{k+1}_{A_{v2}}-(1-\omega_v)v^{k+1}_{A_{v1}}| \\
&\leq |\omega_u u^{k+1}_{A_{u2}}-\omega_u v^{k+1}_{A_{v2}}
+(1-\omega_u)u^{k+1}_{A_{u1}}-(1-\omega_u)v^{k+1}_{A_{v1}}| \\
&\quad +|\omega_u v^{k+1}_{A_{v2}}-\omega_v v^{k+1}_{A_{v2}}
+(1-\omega_u)v^{k+1}_{A_{v1}}-(1-\omega_v)v^{k+1}_{A_{v1}}| \\
&\leq \omega_u |u^{k+1}_{A_{u2}}-v^{k+1}_{A_{v2}}|
+(1-\omega_u)|u^{k+1}_{A_{u1}}-v^{k+1}_{A_{v1}}|
+|\omega_u-\omega_v\| v^{k+1}_{A_{v2}}-v^{k+1}_{A_{v1}}|.
\end{aligned}
\end{equation}
For each departure point, for example $x^k_{D1}$, it reads from \eqref{eq:tdiff},
\[
u^{k+1}_{A_{u1}} = u^k_{D1}+\tau \phi(t_k,x^k_{D1}, u^k_{D1}, \tau, h)
\text{ and }
v^{k+1}_{A_{v1}} = v^k_{D1}+\tau \phi(t_k,x^k_{D1}, v^k_{D1}, \tau, h).
\]
Then
\begin{equation} \label{e8}
\begin{aligned}
|u^{k+1}_{A_{u1}}-v^{k+1}_{A_{v1}}|
&\leq |u^k_{D1}-v^k_{D1}|+\tau |\phi(t_k,x^k_{D1}, u^k_{D1}, \tau, h)
-\phi(t_k,x^k_{D1}, v^k_{D1}, \tau, h)| \\
&\leq |u^k_{D1}-v^k_{D1}|+\tau L |u^k_{D1}-v^k_{D1}|
\leq (1+\tau L) |u^k_{D1}-v^k_{D1}| \\
&\leq (1+\tau L) |E^k_{D1}| \leq (1+\tau L ) E^k.
\end{aligned}
\end{equation}
Similarly,
\begin{equation} \label{e9}
|u^{k+1}_{A_{u2}}-v^{k+1}_{A_{v2}}| \leq (1+\tau L ) E^k.
\end{equation}
From the definition of $\omega_u$ and $\omega_v$ and \eqref{eq:tdiff}, we have
\begin{align*}
|\omega_u - \omega_v|
& = \big|\frac{x^{k+1}_{D2}-x^{k+1}_{A_{u1}}}{x^{k+1}_{A_{u2}}-x^{k+1}_{A_{u1}}}
- \frac{x^{k+1}_{D2}-x^{k+1}_{A_{v1}}}{x^{k+1}_{A_{v2}}-x^{k+1}_{A_{v1}}}\big| \\
& = \big|\frac{h-\tau u^k_{D1}}{h+\tau (u^k_{D2}-u^k_{D1})}
- \frac{h-\tau v^k_{D1}}{h+\tau (v^k_{D2}-v^k_{D1})}\big| \\
& = \frac{|h\tau (u^k_{D2}-v^k_{D2})+\tau^2 ( u^k_{D1}v^k_{D2}
-v^k_{D1}u^k_{D2})|}{|h+\tau (u^k_{D2}-u^k_{D1})\|h+\tau (v^k_{D2}-v^k_{D1})|} .
\end{align*}
From \eqref{e4}, it follows that
\[
h(1-\tau K_1) \leq |h+\tau (u^k_{D2}-u^k_{D1})| \leq h (1+\tau K_1).
\]
Similarly,
\[
h(1-\tau K_1) \leq |h+\tau (v^k_{D2}-v^k_{D1})| \leq h (1+\tau K_1).
\]
Then with help of \eqref{e4} and \eqref{e5},
\begin{equation} \label{e10}
\begin{aligned}
&|\omega_u - \omega_v| \\
&\leq \frac{\tau}{h(1-\tau K_1)^2} \big[ |u^k_{D2}-v^k_{D2}|
+\frac{\tau}{h}(|u^k_{D1}-v^k_{D1}\|v^k_{D2}|+|u^k_{D2}-v^k_{D2}\|v^k_{D1}|)\big] \\
& \leq \frac{\text{$\tau$}}{\text{$h(1-\tau K_1)^2$}}(1+2 \alpha M) E^k,
\quad \text{(where $\tau = \alpha h$)} \\
&\leq 4 \frac{\text{$\tau$}}{\text{$h$}} (1+ 2\alpha M) E^k.
\end{aligned}
\end{equation}
Putting \eqref{e8}, \eqref{e9}, and \eqref{e10} in \eqref{e7}, with the help of
\eqref{e4}, yields
\begin{equation} \label{e11}
\begin{aligned}
&|u^{k+1}_{D2} - v^{k+1}_{D2}| \\
& \leq \omega_u (1+\tau L ) E^k+(1-\omega_u)(1+\tau L ) E^k
+4 \frac{\text{$\tau$}}{\text{$h$}} (1+2 \alpha M) E^k h K_1 \\
&\leq [1+\tau L + 4\tau K_1 (1+2 \alpha M)] E^k
\end{aligned}
\end{equation}
i.e.
\begin{gather*}
|E^{k+1}_D| \leq [1+\tau L + 4\tau K_1 (1+2 \alpha M)] E^k , \\
E^{k+1} \leq (1+\tau L_1) E^k
\end{gather*}
where $L_1 = L+4K_1(1+2\alpha M)$.
Therefore, from Lemmas \ref{lem1} and \ref{lem2},
\[
{ E}^k \le (1+\tau L_1)^k { E}^0 \le e^{k\tau L_1} { E}^0 \le e^{TL_1} { E}^0,
\]
where $k \le N = T/\tau$.
So, the difference method \eqref{eq:tdiff} is stable.
\smallskip
(ii) Let $\phi(t,x,w,0,0) = g(t,x,w)$. From the assumptions of $\phi$, $g$
satisfies the conditions of Theorem \ref{thm1}, then the differential equation
\begin{gather*}
\frac{d v}{dt} = g(t, x, v), \quad 0 < t < T,\; a \le x \le b, \\
v(0,x) = u_0(x), \quad a \le x \le b,
\end{gather*}
has a solution $v(t,x)$. The numerical solution satisfies
\begin{equation} \label{e12}
\begin{gathered}
z^0_D = u_0(x^0_D), \quad x^0_D \in \{x_i: i = 0, 1, \dots, M\},\\
x^{k+1}_A = x^{k}_D + \tau \, z^k_D, \\
z^{k+1}_A = z^k_D + \tau \phi(t_k,x^k_D,z^k_D,\tau,h), \\
z^{k+1}_D \quad \text{ Calculated from linear interpolation of } z^{k+1}_A
\end{gathered}
\end{equation}
for $k = 0, 1, \dots, N-1$.
By the mean value theorem and referring to the Figure \ref{fig1},
\[
v(t_{k+1},x^{k+1}_{A2})
= v(t_k,x^k_{D2}) + \tau g(t_k+\xi \tau,x^k_{D2}+\eta h,v(t_k+\xi \tau,x^k_{D2}+\eta h)),
\]
for some constants $\xi$ and $\eta$ ($0 \le \xi, \eta \le 1$).
Similarly, define ${ E}^k_D = z^k_D - v(t_k,x^k_D)$ and ${ E}^k = \max |{ E}^k_D|$
on the regular grid points (departure points) at $t = t_k$.
Subtract this from the numerical solution, we have
\begin{align*}
&{ z}^{k+1}_{A2} -v(t_{k+1},x_{A2}^{k+1}) = { z}^k_{D2} - v(t_k,x^k_{D2}) \\
& + \tau [\phi(t_k,x^k_{D2},z^k_{D2},\tau,h)-g(t_k+\xi \tau,x^k_{D2}
+\eta h,v(t_k+\xi \tau,x^k_{D2}+\eta h))] \\
&= { E}^k_{D2} + \tau [\phi(t_k,x^k_{D2},z^k_{D2},\tau,h)
-\phi(t_k,x^k_{D2},v(t_k,x^k_{D2}),\tau,h) \\
&\quad + \phi(t_k,x^k_{D2},v(t_k,x^k_{D2}),\tau,h)-\phi(t_k,x^k_{D2},v(t_k,x^k_{D2}), \tau,0) \\
&\quad + \phi(t_k,x^k_{D2},v(t_k,x^k_{D2}),\tau,0)-\phi(t_k,x^k_{D2},v(t_k,x^k_{D2}), 0,0) \\
&\quad + \, \phi(t_k,x^k_{D2},v(t_k,x^k_{D2}),0,0)
- g(t_k+\xi \tau,x^k_{D2}+\eta h,v(t_k+\xi \tau,x^k_{D2}+\eta h))],
\end{align*}
Thanks to the Lipschitz condition of $\phi$,
\begin{gather*}
|\phi(t_k,x^k_{D2},z^k_{D2},\tau,h)-\phi(t_k,x^k_{D2},v(t_k,x^k_{D2}),\tau,h)|
\le L |z^k_{D2} -v(t_k,x^k_{D2})| \le L |{ E}^k_{D2}|,
\\
|\phi(t_k,x^k_{D2},v(t_k,x^k_{D2}),\tau,h)-\phi(t_k,x^k_{D2},v(t_k,x^k_{D2}), \tau,0)|
\le L h, \\
|\phi(t_k,x^k_{D2},v(t_k,x^k_{D2}),\tau,0)-\phi(t_k,x^k_{D2},v(t_k,x^k_{D2}), 0,0)|
\le L \tau,
\end{gather*}
and
\begin{align*}
&|\phi(t_k,x^k_{D2},v(t_k,x^k_{D2}),0,0)- g(t_k+\xi \tau,x^k_{D2}+\eta h,v(t_k
+\xi \tau,x^k_{D2}+\eta h))| \\
&\le |g(t_k,x^k_{D2},v(t_k,x^k_{D2})) - g(t_k+\xi \tau,x^k_{D2},v(t_k,x^k_{D2}))| \\
&\quad + |g(t_k+\xi \tau,x^k_{D2},v(t_k,x^k_{D2}))
- g(t_k+\xi \tau,x^k_{D2}+\eta h,v(t_k,x^k_{D2}))| \\
&\quad + |g(t_k+\xi \tau,x^k_{D2}+\eta h,v(t_k,x^k_{D2})) \\
&\quad - g(t_k+\xi \tau,x^k_{D2}+\eta h,v(t_k+\xi \tau,x^k_{D2}+\eta h))| \\
& \le L\xi \tau + L \eta h +L|v(t_k,x^k_{D2})-v(t_k+\xi \tau,x^k_{D2}+\eta h))| \\
& \le L\xi \tau + L \eta h +L K_1\xi \tau+L K_1 \eta h \\
& \le L\tau + L h +L K_1 \tau+L K_1 h \\
& \le L(1+K_1)(\tau+h).
\end{align*}
Then
\begin{equation} \label{e13}
\begin{aligned}
|{ E}^{k+1}_{A2}|
&\le |{ E}^k_{D2}| + L \tau |{ E}^k_{D2}| + L \tau h +L \tau^2 +L(1+K_1)\tau(\tau + h) \\
&\le (1+\tau L) { E}^k + L(2+K_1) \tau (\tau + h).
\end{aligned}
\end{equation}
From the Figure \ref{fig1}, define a constant $\omega$ as
\[
\omega = \frac{x_{D2}^{k+1}-x_{A1}^{k+1}}{x_{A2}^{k+1}-x_{A1}^{k+1}}.
\]
Then
\[
z_{D2}^{k+1} = \omega z_{A2}^{k+1}+(1-\omega) z_{A1}^{k+1}
\]
and there exists a constant $\eta$ such that
\begin{align*}
v(t_{k+1},x_{D2}^{k+1})
&= \omega v(t_{k+1},x_{A2}^{k+1})+(1-\omega)v(t_{k+1},x_{A1}^{k+1}) \\
&\quad +\frac{1}{2!}\frac{\partial ^2 v(t_{k+1},\eta)}{\partial x^2} (x_{D2}^{k+1}
- x_{A1}^{k+1})(x_{D2}^{k+1} - x_{A2}^{k+1}),
\end{align*}
Then
\begin{equation} \label{e14}
\begin{aligned}
&|{ z}^{k+1}_{D2} -v(t_{k+1},x_{D2}^{k+1})| \\
& \le \omega |{ z}^{k+1}_{A2} -v(t_{k+1},x_{A2}^{k+1})| \\
&\quad +(1-\omega)|{ z}^{k+1}_{A1} -v(t_{k+1},x_{A1}^{k+1})|
+\frac{1}{2}K_2 |x^{k+1}_{D2}-x^{k+1}_{A1}\|x^{k+1}_{D2}-x^{k+1}_{A2}| \\
&\le \omega |E^{k+1}_{A2}|+(1-\omega)|E^{k+1}_{A1}|
+\frac{1}{2}K_2 |h+\tau v(t_k,x^k_{D1})| \tau |v(t_k,x^k_{D2})| \\
&\le \omega |E^{k+1}_{A2}|+(1-\omega)|E^{k+1}_{A1}|+\frac{1}{2} K_2 M(1+\alpha M)\tau h.
\end{aligned}
\end{equation}
Thanks to \eqref{e13} and the similar result for $E^{k+1}_{A1}$, it follows that
\begin{align*}
E^{k+1} & \le (1+\tau L) { E}^k + L(2+K_1) \tau (\tau + h)
+\frac{1}{2}K_2 M (1+\alpha M) \tau h \\
&\le (1+\tau L) { E}^k + L_1 (\tau ^2+\tau h),
\end{align*}
where $L_1 = L(2+K_1)+\frac{1}{2}K_2 M (1+\alpha M)$.
Thanks to Lemma \ref{lem2} we have
\begin{gather*}
E^k \le e^{k\tau L} E^0 + \frac{L_1}{L}(\tau + h)(e^{k\tau L} -1),\\
E^k \le e^{T L} E^0 + \frac{L_1}{L}(\tau + h)(e^{T L} -1).
\end{gather*}
This converges to zero as $\tau$ and $h$ ($E^0 = 0$) go to zero.
So the numerical solution converges to the solution of initial-value boundary
problem \eqref{eq:tpde} if $f=g$. Therefore the numerical solution from the
difference method \eqref{eq:tdiff} is convergent.
On the other hand, if we have convergence, then $u$ and $v$ are the same.
Suppose $f$ and $g$ differ at some point.Then we consider the initial-value boundary
problem starting from that point. Two solutions should be different.
This leads to a contradiction.
\smallskip
(iii) Let ${ E}^k_A$, ${ E}^k_D$, and ${ E}^k$ are defined as in (i).
Then from the definition of the local truncation error, we have
\[
u(t_{k+1},x^{k+1}_A) = u(t_k,x^k_D)+\tau \phi(t_k,x^k_D, u(t_k,x^k_D), \tau, h)
+ \tau \varepsilon_{k+1}(\tau).
\]
Subtracting this from the difference method yields
\[
E^{k+1}_A = E^k_D + \tau [\phi(t_k,x^k_D, w^k_D, \tau, h)
-\phi(t_k,x^k_D, u(t_k,x^k_D), \tau, h)]-\tau \varepsilon_{k+1}(\tau).
\]
therefore,
\[
|E^{k+1}_A| \le |E^k_D| + \tau L |w^k_D - u(t_k,x^k_D)|
+ \tau |\varepsilon_{k+1}(\tau)| \le (1+\tau L) |E^k_D| + \tau |\varepsilon_{k+1}(\tau)|.
\]
Using a similar technique as in \eqref{e13} and \eqref{e14}, with the assumption
for $\varepsilon_k(\tau)$ leads to
\[
E^{k+1} \le (1+\tau L)E^k + \tau \varepsilon(\tau)+\frac{1}{2} K_2 M(1+\alpha M)\tau h.
\]
Then thanks Lemma \ref{lem2}, we have
\[
E^k \le e^{k\tau L}(E^0+\frac{\varepsilon(\tau)+C H}{L})
- \frac{\varepsilon(\tau)+C h}{L},
\]
i.e.,
\[
E^k \le e^{T L}E^0+\frac{\varepsilon(\tau)+C h}{L}(e^{T L}-1),
\]
for $k = 1, 2, \dots, N$ and $C = \frac{1}{2} K_2 M(1+\alpha M)$.
\end{proof}
\section{Semi-Lagrangian forward methods}
Now we study the modified semi-Lagrangian forward Euler method for the initial-value
boundary problem \eqref{eq:tpde},
\begin{equation}
\begin{gathered}
w^0_D = u_0(x^0_D), \quad x^0_D \in \{x_i: i = 0, 1, \dots, M\},\\
k_1 = g(t_k,x_D^k,w^k_D), \\
\bar{w} = w^k_D + \tau k_1, \\
x^{k+1}_A = x^{k}_D + \frac{1}{2}\tau (w^k_D+\bar{w}), \\
k_2 = g(t_{k+1},x^{k+1}_A, \bar{w}), \\
w^{k+1}_A = w^k_D + \frac{1}{2} \tau (k_1 + k_2), \\
w^{k+1}_D \text{ Calculated from second-order interpolation of } w^{k+1}_A \\
\end{gathered}
\label{eq:tdiff2}
\end{equation}
for $k = 0, 1, \dots, N-1$,
for all of arrival points $A$ and $k = 0, 1, 2, \dots, N-1$.
\begin{theorem} \label{thm3}
If the function $g(t,x,u)$ in \eqref{eq:tpde} is continuous and satisfies a Lipschitz
condition in all variables with Lipschitz constant $L$ on the set
\[
D = \{ (t,x,w) : 0 \le t \le T, a \le x \le b, -\infty < u < \infty\}.
\]
Then
\begin{itemize}
\item[(i)] The difference method \eqref{eq:tdiff2} is stable;
\item[(ii)] The difference method \eqref{eq:tdiff2} is convergent.
\end{itemize}
\end{theorem}
The stability and convergence of \eqref{eq:tdiff2} could be proved following the
similar idea with Lipschitz continuity in Theorem \ref{thm2}.
\begin{theorem} \label{thm4}
Suppose $h = O(\tau)$ and the function $g(t,x,u)$ in \eqref{eq:tpde} is continuous
and differentiable. If $\varepsilon_{k}$ is the local truncation error of \eqref{eq:tdiff2},
then
\[
\|\varepsilon_{k}(\tau)\|_{L_{\infty}} = O(\tau^2).
\]
\end{theorem}
\begin{proof}
Let $u(t,x)$ be the unique solution with up to order three continuous bounded partial
derivatives on $[0,T] \times [a,b]$, so that for each $k = 0, 1, \dots, N-1$,
\begin{align*}
& g(t_{k+1},x^{k+1}_A, w^k_D+\tau g(t_k,x^k_D,w^k_D)) \\
&= g(t_k,x_D^k, w_D^k)+\tau \frac{\partial g(t_k,x_D^k, w_D^k)}{\partial t}
+\tau w_D^k \frac{\partial g(t_k,x_D^k, w_D^k)}{\partial x} \\
&\quad +\tau g(t_k,x_D^k, w_D^k) \frac{\partial g(t_k,x_D^k, w_D^k)}{\partial u}+O(\tau^2),
\end{align*}
and
\begin{align*}
u(t_{k+1}, x^{k+1}_A)
& = u(t_k,x^k_D)+\tau \frac{du(t_k,x^k_D)}{dt}
+\frac{1}{2}\tau^2 \frac{d^2u(t_k,x^k_D)}{dt^2}+O(\tau^3) \\
&= u(t_k,x^k_D)+\tau g(t_k,x_D^k, u_D^k)+\frac{1}{2} \tau^2
[ \frac{\partial g(t_k,x_D^k, u_D^k)}{\partial t} \\
&\quad + u_D^k \frac{\partial g(t_k,x_D^k, u_D^k)}{\partial x}
+ g(t_k,x_D^k, u_D^k) \frac{\partial g(t_k,x_D^k, u_D^k)}{\partial u}] +O(\tau^3).
\end{align*}
From the definition of local truncation error and the method (15), the local
truncation error reads
\begin{align*}
\varepsilon_{k+1}(\tau)
&= \frac{u(t_{k+1}, x^{k+1}_A)-u(t_k,x^k_D)}{\tau}
-\frac{1}{2}\big[g(t_k,x_D^k, u(t_k,x^k_D)) \\
&\quad +g(t_{k+1},x^{k+1}_A, u(t_k,x^k_D)+\tau g(t_k,x_D^k, u(t_k,x^k_D))]
\end{align*}
So it is easy to see that the local truncation error at the $(k+1)$th step
for the method \eqref{eq:tdiff2} is
\[
\|\varepsilon_{k+1}(\tau)\|_{L_{\infty}} = O(\tau^2). \qedhere
\]
\end{proof}
Note that the big difference in forward and backward methods is the computing
of arrival points and departure point. In backward method, it needs the iteration
to find better approximation.
Using similar techniques, the semi-Lagrangian method of order-three is constructed
as follows
\begin{equation}
\begin{gathered}
w^0_D = u_0(x^0_D), \quad x^0_D \in \{x_i: i = 0, 1, \dots, M\} \\
k_1 = g(t_k,x_D^k,w^k_D) \\
x^{k+1}_A = x^{k}_D + \tau \, w^k_D +\frac{1}{2} \tau^2 k_1 \\
k_2 = g(t_k+\frac{1}{2} \tau,x_D+\frac{1}{2}\tau w_D^k,w^k_D+\frac{1}{2}\tau k_1) \\
k_3 = g(t_{k+1},x_A^{k+1},w^k_D-\tau k_1+2 \tau k_2) \\
w^{k+1}_A = w^k_D + \frac{1}{6} \tau (k_1+4k_2+k_3) \\
w^{k+1}_D \text{ Calculated from the third-oder interpolation of } w^{k+1}_A \\
\end{gathered} \label{eq:tdiff3}
\end{equation}
for $k = 0, 1, \dots, N-1$.
The numerical analysis is much more complicated because of the need of higher
order methods for computing of departure points and interpolation at the departure points.
\section{Numerical experiments}
The main computational effort in applying the Semi-Lagrangian forward methods (SLFW)
is the evaluation of $g$ and the interpolation of approximation on departure points.
In the second-order method, the local truncation error is $O(\tau^2)$, and the cost is
two functional evaluations per step and the computation of the arrival points.
The Semi-Lagrangian method of order three requires three evaluations per step and
the computation of the arrival points, and the local truncation error is $O(\tau^3)$.
\begin{figure}[htb]
\centering
\includegraphics[width=0.7\textwidth]{fig3} % sul01.eps
\caption{Approximated Solution of the initial-value problem.}
\label{fig3}
\end{figure}
\begin{example} \label{examp1}\rm
Consider the problem
\begin{equation}
\begin{gathered}
u_t+u \, u_x = t+\sin(2 \pi x) +2 \pi u\, t \cos(2 \pi x), \\
u(0,x) = 0, \quad 0 < t \le 1,\; 0 < x < 1,
\end{gathered}\label{eq:example1}
\end{equation}
with periodic boundary condition for $x$.
\begin{figure}[htb]
\centering
\includegraphics[width=0.7\textwidth]{fig4} % contn.eps
\caption{Contour lines of the difference between the solution and its approximation.}
\label{fig4}
\end{figure}
The exact solution for this problem is $u(t,x)=\frac{1}{2}t^2+t \sin(2 \pi x)$.
Figure \ref{fig3} shows the graph of the approximate solution of $u(t,x)$ by
using the Semi-Lagrangian method of order three with $\tau = 0.01$ and $h = 0.01$.
To focus on the truncation error from the proposed Semi-Lagrangian forward methods,
the periodic boundary condition is supplied. For other boundary conditions,
the treatment near boundary is more complicated when the trajectories of the
departure points hit the boundary.
\begin{figure}[ht]
\centering
\includegraphics[width=3.50in]{fig5} % error.eps
\caption{Truncation errors: $-*-$ first-order Semi-Lagrangian forward method (SLFW1);
$-+-$ second-order method (SLFW2); $-\circ -$ third-order method (SLFW3).}
\label{fig5}
\end{figure}
\begin{table}
\caption{Numerical orders of the Semi-Lagrangian forward method (SLFW1), SLFW2,
and SLFW3 as applied to the problem \eqref{eq:example1}} \label{table1}
\centering
\renewcommand{\arraystretch}{1.2}
\begin{tabular}{|c|l|l|l|l|}\hline
$\tau$ & 0.1 & 0.05 & 0.025 & 0.02 \\ \hline
SLFW1 $(\frac{\|u-u_0\|_{\infty}}{\tau})$ & 1.2189 & 1.1656 & 1.1347 & 1.1270 \\ \hline
SLFW2 $(\frac{\|u-u_0\|_{\infty}}{\tau^2})$ & 2.0234 & 2.0242 & 2.0263 & 2.0262 \\ \hline
SLFW3 $(\frac{\|u-u_0\|_{\infty}}{\tau^3})$ & 2.7279 & 2.60209 & 2.5076 & 2.4823 \\ \hline
\end{tabular}
\end{table}
Figure \ref{fig4} shows the contour lines of the difference between the exact solution
and the approximation from the Semi-Lagrangian forward method of order three at the common
grid points. The truncation error is increasing as $t$ goes from $0$ to $0.5$.
The contour lines are matching the changes of the function. The larger errors occur
when the function has the bigger changes. The maximum error reached at about $t=0.5$ for
some points of $x$ as we expected.
\end{example}
The first-order Semi-Lagrangian forward method, the second-order Semi-Lagrangian
forward method, and the Semi-Lagrangian forward method of order three are compared with
$\tau = 0.01$ and $h = 0.005$ in the Figure \ref{fig5}. The maximum truncation error presented
at the each step of $t$ for those three methods.
\begin{figure}[htb]
\centering
\includegraphics[width=0.7\textwidth]{fig6} % diffn.eps
\caption{Scaling of the error of Semi-Lagrangian methods of order one, two and three to the problem \eqref{eq:example1}}
\label{fig6}
\end{figure}
In the table \ref{table1} and the Figure \ref{fig6}, we studied the numerical
orders of convergence with $ h = (1/2)\tau$. It is clear that the first-order
Semi-Lagrangian forward method and the second-order Semi-Lagrangian forward method
have the order as we expected. However, the Semi-Lagrangian forward method of order
three is not exactly order three, but is close to 2.5. The main difficulty for
Semi-Lagrangian method of order three are the calculation of the departure points
and the third-order interpolation from arrival points(non-uniform grids) to departure
points (uniform grids).
\begin{example} \label{examp2} \rm
Consider the Kortweg-deVries (KdV) equation \cite{TA}
\[
u_t + u u_x + \delta^2 u_{xxx} = 0
\]
where $\delta = 0.022$ with the initial conditions
$$
u(x,0) = \cos(\pi x), \qquad 0 \le x \le 2
$$
and $u_x$, $u_{xx}$, and $u_{xxx}$ are periodic on $[0,2]$ for all $t$.
\begin{figure}[htb]
\centering
\includegraphics[width=0.7\textwidth]{fig7} % kdv1.eps
\caption{First-order SLFW method for $\delta = 0.022$: initial values and solutions
at $t= 1/\pi$ and $t = 3.6/\pi$.} \label{fig7}
\end{figure}
\begin{figure}[htb]
\centering
\includegraphics[width=0.7\textwidth]{fig8} % kdv2.eps
\caption{Second-order SLFW method for $\delta = 0.022$: initial values and solutions
at $t= 1/\pi$ and $t = 3.6/\pi$.} \label{fig8}
\end{figure}
The KdV equation describes the long time asymptotic behaviour of small but finite
amplitude shallow water waves in one dimension. In the equation, $u=u(x,t)$ measures
the elevation (the height of water above the equilibrium level) at time $t$ and position $x$.
Two different mechanisms are involved: Non-linearity ($uu_x$) and Dispersion ($u_{xxx}$).
The delicate balance between these two effects leads to travelling waves of permanent
form, which is called Solitary Waves, or Solitons.
Various numerical methods are used to approximate the KdV equations. The comparison of
efficiencies of different methods was studied in \cite{SCN}. The aim here is to show the
efficiency of SLFW as an explicit method applying to KdV equations.
It is possible to apply higher order (higher than order two) to the KdV equations.
By using SLFW, the nonlinear term was perfectly combined in the direction of moving waves.
In Figure \ref{fig7}, the first-order SLFW is applied to the KdV equations with the initial
condition of $u(x,t) = \cos(\pi x)$ on the interval $[0,2]$ with $\Delta x = 0.01$
and $\Delta t = 0.00002$. The wave profile at $t=1/\pi$ does not produce a shock,
but close. When computed at time $t=3.6/\pi$, another shock is also present.
The trains of eight solitons are detected, but not well-defined.
The results present in Figure \ref{fig8} are from applying the second-order SLFW
to the KdV equations with the initial condition of $u(x,t) = \cos(\pi x)$ on the
interval $[0,2]$ with $\Delta x = 0.01$ and $\Delta t = 0.00002$.
The dot line is the initial wave profile. The dash line presents the wave profile
at $t=1/\pi$. The wave almost produces a shock and has a noticeable oscillation for $x<1/2$.
At time $t=3.6/\pi$ the profile shows a train of eight well-defined waves which have
developed since the early stage of evolution.
\end{example}
\section{Conclusions}
In this article, one-step Semi-Lagrangian forward methods were proposed and studied for
the first-order time-dependent nonlinear partial differential equations.
The second-order Semi-Lagrangian forward method was proved to be consistent and stable.
When $\tau$ goes to zero, the numerical solution converges to the solution of the
corresponding initial-boundary problem.
Theorem 2 in the section 4 and the examples in the section 6 have shown that the
Semi-Lagrangian forward methods can be applied to some time-dependent partial differential
equations. They are very convenient for implementation. The only disadvantage is that at
each time step, spatial linear or higher order interpolation is needed to find the
solutions on the departure points.
In the example, only one spatial variable was used. However, we believe that Semi-Lagrangian
forward methods should work for more spatial variables, but more complicated for numerical
analysis. The computing cost will increase a lot, like interpolation for more variables
at each step. We will report the similar results for the systems with more spatial
variables. Also more terms involving higher order derivative only in spatial variables
will be added to the right side function later. More work of applications to
Shallow-water equations and Navier-Stokes equations will be reported later.
\begin{thebibliography}{00}
\bibitem{BB} A. Bourchtein, L. Bourchtein;
\emph{Semi-Lagrangian semi-implicit time-splitting scheme for the shallow water equations},
Int. J. Numer. Meth. Fluids \textbf{54} (2007) 453--471.
\bibitem{BS} S. Bak;
\emph{High-order characteristic-tracking strategy for simulation of a nonlinear
advection--diffusion equation}, Numer Methods Partial Differential Eq. \textbf{35} (2019) 1756--1776.
\bibitem{BKK} S. Bak, P. Kim, D. Kim;
\emph{A semi-Lagrangian approach for numerical simulation of coupled Burgersâ€™ equations},
Commun Nonlinear Sci Numer Simulat \textbf{69} (2019) 31--44.
\bibitem{BM} J. R. Bates, A. McDonald;
\emph{Multiply-upstream, semi-Lagrangian advective schemes: analysis and application to
a multi-level primitive equation model},
Monthly Weather Review, vol \textbf{110}, (1982) 1831--1842.
\bibitem{RBJC} R. Bermejo, J. Carpio;
\emph{An adaptive finite element semi-Lagrangian implicit-explicit Runge-Kutta-Chebyshev
method for convection dominated reaction-diffusion problems},
Applied Numerical Mathematics \textbf{58} (2009) 16--39.
\bibitem{BF} R. L. Burden, J. D. Faires;
\emph{Numerical Analysis}, Thomason Brooks/Cole, 2005.
\bibitem{BFR} L. Bonaventura, R. Ferretti, L. Rocchi;
\emph{A fully semi-Lagrangian discretization for the 2D incompressible Navier-Stokes
equations in the vorticity-streamfunction formulation},
Applied Mathematics and Computation \textbf{323} (2018) 132--144.
\bibitem{CFN} M. P. Calvo, J. de Frutos, J. Novo;
\emph{Linearly implicit Runge-Kutta methods for advection-reaction-diffusion equations},
Applied Numerical Mathematics \textbf{37} (2001) 535--549.
\bibitem{CAR} M. F. Carfora;
\emph{An unconditionally stable semi-Lagrangian method for
the spherical atmospherical shallow water equations},
Int. J. Numer. Meth. Fluids \textbf{34} (2000) 527--558.
\bibitem{CRS} N. Crouseilles, T. Respaud, E. SonnendrĂĽcker;
\emph{A forward semi-Lagrangian method for the numerical solution of the Vlasov equation},
Computer Physics Communications \textbf{180} (2009) 1730--1745.
\bibitem{EKSV} A. Efremov, E. Karepova, V. Shaydurov, A. Vyatkin;
\emph{A Computational Realization of a Semi-Lagrangian Method for Solving the Advection
Equation}, Journal of Applied Mathematics
Volume 2014, Article ID 610398, 12 pages, https://doi.org/10.1155/2014/610398.
\bibitem{GEAR} C. W. Gear;
\emph{Numerical Initial Value Problems in Ordinary Differential Equations},
Prentice-Hall, Englewood Cliffs, NJ, 1971.
\bibitem{GH} D. F. Griffiths, D. J. Higham;
\emph{Numerical Methods for Ordinary Differential Equations: Initial Value Problems},
Springer Undergrauate Mathematics Series, Springer-Verlag London Limited, 2010.
\bibitem{DGJD1} D. X. Guo, J. B. Drake;
\emph{A global semi-Lagrangian spectral model for the reformulated shallow water equations},
Dynamical systems and differential equations, Discrete Contin. Dyn. Syst. suppl., (2003),
375--385.
\bibitem{DGJD2} D. X. Guo, J. B. Drake;
\emph{A global semi-Lagrangian spectral model of the shallow water equations with variable
resolution}, J. Comput. Phys. \textbf{206} (2005) 559--577.
\bibitem{DG3} D. X. Guo;
\emph{A semi-Lagrangian Runge-Kutta Method for time-dependent partial differential equations},
J. of Applied Analysis \& Computation, \textbf{3} (3) (2013) 251--263.
\bibitem{DG4} D. X. Guo;
\emph{On stability and convergence of semi-Lagrangian methods for the first-order
time-dependent nonlinear partial differential equations in 1D}, J. of Computational
and Applied Mathematics, \textbf{324} (2017) 72--84.
\bibitem{LKM} P. H. Lauritzen, E. Kaas, B. Machenhauer;
\emph{A massive-conservative semi-implicit semi-Lagrangian limite-area shallow-water model
on the sphere}, Monthly Weather Review, vol \textbf{134}, (2006) 1205--1221.
\bibitem{MWL} M. Mortezazadeh, L. Wang;
\emph{A high-order backward forward sweep interpolating algorithm for semi-Lagrangian method},
Int. J. Numer. Meth. Fluids, \textbf{84} (2017) 584--597
\bibitem{NSS} R. D. Nair, J. S. Scroggs, F. H. Semazzi;
\emph{A forward-trajectory global semi-Lagrangian transport scheme},
J. Comput. Phys., vol \textbf{190} (2003), 275--294.
\bibitem{ROB} A. Robert;
\emph{A stable numerical integration scheme for the primitive meteorological equations},
Atmos. Ocean., vol \textbf{19} (1981), 35--46.
\bibitem{MS} M. Seaid;
\emph{An Eulerian-Lagrangian method for parabolic-hyperbolic equations},
Applied Numerical Mathematics 59 (2009), 745--768.
\bibitem{SCN} M. Shahrill, M. Chong, \& H. Nor;
\emph{Applying Explicit Schemes to the Korteweg-de Vries Equation},
Modern Applied Science, Vol. \textbf{9}, No. 4 (2015), 200--224.
\bibitem{SCT} A. Staniforth, J. Cote;
\emph{Semi-Lagrangian integration schemes for atmospheric models - A review},
Monthly Weather Review, vol \textbf{119} (1991), 2206--2223.
\bibitem{TA} T. R. Taha, M. J. Ablowitz;
\emph{Analytical and Numerical Aspects of Certain Nonlinear Ecolution Equations. III.
Numerical, Korteweg-de Vries Equation}, J. Comput. Phys., Vol. \textbf{55} (1984), 231--253.
\bibitem{RT} R. Temam;
\emph{Navier-Stokes equations -- Theory and numerical analsys},
North-Holland Publishing Company, 1977.
\bibitem{JGV} J. G. Verwer;
\emph{Explicit Runge-Kutta methods for parabolic partial differential equations},
Applied Numerical Mathematics \textbf{22} (1996), 359--379.
\end{thebibliography}
\end{document}