# General Formulation of Second-Order Semi-Lagrangian Methods for Convection-Diffusion Problems

Authors: Xiaohan Long; Chuanjun Chen
Journal: Abstract and Applied Analysis
 (2013)

---

## Abstract

The general formulation of the second-order semi-Lagrangian methods was presented for convection-dominated diffusion problems. In view of the method of lines, this formulation is in a sufficiently general fashion as to include two-step backward difference formula and Crank-Nicolson type semi-Lagrangian schemes as particular ones. And it is easy to be extended to higher-order schemes. We show that it maintains second-order accuracy even if the involved numerical characteristic lines are first-order accurate. The relationship between semi-Lagrangian methods and the modified method of characteristic is also addressed. Finally convergence properties of the semi-Lagrangian finite difference schemes are tested.


---

## Body

## 1. Introduction

There have been many numerical methods developed to deal with convection-dominated diffusion problems, and among them characteristic-based methods are the most popular. Numerical methods that follow characteristic lines backwards in time and then interpolate at their feet date back to the work of Courant et al. [1]. As one kind of characteristics-based method, the semi-Lagrangian (SL) method was introduced in the beginning of the 1980s by Robert [2]. Its basic idea is to discretize the Lagrangian derivative of the solution instead of the Eulerian derivative. This technique can increase significantly the maximum allowable time step while maintaining the efficiency of symmetric solvers. The SL methods have been extensively applied in numerical simulations of models for weather forecast and oceanography (see Staniforth and Côté [3] and Smolarkiewicz and Pudykiewicz [4] for review). As another kind of characteristic-based methods, the modified method of characteristic (MMOC) was introduced by Douglas Jr. and Russell [5] at roughly the same time as the SL method and has been extensively implemented in numerical simulations of fluid flows in porous media (see [6–8]) and many other transport problems (see [9–12] for review).Strang [13] pointed out that the first-order methods were often too crude and the third-order methods too complicated. The computations are thus made expensive either by the fine mesh required by a first order scheme in order to provide enough detail or else by the delicate differencing which maintains a high-order accuracy. Second-order schemes are the obvious compromise.Based on MMOC, in the beginning of 1980s Ewing and Russell [14] introduced the backward difference formula (BDF) of characteristic for linear constant-coefficient convection-diffusion problems. Afterwards characteristic schemes of Crank-Nicolson (CN) type for convection-diffusion equations and the Navier-Stokes (NS) equations were studied by Rui and Tabata [15] and Notsu and Tabata [16], respectively. The SL-BDF schemes for incompressible NS problems were proposed by Boukir et al. [17, 18]. The SL-CN methods were proposed for convection-diffusion equations and/or the incompressible NS equations by Bermúdez et al. [19, 20], Fourestey and Piperno [21], Xiu and Karniadakis [22], and Xiu et al. [23]. Al-Lawatia et al. [24] and Falcone and Ferretti [25] presented and analyzed the single-step high-order semi-Lagrangian schemes of the Runge-Kutta type. More research was given by Bermejo and Conde [26], Xiao and Yabe [27], and Toda et al. [28].Up to now we have not seen the SL method to be formulated in a sufficiently general fashion. In this paper we present a general second-order SL formulation for convection-dominated diffusion problems. By the method of lines approach, this formulation includes most of the second-order SL schemes mentioned previously. Extension of it to higher-order schemes is also addressed. Our development of the general SL formulation is mainly motivated by the work of Rui and Tabata [15] and Notsu and Tabata [16]. They treated convection-diffusion and the NS problem, respectively, by MMOC. They emphasized that an “additional correction term” was indispensable to maintain the second-order accuracy of the MMOC schemes. We will see that the correction term of the MMOC schemes is actually a natural term of the SL schemes.First, we show that by the method of lines (MOL) approach, the general SL formulation includes the SL-BDF2 and the SL-CN schemes as specific ones. Second, we verify that the formulation maintains second-order accuracy even if the involved numerical characteristic lines are first-order accurate. Third, we show that MMOC can be considered as a special version of the SL method. Finally, combining finite difference discretization in spaces, a fully discretized SL scheme is presented. Numerical tests demonstrate that the SL finite difference scheme is second-order convergent.The outline of the rest of this paper is as follows. In Section2 the forming of the SL schemes is recalled and the relationship between the first-order SL methods and MMOC is addressed. In Section 3 the general SL formulation is derived and its second-order consistency is proved. The relation of the second-order SL method and MMOC is also addressed. In Section 4 the finite difference SL schemes are derived. In Section 5 several SL finite difference schemes are applied to the Gaussian hill problem. In Section 6 conclusions of this paper are drawn.


