# Time-Splitting Procedures for the Numerical Solution of the 2D Advection-Diffusion Equation

Authors: A. R. Appadu; H. H. Gidey
Journal: Mathematical Problems in Engineering
 (2013)

---

## Abstract

We perform a spectral analysis of the dispersive and dissipative properties of two time-splitting procedures, namely, locally one-dimensional (LOD) Lax-Wendroff and LOD (1, 5) [9] for the numerical solution of the 2D advection-diffusion equation. We solve a 2D numerical experiment described by an advection-diffusion partial differential equation with specified initial and boundary conditions for which the exact solution is known. Some errors are computed, namely, the error rate with respect to theL1 norm, dispersion and dissipation errors. Lastly, an optimization technique is implemented to find the optimal value of temporal step size that minimizes the dispersion error for both schemes when the spatial step is chosen as 0.025, and this is validated by numerical experiments.


---

## Body

## 1. Introduction

The advection-diffusion equation is a parabolic partial differential equation combining the diffusion and advection (convection) equations, which describes physical phenomena where particles, energy, or other physical quantities are transferred inside a physical system due to two processes: diffusion and advection [1]. The numerical solution of advection-diffusion equation plays an important role in many fields of science and engineering. These include the transport of air and groundwater pollutants, oil reservoir flow [2], heat transfer in draining film, flow through porous media, the dispersion of pollutants in rivers and streams, long range transport of pollutants in the atmosphere, thermal pollution in river systems, and dispersion of dissolved salts in groundwater [3].The 3D advection-diffusion equation is given by(1)∂u∂t+βx∂u∂x+βy∂u∂y+βz∂u∂z=αx∂2u∂x2+αy∂2u∂y2+αz∂2u∂z2,
where βx, βy, and βz are the velocity components of advection in the directions of x, y, and z, respectively, and αx, αy, and αz are the coefficients of diffusivity in the x-, y-, and z-directions, respectively.This study deals with the 2D advection-diffusion equation,(2)∂u∂t+βx∂u∂x+βy∂u∂y=αx∂2u∂x2+αy∂2u∂y2,
where 0<t≤T, in the domain 0≤x≤1, 0≤y≤1, with specified initial and boundary conditions.Dehghan [3] proposed two time-splitting procedures for the solution of the two-dimensional transport equation. The time-splitting procedure used is the locally-one dimensional (LOD) in proceeding from one time step to the next step. LOD replaces the complicated multidimensional partial differential equations by a sequence of solutions of simpler one-dimensional partial differential equations. The originality in this work is that we perform a spectral analysis of the dispersion and dissipation properties of the two schemes at some values of the temporal and spatial step sizes. 3D plots of the relative phase error per unit time step (RPE) and the modulus of the amplification factor (AFM) versus phase angles in x- and y-directions are obtained. We then use optimization strategies to compute the optimal values of the temporal step size for the two schemes when the spatial step size is 0.025. We then validate this result by performing a 2D numerical experiment with specified initial and boundary conditions. Lastly, we have the conclusion and references.


## 2. Numerical Dispersion and Dissipation

Dissipation reduces the amplitude of sinusoids in a Fourier series. This is caused by the presence of derivatives likeuxx and -uxxxx in the modified equation [4]. On the other hand, the amplitude of sinusoids in a Fourier series is increased by antidissipation. Derivatives like -uxx and uxxxx in the modified equation are generally antidissipative. Dispersion affects the speed of sinusoids in a Fourier series causing phase lag or phase lead and is caused due to the presence of odd-order derivatives in the modified equation [4].The modulus of the amplification factor (AFM) is a measure of the stability of a scheme and it is also used to measure the dissipative characteristics of the scheme. If the modulus of the amplification factor is equal to one, a disturbance neither grows nor damps [5]. If the modulus of the amplification factor is greater than one, then the scheme is unstable [6]; if it is less than one, damping occurs [7]. The partial differential equation given by (2) is dissipative in nature due to the terms uxx and uyy.The relative phase error (RPE) is a measure of the dispersive characteristics of a scheme. The relative phase error of a scheme approximating the 1D advection-diffusion equation is given by(3)RPE=-1cωarctan(ℑ(ξ)ℜ(ξ)),
where c is the Courant number, ω is phase angle, ξ is the amplification factor of the numerical scheme approximating the 1D advection-diffusion equation, and ℜ(ξ) and ℑ(ξ) are the real and imaginary parts of ξ, respectively [8].We extend the work on the relative phase error in [8] for the case of the 2D advection-diffusion equation. The relative phase error for a numerical scheme approximating (2) is obtained on substituting u by exp(I(ω1t-θ1x-θ2y)) [9], where I=-1, t is the time, θ1 and θ2 are wavenumbers in the directions of x and y, respectively, and ω1 is the dispersion relation. Thus, we get
(4)∂u∂t=Iω1exp(I(ω1t-θ1x-θ2y)),∂u∂x=-Iθ1exp(I(ω1t-θ1x-θ2y)),∂u∂y=-Iθ2exp(I(ω1t-θ1x-θ2y)),∂2u∂x2=-(θ1)2exp(I(ω1t-θ1x-θ2y)),∂2u∂y2=-(θ2)2exp(I(ω1t-θ1x-θ2y)).
On substituting (4) into (2), we obtain
(5)ω1=βxθ1+βyθ2+(αx(θ1)2)I+(αy(θ2)2)I.The exact phase velocity isℜ(ω1)/wavenumber and we obtain
(6)exactphasevelocity=βxθ1+βyθ2wavenumber.The amplification factor can be written asξ=ξ1+Iξ2, where ξ1 and ξ2 are the real and imaginary parts of ξ, respectively. Also, we can express ξ as ξ=exp(-bk), where k is time step and b is exponential growth rate [9]. Thus, we obtain
(7)b=1klog(ξ1-Iξ2ξ12+ξ22).The numerical phase velocity is calculated asℑ(b)/wavenumber and we get
(8)numericalphasevelocity=-(1/k)arctan(ξ2/ξ1)wavenumber.
The relative phase error is the ratio of the numerical phase velocity to the exact phase velocity [10]. It is calculated as
(9)RPE=-1k(βxθ1+βyθ2)arctan(ξ2ξ1).
The phase angles in the directions of x and y are given by ωx=Δxθ1 and ωy=Δyθ2, respectively. Hence, the relative phase error is
(10)RPE=-1cxωx+cyωyarctan(ξ2ξ1),
where cx=βxk/Δx and cy=βyk/Δy are the Courant numbers in the directions of x and y, respectively.


## 3. Time-splitting Procedures and Numerical Experiments

The domain we consider isx,y∈[0,1]. We divide the spatial interval [0,1] along x- and y-directions into N and M nodes, respectively, such that (N-1)Δx=1 and (M-1)Δy=1 and also divide the time interval [0,T] into P grid points such that (P-1)Δt=T. Then the grid points (xi,yj,tn) are defined by
(11)xi=iΔx,i=1,2,…,N,yj=jΔy,j=1,2,…,M,tn=nΔt,n=1,2,…,P.
For simplification we take Δx=Δy=h. Let Δt=k; then the parameters h and k represent the spatial and temporal grid spacing, respectively. We denote the approximated value of u at the grid point (i,j,n) by ui,jn.


### 3.1. Time-Splitting Procedures

Since one-dimensional schemes are easier to use than two-dimensional schemes, (2) is split into the following two one-dimensional equations:
(12)12∂u∂t+βx∂u∂x=αx∂2u∂x2,(13)12∂u∂t+βy∂u∂y=αy∂2u∂y2.
Each of (12) and (13) can be solved over half of a time step to be used for the complete 2D advection-diffusion equation, using the procedures developed for the 1D advection-diffusion equation.Some work on time-splitting procedures can be found in [3, 11]. In this paper, we refer to [3] on how a 2D advection-diffusion equation is converted into two 1D advection-diffusion equations using the locally one-dimensional (LOD) time-splitting procedure. Solving (12) and (13) in each half time step is equivalent to solving the following equations over a full-time step:
(14)∂u∂t+βx∂u∂x=αx∂2u∂x2,(15)∂u∂t+βy∂u∂y=αy∂2u∂y2.
Then we can use the schemes used for solving the 1D advection-diffusion equation to solve (14) and (15).


### 3.2. Numerical Experiments

In [12], different explicit and implicit finite differences schemes are used to solve the 1D advection-diffusion equation
(16)∂u∂t+∂u∂x=0.01∂2u∂x2.
Three values of the cell Reynolds number, namely, RΔ=2,4, and 8, are used for the numerical experiments. In [3], the time-splitting procedure is used to solve the two-dimensional transport equation in the region 0≤x≤1, 0≤y≤1, with
(17)αx=αy=α=0.01,βx=βy=β=0.8
and with the following boundary and initial conditions:
(18)g0(x,t)=14t+1exp[-(x-0.8t-0.5)20.01(4t+1)-(0.8t+0.5)20.01(4t+1)],(19)g1(x,t)=14t+1exp[-(x-0.8t-0.5)20.01(4t+1)-(0.5-0.8t)20.01(4t+1)],(20)h0(y,t)=14t+1exp[-(0.8t+0.5)20.01(4t+1)-(y-0.8t-0.5)20.01(4t+1)],(21)h1(y,t)=14t+1exp[-(0.5-0.8t)20.01(4t+1)-(y-0.8t-0.5)20.01(4t+1)],(22)f(x,y)=exp[-(x-0.5)20.01-(y-0.5)20.01],
for which the exact solution is
(23)u(x,y,t)=14t+1exp[-(x-0.8t-0.5)20.01(4t+1)u(x,y,t)=14t+1exp1-(y-0.8t-0.5)20.01(4t+1)].
We consider two time-splitting procedures LOD Lax-Wendroff and LOD (1,5) to solve (2) with αx=αy=α=0.01 and βx=βy=β=0.8, subject to the boundary conditions (18)–(22) at time T=0.3.We consider two values ofRΔ; RΔ=2 and RΔ=4. Since c=0.8k/h and s=0.01k/h2, we have RΔ=80h. Thus for RΔ=2 and 4, we have h=0.025 and h=0.05, respectively. When h=0.025, we have c=32k and s=16k and for h=0.05 we have c=16k and s=4k.


