Semi-Lagrangian one-step methods for two classes of time-dependent partial differential systems

TL;DR

This paper introduces semi-Lagrangian one-step methods utilizing Runge-Kutta schemes for efficiently solving time-dependent PDE systems, with analysis of stability, convergence, and error bounds.

Contribution

It develops a novel semi-Lagrangian approach for systems of PDEs, combining ODE solvers along characteristics with high-order accuracy and rigorous numerical analysis.

Findings

Achieves high order local truncation error using Runge-Kutta methods

Provides stability and convergence analysis for the proposed methods

Establishes maximum error bounds for the numerical solutions

Abstract

Semi-Lagrangian methods are numerical methods designed to find approximate solutions to particular time-dependent partial differential equations (PDEs) that describe the advection process. We propose semi-Lagrangian one-step methods for numerically solving initial value problems for two general systems of partial differential equations. Along the characteristic lines of the PDEs, we use ordinary differential equation (ODE) numerical methods to solve the PDEs. The main benefit of our methods is the efficient achievement of high order local truncation error through the use of Runge-Kutta methods along the characteristics. In addition, we investigate the numerical analysis of semi-Lagrangian methods applied to systems of PDEs: stability, convergence, and maximum error bounds.

Equations

$$\left\{\begin{array}[]{cccc}\dfrac{\partial y_{1}}{\partial t}+\omega\dfrac{\partial y_{1}}{\partial x}&=&f_{1}(t,x,y_{1},\dots,y_{n}),&\ \ y_{1_{0}}=y_{1}(0,x),\\ \vdots&\vdots&\vdots&\vdots\\ \dfrac{\partial y_{n}}{\partial t}+\omega\dfrac{\partial y_{n}}{\partial x}&=&f_{n}(t,x,y_{1},\dots,y_{n}),&\ \ y_{n_{0}}=y_{n}(0,x),\end{array}\right.$$

$$\left\{\begin{array}[]{cccc}\dfrac{\partial u}{\partial t}+u\dfrac{\partial u}{\partial x}+v\dfrac{\partial u}{\partial y}&=&f(t,x,y,u,v),&\ \ u_{0}=u(0,x,y),\\ &&&\\ \dfrac{\partial v}{\partial t}+u\dfrac{\partial v}{\partial x}+v\dfrac{\partial v}{\partial y}&=&g(t,x,y,u,v),&\ \ v_{0}=u(0,x,y),\\ \end{array}\right.$$

$$\frac{dF}{dt}=\frac{\partial F}{\partial t}+\frac{dx}{dt}\frac{\partial F}{\partial x}=0,$$

$$\frac{dx}{dt}=U(t,x)$$

$$\frac{F(t_k+\Delta t,x_A)-F(t_k-\Delta t,x_A-2\Delta x)}{2\Delta t}=0.$$

$$\Delta x=U(t_k,x_A-\Delta x)\Delta t,$$

$$\begin{array}[]{rll}\int_{x_{D}^{k}}^{x_{A}^{k+1}}dx&=&\int_{t_{k}}^{t_{k+1}}\omega dt\vspace{2mm}\\ x_{D}^{k}&=&x_{A}^{k+1}-\int_{t_{k}}^{t_{k+1}}\omega dt\vspace{2mm}\\ &\approx&x_{A}^{k+1}-\omega\tau,\end{array}$$

$${\mathbf{y}}_t+\omega {\mathbf{y}}_x=\mathbf{f}(t,x,\mathbf{y}),0\le t\le T,a\le x\le b,\Vert \mathbf{y}{\Vert }_{\infty }<\infty .$$

$$\frac{d\mathbf{y}}{dt}=\mathbf{f}(t,x,\mathbf{y}),\mathbf{y}(0,x)={\mathbf{y}}_0(x),(t,x)\in [0,T]\times [a,b],\Vert \mathbf{y}{\Vert }_{\infty }<\infty .$$

$$\mathbf{y}(t_{k+1},x_A^{k+1})=\mathbf{y}(t_k,x_D^k)+{\int }_{t_k}^{t_{k+1}}\mathbf{f}(t,x,\mathbf{y})dt.$$

$${\int }_{t_k}^{t_{k+1}}\omega dt\approx \frac{\tau }{2}({\omega }_D^k+{\omega }_A^{k+1}),$$

$${\omega }_D^k=\omega (t_k,x_D^k,{\mathbf{y}}_D^k)$$

$${\omega }_A^{k+1}=\omega (t_{k+1},x_A^{k+1},{\mathbf{y}}_A^{k+1}).$$

$${\int }_{t_k}^{t_{k+1}}\omega dt\approx \frac{\tau }{6}({\omega }_D^k+4{\tilde{\omega }}_I+{\omega }_A^{k+1}),$$

$${\omega }_D^k=\omega (t_k,x_D^k,{\mathbf{y}}_D^k),$$

$${\omega }_A^{k+1}=\omega (t_{k+1},x_A^{k+1},{\mathbf{y}}_A^{k+1}),$$

$$\tilde{\omega}_{I}=\omega\Big{(}t_{k+1}+\frac{\tau}{2},x_{D}^{k}+\frac{\tau}{2}\omega_{D}^{k},\tilde{\boldsymbol{y}}_{I}\Big{)}.$$

$$a_{k+1}\le e^{(k+1)s}(a_0+\frac{t}{s})-\frac{t}{s}.$$

$$|p(x)-q(x)|\le \underset{0\le i\le M}{\max }|u_i-v_i|.$$

$$\left\{\begin{array}[]{lr}\boldsymbol{y}_{A}^{0}=\boldsymbol{y}_{0}(x_{A}^{0}),&x_{A}^{0}=x_{A}^{k}=x_{i},\ \text{for}\ i=0,\dots,M,\vspace{2mm}\\ x_{D}^{k}=x_{A}^{k+1}-\omega_{D}^{k}\tau,&\vspace{2mm}\hfil\\ \boldsymbol{y}_{A}^{k+1}=\boldsymbol{y}_{D}^{k}+\tau\boldsymbol{\phi}(t_{k},x_{D}^{k},\boldsymbol{y}_{D}^{k},\tau,h)&\text{for each}\ k=0,\dots,N-1.\end{array}\right.$$

$$D=\{(t,x,\mathbf{y},\tau ,h):(t,x)\in [0,T]\times [a,b],\Vert \mathbf{y}{\Vert }_{\infty }<\infty ,0\le \tau \le {\tau }_0,0\le h\le \alpha {\tau }_0\}.$$

$$ma{x}_{x_A}\Vert \mathbf{y}(t_k,x_A^k)-{\mathbf{y}}_A^k{\Vert }_{\infty }\le e^{TL}\underset{x_A}{\max }\Vert {\mathbf{y}}_0(x_A)-{\mathbf{y}}_A^0{\Vert }_{\infty }+\frac{\epsilon (\tau )}{L}(e^{TL}-1)$$

$${\mathbf{E}}_A^k={\mathbf{u}}_A^k-{\mathbf{v}}_A^k$$

$${\mathbf{E}}_D^k={\mathbf{u}}_D^k-{\mathbf{v}}_D^k$$

$${\mathbf{E}}^k=\underset{x}{\max }\Vert {\mathbf{E}}_A^k{\Vert }_{\infty }$$

$${\mathbf{u}}_A^{k+1}-{\mathbf{v}}_A^{k+1}={\mathbf{u}}_D^k-{\mathbf{v}}_D^k+\tau [\mathbf{\varphi }(t_k,x_D^k,{\mathbf{u}}_D^k,\tau ,h)-\mathbf{\varphi }(t_k,x_D^k,{\mathbf{v}}_D^k,\tau ,h)]$$

$$\begin{array}[]{lll}\|\boldsymbol{u}_{A}^{k+1}-\boldsymbol{v}_{A}^{k+1}\|_{\infty}&=&\|\boldsymbol{u}_{D}^{k}-\boldsymbol{v}_{D}^{k}+\tau[\boldsymbol{\phi}(t_{k},x_{D}^{k},\boldsymbol{u}_{D}^{k},\tau,h)-\boldsymbol{\phi}(t_{k},x_{D}^{k},\boldsymbol{v}_{D}^{k},\tau,h)]\|_{\infty}\vspace{2mm}\\ &\leq&\|\boldsymbol{u}_{D}^{k}-\boldsymbol{v}_{D}^{k}\|_{\infty}+\tau L\|\boldsymbol{u}_{D}^{k}-\boldsymbol{v}_{D}^{k}\|_{\infty}\vspace{2mm}\\ &=&(1+\tau L)\|\boldsymbol{u}_{D}^{k}-\boldsymbol{v}_{D}^{k}\|_{\infty}\end{array}$$

$$\Vert {\mathbf{E}}_A^{k+1}{\Vert }_{\infty }\le (1+\tau L)\Vert {\mathbf{E}}_D^k{\Vert }_{\infty }.$$

$$\Vert {\mathbf{E}}_A^{k+1}{\Vert }_{\infty }\le (1+\tau L)\Vert {\mathbf{E}}^k{\Vert }_{\infty }$$

$$\Vert {\mathbf{E}}^{k+1}{\Vert }_{\infty }\le (1+\tau L)\Vert {\mathbf{E}}^k{\Vert }_{\infty }.$$

Full text

Nikolai D. Lipscomb and Daniel X. Guo

Department of Mathematics and Statistics

University of North Carolina Wilmington

Wilmington, North Carolina, USA

2010 Mathematics Subject Classification. Primary: 65M25; Secondary: 35Q35, 65M12.

1 Introduction

Time-dependent partial differential equations are at the core of particle physics. Due to the difficult and often analytically unsolvable nature of most of these equations, the next best approach is to use an algorithm to approximate the solution. Research on numerical computation for approximating solutions to advection equations goes back to the 1950s with finite difference method approaches to nonlinear hyperbolic partial differential equations by Courant et al. [3] and numerical integration of the barotropic vorticity equation by Fjørtoft [6]. In the fields of weather forecasting and climate modelling, particle trajectory methods were proposed by Wiin-Nielsen [14] which led to the only reliable forecasting model of its time. Meanwhile, researchers in plasma physics saw the promise of semi-Lagrangian approaches; for example, Cheng and Knorr [2] produced an efficient numerical splitting scheme for solving the Vlasov-Maxwell equations. While there was much research on advection processes and weather prediction over the following decades, the 1980s produced a slew of research that brought characteristic-based methods into different numerical approaches: Douglas Jr. and Russell [5] brought the method of characteristics to finite difference and finite element methods, André Robert’s meterological contributions produced stable numerical solutions to the shallow-water equations [11], and many more. Today, semi-Lagrangian models are frequently used by organisations that focus on atmospheric modelling such as the European Centre for Medium-Range Weather Forecasts (ECMWF) [4], the National Oceanic and Atmospheric Administration (NOAA), the National Center for Atmospheric Research [9], and the High Resolution Local Area Modelling (HIRLAM) programme.

The popularity of semi-Lagrangian methods in today’s atmospheric models lies in the resolution problem. Eulerian schemes’ accuracy is dependent on the resolution of the solution grid, a function of the problem domain’s discretisation–specifically the temporal and spatial discretisation. For earlier Eulerian schemes, resolution was greatly dependent on stability, which demanded a very small time discretisation relative to the spatial discretisation [13]. The gradual development of methods to avoid such constraints brought attention to semi-Lagrangian schemes: a pairing of the equal spacing of solutions from an Eulerian approach and the particle-tracing of a Lagrangian approach. Modern semi-Lagrangian numerical schemes perform with great numerical stability under a wide range of resolutions and produce little numerical dispersion [4]. In order to maintain competitiveness in the near future, semi-Lagrangian-dependent numerical schemes must continue to improve: reduction of error while maintaining sufficiently fast computation time considers not just the resolution of the Eulerian grid, but the order of the error. Most numerical schemes in weather applications achieve second order results with respect to the spatial and time discretisations. Further, while developed from physical laws, a complete numerical analysis of semi-Lagrangian theory is still underway.

We will examine numerical methods for solving initial value problems involving systems of time-dependent partial differential equations (PDEs). The particular types of PDEs we will examine are PDEs that describe the advection process, mostly found in fluid dynamics and atmospheric modelling. We will consider two general cases for systems of time-dependent PDEs. The first system, general advection in one dimension, is

[TABLE]

where $y_{i}=y_{i}(t,x)$, $(t,x)\in[0,\infty)\times[a,b]$ and $\omega$ belongs to one of three cases:

$\omega\in\mathbb{R}$, 2. 2.

$\omega=\omega(t,x)$, $|\omega|<\infty$, 3. 3.

$\omega=\omega(t,x,\boldsymbol{y})$, $|\omega|<\infty$.

The second system, nonlinear advection in two dimensions, is of the form

[TABLE]

where $u=u(t,x,y)$, $v=v(t,x,y)$, $(t,x,y)\in[0,\infty)\times[a,b]\times[c,d]$.

An Eulerian scheme examines a prescribed set of points and examines how the solution to the PDEs change at these points as time goes on. Since these time-dependent PDEs model physical processes, such as the movement of particles, a Lagrangian approach would examine individual particles (or parcels) with respect to the solution and trace their trajectory, examining how the solution updates at the new arrival points as time goes on. A semi-Lagrangian method is a marriage of the two concepts: we preserve an Eulerian framework by constructing a grid that keeps the analysis of the solution spread evenly throughout a region of interest; however, we also examine the parcels that pass through these grid points, tracing their trajectory and using that information to update the solution at later times. By considering the system of ordinary differential equations (ODEs) along the characteristic lines of the PDE, we can take advantage of numerical methods that solve ODEs [7, 8].

Figure 1 demonstrates the concept in one spatial dimension: the $t,x$-plane. Each $x_{A}^{k+1}$ refers to the grid point a particle will arrive on at time $t_{k+1}$. Each $x_{D}^{k}$ refers to the departure point for the particle at time $t_{k}$. Departure points and arrival points are paired as they represent the position of the same particle at two different times. We note that for a backward trace, the calculated departure point may not necessarily be on a grid intersection.

Early semi-Lagrangian theory began with the development around a finite difference framework [13]. To understand the elementary theory, consider the following advection equation in 1D,

[TABLE]

where

[TABLE]

determines our characteristic lines. By using a central difference approximation for $F$ about the point $(t_{k},x_{A}-\Delta x)$ along a characteristic line and plugging into Equation (3), we have the finite difference approach,

[TABLE]

Point $(t_{k}+\Delta t,x_{A})$ is a regular point on the Eulerian grid and $t_{k}-\Delta t=t_{k-1}$. However, due to varying values of $\Delta x$, we rarely find $(t_{k}-\Delta t,x_{A}-2\Delta x)$ to be a grid point. This requires us to use interpolation to determine values of $F$ at such a point. By using the above characteristic equation, we can approximate $\Delta x$ from the implicit formula,

[TABLE]

which may be solved via iteration.

We recall that most of the present semi-Lagrangian-based methods take advantage of second order error with respect to the temporal discretisation, $\Delta t=t_{k+1}-t_{k}$. Our proposed algorithms achieve not just first and second order error with respect to the temporal discretisation, but also third and fourth order error for two general classes of time-dependent partial differential systems.

The rest of this paper is organised into five more sections. Section 2 will deal with the construction and implementation of our semi-Lagrangian algorithms in the context of the 1D initial value problem from System (1). Section 3 will consider the stability and convergence of the methods from section 2. Section 4 will explain how these methods can be developed for higher dimensional initial value problems, such as System (2). Section 5 consists of numerical results for two test cases, one for System (1) and another for System (2). Section 6 contains our conclusions and considerations for future research.

2 Semi-Lagrangian Methods

Recall the general advection in one dimension problem from System (1). Many semi-Lagrangian methods are developed in [7, 8] to numerically solve a single PDE. We seek to adapt these methods to solve the system of interest. Since $\omega=\frac{dx}{dt}$, we have

[TABLE]

where $\tau=t_{k+1}-t_{k}=\frac{T}{N}$. Also, let $h=\frac{b-a}{M}$.

Let $\boldsymbol{y}=\boldsymbol{(}y_{1},\dots,y_{n}\boldsymbol{)}^{\intercal}$ and $\boldsymbol{f}=\boldsymbol{(}f_{1},\dots,f_{n}\boldsymbol{)}^{\intercal}$, then the reformulation of PDE System (1) is

[TABLE]

Since $\omega=\dfrac{dx}{dt}$, we also have, from the chain rule, $\dfrac{d}{dt}=\dfrac{\partial}{\partial t}+\omega\dfrac{\partial}{\partial x}$. This allows us to write System (1) in the compact form

[TABLE]

Integrating from $t_{k}$ to $t_{k+1}$, we get

[TABLE]

Let $\boldsymbol{y}_{A}^{k+1}$ be the numerical approximation of $\boldsymbol{y}(t_{k+1},x_{A}^{k+1})$ and $\boldsymbol{y}_{D}^{k}$ be the numerical approximation of $\boldsymbol{y}(t_{k},x_{D}^{k})$. We can now construct the semi-Lagrangian Euler method for finding $\boldsymbol{y}_{A}^{k+1}$ at all grid points $x_{A}^{k+1}$ given the values for $\boldsymbol{y}_{A}^{k}$ at all grid points $x_{A}^{k}$ as seen in Algorithm 1.

Lipscomb [10] shows that the semi-Lagrangian Euler method in Algorithm 1 has a local truncation error of $O(\tau)$. The integer $n$ is the number of implicit calculations used to approximate each departure point. While $n$ can vary, Williamson & Olson [15] and Simmons [12] have demonstrated in most cases that just a few iterations provide rapid convergence to the correct departure point, a common choice being $n=5$.

We will also consider the departure point calculations by case of $\omega$.

If $\omega\in\mathbb{R}$, then the calculation of the departure requires neither iteration nor interpolation. 2. 2.

If $\omega=\omega(t,x)$, then the calculation of the departure points is defined implicitly. We employ an iterative method; however, interpolation is not required to determine the departure points. 3. 3.

If $\omega=\omega(t,x,\boldsymbol{y})$, then the calculation of the departure points is again defined implicitly. However, $\omega_{D}^{k}$ is dependent on unknown values of $\boldsymbol{y}_{D}^{k}$ since the numerical solution is only known at the arrival (grid) points, $\boldsymbol{y}_{A}^{k}$. Therefore, an interpolation scheme is employed. Linear interpolation is sufficient for the semi-Lagrangian Euler method.

The modified semi-Lagrangian approach is developed from a Runge-Kutta Order-2 method for solving systems of ODEs, much like the semi-Lagrangian Euler method was based on Euler’s method for solving systems of ODEs. The development of Runge-Kutta methods can be found in many numerical analysis texts such as Burden and Faires [1]. The modified semi-Lagrangian Euler method is shown in Algorithm 2.

We notice that the algorithm utilises a higher order recalculation of the departure points. In the first and second order methods, we used an Euler approximation to numerically integrate along the characteristic lines. Algorithm (2), on the other hand, uses a trapezoid rule approximation:

[TABLE]

where, under case 3 assumptions for $\omega$,

[TABLE]

and

[TABLE]

In order to determine $\omega_{A}^{k+1}$, we may need to evaluate $\boldsymbol{y}_{A}^{k+1}$ at all the arrival points. This is dependent on which case of $\omega$ we have to deal with. We will consider what elements of the algorithm can be omitted by case for $\omega$:

If $\omega\in\mathbb{R}$, then the calculation of the departure points reduces to a simple subtraction formula. Neither iteration nor interpolation is required for computing $x_{D}^{k}$. 2. 2.

If $\omega=\omega(t,x)$, then $\omega_{A}^{k+1}$ is not dependent on $\boldsymbol{y}_{A}^{k+1}$. Only implicit iteration will be required to compute $x_{D}^{k}$. 3. 3.

If $\omega=\omega(t,x,\boldsymbol{y})$, then $\omega_{A}^{k+1}$ is dependent on $\boldsymbol{y}_{A}^{k+1}$. Both implicit iteration and interpolation are required for calculating departure points. Second-order or higher interpolation is used for the modified semi-Lagrangian Euler method.

The modified semi-Lagrangian Euler method has a local truncation error of $O(\tau^{2})$ [10].

Much like the modified semi-Lagrangian Euler method, we can derive third and fourth order Runge-Kutta methods for semi-Lagrangian schemes. The semi-Lagrangian Runge-Kutta method of order-3 is shown in Algorithm 3.

Algorithm 3 uses an even higher order recalculation of the departure points. In the first and second order methods, we used left-side-rule and trapezoid rule approximations to numerically integrate along the characteristic lines. Algorithm (3), on the other hand, uses a Simpson’s rule approximation:

[TABLE]

where, under case 3 assumptions for $\omega$,

[TABLE]

[TABLE]

and

[TABLE]

In order to determine $\omega_{A}^{k+1}$, we may need to evaluate $\boldsymbol{y}_{A}^{k+1}$ at all the arrival points; further, a half-way Runge-Kutta computation is required for $\tilde{\omega}_{I}$. This is dependent on which case of $\omega$ we have to deal with. We will consider what elements of the algorithm can be omitted by case for $\omega$:

If $\omega\in\mathbb{R}$, then calculation of the departure points reduces to simple subtraction. Further, neither iteration nor interpolation is required for computing $x_{D}^{k}$. 2. 2.

If $\omega=\omega(t,x)$, then $\omega_{A}^{k+1}$ is not dependent on $\boldsymbol{y}_{A}^{k+1}$. Also, $\tilde{\omega}_{I}$ is not dependent on $\tilde{\boldsymbol{y}}_{I}$. Implicit iteration will be required to compute $x_{D}^{k}$. 3. 3.

If $\omega=\omega(t,x,\boldsymbol{y})$, then $\omega_{A}^{k+1}$ is dependent on $\boldsymbol{y}_{A}^{k+1}$. We also need to find $\tilde{\boldsymbol{y}}_{I}$ in order to determine $\tilde{\omega}_{I}$. Both implicit iteration and interpolation are required for calculating departure points. Third-order or higher interpolation is used for the semi-Lagrangian Runge-Kutta order-3 method.

The semi-Lagrangian Runge-Kutta order-3 method has a local truncation error of $O(\tau^{3})$ given that we use at least third-order interpolations within the method [10].

We will now introduced the semi-Lagrangian Runge-Kutta method of order-4 as seen in Algorithm 4.

Once again, we are using a Simpson’s rule approximation; therefore, an intermediate evaluation of $\tilde{\omega}_{I}$ is often needed. Depending on the case of $\omega$, this may require a half step-size evaluation of the solution, $\tilde{\boldsymbol{y}}_{I}$. The algorithmic alterations by case of $\omega$ are similar to the semi-Lagrangian order-3 method in Algorithm 3 with the exception that we use at least fourth order interpolation.

The semi-Lagrangian Runge-Kutta order-4 method has a local truncation error of $O(\tau^{4})$ given that we use at least fourth order interpolations [10].

3 Stability and Convergence

An important consideration when developing numerical algorithms for solving differential equations is whether the algorithm is numerically stable, convergent, and whether an upper bound on the error can be determined. We will consider this for the initial value problem posed by System (1). Before we introduce our stability and convergence theorem, we need a few lemmas. The first two lemmas and their proofs can be found in Burden and Faires [1].

Lemma 1**.**

For all $t\geq-1$ and any positive $m$, we have $0\leq(1+t)^{m}\leq e^{mt}$.

Lemma 2**.**

If $s$ and $t$ are positive real numbers and $\{a_{k}\}_{k=0}^{n}$ is a sequence satisfying $a_{0}\geq\frac{-t}{s}$ and $a_{k+1}\leq(1+s)a_{k}+t$, for each $k=0,1,\dots,n-1$, then

[TABLE]

We will also need an interpolation lemma from Guo [8].

Lemma 3**.**

Suppose $x_{i}=a+ih$ for $i=0,1,\dots,M$ and $h=(b-a)/M$. Then, for each $x\in[a,b]$, if $p(x)$ and $q(x)$ are two piece-wise linear interpolations with $p(x_{i})=u_{i}$ and $q(x_{i})=v_{i}$ for $i=0,1,\dots,M,$ then

[TABLE]

Last, we need a theorem regarding the existence and uniqueness of solutions to initial value problems in the form of System (1).

Theorem 1**.**

If $\boldsymbol{f}(t,x,\boldsymbol{y})$ from (1) is a continuous function of $t$ and $x$ and satisfies a Lipschitz condition in $\boldsymbol{y}$ for $(t,x)\in[0,T]\times[a,b]$ with $\|\boldsymbol{y}\|_{\infty}<\infty$, then there exists a unique differentiable solution, $\boldsymbol{y}(t,x)$ for the initial value problem in System (1).

A proof of Theorem 1 can be derived from most textbooks in partial differential equations. We will now move on to our stability and convergence theorem which applies to the initial value problem in System (1). It is a generalisation of the theorem in Guo [8].

Theorem 2**.**

Let the initial value problem in (5) be approximated by the one-step difference method

[TABLE]

If there exists $\tau_{0},\alpha>0$ and $\boldsymbol{\phi}(t,x,\boldsymbol{y},\tau,h)$ is continuous and satisfies a Lipschitz condition in the variable $\boldsymbol{y}$ with Lipschitz constant $L$ on

[TABLE]

Then,

The one-step method in (6) is stable. 2. 2.

If $\boldsymbol{\phi}$ additionally satisfies a Lipschitz condition in all of its variables, then the one-step method in (6) is convergent if and only if it is consistent; that is $\boldsymbol{\phi}(t,x,\boldsymbol{y},0,0)=\boldsymbol{f}(t,x,\boldsymbol{y})$. 3. 3.

If a function $\varepsilon(\tau)$ exists and, for each $k=0,1,\dots,N$, the local truncation error $\boldsymbol{\varepsilon}_{k}(\tau)$ satisfies $\|\boldsymbol{\varepsilon}_{k}(\tau)\|_{\infty}\leq\varepsilon(\tau)$ for $\tau\in[0,\tau_{0}]$, then

[TABLE]

for $k=1,2,\dots,N$.

Proof.

(1.) Let $\{\boldsymbol{u}_{A}^{k}\}_{k=1}^{N}$ and $\{\boldsymbol{v}_{A}^{k}\}_{k=1}^{N}$ each satisfy the difference equation in (6) with initial conditions $\boldsymbol{u}_{0}$ and $\boldsymbol{v}_{0}$ respectively. Let

[TABLE]

for point $(t_{k},x_{A}^{k})$ and

[TABLE]

for point $(t_{k},x_{D}^{k})$. Also, let

[TABLE]

on the regular grid points (arrival points). Now,

[TABLE]

from the difference formula. So,

[TABLE]

from our Lipschitz condition of $\boldsymbol{\phi}$. Thus,

[TABLE]

By Lemma 3, we have

[TABLE]

and

[TABLE]

From Lemmas 1 and 2,

[TABLE]

where $k\leq\frac{T}{\tau}=N$. We notice that $e^{TL}$ is our stability constant. So, the method in (6) is stable. ∎

Proof.

(2.) Let $\boldsymbol{\phi}(t,x,\boldsymbol{y},0,0)=\boldsymbol{g}(t,x,\boldsymbol{y})$. $\boldsymbol{g}$ satisfies the conditions in Theorem 1; so,

[TABLE]

has a unique, differentiable solution $\boldsymbol{v}(t,x)$. The numerical solution $\boldsymbol{z}$ satisfies

[TABLE]

By the Mean Value Theorem,

[TABLE]

for some $\xi,\eta\in[0,1]$. Let

[TABLE]

[TABLE]

on the regular grid points (arrival points), and

[TABLE]

Then,

[TABLE]

By the assumption of $\boldsymbol{\phi}$ satisfying Lipschitz in all of its variables, with Lipschitz constant $L$ that is sufficiently large enough to satisfy the Lipschitz conditions in all the variables, we have

[TABLE]

[TABLE]

[TABLE]

and

[TABLE]

where $L_{1}$ is a positive constant such that

[TABLE]

for all $(t,x)\in[0,T]\times[a,b]$. So,

[TABLE]

and

[TABLE]

Using Lemma 2,

[TABLE]

Clearly, the right hand side of the inequality goes to zero as $\tau,h\to 0$ with $E^{0}=0$. So, $\boldsymbol{z}^{k}\to\boldsymbol{v}^{k}$; that is, our numerical solution converges to the solution of System (5). Thus, given $\boldsymbol{\phi}(t,x,\boldsymbol{y},0,0)=\boldsymbol{g}(t,x,\boldsymbol{y})$, the one-step difference method in Equation (6) converges.

Now we will assume convergence of the method in Equation (6). $\boldsymbol{v}$, our unique solution to System (7), matches $\boldsymbol{y}$, our unique solution to System (5). Let $\boldsymbol{\phi}$ and $\boldsymbol{g}$ differ at some point; now consider the initial value problem starting at that point. Obviously, the two solutions are different, which leads to a contradiction. Thus, convergence requires $\boldsymbol{\phi}(t,x,\boldsymbol{y},0,0)=\boldsymbol{g}(t,x,\boldsymbol{y})$, the condition of consistency. ∎

Proof.

(3.) Let $\boldsymbol{E}_{A}^{k}=\boldsymbol{y}_{A}^{k}-\boldsymbol{y}(t_{k},x_{A}^{k})$ and $E^{k}=\max_{x_{A}}\|\boldsymbol{E}_{A}^{k}\|_{\infty}$. From the definition of local truncation error,

[TABLE]

Subtracting this from the difference method yields

[TABLE]

Then,

[TABLE]

and

[TABLE]

From Lemma 2,

[TABLE]

for $k=1,2,\dots,N$. ∎

This theorem can be adapted for the different methods by determining the correct form of $\boldsymbol{\phi}(t_{k},x_{D}^{k},\boldsymbol{y}_{D}^{k},\tau,h)$, higher order interpolations, and also considering the higher order approximations of the departure points.

4 Higher Dimensions

We have also developed semi-Lagrangian methods for solving initial value problems in the form of nonlinear advection in two dimensions seen in System (1). For example, the semi-Lagrangian Euler method in two dimensions is seen in Algorithm 5.

$\boldsymbol{u}_{A}^{k}$ is the numerical approximation of $\boldsymbol{u}(t_{k},x_{A},y_{A})=\boldsymbol{(}u(t_{k},x_{A},y_{A}),v(t_{k},x_{A},y_{A})\boldsymbol{)}^{\intercal}$. This method has a local truncation error of $O(\tau)$ for $h_{x},h_{y}=O(\tau)$, $h_{x}=\frac{(b-a)}{M_{x}}$, $h_{y}=\frac{d-c}{M_{y}}$. We see that the method is similar to the one dimensional case. However, we have no longer have a generalised advection term and we must consider departure points with respect to an extra dimension. The 2nd–4th order methods are similarly constructed; further, similar results to Theorem 2 also hold for these algorithms as seen in Lipscomb [10].

5 Numerical Results

We will now consider some numerical results. In order to do so, we require a measure of error.

Definition 1**.**

Let $y_{\boldsymbol{x}}^{t}$ be the value of a numerical solution to a PDE at point $(t,\boldsymbol{x})$ on an Eulerian grid. Also, let $y(t,\boldsymbol{x})$ be the value of the exact solution to the same PDE at the same point. We define the residual at point $(t,\boldsymbol{x})$ as

[TABLE]

We further define the maximum residual at time $t$ as

[TABLE]

Clearly, residuals capture the absolute error between the exact solution and numerical solution at specified points in time and space. We will use the maximum residual at time $t$ to create residual plots against time.

Our first initial value problem is in the form of general advection in one dimension (1). Specifically, we will work with $\omega(t,x,u,v)=u+v$, where $u=u(t,x)$ and $v=v(t,x)$ are the solutions to the system:

[TABLE]

on the domain $(t,x)\in[0,1]\times[0,1]$ with initial conditions

[TABLE]

We can verify that the exact solutions are

[TABLE]

In this case, the nonlinear system is specified strictly in terms of the unknown functions and their partial derivatives. We will use 5 iterations for calculating the departure points implicitly. However, we also note that interpolation is required for determining $\omega_{D}^{k}$ for each iteration.

The numerical solutions from the semi-Lagrangian Runge-Kutta Order-4 method can be seen in Figure 2 using 50 spatial steps and 50 time steps; that is, $\tau=h=0.02$. The solutions look very similar as they are simply translations of each other. However, they have different starting points and a full period is completed for both from $t=0$ to $t=1$.

Examining the semilog max residual plots in Figure 3, it is clear we have achieved very strong results for the higher order methods. For all four methods, we chose $h=\tau=0.005$.

We will now consider a nonlinear advection in two dimensions (2) problem. Consider

[TABLE]

on the domain $(t,x,y)\in[0,1]\times[0,1]\times[0,1]$ with initial conditions

[TABLE]

We can verify that the exact solutions are

[TABLE]

The final time evolution of the numerical solutions from the fourth order method can be seen in Figure 4. The step-sizes used were $\tau=h_{x}=h_{y}=0.02$. The initial conditions for both solutions are sine and cosine-based formations whose magnitude increases to what is seen in the picture.

Figure 5 presents the max semilog residual plots for comparative analysis. The spatial step-size was $h_{x}=h_{y}=0.04$ with a time step-size of $\tau=0.02$. Clearly, the fourth order method is superior, followed by the third order, second order, then first order.

We will now provide some numerical confirmation of the order of these algorithms by examining the order of the absolute errors at $t=1$. The best way to examine the order of error is to produce a $\log$-$\log$ plot of the max residuals at $t=1$ as a function of the time step-size $\tau=\Delta t$. This will allow us to observe the slope to determine the order of the absolute error. Figure 6 presents these plots for the first four algorithms in this paper; clearly the results match the order of each algorithm.

6 Conclusion and Future Development

We have constructed several useful algorithms for numerically solving two general systems of time-dependent partial differential equations while examining their convergence, stability, and error. Using numerical simulations for several examples, we have established the effectiveness of the higher order numerical algorithms: the semi-Lagrangian Runge-Kutta order-3 and order-4 methods. The two general systems dealt with two different types of cases:

PDE System 1 is an initial value problem in one spatial dimension; however, the advection term $\omega$ can be generalised to $\omega=\omega(t,x,\boldsymbol{y})$ where $\boldsymbol{y}$ is the unknown solution to the initial value problem. 2. 2.

PDE System 2 is an initial value problem in two spatial dimensions. In this case, we have standard nonlinear advection terms: $u=\dfrac{dx}{dt}$ and $v=\dfrac{dy}{dt}$.

By examining the construction of these algorithms, it is clear that these ideas can be built upon in order to solve PDE initial value problems in $n$ spatial dimensions as long as there are $n$ advection terms. These advection terms can be generalised as in PDE System 1; however, the performance of semi-Lagrangian methods in these more complicated cases will have to be tested.