## 2. The SL Methods and the MMOC

In this section we recall some details of the construction of the SL method and MMOC. Also we study the relationship between them.LetΩ be a bounded domain in ℝ2 with Lipschitz boundary ∂Ω, and let T be a positive constant. Without loss of generality, we consider the initial boundary value problem: find ϕ:Ω×(0,T)→ℝ such that(1a)∂ϕ∂t+u·∇ϕ-νΔϕ=finΩ×(0,T),(1b)ϕ(x,t)=0on∂Ω×(0,T),(1c)ϕ(x,0)=ϕ0(x)inΩ,where ν is a positive constant, u:Ω×(0,T)→ℝ2, and f:Ω×(0,T)→ℝ are given functions. For simplicity, we assume that u is a divergence-free velocity field, that is, ∇·u=0, and vanishes on the boundary ∂Ω.Characteristic line is the trajectory of a fluid particle. The travel of the particle is associated with the velocity fieldu. For a given point (x,s)∈Ω-×[0,T], the characteristic line through (x,s) is represented by the vector function X(x,s;t), which solves the initial value problem(2a)dX(x,s;t)dt=u(X(x,s;t),t),(2b)X(x,s;s)=x.Here X(x,s;t) is the position of the particle on the characteristic line at time t. The particle is located at x at time s. We assume that u∈C0(Ω-×[0,T]) is Lipschitz continuous on ∂Ω with respect to the first variable. By the theory of ODEs, the characteristic line is well defined. By the chain rule we have
(3)dϕdt(X(x,s;t),t)=∂ϕ∂t(X(x,s;t),t)+u(X(x,s;t),t)·∇ϕ(·,t)∘X(x,s;t),
where ϕ(·,t)∘X(x,s;t) is the composition of functions with respect to the first argument of ϕ. Similar to the deriving of (4.1) and (4.2) in [19], (1a) can be written to Lagrangian form
(4)dϕdt(X(x,s;t),t)-νΔϕ(·,t)∘X(x,s;t)=f(X(x,s;t),t).For a positive integerN, let Δt=T/N be time step length and tn=nΔt, for n=1,2,…,N. Let X(x,tn+1;t) denote the characteristic line on [tn,tn+1] (or [tn-1,tn+1] for two-step methods) through (x,tn+1). Thus (4) can locally be written to
(5)dϕdt(X(x,tn+1;t),t)-νΔϕ(·,t)∘X(x,tn+1;t)=f(X(x,tn+1;t),t).Tracking the particle backward fromx to X(x,tn+1;tn) along the characteristic line on [tn,tn+1] (or from x to X(x,tn+1;tn-1) on [tn-1,tn+1] for two-step methods) corresponds to the backward solution of the following Cauchy problem:
(6)dX(x,tn+1;t)dt=u(X(x,tn+1;t),t),X(x,tn+1;tn+1)=x.
For function w:Ω×[tn,tn+1]→ℝ, denote wi=w(x,ti) and Xi=X(x,tn+1,ti) for i=n,n+1. By applying the backward Euler’s method to (5), it follows that
(7)φn+1-φn(Xn)Δt-νΔφn+1=fn+1,
where φn represents the approximate of ϕn. Note that in practice, (6) can usually be solved approximately. By applying the backward Euler’s method to (6), it follows that
(8)x-XEnΔt=un+1(x),
and XEn=x-un+1Δt (similarly we see that XEn-1=x-2un+1Δt for the two-step methods). In (7) with Xn being replaced by XEn, we have
(9)φn+1-φn(x-un+1Δt)Δt-νΔφn+1=fn+1.On the other hand, (9) can be derived by MMOC [5]. In fact, with u=(u1,u2), let s denote the direction vector (1,u1,u2), and define the operator
(10)dds:=1θ(∂∂t+u·∇),
with θ(x,t):=[1+|u(x,t)|2]1/2=[1+|u1(x,t)|2+|u2(x,t)|2]1/2. So (1a) can be written to the form
(11)θdϕds-νΔϕ=f.
Using the backward difference quotient, we have
(12)dϕn+1ds=ϕn+1-ϕn(x-Δtun+1/θn+1)Δt+Δt2d2ϕn+1d2s+O(Δt2).
Substituting (12) into (11), we obtain (9). Thus MMOC can be considered as a special version of SL. In the next section we will see that similar relation exists for second-order case.