### 3.3. Quantification of Errors from Numerical Results

In this subsection, we describe how errors from numerical results can be quantified into dispersion and dissipation by a technique devised by Takacs [13].The total mean square error in the 1D case [13] is calculated as
(24)ETMS=1N∑i=1N‍(ui-vi)2,
where ui is the exact solution, vi is the numerical solution at a grid point i, and N is the number of grid points. The total mean square error can be expressed as
(25)1N∑i=1N‍(ui-vi)2=1N∑i=1N‍(ui-u-)2+1N∑i=1N‍(vi-v-)21N∑i=1N‍(ui-vi)2=+2N∑i=1N‍uiu-+2N∑i=1N‍viv--1N∑i=1N‍(u-)21N∑i=1N‍(ui-vi)2=-1N∑i=1N‍(v-)2-2N∑i=1N‍uivi.The right hand side of (25) can be rewritten as
(26)σ2(u)+σ2(v)+2(u-)2+2(v-)2-(u-)2-(v-)2-2N∑i=1N‍uivi,
where σ2(u) and σ2(v) denote the variance of u and v, respectively, and u- and v- denote the mean values of u and v, respectively. Then we have
(27)ETMS=σ2(u)+σ2(v)+((u-)2-2u-v-+(v-)2)ETOT=+(2u-v--2N∑i=1N‍uivi)ETOT=σ2(u)+σ2(v)+(u--v-)2-2(1N∑i=1N‍uivi-u-v-).
Therefore,
(28)1N∑i=1N‍(ui-vi)2=σ2(u)+σ2(v)+(u--v-)2-2Cov(u,v),
where Cov(u,v)=(1/N)(∑i=1N‍uivi-u-v-).The total mean square error can be expressed as(29)1N∑i=1N‍(ui-vi)2=(σ(u)-σ(v))21N∑i=1N‍(ui-vi)2=+(u--v-)2+2(1-ρ)σ(u)σ(v),
where ρ=Cov(u,v)/(σ(u)σ(v)) is the coefficient of correlation.The term2(1-ρ)σ(u)σ(v) measures the dispersion error and the term (σ(u)-σ(v))2+(u--v-)2 measures the dissipation error.We extend the work on quantification of errors in [14, 15] for the 2D case. The total mean square error for the 2D case is calculated as
(30)ETMS=1NM∑j=1M‍∑i=1N‍(ui,j-vi,j)2,
where ui,j and vi,j are the exact and numerical solutions at a grid point (i,j), respectively, and
(31)1NM∑j=1M‍∑i=1N‍(ui,j-vi,j)2=1NM∑j=1M‍∑i=1N‍(ui,j2-2ui,jvi,j+vi,j2)=1NM(∑j=1M‍∑i=1N‍ui,j2+∑j=1M‍∑i=1N‍vi,j2-2∑j=1M‍∑i=1N‍ui,jvi,j).
Since
(32)∑j=1M‍∑i=1N‍(ui,j-u-)2=∑j=1M‍∑i=1N‍(ui,j2-2ui,ju-+u-2),∑j=1M‍∑i=1N‍(vi,j-v-)2=∑j=1M‍∑i=1N‍(vi,j2-2vi,jv-+v-2),
we have
(33)ETMS=1NM[∑j=1M‍∑i=1N‍(ui,j-u-)2+∑j=1M‍∑i=1N‍(vi,j-v-)2ETMS=1NM1+2∑j=1M‍∑i=1N‍ui,ju-+2∑j=1M‍∑i=1N‍vi,jv-ETMS=1NM1-2∑j=1M‍∑i=1N‍ui,jvi,j-∑j=1M‍∑i=1N‍u-2-∑j=1M‍∑i=1N‍v-2]ETMS=σ2(u)+σ2(v)+2u-2+2v-2ETMS=-2NM∑j=1M‍∑i=1N‍ui,jvi,j-u-2-v-2ETMS=σ2(u)+σ2(v)+(u-2-2u-v-+v-2)ETMS=+2u-v--2NM∑j=1M‍∑i=1N‍ui,jvi,jETMS=σ2(u)+σ2(v)+(u--v-)2ETMS=-2[1NM∑j=1M‍∑i=1N‍ui,jvi,j-u-v-]ETMS=σ2(u)+σ2(v)+(u--v-)2-2Cov(u,v).Hence,(34)ETMS=(σ(u)-σ(v))2+(u--v-)2+2(1-ρ)σ(u)σ(v),
where ρ=Cov(u,v)/σ(u)σ(v). The dissipation error is (σ(u)-σ(v))2+(u--v-)2 and the dispersion error is 2(1-ρ)σ(u)σ(v).The error rate with respect toL1 norm for x,y∈[0,1] is calculated as
(35)Enum=1(N-1)(M-1)∑j=1M‍∑i=1N‍|ui,j-vi,j|.


## 3.1. Time-Splitting Procedures

Since one-dimensional schemes are easier to use than two-dimensional schemes, (2) is split into the following two one-dimensional equations:
(12)12∂u∂t+βx∂u∂x=αx∂2u∂x2,(13)12∂u∂t+βy∂u∂y=αy∂2u∂y2.
Each of (12) and (13) can be solved over half of a time step to be used for the complete 2D advection-diffusion equation, using the procedures developed for the 1D advection-diffusion equation.Some work on time-splitting procedures can be found in [3, 11]. In this paper, we refer to [3] on how a 2D advection-diffusion equation is converted into two 1D advection-diffusion equations using the locally one-dimensional (LOD) time-splitting procedure. Solving (12) and (13) in each half time step is equivalent to solving the following equations over a full-time step:
(14)∂u∂t+βx∂u∂x=αx∂2u∂x2,(15)∂u∂t+βy∂u∂y=αy∂2u∂y2.
Then we can use the schemes used for solving the 1D advection-diffusion equation to solve (14) and (15).


## 3.2. Numerical Experiments

In [12], different explicit and implicit finite differences schemes are used to solve the 1D advection-diffusion equation
(16)∂u∂t+∂u∂x=0.01∂2u∂x2.
Three values of the cell Reynolds number, namely, RΔ=2,4, and 8, are used for the numerical experiments. In [3], the time-splitting procedure is used to solve the two-dimensional transport equation in the region 0≤x≤1, 0≤y≤1, with
(17)αx=αy=α=0.01,βx=βy=β=0.8
and with the following boundary and initial conditions:
(18)g0(x,t)=14t+1exp[-(x-0.8t-0.5)20.01(4t+1)-(0.8t+0.5)20.01(4t+1)],(19)g1(x,t)=14t+1exp[-(x-0.8t-0.5)20.01(4t+1)-(0.5-0.8t)20.01(4t+1)],(20)h0(y,t)=14t+1exp[-(0.8t+0.5)20.01(4t+1)-(y-0.8t-0.5)20.01(4t+1)],(21)h1(y,t)=14t+1exp[-(0.5-0.8t)20.01(4t+1)-(y-0.8t-0.5)20.01(4t+1)],(22)f(x,y)=exp[-(x-0.5)20.01-(y-0.5)20.01],
for which the exact solution is
(23)u(x,y,t)=14t+1exp[-(x-0.8t-0.5)20.01(4t+1)u(x,y,t)=14t+1exp1-(y-0.8t-0.5)20.01(4t+1)].
We consider two time-splitting procedures LOD Lax-Wendroff and LOD (1,5) to solve (2) with αx=αy=α=0.01 and βx=βy=β=0.8, subject to the boundary conditions (18)–(22) at time T=0.3.We consider two values ofRΔ; RΔ=2 and RΔ=4. Since c=0.8k/h and s=0.01k/h2, we have RΔ=80h. Thus for RΔ=2 and 4, we have h=0.025 and h=0.05, respectively. When h=0.025, we have c=32k and s=16k and for h=0.05 we have c=16k and s=4k.


## 3.3. Quantification of Errors from Numerical Results

In this subsection, we describe how errors from numerical results can be quantified into dispersion and dissipation by a technique devised by Takacs [13].The total mean square error in the 1D case [13] is calculated as
(24)ETMS=1N∑i=1N‍(ui-vi)2,
where ui is the exact solution, vi is the numerical solution at a grid point i, and N is the number of grid points. The total mean square error can be expressed as
(25)1N∑i=1N‍(ui-vi)2=1N∑i=1N‍(ui-u-)2+1N∑i=1N‍(vi-v-)21N∑i=1N‍(ui-vi)2=+2N∑i=1N‍uiu-+2N∑i=1N‍viv--1N∑i=1N‍(u-)21N∑i=1N‍(ui-vi)2=-1N∑i=1N‍(v-)2-2N∑i=1N‍uivi.The right hand side of (25) can be rewritten as
(26)σ2(u)+σ2(v)+2(u-)2+2(v-)2-(u-)2-(v-)2-2N∑i=1N‍uivi,
where σ2(u) and σ2(v) denote the variance of u and v, respectively, and u- and v- denote the mean values of u and v, respectively. Then we have
(27)ETMS=σ2(u)+σ2(v)+((u-)2-2u-v-+(v-)2)ETOT=+(2u-v--2N∑i=1N‍uivi)ETOT=σ2(u)+σ2(v)+(u--v-)2-2(1N∑i=1N‍uivi-u-v-).
Therefore,
(28)1N∑i=1N‍(ui-vi)2=σ2(u)+σ2(v)+(u--v-)2-2Cov(u,v),
where Cov(u,v)=(1/N)(∑i=1N‍uivi-u-v-).The total mean square error can be expressed as(29)1N∑i=1N‍(ui-vi)2=(σ(u)-σ(v))21N∑i=1N‍(ui-vi)2=+(u--v-)2+2(1-ρ)σ(u)σ(v),
where ρ=Cov(u,v)/(σ(u)σ(v)) is the coefficient of correlation.The term2(1-ρ)σ(u)σ(v) measures the dispersion error and the term (σ(u)-σ(v))2+(u--v-)2 measures the dissipation error.We extend the work on quantification of errors in [14, 15] for the 2D case. The total mean square error for the 2D case is calculated as
(30)ETMS=1NM∑j=1M‍∑i=1N‍(ui,j-vi,j)2,
where ui,j and vi,j are the exact and numerical solutions at a grid point (i,j), respectively, and
(31)1NM∑j=1M‍∑i=1N‍(ui,j-vi,j)2=1NM∑j=1M‍∑i=1N‍(ui,j2-2ui,jvi,j+vi,j2)=1NM(∑j=1M‍∑i=1N‍ui,j2+∑j=1M‍∑i=1N‍vi,j2-2∑j=1M‍∑i=1N‍ui,jvi,j).
Since
(32)∑j=1M‍∑i=1N‍(ui,j-u-)2=∑j=1M‍∑i=1N‍(ui,j2-2ui,ju-+u-2),∑j=1M‍∑i=1N‍(vi,j-v-)2=∑j=1M‍∑i=1N‍(vi,j2-2vi,jv-+v-2),
we have
(33)ETMS=1NM[∑j=1M‍∑i=1N‍(ui,j-u-)2+∑j=1M‍∑i=1N‍(vi,j-v-)2ETMS=1NM1+2∑j=1M‍∑i=1N‍ui,ju-+2∑j=1M‍∑i=1N‍vi,jv-ETMS=1NM1-2∑j=1M‍∑i=1N‍ui,jvi,j-∑j=1M‍∑i=1N‍u-2-∑j=1M‍∑i=1N‍v-2]ETMS=σ2(u)+σ2(v)+2u-2+2v-2ETMS=-2NM∑j=1M‍∑i=1N‍ui,jvi,j-u-2-v-2ETMS=σ2(u)+σ2(v)+(u-2-2u-v-+v-2)ETMS=+2u-v--2NM∑j=1M‍∑i=1N‍ui,jvi,jETMS=σ2(u)+σ2(v)+(u--v-)2ETMS=-2[1NM∑j=1M‍∑i=1N‍ui,jvi,j-u-v-]ETMS=σ2(u)+σ2(v)+(u--v-)2-2Cov(u,v).Hence,(34)ETMS=(σ(u)-σ(v))2+(u--v-)2+2(1-ρ)σ(u)σ(v),
where ρ=Cov(u,v)/σ(u)σ(v). The dissipation error is (σ(u)-σ(v))2+(u--v-)2 and the dispersion error is 2(1-ρ)σ(u)σ(v).The error rate with respect toL1 norm for x,y∈[0,1] is calculated as
(35)Enum=1(N-1)(M-1)∑j=1M‍∑i=1N‍|ui,j-vi,j|.


## 4. Construction of the LOD Lax-Wendroff Procedure