## 3. General Second-Order SL Formulation

In this section we present a general second-order SL formulation and show that by MOL this formulation includes several widely used schemes as specific ones.Let us first consider the abstract ODEs of the following form: given the Hilbert spaceℋ and y0∈ℋ, find y:(0,T]→ℋ such that
(13)y′(t)=g(t,y(t)),t∈(0,T],y(0)=y0,
where g:(0,T]×ℋ→ℋ. A general second-order scheme for (13) can be of the form (see [29, 30])
(14)∑ℓ=02‍αℓYn+ℓ-1=Δt∑ℓ=02‍βℓgn+ℓ-1(Yn+ℓ-1),n=1,2,…,N-1.
As usual we denote scheme (14) as (ρ,σ), where ρ and σ are the characteristic polynomials of the scheme, with
(15)ρ(ζ):=∑ℓ=02αℓζℓ,σ(ζ):=∑ℓ=02βℓζℓ.
We further assume that (ρ,σ) is normalized by ∑ℓ=02‍βℓ=1 and satisfies the following conditions:
(16)α1=1-2α2,α0=-1+α2,β1=12+α2-2β2,β0=12-α2+β2,
for α2≥1/2, β2≥α2/2, such that (ρ,σ) is Dahlquist and A-stable [31].For a fixedx∈Ω and t∈[tn-1,tn+1], let y(t)=ϕ(X(x,tn+1;t),t). Then the (ρ,σ) scheme for (5) is of the form
(17)∑ℓ=02αℓφn+ℓ-1(Xn+ℓ-1)=Δt∑ℓ=02βℓ[νΔφn+ℓ-1∘Xn+ℓ-1+fn+ℓ-1(Xn+ℓ-1)],n=1,2,…,N-1.
With X being replaced by the approximate characteristic line X¯, we have the analogue of (17):
(18)∑ℓ=02‍αℓφn+ℓ-1(X¯n+ℓ-1)=Δt∑ℓ=02βℓ[νΔφn+ℓ-1∘X¯n+ℓ-1+fn+ℓ-1(X¯n+ℓ-1)],n=1,2,…,N-1.
Scheme (18) is the general second-order SL formulation which by MOL includes all the previously introduced second-order SL schemes. Next we will prove that it includes the SL-BDF2 schemes and the SL-CN schemes as specific cases.


### 3.1. The BDF2 Scheme

In (14), let β0=β1=0; then from (16), we have α2=3/2, α1=-2, α0=1/2. Substituting these coefficients into (14), we obtain the BDF2 scheme:
(19)3Yn+1-4Yn+Yn-12Δt=gn+1(Yn+1).
Due to its stability and damping properties, (19) is one of the most popular second-order schemes [29]. Substituting the previously coefficients into (17), we get the SL analogue of (19):
(20)3φn+1-4φn(Xn)+φn-1(Xn-1)2Δt-νΔφn+1=fn+1.
Furthermore, in (20) with Xn and Xn-1 being, respectively, approximated by XEn=x-un+1Δt, XEn-1=x-2un+1Δt, we have
(21)3φn+1-4φn(x-un+1Δt)+φn-1(x-2un+1Δt)2Δt-νΔφn+1=fn+1.
This is just the multistep characteristic scheme derived based on MMOC in [14]. In fact, using MMOC (analogous to (12)) we see that
(22)dϕn+1ds=(3ϕn+1-4ϕn(x-Δtun+1/θn+1)+ϕn-1(x-2Δtun+1/θn+1))(2Δt)-1+(Δt)23d3ϕn+1ds3+O(Δt3).
Substituting (22) into (11), we obtain (21).The SL-BDF2 schemes were presented and analyzed in [14, 17, 18, 21]. It is deserved to note that though the involved approximate characteristic line is first-order accurate, the resulting SL-BDF2 scheme (21) maintains second-order accurate.


### 3.2. The CN Scheme