We use the following approximations in the first time step of the LOD procedure [3]:
(36)∂u∂t|i,jn≃ui,jn+(1/2)-ui,jnΔt,∂u∂x|i,jn≃(1-γx)(ui,jn-ui-1,jn)+γx(ui+1,jn-ui,jn)Δx,∂2u∂x2|i,jn≃ui+1,jn-2ui,jn-ui-1,jn(Δx)2,
where γx is the spatial weighting factor in the direction of x. The following approximations are used for the second step of the LOD procedure:
(37)∂u∂t|i,jn+(1/2)≃ui,jn+1-ui,jn+(1/2)Δt,∂u∂y|i,jn+(1/2)≃(1-γy)(ui,jn+(1/2)-ui,j-1n+(1/2))+γy(ui,j+1n+(1/2)-ui,jn+(1/2))Δy,∂2u∂y2|i,jn+(1/2)≃ui,j+1n+(1/2)-2ui,jn+(1/2)-ui,j-1n+(1/2)(Δy)2,
where γy is the spatial weighting factor in the direction of y. By using the following relationships
(38)cx=βxkh,sx=αxkh2,γx=1-cx2,cy=βykh,sy=αykh2,γy=1-cy2,
the finite difference formula for the first half time step of the LOD Lax-Wendroff procedure is given by
(39)ui,jn+(1/2)=12(2sx+cx+cx2)ui-1,jn+(1-2sx-cx2)ui,jnui,jn+(1/2)=+12(2sx-cx+cx2)ui+1,jn.
At the second half time step the finite difference is given by
(40)ui,jn+1=12(2sy+cy+cy2)ui,j-1n+(1/2)+(1-2sy-cy2)ui,jn+(1/2)ui,jn+1=+12(2sy-cy+cy2)ui,j+1n+(1/2),
where
(41)ui,j-1n+(1/2)=12(2sx+cx+cx2)ui-1,j-1n+(1-2sx-cx2)ui,j-1nui,j-1n+(1/2)=+12(2sx-cx+cx2)ui+1,j-1n,ui,j+1n+(1/2)=12(2sx+cx+cx2)ui-1,j+1n+(1-2sx-cx2)ui,j+1nui,j-1n+(1/2)=+12(2sx-cx+cx2)ui+1,j+1n.
We obtain a single expression for ui,jn+1 in a complete time step as follows:
(42)ui,jn+1=14(2sy+cy+cy2)(2sx+cx+cx2)ui-1,j-1nui,jn+1=+12(2sy+cy+cy2)(1-2sx-cx2)ui,j-1nui,jn+1=+14(2sy+cy+cy2)(2sx-cx+cx2)ui+1,j-1nui,jn+1=+12(1-2sy-cy2)(2sx+cx+cx2)ui-1,jnui,jn+1=+(1-2sy-cy2)(1-2sx-cx2)ui,jnui,jn+1=+12(1-2sy-cy2)(2sx-cx+cx2)ui+1,jnui,jn+1=+14(2sy-cy+cy2)(2sx+cx+cx2)ui-1,j+1nui,jn+1=+12(2sy-cy+cy2)(1-2sx-cx2)ui,j+1nui,jn+1=+14(2sy-cy+cy2)(2sx-cx+cx2)ui+1,j+1n.
To find the modified equation of the scheme, we first find Taylor’s expansion of each term in (42) about ui,jn. The Taylor series expansion of ui,jn+1 is given by
(43)ui,jn+1=u+kut+k22utt+k36uttt+k44!utttt+⋯.The Taylor series expansions for some grid points aboutui,jn are given as follows:
(44)ui-1,j-1n=u-hux-huy+h22[uxx+2uxy+uyy]ui-1,j-1n=-h36[uxxx+3uxxy+3uxyy+uyyy]ui-1,j-1n=+h44![uxxxx+4uxxxy+6uxxyyui-1,j-1n=+h44!1+4uxyyy+uyyyy]+⋯,ui-1,jn=u-hux+h22uxx-h36uxxx+h44!uxxxx+⋯,ui,j-1n=u-huy+h22uyy-h36uyyy+h44!uyyyy+⋯,ui-1,j+1n=u-hux+huy+h22[uxx-2uxy+uyy]ui-1,j+1n=-h36[uxxx-3uxxy+3uxyy-uyyy]ui-1,j+1n=+h44![uxxxx-4uxxxy+6uxxyyui-1,j+1n=+h44!1-4uxyyy+uyyyy]+⋯,ui+1,j-1n=u+hux-huy+h22[uxx-2uxy+uyy]ui+1,j-1n=-h36[-uxxx+3uxxy-3uxyy+uyyy]ui+1,j-1n=+h44![uxxxx-4uxxxyui+1,j-1n=+h44!1+6uxxyy-4uxyyy+uyyyy]+⋯,ui+1,jn=u+hux+h22uxx+h36uxxx+h44!uxxxx+⋯,ui,j+1n=u+huy+h22uyy+h36uyyy+h44!uyyyy+⋯,ui+1,j+1n=u+hux+huy+h22[uxx+2uxy+uyy]ui+1,j+1n=+h36[uxxx+3uxxy+3uxyy+uyyy]ui+1,j+1n=+h44![uxxxx+4uxxxy+6uxxyyui+1,j+1n=+h44!1+4uxyyy+uyyyy]+⋯.
Using (2), we convert the temporal derivatives utt, uttt, and utttt into spatial derivatives. We thus have
(45)utt=αx2uxxxx+αy2uyyyy+αxαy[uxxyy+uyyxx]utt=-2αxβxuxxx-2αyβyuyyy-αxβy[uxxy+uyxx]utt=-αyβx[uxyy+uyyx]+βx2uxx+βy2uyyutt=+βxβy[uxy+uyx]+⋯,uttt=3αxβx2uxxxx+3αyβy2uyyyyuttt=+αxβy2[uyyxx+uyxxy+uxxyy]uttt=+αyβx2[uxxyy+uyyxx+uxyyx]uttt=+αyβxβy[2uyyyx+uyyxy+2uxyyy+2uyxyy]uttt=+αxβxβy[2uyxxx+2uxxxy+uxxyx+2uxyxx]uttt=-βx3uxxx-βxβy2[uyyx+uyxy+uxyy]uttt=-βy3uyyy-βx2βy[uxxy+uxyx+uyxx]+⋯,utttt=βx[βx3uxxxx+βx2βy[uxxyx+uxyxx+uyxxx]utttt=βx+βxβy2[uyyxx+uyxyx+xyyx]+βy3uyyyx]utttt=+βy[βy3uyyyy+βy2βx[uyyxy+uyxyy+uxyyy]utttt=+βy+βyβx2[uxxyy+uxyxy+uyxxy]+βx3uxxxy]+⋯.
Then on substituting the Taylor series expansions of each term of the difference scheme we get the following modified equation for the LOD Lax-Wendroff
(46)ut+βxux+βyuy=αxuxx+αyuyy-16βxh2(1-cx2-6sx)uxxx=-16βyh2(1-cy2-6sy)uyyy=-h324[-2sx-cx2+12sxcx2+12sx2+cx4]uxxxx=-h324[-2sy-cy2+12sycy2+12sy2+cy4]uyyyy+⋯.
The scheme is second order accurate in space and the leading error terms are dispersive in nature (presence of odd-order derivatives uxxx and uyyy). As the time and spatial increments go to zero, the modified equation (46) reduces to its original equation, that is, (2). Hence LOD Lax-Wendroff is consistent.We now study the spectral analysis of the dispersive and dissipative properties of the scheme for the casecx=cy=c and sx=sy=s. To obtain the amplification factor we use the Von Neumann stability analysis by substituting ui,jn by ξnexp(I(iωx+jωy)) in (42), where I=-1. We thus have
(47)ξ=14(2sy+cy+cy2)(2sx+cx+cx2)exp(-I(ωx+ωy))ξ=+12(2sy+cy+cy2)(1-2sx-cx2)exp(-Iωy)ξ=+14(2sy+cy+cy2)(2sx-cx+cx2)exp(I(ωx-ωy))ξ=+12(1-2sy-cy2)(2sx+cx+cx2)exp(-Iωx)ξ=+(1-2sy-cy2)(1-2sx-cx2)ξ=+12(1-2sy-cy2)(2sx-cx+cx2)exp(Iωx)ξ=+14(2sy-cy+cy2)(2sx+cx+cx2)exp(-I(ωx-ωy))ξ=+12(2sy-cy+cy2)(1-2sx-cx2)exp(Iωy)ξ=+14(2sy-cy+cy2)(2sx-cx+cx2)exp(I(ωx+ωy)).
The real and imaginary parts of the amplification factor are given by
(48)Real(ξ)=1+4sc2cos(ωx)cos(ωy)+c4cos(ωx)cos(ωy)Real(ξ)=+4s2cos(ωx)cos(ωy)-4sc2cos(ωy)Real(ξ)=-c2sin(ωx)sin(ωy)+4s2+4sc2-4sc2cos(ωx)Real(ξ)=+c2cos(ωx)-c4cos(ωx)+4scos(ωx)Real(ξ)=+2scos(ωy)+2scos(ωx)-4s2cos(ωy)Real(ξ)=-4s2cos(ωx)-2c2+c2cos(ωy)Real(ξ)=-c4cos(ωy)+c4-4s,Imag(ξ)=-c3cos(ωx)sin(ωy)+2scsin(ωx)Imag(ξ)=+2scsin(ωy)-2sccos(ωx)sin(ωy)Imag(ξ)=-2scsin(ωx)cos(ωy)-csin(ωx)+c3sin(ωx)Imag(ξ)=-csin(ωy)+c3sin(ωy)-c3sin(ωx)cos(ωy),
respectively. The modulus of the amplification factor, AFM, is obtained as
(49)AFM=(ℜ(ξ))2+(ℑ(ξ))2.
We find the region of stability using the approach of Hindmarsh et al. [16] and Sousa [17]. We consider the case when ωx=ωy=π and ωx→0 and ωy→0. When ωx=ωy=π, (49) gives
(50)|ξ|=1+16sc2-4c2+4c4+16s2-8s.
From the Von Neumann stability analysis, the scheme is stable if and only if |ξ|≤1. Thus, we get
(51)1+16sc2-4c2+4c4+16s2-8s≤1.On simplification, we get(52)|s+2c2-14|≤14,
which reduces to
(53)s≤1-c22.Whenωx→0 and ωy→0, we use the following approximations:
(54)cos(ωx)≃1-ωx22,cos(ωy)≃1-ωy22,sin(ωx)≃ωx,sin(ωy)≃ωy.
We consider (49) and use the approximations in (54) to obtain
(55)|ξ|2≃1-2s(ωx2+ωy2).Thus LOD Lax-Wendroff is stable if-2s≤0. Therefore, we have
(56)s≥0.On combining (53) and (56), we obtain the region of stability for the LOD Lax-Wendroff procedure as
(57)0≤s≤1-c22.
We choose α=0.01 and β=0.8 [3]. For h=0.025 from (57), we have
(58)0≤0.01k(0.025)2≤1-(0.8k/0.025)22.
On solving for k, we get
(59)0≤k≤0.0193.Therefore, forh=0.025, the stability region for the LOD Lax-Wendroff procedure is 0≤k≤0.0193. We choose k such that 0.3/k is an integer as for our numerical experiments, T=0.3. Then 0.3/k gives the number of time intervals. We choose k=0.005,0.01,3/160, and 1/60.Forh=0.05, we have c=16k and s=4k. When ωx=π and ωy=π, using (51) we get
(60)0≤1-32k+1638k3-768k2+262144k4≤1,
which gives
(61)0≤k≤0.0487985.
Therefore, for h=0.05, LOD Lax-Wendroff is stable if 0≤k≤0.0487985. We then choose some values of k∈(0,0.048) for the numerical experiments. for our numerical experiments.3D plots of the modulus of the amplification factor versus phase angles inx- and y-directions for some different values of h and k are depicted in Figures 1 and 2. 2D plots of the modulus of the amplification factor versus ωx, when ωy=0, are shown in Figure 3. The scheme is in general less dissipative at h=0.05 as compared to h=0.025. Out of the eight combinations of values of h and k, the scheme is least dissipative when h=0.05 and k=0.005 and it is most dissipative when h=0.025 and k=0.01.Plots of modulus of amplification factor versus phase angle inx-direction, ωx, versus phase angle in y-direction, ωy, at h=0.025 with some different values of k for the LOD Lax-Wendroff scheme.
(a)
h
=
0.025
,
k
=
0.005
(b)
h
=
0.025
,
k
=
0.01
(c)
h
=
0.025
,
k
=
1
/
60
(d)
h
=
0.025
,
k
=
3
/
160Plots of modulus of amplification factor versus phase angle inx-direction, ωx, versus phase angle in y-direction, ωy, at h=0.05 with some different values of k for the LOD Lax-Wendroff scheme.
(a)
h
=
0.05
,
k
=
0.005
(b)
h
=
0.05
,
k
=
0.01
(c)
h
=
0.05
,
k
=
0.02
(d)
h
=
0.05
,
k
=
0.03Plots of modulus of amplification factor versusωx when ωy=0 at h=0.025 and h=0.05 with some different values of k for the LOD Lax-Wendroff scheme.
(a)
h
=
0.025
(b)
h
=
0.05Figures4 and 5 show the 3D plots of the relative phase error versus ωx versus ωy for some different values of h and k. Figure 6 shows the 2D plots of the relative phase error versus ωx, for the case ωy=0 at h=0.025 and h=0.05. For h=0.025, we observe phase lag behaviour at k=0.01 and k=0.005 and phase lead behaviour at k=1/60 and k=3/160. For h=0.05, we have phase lag behaviour at k=0.005,0.01,0.02, and 0.03. The scheme is least dispersive when h=0.05; k=0.03 and h=0.025; k=0.01. In the following section, we consider the LOD (1,5) scheme.Plots of relative phase error versus phase angle inx-direction (ωx) versus phase angle in y-direction (ω-y), at h=0.025 with some different values of k for the LOD Lax-Wendroff scheme.
(a)
h
=
0.025
,
k
=
0.005
(b)
h
=
0.025
,
k
=
0.01
(c)
h
=
0.025
,
k
=
1
/
60
(d)
h
=
0.025
,
k
=
3
/
160Plots of relative phase error versusωx versus ωy, at h=0.05 with some different values of k for the LOD Lax-Wendroff scheme.
(a)
h
=
0.05
,
k
=
0.005
(b)
h
=
0.05
,
k
=
0.01
(c)
h
=
0.05
,
k
=
0.02
(d)
h
=
0.05
,
k
=
0.03Plots of relative phase error versusωx when ωy=0 at h=0.025 and h=0.05 with some different values of k for the LOD Lax-Wendroff scheme.
(a)
h
=
0.025
(b)
h
=
0.05