In (14), let
(23)α2=1,α1=-1,α0=0;β2=12,β1=β0=0.
From (14) we obtain the CN scheme
(24)Yn+1-YnΔt=12(gn(Yn)+gn+1(Yn+1)).
From (17) we get the SL analogue of (24):
(25)φn+1-φn(Xn)Δt-12νΔ(φn∘Xn+φn+1)=12(fn+1+fn(Xn)),
and with Xn being replaced by XE, we have
(26)φn+1-φn(x-un+1Δt)Δt-12νΔ(φn∘(x-un+1Δt)+φn+1)=12(fn(x-un+1Δt)+fn+1).Remark 1.
It is deserved to note that Rui and Tabata [15] and Notsu and Tabata [16] called (1/2)νΔφn∘(x-un+1Δt)-term in (26) the “additional corrected term,” since it is introduced from “outside” to recover the second-order consistency of the MMOC schemes. But in view of the previous discussion, the “additional corrected term” is a natural one in the SL schemes. Thus we think the SL method is more general than MMOC.


### 3.3. The Consistency of the General Formulation

Now we show that both (17) and (18) with first-order approximate characteristic lines are second-order accurate. Let us denote
(27)ℒ[w(t);Δt]:=∑ℓ=02[αℓwn+ℓ-1-Δtβℓgn+ℓ-1(wn+ℓ-1)],
for w:(0,T]→ℋ. By the theory of ODEs (see [30]), using the Taylor’s expansion, from (13) and (14) we have
(28)ℒ[y(t);Δt]=(α0+α1+α2)y(t)+Δt[(α1+2α2)-(β0+β1+β2)]y′(t)+Δt2[12(α1+4α2)-(β1+2β2)]y′′(t)+O(Δt3).
If the following conditions hold
(29)α0+α1+α2=0,(α1+2α2)-(β0+β1+β2)=0,12(α1+4α2)-(β1+2β2)=0,
then ℒ[y(t);Δt]=O(Δt3). It is easy to see that conditions (29) and (16) are equivalent if ∑ℓ=02‍βℓ=1. In (13), let
(30)y(t)=φ(X(·,tn+1;t),t),g(y(t),t)=νΔφ(·,t)∘X(·,tn+1;t)+f(X(·,tn+1;t),t).
From (28)–(30) it follows that
(31)ℒ[ϕ(X(·,tn+1;t),t);Δt]=O(Δt3),
where ϕ is the exact solution of (5). Thus we have confirmed the following proposition.Proposition 2.
Suppose thatu and f in (1a)–(1c) are smooth functions in Ω-×(0,T), and X(t) is the exact characteristic line which solves (6). Then (17) is second-order consistent with (5).Analogous to Proposition2, if first-order characteristic lines are involved, then the following proposition holds.Proposition 3.
Suppose thatu and f in (1a)–(1c) are smooth functions in Ω-×(0,T). In (18), if X¯n=XEn, X¯n-1=XEn-1, then the scheme
(32)∑ℓ=02αℓφn+ℓ-1(XEn+ℓ-1)-νΔt∑ℓ=02βℓΔφn+ℓ-1∘XEn+ℓ-1=Δt∑ℓ=02βℓfn+ℓ-1(XEn+ℓ-1)
is second-order consistent with (5).Proof.
Analogous to the deriving of(4.1) and (4.2) in [19], for t∈[tn-1,tn+1] (5) can be written to the form
(33)∂ϕ(·,t)∂t∘X(x,tn+1;t)+u(X(x,tn+1;t),t)·∇ϕ(·,t)∘X(x,tn+1;t)-νΔϕ(·,t)∘X(x,tn+1;t)=f(X(x,tn+1;t),t).
Noting that
(34)X(x,tn+1;t)=x-Δtun+1(x)+O(Δt2),X(x,tn+1;t)=x-2Δtun+1(x)+O(Δt2),
we denote
(35)XE(x,tn+1;t)=x-Δtun+1(x),XE(x,tn+1;t)=x-2Δtun+1(x).
With X(x,tn+1;t) being replaced by XE(x,tn+1;t) and φ being the approximate of ϕ, we obtain O(Δt2)-perturbation of (33) of the form
(36)∂φ(·,t)∂t∘XE(x,tn+1;t)+u(XE(x,tn+1;t),t)·∇φ(·,t)∘X(x,tn+1;t)-νΔφ(·,t)∘XE(x,tn+1;t)=f(XE(x,tn+1;t),t).
Rewrite (36) to the Lagrangian form
(37)dφ(·,t)dt(XE(x,tn+1;t),t)-νΔφ(·,t)∘XE(x,tn+1;t)=f(XE(x,tn+1;t),t).
Thus (37) is the second-order approximation to (5).
Corresponding to (30), let
(38)y(t)=φ(XE(·,tn+1;t),t),g(y(t),t)=νΔφ(·,t)∘XE(x,tn+1;t)+f(XE(·,tn+1;t),t).
Similar to the proof of Proposition 2, we can see that (32) is second-order consistent with (37). We deduce that (32) is second-order consistent with (5).Remark 4.
In (32) higher-order numerical characteristic lines are usually preferred, though the characteristic line computed by the backward Euler’s method can retain the second-order accuracy. For example, Bermúdez et al. [19, 20], Rui and Tabata [15], and Notsu and Tabata [16] computed the numerical characteristic by the higher-order Runge-Kutta methods. In the following numerical tests we will see that the first-order characteristic line is too coarse to ensure reasonable convergence of the SL-CN schemes.


## 3.1. The BDF2 Scheme

In (14), let β0=β1=0; then from (16), we have α2=3/2, α1=-2, α0=1/2. Substituting these coefficients into (14), we obtain the BDF2 scheme:
(19)3Yn+1-4Yn+Yn-12Δt=gn+1(Yn+1).
Due to its stability and damping properties, (19) is one of the most popular second-order schemes [29]. Substituting the previously coefficients into (17), we get the SL analogue of (19):
(20)3φn+1-4φn(Xn)+φn-1(Xn-1)2Δt-νΔφn+1=fn+1.
Furthermore, in (20) with Xn and Xn-1 being, respectively, approximated by XEn=x-un+1Δt, XEn-1=x-2un+1Δt, we have
(21)3φn+1-4φn(x-un+1Δt)+φn-1(x-2un+1Δt)2Δt-νΔφn+1=fn+1.
This is just the multistep characteristic scheme derived based on MMOC in [14]. In fact, using MMOC (analogous to (12)) we see that
(22)dϕn+1ds=(3ϕn+1-4ϕn(x-Δtun+1/θn+1)+ϕn-1(x-2Δtun+1/θn+1))(2Δt)-1+(Δt)23d3ϕn+1ds3+O(Δt3).
Substituting (22) into (11), we obtain (21).The SL-BDF2 schemes were presented and analyzed in [14, 17, 18, 21]. It is deserved to note that though the involved approximate characteristic line is first-order accurate, the resulting SL-BDF2 scheme (21) maintains second-order accurate.


## 3.2. The CN Scheme

In (14), let
(23)α2=1,α1=-1,α0=0;β2=12,β1=β0=0.
From (14) we obtain the CN scheme
(24)Yn+1-YnΔt=12(gn(Yn)+gn+1(Yn+1)).
From (17) we get the SL analogue of (24):
(25)φn+1-φn(Xn)Δt-12νΔ(φn∘Xn+φn+1)=12(fn+1+fn(Xn)),
and with Xn being replaced by XE, we have
(26)φn+1-φn(x-un+1Δt)Δt-12νΔ(φn∘(x-un+1Δt)+φn+1)=12(fn(x-un+1Δt)+fn+1).Remark 1.
It is deserved to note that Rui and Tabata [15] and Notsu and Tabata [16] called (1/2)νΔφn∘(x-un+1Δt)-term in (26) the “additional corrected term,” since it is introduced from “outside” to recover the second-order consistency of the MMOC schemes. But in view of the previous discussion, the “additional corrected term” is a natural one in the SL schemes. Thus we think the SL method is more general than MMOC.


## 3.3. The Consistency of the General Formulation