## 5. LOD(1,5) Explicit Procedure

In this procedure the following approximations are used in the first half time step [3]:
(62)∂u∂t|i,jn≃ui,jn+(1/2)-ui,jnΔt,∂u∂x|i,jn≃(12sx+2cx2-3cx-212)(ui+2,jn-ui,jn2Δx)∂u∂x|i,jn≃+(12sx+2cx2+3cx-212)(ui,jn-ui-2,jn2Δx)∂u∂x|i,jn≃-(cx2+6sx-412)(ui+1,jn-ui-1,jn2Δx),∂2u∂x2|i,jn≃(-cx4+4cx2-12sx2-12sxcx2+8sx6sx)∂2u∂x2|i,jn≃×(ui+1,jn-2ui,jn+ui-1,jn(Δx)2)∂2u∂x2|i,jn≃+(cx4-4cx2+12sx2+12sxcx2-2sx6sx)∂2u∂x2|i,jn≃×(ui+2,jn-2ui,jn+ui-2,jn4(Δx)2).
On substituting (62) into (14), we get the following finite difference equation for the first half time step:
(63)ui,jn+1/2=Axui-2,jn+Bxui-1,jn+Cxui,jnui,jn+1/2=+Dxui+1,jn+Exui+2,jn,
where
(64)Ax=124(12sx(sx+cx2)+2sx(6cx-1)Ax=1241+cx(cx-1)(cx+1)(cx+2)(sx+cx2)),Bx=-16(12sx(sx+cx2)+2sx(3cx-4)Bx=-161+cx(cx-2)(cx+1)(cx+2)(sx+cx2)),Cx=14(12sx(sx+cx2)-10sxCx=141+(cx-1)(cx-2)(cx+1)(cx+2)(sx+cx2)),Dx=-16(12sx(sx+cx2)-2sx(3cx+4)Dx=-161+cx(cx-2)(cx-1)(cx+2)(sx+cx2)),Ex=124(12sx(sx+cx2)-2sx(6cx+1)Ex=1241+cx(cx-1)(cx+1)(cx-2)(sx+cx2)).The following approximations are used in the second half time step:(65)∂u∂t|i,jn+(1/2)≃ui,jn+1-ui,jn+(1/2)Δt,∂u∂x|i,jn+(1/2)≃(12sy+2cy2-3cy-212)(ui,j+2n+(1/2)-ui,jn+(1/2)2Δy)+(12sy+2cy2+3cy-212)(ui,jn+(1/2)-ui,j-2n+(1/2)2Δy)-(cy2+6sy-412)(ui,j+1n+(1/2)-ui,j-1n+(1/2)2Δy),∂2u∂y2|i,jn+(1/2)≃(-cy4+4cy2-12sy2-12sycy2+8sy6sy)∂2u∂y2|i,jn+(1/2)≃×(ui,j+1n+(1/2)-2ui,jn+(1/2)+ui,j-1n+(1/2)(Δy)2)∂2u∂y2|i,jn+(1/2)≃+(cy4-4cy2+12sy2+12sycy2-2sy6sy)∂2u∂y2|i,jn+(1/2)≃×(ui,j+2n+(1/2)-2ui,jn+(1/2)+ui,j-2n+(1/2)4(Δy)2).
Then on substituting (65) into (15), we get the following difference equation for the second half time step:
(66)ui,jn+1=Ayui,j-2n+(1/2)+Byui,j-1n+(1/2)+Cyui,jn+(1/2)ui,jn+1=+Dyui,j+1n+(1/2)+Eyui,j+2n+(1/2),
where
(67)Ay=124(12sy(sy+cy2)+2sy(6cy-1)Ay=1241+cy(cy-1)(cy+1)(cy+2)),By=-16(12sy(sy+cy2)+2sy(3cy-4)By=-161+cy(cy-2)(cy+1)(cy+2)),Cy=14(12sy(sy+cy2)-10syCy=141+(cy-1)(cy-2)(cy+1)(cy+2)),Dy=-16(12sy(sy+cy2)-2sy(3cy+4)Dy=-161+cy(cy-2)(cy-1)(cy+2)),Ey=124(12sy(sy+cy2)-2sy(6cy+1)Ey=1241+cy(cy-1)(cy+1)(cy-2)).The complete LOD(1,5) scheme is given by
(68)ui,jn+1=AyAxui-2,j-2n+AyBxui-1,j-2n+AyCxui,j-2nui,jn+1=+AyDxui+1,j-2n+AyExui+2,j-2nui,jn+1=+ByAxui-2,j-1n+ByBxui-1,j-1n+ByCxui,j-1nui,jn+1=+ByDxui+1,j-1n+ByExui+2,j-1nui,jn+1=+CyAxui-2,jn+CyBxui-1,jn+CyCxui,jnui,jn+1=+CyDxui+1,jn+CyExui+2,jnui,jn+1=+DyAxui-2,j+1n+DyBxui-1,j+1n+DyCxui,j+1nui,jn+1=+DyDxui+1,j+1n+DyExui+2,j+1nui,jn+1=+EyAxui-2,j+2n+EyBxui-1,j+2n+EyCxui,j+2nui,jn+1=+EyDxui+1,j+2n+EyExui+2,j+2n.
The Taylor series expansion of the terms on the right hand side of (68) about ui,jn is given as follows:
(69)ui+a,j+bn=u+ahux+bhuy+h22(a2uxx+2abuxy+b2uyy)ui+a,j+bn=+h36(a3uxxx+3a2buxxy+3ab2uxyy+b3uyyy)ui+a,j+bn=+h424(a4uxxxx+4a3buxxxy+6a2b2uxxyyui+a,j+bn=+h4241+4ab3uxyyy+b4uyyyy)ui+a,j+bn=+h5120(a5uxxxxx+5a4buxxxxyui+a,j+bn=+h51201+10a3b2uxxxyy+10a2b3uxxyyyui+a,j+bn=+h51201+5ab4uxyyyy+b5uyyyyy)+⋯,
for a=-2,-1,0,1,2 and b=-2,-1,0,1,2.We obtain the following modified equation:(70)ut+βxux+βyuy=αxuxx+αyuyyut+βxux+βyuy=+βxh4120(60sx2+20sxcx2+cx4ut+βxux+βyuy=+βxh41201-5cx2+4-30sx)uxxxxxut+βxux+βyuy=+βyh4120(60sy2+20sycy2+cy4ut+βxux+βyuy=+βxh41201-5cy2+4-30sy)uyyyyy+⋯.
The scheme is essentially dispersive as the leading error terms are dispersive in nature due to the presence of the odd-order derivative terms uxxxxx and uyyyyy. Also, the scheme is consistent and it is fourth order accurate in space.The amplification factor of the LOD(1,5) scheme is given by
(71)ξ=Ay[Axexp(I(-2ωx-2ωy))+Bxexp(I(-ωx-2ωy))ξ=Ay1+Cxexp(-2Iωy)ξ=Ay1+Dxexp(I(ωx-2ωy))+Exexp(I(2ωx-2ωy))]ξ=+By[Axexp(I(-2ωx-ωy))+Bxexp(I(-ωx-ωy))ξ=+By1+Cxexp(-Iωy)+Dxexp(I(ωx-ωy))ξ=+By1+Exexp(I(2ωx-ωy))]ξ=+Cy[Axexp(-2Iωx)+Bxexp(-Iωx)ξ=+Cy1+Cx+Dxexp(Iωx)+Exexp(2Iωx)]ξ=+Dy[Axexp(I(-2ωx+ωy))+Bxexp(-I(ωx+ωy))ξ=+Dy1+Cxexp(Iωy)+Dxexp(I(ωx+ωy))ξ=+Dy1+Exexp(I(2ωx+ωy))]ξ=+Ey[Axexp(I(-2ωx+2ωy))ξ=+Ey1+Bxexp(I(-ωx+2ωy))+Cxexp(2Iωy)ξ=+Ey1+Dxexp(I(ωx+2ωy))ξ=+Ey1+Exexp(I(2ωx+2ωy))].
We consider cx=cy=c and sx=sy=s. We use the Von Neumann stability analysis and the approach of Hindmarsh et al. [16] to obtain the stability region. When ωx=ωy=π, on simplification of (71), we get
(72)|ξ|=19(3-8c2+2c4-16s+24sc2+24s2)2.
For stability, we must have
(73)19(3-8c2+2c4-16s+24sc2+24s2)2≤1.
We consider (71) and for ωx→0 and ωy→0 we use the following approximations:
(74)cos(ωx)≃1-ωx22,cos(ωy)≃1-ωy22,cos(2ωx)≃1-2ωx2,cos(2ωy)≃1-2ωy2,sin(ωx)≃ωx,sin(ωy)≃ωy,sin(2ωx)≃2ωx,sin(2ωy)≃2ωy.
We thus have
(75)|ξ|2≃1-2s(ωx2+ωy2)+4s2ωx2ωy2|ξ|2≃+(-2s3ωy2+sc2+14c4-2c2s2ωy2-12sc4ωy2+s2)ωx4|ξ|2≃+(s4+116c8+32s2c4+12sc6+2s3c2)ωx4ωy4|ξ|2≃+(-12c4sωx2-2c2s2ωx2+14c4+s2-2s3ωx2+c2s)ωy4.
On neglecting higher order terms, we have
(76)|ξ|2≃1-2s(ωx2+ωy2).
Thus, the numerical method is stable if -2s≤0. Therefore, we have
(77)s≥0.
When h=0.025, we have c=32k and s=16k. Using (73) and (77), we get
(78)19(2097152k4+393216k3-2048k2-256k+3)2≤1,(79)k≥0,
respectively. On combining (78) and (79), we get the stability region when h=0.025 for the LOD (1,5) procedure as
(80)0≤k≤0.026288.
When h=0.05, we have c=16k and s=4k and the modulus of the amplification factor is given by
(81)|ξ|=19(131072k4+24576k3-1664k2-64k+3)2.
From (77) and (81), we obtain the stability region for the LOD (1,5) procedure when h=0.05 as
(82)0≤k≤0.073865.
We analyse the spectral analysis; for h=0.025, we choose k=0.005,0.01,0.02,0.025 and for h=0.05, we choose k=0.01,0.02,0.03 and 0.05.3D plots of the modulus of the amplification factor versus phase angle inx-direction versus phase angle in y-direction at two values of h, namely, 0.025 and 0.05, at some different values of k are shown in Figures 7 and 8. 2D plots of the modulus of the amplification factor versus ωx when ωy=0 are illustrated in Figure 9. The scheme is in general less dissipative at h=0.05 as compared to h=0.025. Out of the eight combinations of h and k values, the scheme is least dissipative when h=0.05 and k=0.01.Plots of modulus of amplification factor versus phase angle inx-direction (ωx) versus phase angle in y-direction (ωy), at h=0.025 for some different values of k for the LOD (1,5) procedure.
(a)
h
=
0.025
,
k
=
0.005
(b)
h
=
0.025
,
k
=
0.01
(c)
h
=
0.025
,
k
=
0.02
(d)
h
=
0.025
,
k
=
0.025Plots of modulus of amplification factor versus phase angle inx-direction (ωx) versus phase angle in y-direction (ωy) at h=0.05 for some different values of k for the LOD (1,5) procedure.
(a)
h
=
0.05
,
k
=
0.01
(b)
h
=
0.05
,
k
=
0.02
(c)
h
=
0.05
,
k
=
0.03
(d)
h
=
0.05
,
k
=
0.05Plots of modulus of amplification factor versusωx when ωy=0 at h=0.025 and h=0.05 with some different values of k for the LOD (1,5) procedure.
(a)
h
=
0.025
(b)
h
=
0.05Figures10 and 11 show the 3D plots of the relative phase error versus ωx versus ωy for some different values of h and k. Figure 12 shows the 2D plots of the relative phase error versus ωx, for the case ωy=0 at h=0.025 and h=0.05, respectively. We observe that the scheme is least dispersive when h=0.05 and k=0.05 and most dispersive when h=0.025 and k=0.025. We have phase lead when h=0.05 and k=0.05 and there is phase lag for the other seven combinations of h and k. In the following section, we present the results of the numerical experiment described in Section 3.2.3D plots of relative phase error versusωx versus ωy at h=0.025 at some different values of k for the LOD (1,5) procedure.
(a)
h
=
0.025
,
k
=
0.005
(b)
h
=
0.025
,
k
=
0.01
(c)
h
=
0.025
,
k
=
0.02
(d)
h
=
0.025
,
k
=
0.0253D plots of relative phase error versusωx versus ωy at h=0.05 at some different values of k for the LOD (1,5) procedure.
(a)
h
=
0.05
,
k
=
0.01
(b)
h
=
0.05
,
k
=
0.02
(c)
h
=
0.05
,
k
=
0.03
(d)
h
=
0.05
,
k
=
0.05Plots of relative phase error versusωx when ωy=0 at h=0.025 and h=0.05 with some different values of k for the LOD (1,5) procedure.
(a)
h
=
0.025
(b)
h
=
0.05