Now we show that both (17) and (18) with first-order approximate characteristic lines are second-order accurate. Let us denote
(27)ℒ[w(t);Δt]:=∑ℓ=02[αℓwn+ℓ-1-Δtβℓgn+ℓ-1(wn+ℓ-1)],
for w:(0,T]→ℋ. By the theory of ODEs (see [30]), using the Taylor’s expansion, from (13) and (14) we have
(28)ℒ[y(t);Δt]=(α0+α1+α2)y(t)+Δt[(α1+2α2)-(β0+β1+β2)]y′(t)+Δt2[12(α1+4α2)-(β1+2β2)]y′′(t)+O(Δt3).
If the following conditions hold
(29)α0+α1+α2=0,(α1+2α2)-(β0+β1+β2)=0,12(α1+4α2)-(β1+2β2)=0,
then ℒ[y(t);Δt]=O(Δt3). It is easy to see that conditions (29) and (16) are equivalent if ∑ℓ=02‍βℓ=1. In (13), let
(30)y(t)=φ(X(·,tn+1;t),t),g(y(t),t)=νΔφ(·,t)∘X(·,tn+1;t)+f(X(·,tn+1;t),t).
From (28)–(30) it follows that
(31)ℒ[ϕ(X(·,tn+1;t),t);Δt]=O(Δt3),
where ϕ is the exact solution of (5). Thus we have confirmed the following proposition.Proposition 2.
Suppose thatu and f in (1a)–(1c) are smooth functions in Ω-×(0,T), and X(t) is the exact characteristic line which solves (6). Then (17) is second-order consistent with (5).Analogous to Proposition2, if first-order characteristic lines are involved, then the following proposition holds.Proposition 3.
Suppose thatu and f in (1a)–(1c) are smooth functions in Ω-×(0,T). In (18), if X¯n=XEn, X¯n-1=XEn-1, then the scheme
(32)∑ℓ=02αℓφn+ℓ-1(XEn+ℓ-1)-νΔt∑ℓ=02βℓΔφn+ℓ-1∘XEn+ℓ-1=Δt∑ℓ=02βℓfn+ℓ-1(XEn+ℓ-1)
is second-order consistent with (5).Proof.
Analogous to the deriving of(4.1) and (4.2) in [19], for t∈[tn-1,tn+1] (5) can be written to the form
(33)∂ϕ(·,t)∂t∘X(x,tn+1;t)+u(X(x,tn+1;t),t)·∇ϕ(·,t)∘X(x,tn+1;t)-νΔϕ(·,t)∘X(x,tn+1;t)=f(X(x,tn+1;t),t).
Noting that
(34)X(x,tn+1;t)=x-Δtun+1(x)+O(Δt2),X(x,tn+1;t)=x-2Δtun+1(x)+O(Δt2),
we denote
(35)XE(x,tn+1;t)=x-Δtun+1(x),XE(x,tn+1;t)=x-2Δtun+1(x).
With X(x,tn+1;t) being replaced by XE(x,tn+1;t) and φ being the approximate of ϕ, we obtain O(Δt2)-perturbation of (33) of the form
(36)∂φ(·,t)∂t∘XE(x,tn+1;t)+u(XE(x,tn+1;t),t)·∇φ(·,t)∘X(x,tn+1;t)-νΔφ(·,t)∘XE(x,tn+1;t)=f(XE(x,tn+1;t),t).
Rewrite (36) to the Lagrangian form
(37)dφ(·,t)dt(XE(x,tn+1;t),t)-νΔφ(·,t)∘XE(x,tn+1;t)=f(XE(x,tn+1;t),t).
Thus (37) is the second-order approximation to (5).
Corresponding to (30), let
(38)y(t)=φ(XE(·,tn+1;t),t),g(y(t),t)=νΔφ(·,t)∘XE(x,tn+1;t)+f(XE(·,tn+1;t),t).
Similar to the proof of Proposition 2, we can see that (32) is second-order consistent with (37). We deduce that (32) is second-order consistent with (5).Remark 4.
In (32) higher-order numerical characteristic lines are usually preferred, though the characteristic line computed by the backward Euler’s method can retain the second-order accuracy. For example, Bermúdez et al. [19, 20], Rui and Tabata [15], and Notsu and Tabata [16] computed the numerical characteristic by the higher-order Runge-Kutta methods. In the following numerical tests we will see that the first-order characteristic line is too coarse to ensure reasonable convergence of the SL-CN schemes.


## 4. The SL Finite Difference Method