## 6. Numerical Results

### 6.1. LOD Lax-Wendroff Procedure

The results of the numerical experiment at some values ofk using the LOD Lax-Wendroff procedure at T=0.3 at h=0.025 and h=0.05 are shown in Tables 1 and 2, respectively. We observe that the dispersion error is significantly greater than the dissipation error. Out of the five combinations of values of k and h, we observe that the dispersion error and total mean square error are both least when h=0.05 and k=0.03 and are both greatest when h=0.05 and k=0.0025.Table 1
Errors obtained from LOD Lax-Wendroff ath=0.025.
k
CFL
E
num
Total mean square error
DISP. ERROR
DISS. ERROR
0.0025
0.08
1.2043
×
1
0
-
3
1.0930
×
1
0
-
5
1.0602
×
1
0
-
5
2.5787
×
1
0
-
7
0.005
0.16
9.4236
×
1
0
-
4
7.4519
×
1
0
-
6
7.2355
×
1
0
-
6
2.4370
×
1
0
-
7
0.01
0.32
4.1557
×
1
0
-
4
4.1970
×
1
0
-
6
4.0222
×
1
0
-
6
1.6465
×
1
0
-
7
3
/
160
0.6
1.7126
×
1
0
-
3
2.0913
×
1
0
-
5
2.0819
×
1
0
-
5
5.5626
×
1
0
-
8
1
/
60
8
/
15
1.3626
×
1
0
-
3
1.3761
×
1
0
-
5
1.3675
×
1
0
-
5
8.3438
×
1
0
-
8Table 2
Errors obtained from LOD Lax-Wendroff ath=0.05.
k
CFL
E
num
Total mean square error
DISP. ERROR
DISS. ERROR
0.0025
0.04
1.3398
×
1
0
-
3
4.1073
×
1
0
-
5
4.1056
×
1
0
-
5
2.1203
×
1
0
-
8
0.005
0.08
1.2483
×
1
0
-
3
3.5560
×
1
0
-
5
3.5562
×
1
0
-
5
3.5412
×
1
0
-
9
0.01
0.16
1.0547
×
1
0
-
3
2.5261
×
1
0
-
5
2.5255
×
1
0
-
5
8.1665
×
1
0
-
9
1
/
60
4
/
15
7.7546
×
1
0
-
4
1.3488
×
1
0
-
5
1.3433
×
1
0
-
5
5.4469
×
1
0
-
8
3
/
160
0.3
6.8350
×
1
0
-
4
1.0388
×
1
0
-
5
1.0323
×
1
0
-
5
6.9368
×
1
0
-
8
0.02
0.32
6.2698
×
1
0
-
4
8.6851
×
1
0
-
6
8.6112
×
1
0
-
6
7.7093
×
1
0
-
8
0.025
0.4
3.9085
×
1
0
-
4
3.2351
×
1
0
-
6
3.1467
×
1
0
-
6
9.3072
×
1
0
-
8
1
/
34
8
/
17
4.5965
×
1
0
-
4
5.2537
×
1
0
-
6
5.2510
×
1
0
-
6
3.708
×
1
0
-
9
0.03
0.48
1.4596
×
1
0
-
4
4.6945
×
1
0
-
7
3.9599
×
1
0
-
7
7.7350
×
1
0
-
8
0.04
0.64
1.1718
×
1
0
-
3
2.4062
×
1
0
-
5
2.3462
×
1
0
-
5
6.0153
×
1
0
-
7


### 6.2. LOD(1,5) Procedure

The results of the numerical experiment using the LOD(1,5) procedure at T=0.3 for h=0.025 and h=0.05 are shown in Tables 3 and 4, respectively. Out of the 12 combinations of values of k and h, we observe that the dissipation, dispersion error, and error rate are all least when h=0.05 and k=0.05. Also, the dissipation and dispersion errors and error rate are greatest when h=0.05 and k=0.0625. We observe that at h=0.025 the total mean square error, dispersion error, and dissipation error are not much affected by the values of k. Also, the dispersion error is greater than the at all values of h and k considered.Table 3
Errors obtained from the LOD (1,5) procedure at h=0.025.
k
CFL
E
num
Total mean square error
DISP. ERROR
DISS. ERROR
0.005
0.16
4.0539
×
1
0
-
4
4.1593
×
1
0
-
6
4.0009
×
1
0
-
6
1.6455
×
1
0
-
7
0.01
0.32
4.0072
×
1
0
-
4
4.1573
×
1
0
-
6
3.9874
×
1
0
-
6
1.6540
×
1
0
-
7
1
/
60
8
/
15
4.0565
×
1
0
-
4
4.1577
×
1
0
-
6
3.9797
×
1
0
-
6
1.6390
×
1
0
-
7
3
/
160
0.60
3.9839
×
1
0
-
4
4.1494
×
1
0
-
6
3.9782
×
1
0
-
6
1.6268
×
1
0
-
7
0.02
0.64
3.9878
×
1
0
-
4
4.1481
×
1
0
-
6
3.9768
×
1
0
-
6
1.6174
×
1
0
-
7
0.025
0.80
4.0890
×
1
0
-
4
4.1640
×
1
0
-
6
4.0430
×
1
0
-
6
1.5618
×
1
0
-
7Table 4
Errors obtained from the LOD (1,5) procedure at h=0.05.
k
CFL
E
num
Total mean square error
DISP. ERROR
DISS. ERROR
0.01
0.16
2.0000
×
1
0
-
4
7.8149
×
1
0
-
7
7.7377
×
1
0
-
7
1.1297
×
1
0
-
8
0.02
0.32
1.3905
×
1
0
-
4
3.9277
×
1
0
-
7
3.8065
×
1
0
-
7
1.4712
×
1
0
-
8
1
/
34
8/17
4.1603
×
1
0
-
4
4.0763
×
1
0
-
6
4.0663
×
1
0
-
6
1.1482
×
1
0
-
8
0.03
0.48
9.5979
×
1
0
-
5
2.1598
×
1
0
-
7
2.0584
×
1
0
-
7
1.4193
×
1
0
-
8
0.05
0.8
7.9301
×
1
0
-
5
1.7934
×
1
0
-
7
1.7059
×
1
0
-
7
1.0449
×
1
0
-
8
0.0625
1
8.5551
×
1
0
-
4
1.4674
×
1
0
-
5
1.4310
×
1
0
-
5
3.6790
×
1
0
-
7


## 6.1. LOD Lax-Wendroff Procedure

The results of the numerical experiment at some values ofk using the LOD Lax-Wendroff procedure at T=0.3 at h=0.025 and h=0.05 are shown in Tables 1 and 2, respectively. We observe that the dispersion error is significantly greater than the dissipation error. Out of the five combinations of values of k and h, we observe that the dispersion error and total mean square error are both least when h=0.05 and k=0.03 and are both greatest when h=0.05 and k=0.0025.Table 1
Errors obtained from LOD Lax-Wendroff ath=0.025.
k
CFL
E
num
Total mean square error
DISP. ERROR
DISS. ERROR
0.0025
0.08
1.2043
×
1
0
-
3
1.0930
×
1
0
-
5
1.0602
×
1
0
-
5
2.5787
×
1
0
-
7
0.005
0.16
9.4236
×
1
0
-
4
7.4519
×
1
0
-
6
7.2355
×
1
0
-
6
2.4370
×
1
0
-
7
0.01
0.32
4.1557
×
1
0
-
4
4.1970
×
1
0
-
6
4.0222
×
1
0
-
6
1.6465
×
1
0
-
7
3
/
160
0.6
1.7126
×
1
0
-
3
2.0913
×
1
0
-
5
2.0819
×
1
0
-
5
5.5626
×
1
0
-
8
1
/
60
8
/
15
1.3626
×
1
0
-
3
1.3761
×
1
0
-
5
1.3675
×
1
0
-
5
8.3438
×
1
0
-
8Table 2
Errors obtained from LOD Lax-Wendroff ath=0.05.
k
CFL
E
num
Total mean square error
DISP. ERROR
DISS. ERROR
0.0025
0.04
1.3398
×
1
0
-
3
4.1073
×
1
0
-
5
4.1056
×
1
0
-
5
2.1203
×
1
0
-
8
0.005
0.08
1.2483
×
1
0
-
3
3.5560
×
1
0
-
5
3.5562
×
1
0
-
5
3.5412
×
1
0
-
9
0.01
0.16
1.0547
×
1
0
-
3
2.5261
×
1
0
-
5
2.5255
×
1
0
-
5
8.1665
×
1
0
-
9
1
/
60
4
/
15
7.7546
×
1
0
-
4
1.3488
×
1
0
-
5
1.3433
×
1
0
-
5
5.4469
×
1
0
-
8
3
/
160
0.3
6.8350
×
1
0
-
4
1.0388
×
1
0
-
5
1.0323
×
1
0
-
5
6.9368
×
1
0
-
8
0.02
0.32
6.2698
×
1
0
-
4
8.6851
×
1
0
-
6
8.6112
×
1
0
-
6
7.7093
×
1
0
-
8
0.025
0.4
3.9085
×
1
0
-
4
3.2351
×
1
0
-
6
3.1467
×
1
0
-
6
9.3072
×
1
0
-
8
1
/
34
8
/
17
4.5965
×
1
0
-
4
5.2537
×
1
0
-
6
5.2510
×
1
0
-
6
3.708
×
1
0
-
9
0.03
0.48
1.4596
×
1
0
-
4
4.6945
×
1
0
-
7
3.9599
×
1
0
-
7
7.7350
×
1
0
-
8
0.04
0.64
1.1718
×
1
0
-
3
2.4062
×
1
0
-
5
2.3462
×
1
0
-
5
6.0153
×
1
0
-
7


## 6.2. LOD(1,5) Procedure

The results of the numerical experiment using the LOD(1,5) procedure at T=0.3 for h=0.025 and h=0.05 are shown in Tables 3 and 4, respectively. Out of the 12 combinations of values of k and h, we observe that the dissipation, dispersion error, and error rate are all least when h=0.05 and k=0.05. Also, the dissipation and dispersion errors and error rate are greatest when h=0.05 and k=0.0625. We observe that at h=0.025 the total mean square error, dispersion error, and dissipation error are not much affected by the values of k. Also, the dispersion error is greater than the at all values of h and k considered.Table 3
Errors obtained from the LOD (1,5) procedure at h=0.025.
k
CFL
E
num
Total mean square error
DISP. ERROR
DISS. ERROR
0.005
0.16
4.0539
×
1
0
-
4
4.1593
×
1
0
-
6
4.0009
×
1
0
-
6
1.6455
×
1
0
-
7
0.01
0.32
4.0072
×
1
0
-
4
4.1573
×
1
0
-
6
3.9874
×
1
0
-
6
1.6540
×
1
0
-
7
1
/
60
8
/
15
4.0565
×
1
0
-
4
4.1577
×
1
0
-
6
3.9797
×
1
0
-
6
1.6390
×
1
0
-
7
3
/
160
0.60
3.9839
×
1
0
-
4
4.1494
×
1
0
-
6
3.9782
×
1
0
-
6
1.6268
×
1
0
-
7
0.02
0.64
3.9878
×
1
0
-
4
4.1481
×
1
0
-
6
3.9768
×
1
0
-
6
1.6174
×
1
0
-
7
0.025
0.80
4.0890
×
1
0
-
4
4.1640
×
1
0
-
6
4.0430
×
1
0
-
6
1.5618
×
1
0
-
7Table 4
Errors obtained from the LOD (1,5) procedure at h=0.05.
k
CFL
E
num
Total mean square error
DISP. ERROR
DISS. ERROR
0.01
0.16
2.0000
×
1
0
-
4
7.8149
×
1
0
-
7
7.7377
×
1
0
-
7
1.1297
×
1
0
-
8
0.02
0.32
1.3905
×
1
0
-
4
3.9277
×
1
0
-
7
3.8065
×
1
0
-
7
1.4712
×
1
0
-
8
1
/
34
8/17
4.1603
×
1
0
-
4
4.0763
×
1
0
-
6
4.0663
×
1
0
-
6
1.1482
×
1
0
-
8
0.03
0.48
9.5979
×
1
0
-
5
2.1598
×
1
0
-
7
2.0584
×
1
0
-
7
1.4193
×
1
0
-
8
0.05
0.8
7.9301
×
1
0
-
5
1.7934
×
1
0
-
7
1.7059
×
1
0
-
7
1.0449
×
1
0
-
8
0.0625
1
8.5551
×
1
0
-
4
1.4674
×
1
0
-
5
1.4310
×
1
0
-
5
3.6790
×
1
0
-
7


## 7. Optimization

In this section, we obtain the optimal value ofk at h=0.025 that minimizes the dispersion error for the two time-splitting procedures.Since the partial differential equation we consider is slightly dissipative and also we observe from the numerical experiments carried out that the dissipative errors are much less than the dispersive errors, we choose to minimize the square of the dispersion error of the two splitting schemes.


### 7.1. Proposed Techniques of Optimization