In this section we present a full-discretized SL formulation which combines finite difference for spatial discretizations. We also numerically verify the convergence of the formulation.First, we build the finite difference scheme for (5). Assume that Ω is a unit square, with boundary ∂Ω. Denote ∂Ω=Γ1∪Γ2∪Γ3∪Γ4, with
(39)Γ1:={(x,y)∣0≤x≤1,y=0},Γ2:={(x,y)∣x=1,0≤y≤1},Γ3:={(x,y)∣0≤x≤1,y=1},Γ4:={(x,y)∣x=0,0≤y≤1}.
For the partition of Ω, we denote h:=M-1,xij:=(ih,jh) and w(xij,t):=wij(t). Let Ω-h:=ω-1×ω-2, with
(40)ω-1:={xi∣0≤i≤M},ω-2:={yj∣0≤j≤M}.
Let Ωh:=ω1×ω2, with
(41)ω1:={xi∣1≤i≤M-1},ω2:={yj∣1≤j≤M-1}.
Let x:=(x,y), X:=(X,Y) with X(x,y,t)=(X(x,y,t),Y(x,y,t)). By the transformation given in Appendix section, term Δϕ(·,t)∘X(x,tn+1;t) in (5) changes to expression that consists of ∂2ϕ/∂x2, ∂2ϕ/∂y2, ∂2ϕ/∂x∂y, ∂ϕ/∂x, ∂ϕ/∂y, ∂2X/∂x2, ∂2X/∂y2, ∂2X/∂x∂y, ∂X/∂x, ∂X/∂y, and so forth. If X(x,y,t), Y(x,y,t) are polynomials of x,y of degrees not more than one, then it holds that (see Appendix)
(42)Δϕ(·,t)∘X(x,tn+1;t)={∂2ϕ∂x2·[(∂X∂y)2+(∂Y∂y)2]+∂2ϕ∂y2·[(∂X∂x)2+(∂Y∂x)2]-2∂2ϕ∂x∂y[∂X∂x·∂X∂y+∂Y∂x·∂Y∂y]}·[∂Y∂x·∂X∂y-∂Y∂y·∂X∂x]-2.
Using centered difference to discretize these partial derivatives, we obtain the finite difference approximation of (5) as follows: find ϕh:Ωh×(0,T)→ℝ such that
(43)dϕhdt(X(xij,tn+1;t),t)-νΔhϕ(·,t)∘X(xij,tn+1;t)=f(X(xij,tn+1;t),t),
where Δhϕ(·,t)∘X(xij,tn+1;t) is the approximate of Δϕ(·,t)∘X(x,tn+1;t). Applying (17) to (43), with Φ~ representing the approximate of ϕh, we obtain the second-order-in-time finite difference scheme:(44)∑l=02‍αlΦ~n+ℓ-1(Xijn+ℓ-1)-νΔtΔh∑l=02‍βlΦ~n+ℓ-1∘Xijn+ℓ-1=Δt∑l=02‍βlfn+ℓ-1(Xijn+ℓ-1),i,j=1,2,…,M-1,n=1,2,…,N-1.
In (44) with X being replaced by XE, Φ representing the approximate of Φ~, we finally discretize (1a)–(1c) as follows:
(45)∑l=02‍αlΦijn+ℓ-1(XEn+ℓ-1)-νΔtΔh∑l=02‍βlΦijn+ℓ-1∘XEn+ℓ-1=Δt∑l=02‍βlfijn+ℓ-1(XEn+ℓ-1)i,j=1,2,…,M-1,n=1,2,…,N-1.
The computation of the starting values of this multistep scheme is similar to the general Eulerian schemes. Since XEn and XEn-1 in (45) are not grid points of Ωh, interpolations are needed. In the following numerical tests, we will use cubic interpolations (CI) and cubic spline interpolations, respectively (CSI).


## 5. Numerical Results