Tam and Webb [18], Bogey and Bailly [19], and Hixon [20] among others have implemented techniques which enable coefficients to be determined in numerical schemes specifically designed for computational aeroacoustics. We now describe briefly how Tam and Webb [18] and Bogey and Bailly [19] define their measures and consequently their technique of optimization in computational aeroacoustics.The dispersion-relation-preserving (DRP) scheme was designed so that the dispersion relation of the finite difference scheme is formally the same as that of the original partial differential equations. The integrated error is defined as(83)E=∫-ηη‍|θ*h-θh|2d(θh),
where the quantities θ*h and θh represent the numerical and exact wavenumbers, respectively. The dispersion error and dissipation error are calculated as |ℜ(θ*h)-θh| and |Im(θ*h)|, respectively.Tam and Shen [21] set η as 1.1 and optimize the coefficients in the numerical scheme such that the integrated error is minimized.Bogey and Bailly minimize the relative difference between the exact wavenumberθh and the effective/numerical wavenumber θ*h and define their integrated errors as
(84)E=∫(θh)l(θh)h‍|θ*h-θh|θhd(θh)
or
(85)E=∫ln(θh)lln(θh)h‍|θ*h-θh|d(ln(θh)).In computational fluid dynamics for a particular method under consideration, the dispersion error is calculated as(86)|1-RPE|.We have modified the measures used by Tam and Webb and Bogey and Bailly in a computational aeroacoustics framework to suit them in a computational fluid dynamics framework [22] such that the optimal parameter can be obtained. We have defined the following integrated errors integrated error from Tam and Webb (IETAM) integrated error from Bogey and Bailly (IEBOGEY) [22] as follows:
(87)IETAM=∫0w1‍|1-RPE|2dw,IEBOGEY=∫0w1‍|1-RPE|dw.In [8], the integrated error for a scheme discretising the 1D advection-diffusion equation
(88)∂u∂t+∂u∂x=0.01∂2u∂x2
is obtained as
(89)∫01.1‍(RPE-1)2dω.
The value of h was fixed as 0.02 and the range of values of k, was determined. Then the integrated error, which is a function of k, was minimized and the optimal value of k is determined using NLPSolve function in maple.We extend the work on optimization of parameters in [8] which is for the 1D advection-diffusion equation for the case of the 2D advection-diffusion equation. We define the integrated error as
(90)∫01.1‍∫01.1‍(RPE-1)2dωxdωy.
We first obtain an expression for the RPE of the LOD Lax-Wendroff when discretizing the equation
(91)∂u∂t+0.8∂u∂x+0.8∂u∂y=0.01∂2u∂x2+0.01∂2u∂y2.Sincec=0.8k/h and s=0.8k/h2 and we choose h=0.025, we have c=32k and s=16k. Hence, the RPE is a function of k, ωx, and ωy. Since we can have phase wrapping, we make use of Taylor’s expansion to obtain an approximation for the RPE up to the terms
(92)(ωx)3(ωy)2,(ωx)5,(ωx)2(ωy)3,(ωy)5,ωx(ωy)4,(ωx)4ωy.The integrated error∫01.1‍∫01.1‍(RPE-1)2dωxdωy is obtained by using Simpson's method and it is a function of k only. A plot of the integrated error versus k is shown in Figure 13(a) for k∈[0,0.0193]. Using the NLPSolve function in maple, this optimal value of k is found to be 0.009593 and the minimum value of the integrated error is 1.883960×10-7.Plots of integrated errors versusk at h=0.025 for the two time-splitting schemes.
(a)
LOD Lax-Wendroff
(b)
LOD(1,5) procedureTo validate our results, we perform the same numerical experiment described in Section3.2 at h=0.025 and use the optimal value of k or a value of k close to this optimal value, in that case, k=3/310≈0.0096, and compute the errors. The error rate, total mean square error, and dispersion error are 3.9967×10-4, 4.1603×10-6, 3.9878×10-6. These three errors are all least as compared to when other values of k are used as shown in Table 1.We adopt the same procedure to compute the optimal value ofk for the LOD (1,5) scheme when h=0.025. We obtain an approximate expression for the RPE of the LOD (1,5) when c=32k and s=16k. We use Simpson's rule to approximate the integral given by (90) which is a function of k. A plot of the integrated error versus k for k∈[0,0.026288] is shown in Figure 13(b) and using NLPSolve function in maple the optimal value of k is 0.013782 and also the minimum value of the integral is 1.139313×10-6.We perform the numerical experiment described in Section3.2 with k=3/220 which is close to the optimal value of k we have obtained with h=0.025. The error rate, total mean square error, dispersion error, and dissipation error are 3.9948×10-4, 4.1504×10-6, 3.9459×10-6, and 1.6501×10-7, respectively. The total mean square error and dispersion error are both least when k=3/220≈0.01378.


## 7.1. Proposed Techniques of Optimization

Tam and Webb [18], Bogey and Bailly [19], and Hixon [20] among others have implemented techniques which enable coefficients to be determined in numerical schemes specifically designed for computational aeroacoustics. We now describe briefly how Tam and Webb [18] and Bogey and Bailly [19] define their measures and consequently their technique of optimization in computational aeroacoustics.The dispersion-relation-preserving (DRP) scheme was designed so that the dispersion relation of the finite difference scheme is formally the same as that of the original partial differential equations. The integrated error is defined as(83)E=∫-ηη‍|θ*h-θh|2d(θh),
where the quantities θ*h and θh represent the numerical and exact wavenumbers, respectively. The dispersion error and dissipation error are calculated as |ℜ(θ*h)-θh| and |Im(θ*h)|, respectively.Tam and Shen [21] set η as 1.1 and optimize the coefficients in the numerical scheme such that the integrated error is minimized.Bogey and Bailly minimize the relative difference between the exact wavenumberθh and the effective/numerical wavenumber θ*h and define their integrated errors as
(84)E=∫(θh)l(θh)h‍|θ*h-θh|θhd(θh)
or
(85)E=∫ln(θh)lln(θh)h‍|θ*h-θh|d(ln(θh)).In computational fluid dynamics for a particular method under consideration, the dispersion error is calculated as(86)|1-RPE|.We have modified the measures used by Tam and Webb and Bogey and Bailly in a computational aeroacoustics framework to suit them in a computational fluid dynamics framework [22] such that the optimal parameter can be obtained. We have defined the following integrated errors integrated error from Tam and Webb (IETAM) integrated error from Bogey and Bailly (IEBOGEY) [22] as follows:
(87)IETAM=∫0w1‍|1-RPE|2dw,IEBOGEY=∫0w1‍|1-RPE|dw.In [8], the integrated error for a scheme discretising the 1D advection-diffusion equation
(88)∂u∂t+∂u∂x=0.01∂2u∂x2
is obtained as
(89)∫01.1‍(RPE-1)2dω.
The value of h was fixed as 0.02 and the range of values of k, was determined. Then the integrated error, which is a function of k, was minimized and the optimal value of k is determined using NLPSolve function in maple.We extend the work on optimization of parameters in [8] which is for the 1D advection-diffusion equation for the case of the 2D advection-diffusion equation. We define the integrated error as
(90)∫01.1‍∫01.1‍(RPE-1)2dωxdωy.
We first obtain an expression for the RPE of the LOD Lax-Wendroff when discretizing the equation
(91)∂u∂t+0.8∂u∂x+0.8∂u∂y=0.01∂2u∂x2+0.01∂2u∂y2.Sincec=0.8k/h and s=0.8k/h2 and we choose h=0.025, we have c=32k and s=16k. Hence, the RPE is a function of k, ωx, and ωy. Since we can have phase wrapping, we make use of Taylor’s expansion to obtain an approximation for the RPE up to the terms
(92)(ωx)3(ωy)2,(ωx)5,(ωx)2(ωy)3,(ωy)5,ωx(ωy)4,(ωx)4ωy.The integrated error∫01.1‍∫01.1‍(RPE-1)2dωxdωy is obtained by using Simpson's method and it is a function of k only. A plot of the integrated error versus k is shown in Figure 13(a) for k∈[0,0.0193]. Using the NLPSolve function in maple, this optimal value of k is found to be 0.009593 and the minimum value of the integrated error is 1.883960×10-7.Plots of integrated errors versusk at h=0.025 for the two time-splitting schemes.
(a)
LOD Lax-Wendroff
(b)
LOD(1,5) procedureTo validate our results, we perform the same numerical experiment described in Section3.2 at h=0.025 and use the optimal value of k or a value of k close to this optimal value, in that case, k=3/310≈0.0096, and compute the errors. The error rate, total mean square error, and dispersion error are 3.9967×10-4, 4.1603×10-6, 3.9878×10-6. These three errors are all least as compared to when other values of k are used as shown in Table 1.We adopt the same procedure to compute the optimal value ofk for the LOD (1,5) scheme when h=0.025. We obtain an approximate expression for the RPE of the LOD (1,5) when c=32k and s=16k. We use Simpson's rule to approximate the integral given by (90) which is a function of k. A plot of the integrated error versus k for k∈[0,0.026288] is shown in Figure 13(b) and using NLPSolve function in maple the optimal value of k is 0.013782 and also the minimum value of the integral is 1.139313×10-6.We perform the numerical experiment described in Section3.2 with k=3/220 which is close to the optimal value of k we have obtained with h=0.025. The error rate, total mean square error, dispersion error, and dissipation error are 3.9948×10-4, 4.1504×10-6, 3.9459×10-6, and 1.6501×10-7, respectively. The total mean square error and dispersion error are both least when k=3/220≈0.01378.


## 8. Conclusion

In this paper, two time-splitting procedures are used to solve a 2D advection-diffusion equation with constant coefficients when the advection velocity in bothx- and y-directions is 0.8 and also when the coefficient of diffusivity in both x- and y-directions is 0.01. We perform a stability analysis and spectral analysis of the dispersion and dissipation properties of the two schemes at some values of h and k. Numerical experiments are carried out and various errors are computed. These errors are dependent on the values of h and k. It is observed that in general the dispersion error is more affected by the values of k and h for the LOD Lax-Wendroff scheme as compared to that of the LOD (1,5) scheme at a given value of h. We then use an optimization technique based on minimisation of the square of the dispersion error to find the optimal value of k when h is chosen as 0.025 and this is validated by numerical experiments.Future extension of this work to consider other types of advection-diffusion equations when dissipation dominates and to find out which optimization techniques are suitable in these cases. Also, the work can be extended to 2D nonlinear convection-diffusion problems.



---

*Source: 634657-2013-11-17.xml*