In this section, we test some specific cases of (45). The problem of rotating Gaussian hill has been widely used to test numerical schemes for convection-diffusion equations. In (1a)–(1c), let Ω=(-0.5,0.5)×(-0.5,0.5), T=π/4, u=(-y,x), f=0, and ν=1.25×10-4 or ν=1.0×10-3. The Dirichlet boundary conditions and initial condition of (1a)–(1c) are given such that the exact solution is
(46)c(x,y,t)=σσ+4νtexp{-(x-(t)-xc)2+(y-(t)-yc)2σ+4νt},
where x-(t)=xcost+ysint, y-(t)=-xsint+ycost, (xc,yc)=(0.25,0) and σ=0.01. The exact characteristic is given by
(47)X(x,y,t)=xcost-ysint,Y(x,y,t)=xsint+ycost.
By the Euler’s method the approximate characteristic is given by
(48)XEn(x,y)=x+yΔt,XEn-1(x,y)=x+2yΔt,YEn(x,y)=y-xΔt,YEn-1(x,y)=y-2xΔt.Note that since the velocityu does not satisfy the nonflow boundary condition (1b), the characteristics and its approximations are not necessarily contained in Ω. We assume that ϕ=0 outside of Ω (just as in Notsu and Tabata [16], Rui and Tabata [15], and Long and Yuan [32]).Using the standard central difference along the numerical characteristic and the transformation given in Appendix, we have(49)ΔhΦij∘XEn+1=ΔhΦij=Φi+1,jn+1+Φi-1,jn+1+Φi,j+1n+1+Φi,j-1n+1-4Φi,jn+1h2,ΔhΦij∘XEn=(Φi+1,jn+Φi-1,jn+Φi,j+1n+Φi,j-1n-4Φi,jn)(XEn)(1+Δt2)h2,ΔhΦij∘XEn-1=(Φi+1,jn-1+Φi-1,jn-1+Φi,j+1n-1+Φi,j-1n-1-4Φi,jn-1)(XEn-1)(1+4Δt2)h2.
Therefore the BDF2 scheme takes the form
(50)3Φijn+1-4Φijn(XEn)+Φijn-1(XEn-1)-2νΔtΔhΦijn+1=0,
and the CN scheme takes the form
(51)Φijn+1-Φijn(Xn)-(νΔt2)ΔhΦijn+1-(νΔt2)ΔhΦijn∘XEn=0.
To compare the results, we need the first-step backward difference scheme (BDF1):
(52)Φijn+1-Φijn(XEn)-νΔtΔhΦijn+1=0.
Let El2 denote the l2-error,
(53)El2=(h2∑ij|ϕijN-ΦijN|2)1/2,
where xij∈Ωh and tN=π/4.Figure1 illustrates convergence of BDF1 and BDF2 using numerical characteristics and using CI (Figure 1(a)) and CSI (Figure 1(b)), respectively. The straight lines are the first-order line and the second-order line, respectively. We can see that BDF2 exhibits higher-order convergence than BDF1.E
l
2-error of BDF1 and BDF2 versus the time step Δt in log-log scale, with T=π/4, ν=1.00×10-3; h=Δt=0.1,0.05,0.04,0.25,0.02,0.005. CI are involved in (a) and CSI in (b).
(a)
(b)Figure2 shows the convergence of CN schemes using exact characteristics (Figure 2(a)) and using numerical characteristics (Figure 2(b)), respectively. We can see that when first-order numerical characteristics are involved, the computed result is poor. This confirms the fact that the higher-order numerical characteristic line is necessary to ensure reasonable convergence.E
l
2-error of CN scheme versus Δt in log-log scale, with T=π/4, ν=1.25×10-4; h=Δt=0.1,0.05,0.04,0.25,0.02. Exact characteristic lines are involved in (a) and numerical ones in (b).
(a)
(b)Figure3 exhibits convergence of BDF2 with the smaller diffusion coefficient ν=1.25×10-4 using numerical characteristics. We see that either using CI or CSI, convergence is approximately second-order.Figure 3
E
l
2-errors of BDF2 using CI and CSI with ν=1.25×10-4; T=π/4, h=Δt=0.1,0.08,0.05,0.02,0.008,0.005.Remark 5.
As the particular cases of the general formulation, the following schemes are also tested:(54)α2=2,α1=-3,α0=1;β2=3/2,β1=-1/2,β0=0,α2=3/4,α1=-1/2,α0=-1/4;β2=3/2,β1=-7/4,β0=5/4.
Both of them exhibit similar convergence properties as the BDF2 schemes.


## 6. Conclusions

In this paper we formulate a general SL scheme of second-order. By the MOL approach, this formulation includes all previously introduced schemes. We prove that this general scheme is second-order accurate even if the first-order characteristic line is involved. This formulation is very easy to extend to higher order cases. We see that MMOC can be considered a special version of SL method. Convergence properties of SL finite difference schemes are numerically tested.



---

*Source: 763630-2013-01-31.xml*