**A Fourth-Order Compact Finite Difference Scheme for Solving Unsteady Convection-Diffusion Equations**

Wenyuan Liao<sup>1</sup> and Jianping Zhu<sup>2</sup>

<sup>1</sup>*University of Calgary* <sup>2</sup>*University of Texas at Arlington* <sup>1</sup>*Canada* <sup>2</sup>*USA*

#### **1. Introduction**

80 Computational Simulations and Applications

Akbari, M. H. & Price, S. J. (2005). Numerical investigation of flow patterns for staggered cylinder pairs in cross-flow, *Journal of Fluids and Structures,* Vol.20, pp. 533-554 Cheng, M.; Liu, G. R. & Lam, K. Y. (2001a). Numerical simulation of flow past a rotationally

Cheng, M.; Chew, Y. T. & Luo, S. C. (2001b). Numerical investigation of a rotationally oscillating cylinder in mean flow*, Journal of Fluids and Structures*, Vol.15, pp. 981-1007 Chorin, A. (1968). Numerical solution of the Navier-Stokes equations*, Math. Comp.* , Vol.22,

Deng, J.; Ren, A-L.; Zou, J-F. & Shao, X-M. (2006). Three-dimensional flow around two circular cylinders in tandem arrangement, *Fluid Dynamics Research*, Vol.38, pp. 386-404 Ferziger, J. H. & Peric, M. (2002). *Computational Methods for Fluid Dynamics, 3.* Ed. Berlin

Kang, S.; Choi, H. & Lee, S. (1999). Laminar Flow past a Rotating Circular Cylinder*, Physics* 

Lima e Silva, A. L. F.; Silva, A. R. & Silveira-Neto, A. (2007). Numerical Simulation of Two-

Naudascher, E. & Rockwell, D. (1994). *Flow-Induced Vibrations: An Engineering Guide,* Dover

Peskin, C. S. (1977). Numerical Analysis of Blood Flow in the Heart*, Journal of Computational* 

Peskin, C. S. & McQueen, D. M. (1994). A General Method for the Computer Simulation of

Reynolds, O. (1894). On the dynamical theory of incompressible viscous fluids and the

Schneider, G. E. & Zedan, M. (1981). A Modified Strongly Implicit Procedure for the Numerical Solution of Field Problems*, Numerical Heat Transfer*, Vol.4, No.1 Silva, A. R.; Silveira-Neto, A.; Rade, D. A.; Francis, R. & Santos, E. A. (2009). Numerical

Smagorinsky, J. (1963). General Circulation Experiments with Primitive Equations*, Mon.* 

Sumner, D.; Price, S. J. & Païdoussis, M. P. (1999). Tandem cylinders in impulsively started

Sumner, D.; Richards, M. D. & Okasile, O. O. (2005). Two staggered circular cylinders of equal diameter in cross-flow*, J. of Fluids and Structures,* Vol. 20, pp. 255-276 Sumner, D.; Richards, M. D. & Okasile, O. O. (2008). Strouhal number data for two

Tuszynski, J. & Löhner, R. (1998). Control of a Kármán Vortex Flow by Rotational Oscillations of a Cylinder*, George Mason University, USA,* Vol.15, pp. 1-12 Williamson, C. H. K. & Jauvtis, N. (2004). A high-amplitude 2T mode of vortex-induced

staggered circular cylinders*, Journal of Wind Engineering and Industrial Aerodynamics,*

vibration for a light body in X-Y motion*, European Journal of Mechanics B/Fluids,* 

flow*, Journal of Fluids and Structures,* Vol.13, pp. 955-965

Method*, J. of the Braz. Soc. Of Mech. Sci. & Eng.*, Vol.4, pp. 378-386

Publications, Inc. Mineola, New York*,* pp. 414

Dimensional Complex Flows around Bluff Bodies using the Immersed Boundary

Biological Systems Interacting with Fluids*, SEB Symphosium on Biological Fluid* 

determination of the criterion, *Philosophical Transaction of the Royal Society of London,*

Simulations of Flows over a Pair of Cylinders at Different Arrangements using the Immersed Boundary Method*, CMES: Computer Modeling in Engineering & Sciences,* 

oscillating cylinder*, Computer & Fluids*, Vol.30, pp. 365-392

**8. References** 

pp. 745-762

Springer-Verlag, pp. 423

*of Fluid*, Vol.11, pp. 3312-3320

*Physics*, Vol.25, pp. 220-252

Vol.186, Part I, pp. 122-164

Vol.50, No.3, pp. 285-303

Vol. 96, pp. 859-871

Vol.23, pp. 107-114

*Weather Rev.*, Vol.91, pp. 99-164

*Dynamics*, Leeds, England, (July 5-8)

Convection-diffusion equations are widely used for modeling and simulations of various complex phenomena in science and engineering (Hundsdorfer & Verwer, 2003; Morton, 1996). Since for most application problems it is impossible to solve convection-diffusion equations analytically, efficient numerical algorithms are becoming increasingly important to numerical simulations involving convection-diffusion equations.

Recently a great deal of efforts have been devoted to developing high-order compact schemes, which utilize only the grid nodes directly adjacent to the central node. In (Noye & Tan, 1989), Noye and Tan derived a class of high-order implicit schemes for solving the one-dimensional unsteady convection-diffusion equations. This method is very stable and accurate (third-order in space and second-order in time). In (Gupta et al., 1984), a fourth-order finite difference scheme for a steady convection-diffusion equation with variable coefficients was proposed. The scheme is defined on a single square cell of size 2Δ*x* over a nine-point stencil. In (Rigal, 1994), Rigal provided an extensive analysis of the properties of a class of two- and three-level second-order difference schemes which have been proposed in (Rigal, 1989; 1990). In (Spotz & Carey, 2001), the two-dimensional HOC (High Order Compact) scheme proposed in (Gupta et al., 1984) was extended to solve unsteady one-dimensional convection-diffusion equations with variable coefficients and two-dimensional diffusion equations. This method was further extended by Kalita et al. in (Kalita et al., 2002) to a class of HOC schemes with weighted time discretization, and successfully used to solve unsteady two-dimensional convection-diffusion equations. In (Karaa & Zhang, 2004), Karaa and Zhang proposed a novel high-order alternating direction implicit method, based on the technique in (Zhang et al., 2002), for solving unsteady two-dimensional convection-diffusion problems. This new method is second-order in time and fourth-order in space, and is computationally efficient. In (Tian & Dai, 2007), Tian and Dai proposed a class of high-order compact exponential finite difference methods for solving one- and two-dimensional steady-state convection-diffusion problems. This method is nonoscillatory, fourth-order in space, and easy to implement. Some more recent high-order ADI methods for unsteady convection-diffusion equations can be found in (Tian & Ge, 2007; You, 2006).

of convection-diffusion equations. It is unconditionally stable, and is particularly suitable for problems that require the solution of both *u* and *ux*. For example, when solving the

<sup>83</sup> A Fourth-Order Compact Finite Difference Scheme

*<sup>∂</sup>S*<sup>2</sup> <sup>+</sup> *rS <sup>∂</sup><sup>V</sup>*

both the solution *V* and its derivative *VS* are desired. The solution *V* is the price of an option and its derivative *VS* is called the hedge delta that represents the sensitivity of the option value

The rest of this chapter is organized as follows: The description of the new method is given in the next section. Proof of unconditional stability of the new method is given in Section 3. Computational complexity of this new method is analyzed and compared with upwind and standard central finite difference schemes in Section 4. Several numerical examples are

In this section, we will first outline the new algorithm that eliminates the convection term in a convection-diffusion equation to facilitate the use of central finite difference schemes and then discuss how the initial and boundary conditions are handled using the new algorithm.

Eq. (5) and (7) now form a system of equations for *u* and *v*. They only involve diffusion term *uxx* and *vxx*, which can be discretized by the standard central finite difference schemes. If the 3-point second-order central difference scheme is used, the discretized equations will form a block tri-diagonal algebraic equation system with 2 × 2 blocks. For nonlinear *f*(*u*), the Newton's method or its variations can be used to solve the nonlinear system of algebraic

> *δ*2*un*+<sup>1</sup> *i* Δ*x*2(1 + *<sup>δ</sup>*<sup>2</sup>

*δ*2*un*+<sup>1</sup> *i* Δ*x*2(1 + *<sup>δ</sup>*<sup>2</sup>

+[(*cx*)*<sup>i</sup>* + *fu*(*un*+<sup>1</sup>

<sup>12</sup> )

12 ) + *a*

If the fourth-order Padé approximation is used for Eqs. (5) and (7), we have

*<sup>∂</sup><sup>S</sup>* <sup>−</sup> *rV* <sup>=</sup> 0, (4)

*ut* = *auxx* + *c*(*x*)*v* + *f*(*u*) (5)

(*ux*)*<sup>t</sup>* = *a*(*uxx*)*<sup>x</sup>* + *cx*(*x*)*ux* + *c*(*x*)(*ux*)*<sup>x</sup>* + *fu*(*u*)*ux* (6)

*vt* = *c*(*x*)*uxx* + *avxx* + *cx*(*x*)*v* + *fu*(*u*)*v* (7)

+ *civi* + *f*(*un*+<sup>1</sup>

*<sup>i</sup>* )]*vn*+<sup>1</sup>

*δ*2*vn*+<sup>1</sup> *i* Δ*x*2(1 + *<sup>δ</sup>*<sup>2</sup>

12 )

*<sup>i</sup>* ), (8)

*<sup>i</sup>* . (9)

*<sup>σ</sup>*2*S*<sup>2</sup> *<sup>∂</sup>*2*<sup>V</sup>*

Black-Scholes option pricing model(Seydel, 2002)

for Solving Unsteady Convection-Diffusion Equations

to the change of the underlining stock price.

**2. The new method**

equations.

**2.1 Description of the new method** Setting *v* = *ux* in Eq. (1), we have

which, considering *v* = *ux*, can be written as

*un*+<sup>1</sup> *<sup>i</sup>* <sup>−</sup> *<sup>u</sup><sup>n</sup> i* <sup>Δ</sup>*<sup>t</sup>* <sup>=</sup> *<sup>a</sup>*

*vn*+<sup>1</sup> *<sup>i</sup>* <sup>−</sup> *<sup>v</sup><sup>n</sup> i* <sup>Δ</sup>*<sup>t</sup>* <sup>=</sup> *ci*

*∂V <sup>∂</sup><sup>t</sup>* <sup>+</sup> 1 2

presented in Section 5, followed by conclusions in Section 6.

Differentiating both sides of Eq. (1) with respect to *x* leads to

For simplicity, we will use the following one-dimensional equation to describe the new method developed in this chapter:

$$u\_t = a u\_{xx} + c(x)u\_x + f(u)\_{\prime} \tag{1}$$

$$0 < x < 1, 0 < t \le T,$$

$$u(0, t) = g\_1(t), \\ u(1, t) = g\_2(t),$$

$$u(x, 0) = h(x),$$

where *a* is a constant. Extensions to more complicated equations in one-dimension(for example when *a* is not a constant) are straightforward; and extensions to higher-dimensional equations will be briefly discussed in Section 6.

The existence of the convection term in Eq. (1) creates several difficulties when the equation is solved numerically using the finite difference schemes. It is well-known (Hundsdorfer & Verwer, 2003; Morton, 1996) that the convection term needs to be discretized using proper upwind finite difference schemes to avoid oscillations in convection dominated problems. If the sign of *c*(*x*) changes over the solution domain, the direction of the upwind scheme must also be changed accordingly. The order of accuracy of the upwind schemes is usually lower than the central difference schemes on the same finite difference stencil.

In addition, the convection term in Eq. (1) also makes it more difficult to use the fourth-order Padé approximations. For reaction-diffusion equations (*c*(*x*) = 0), the Padé approximation can be used to achieve fourth-order accuracy on a 3-point stencil typically used for the standard second-order algorithms(Gu et al., 2003):

$$\frac{u\_i^{n+1} - u\_i^n}{\Delta t} = a \frac{\delta^2 u\_i^{n+1}}{\Delta x^2 (1 + \frac{\delta^2}{12})} + f(u\_i^{n+1}),\tag{2}$$

where *i* = 0, 1, ··· , *M*, and *n* = 0, 1, ··· , *N* are indices for spatial and temporal grid points, respectively, *u<sup>n</sup> <sup>i</sup>* is the numerical approximation to the exact solution *u*(*xi*, *t <sup>n</sup>*), and the central difference operator *δ*<sup>2</sup> is defined as *δ*2*u<sup>n</sup> <sup>i</sup>* <sup>=</sup> *<sup>u</sup><sup>n</sup> <sup>i</sup>*−<sup>1</sup> <sup>−</sup> <sup>2</sup>*u<sup>n</sup> <sup>i</sup>* <sup>+</sup> *<sup>u</sup><sup>n</sup> <sup>i</sup>*+1. Applying the operator (<sup>1</sup> <sup>+</sup> *<sup>δ</sup>*<sup>2</sup> <sup>12</sup> ) to both sides of Eq. (2), we have

$$(1 + \frac{\delta^2}{12})\frac{u\_i^{n+1} - u\_i^n}{\Delta t} = a\frac{\delta^2 u\_i^{n+1}}{\Delta x^2} + (1 + \frac{\delta^2}{12})f(u\_i^{n+1}).\tag{3}$$

Eq. (3) is fourth-order accurate in space but contains only the second-order central difference operator *δ*<sup>2</sup> that requires only a 3-point stencil. This approach, however, does not work when the convection term is present in the equation. This is because the fourth-order Padé approximation for Eq. (1) leads to the following equation

$$\frac{u\_i^{n+1} - u\_i^n}{\Delta t} = a \frac{\delta^2 u\_i^{n+1}}{\Delta x^2 (1 + \frac{\delta^2}{12})} + c\_l \frac{\delta u\_i^{n+1}}{2 \Delta x (1 + \frac{\delta^2}{6})} + f(u\_i^{n+1})\_\lambda$$

where the central difference operator *δ* is defined as *δu<sup>n</sup> <sup>i</sup>* <sup>=</sup> *<sup>u</sup><sup>n</sup> <sup>i</sup>*+<sup>1</sup> <sup>−</sup> *<sup>u</sup><sup>n</sup> <sup>i</sup>*−1. This equation cannot be simplified to an equation that is defined on a 3-point stencil in the same way as Eq. (2) is reduced to Eq. (3).

The new method discussed in this chapter eliminates the convection term *ux* from Eq. (1) and solves *v* = *ux* directly along with *u*. This makes it possible to use central finite difference schemes and higher-order Padé approximations for accurate and efficient numerical solutions of convection-diffusion equations. It is unconditionally stable, and is particularly suitable for problems that require the solution of both *u* and *ux*. For example, when solving the Black-Scholes option pricing model(Seydel, 2002)

$$\frac{\partial V}{\partial t} + \frac{1}{2}\sigma^2 S^2 \frac{\partial^2 V}{\partial S^2} + rS \frac{\partial V}{\partial S} - rV = 0,\tag{4}$$

both the solution *V* and its derivative *VS* are desired. The solution *V* is the price of an option and its derivative *VS* is called the hedge delta that represents the sensitivity of the option value to the change of the underlining stock price.

The rest of this chapter is organized as follows: The description of the new method is given in the next section. Proof of unconditional stability of the new method is given in Section 3. Computational complexity of this new method is analyzed and compared with upwind and standard central finite difference schemes in Section 4. Several numerical examples are presented in Section 5, followed by conclusions in Section 6.

#### **2. The new method**

2 Will-be-set-by-IN-TECH

For simplicity, we will use the following one-dimensional equation to describe the new

*u*(0, *t*) = *g*1(*t*), *u*(1, *t*) = *g*2(*t*),

where *a* is a constant. Extensions to more complicated equations in one-dimension(for example when *a* is not a constant) are straightforward; and extensions to higher-dimensional

The existence of the convection term in Eq. (1) creates several difficulties when the equation is solved numerically using the finite difference schemes. It is well-known (Hundsdorfer & Verwer, 2003; Morton, 1996) that the convection term needs to be discretized using proper upwind finite difference schemes to avoid oscillations in convection dominated problems. If the sign of *c*(*x*) changes over the solution domain, the direction of the upwind scheme must also be changed accordingly. The order of accuracy of the upwind schemes is usually lower than the central difference schemes on the same finite difference stencil. In addition, the convection term in Eq. (1) also makes it more difficult to use the fourth-order Padé approximations. For reaction-diffusion equations (*c*(*x*) = 0), the Padé approximation can be used to achieve fourth-order accuracy on a 3-point stencil typically used for the

> *δ*2*un*+<sup>1</sup> *i* Δ*x*2(1 + *<sup>δ</sup>*<sup>2</sup>

where *i* = 0, 1, ··· , *M*, and *n* = 0, 1, ··· , *N* are indices for spatial and temporal grid points,

*<sup>i</sup>*−<sup>1</sup> <sup>−</sup> <sup>2</sup>*u<sup>n</sup>*

*δ*2*un*+<sup>1</sup> *i* <sup>Δ</sup>*x*<sup>2</sup> + (<sup>1</sup> <sup>+</sup>

Eq. (3) is fourth-order accurate in space but contains only the second-order central difference operator *δ*<sup>2</sup> that requires only a 3-point stencil. This approach, however, does not work when the convection term is present in the equation. This is because the fourth-order Padé

*<sup>i</sup>* is the numerical approximation to the exact solution *u*(*xi*, *t*

*<sup>i</sup>* <sup>=</sup> *<sup>u</sup><sup>n</sup>*

*δ*2*un*+<sup>1</sup> *i* Δ*x*2(1 + *<sup>δ</sup>*<sup>2</sup>

12 ) + *ci*

be simplified to an equation that is defined on a 3-point stencil in the same way as Eq. (2) is

The new method discussed in this chapter eliminates the convection term *ux* from Eq. (1) and solves *v* = *ux* directly along with *u*. This makes it possible to use central finite difference schemes and higher-order Padé approximations for accurate and efficient numerical solutions

<sup>12</sup> )

*<sup>i</sup>* <sup>+</sup> *<sup>u</sup><sup>n</sup>*

*δun*+<sup>1</sup> *i* 2Δ*x*(1 + *<sup>δ</sup>*<sup>2</sup>

*<sup>i</sup>* <sup>=</sup> *<sup>u</sup><sup>n</sup>*

+ *f*(*un*+<sup>1</sup>

*δ*2 <sup>12</sup> )*f*(*un*+<sup>1</sup>

6 )

*<sup>i</sup>*+<sup>1</sup> <sup>−</sup> *<sup>u</sup><sup>n</sup>*

+ *f*(*un*+<sup>1</sup> *<sup>i</sup>* ),

*<sup>i</sup>* ), (2)

*<sup>i</sup>*+1. Applying the operator (<sup>1</sup> <sup>+</sup> *<sup>δ</sup>*<sup>2</sup>

*<sup>n</sup>*), and the central

*<sup>i</sup>* ). (3)

*<sup>i</sup>*−1. This equation cannot

<sup>12</sup> )

0 < *x* < 1, 0 < *t* ≤ *T*,

*u*(*x*, 0) = *h*(*x*),

*ut* = *auxx* + *c*(*x*)*ux* + *f*(*u*), (1)

method developed in this chapter:

equations will be briefly discussed in Section 6.

standard second-order algorithms(Gu et al., 2003):

difference operator *δ*<sup>2</sup> is defined as *δ*2*u<sup>n</sup>*

(1 + *δ*2 12 )

*un*+<sup>1</sup> *<sup>i</sup>* <sup>−</sup> *<sup>u</sup><sup>n</sup> i* <sup>Δ</sup>*<sup>t</sup>* <sup>=</sup> *<sup>a</sup>*

to both sides of Eq. (2), we have

respectively, *u<sup>n</sup>*

reduced to Eq. (3).

*un*+<sup>1</sup> *<sup>i</sup>* <sup>−</sup> *<sup>u</sup><sup>n</sup> i* <sup>Δ</sup>*<sup>t</sup>* <sup>=</sup> *<sup>a</sup>*

> *un*+<sup>1</sup> *<sup>i</sup>* <sup>−</sup> *<sup>u</sup><sup>n</sup> i* <sup>Δ</sup>*<sup>t</sup>* <sup>=</sup> *<sup>a</sup>*

approximation for Eq. (1) leads to the following equation

where the central difference operator *δ* is defined as *δu<sup>n</sup>*

In this section, we will first outline the new algorithm that eliminates the convection term in a convection-diffusion equation to facilitate the use of central finite difference schemes and then discuss how the initial and boundary conditions are handled using the new algorithm.

#### **2.1 Description of the new method**

Setting *v* = *ux* in Eq. (1), we have

$$u\_t = au\_{xx} + c(\mathbf{x})v + f(u) \tag{5}$$

Differentiating both sides of Eq. (1) with respect to *x* leads to

$$(\mathfrak{u}\_{\mathbf{x}})\_{\mathfrak{l}} = a(\mathfrak{u}\_{\mathbf{xx}})\_{\mathfrak{x}} + c\_{\mathfrak{x}}(\mathfrak{x})\mathfrak{u}\_{\mathbf{x}} + c(\mathfrak{x})(\mathfrak{u}\_{\mathbf{x}})\_{\mathfrak{x}} + f\_{\mathfrak{u}}(\mathfrak{u})\mathfrak{u}\_{\mathbf{x}} \tag{6}$$

which, considering *v* = *ux*, can be written as

$$
\omega v\_t = \mathfrak{c}(\mathfrak{x}) u\_{\text{xx}} + a v\_{\text{xx}} + c\_{\mathfrak{x}}(\mathfrak{x}) v + f\_{\mathfrak{u}}(\mathfrak{u}) v \tag{7}
$$

Eq. (5) and (7) now form a system of equations for *u* and *v*. They only involve diffusion term *uxx* and *vxx*, which can be discretized by the standard central finite difference schemes. If the 3-point second-order central difference scheme is used, the discretized equations will form a block tri-diagonal algebraic equation system with 2 × 2 blocks. For nonlinear *f*(*u*), the Newton's method or its variations can be used to solve the nonlinear system of algebraic equations.

If the fourth-order Padé approximation is used for Eqs. (5) and (7), we have

$$\frac{u\_i^{n+1} - u\_i^n}{\Delta t} = a \frac{\delta^2 u\_i^{n+1}}{\Delta x^2 (1 + \frac{\delta^2}{12})} + c\_i v\_i + f(u\_i^{n+1}),\tag{8}$$

$$\frac{v\_i^{n+1} - v\_i^n}{\Delta t} = c\_i \frac{\delta^2 u\_i^{n+1}}{\Delta x^2 (1 + \frac{\delta^2}{12})} + a \frac{\delta^2 v\_i^{n+1}}{\Delta x^2 (1 + \frac{\delta^2}{12})}$$

$$+ [(c\_x)\_i + f\_u(u\_i^{n+1})] v\_i^{n+1}.\tag{9}$$

Direct integration of Eq. (1): This approach is particularly convenient for special cases of Eq. (1). For example, when dealing with steady state equation with a constant convection

<sup>85</sup> A Fourth-Order Compact Finite Difference Scheme

where *a*, *c* and *d* are constants. Integrating this equation from *x*<sup>0</sup> to *x*1, we have

*x*1

This boundary condition has no truncation error since the integrations are carried out exactly. For more general cases when some of the terms in Eq. (1) cannot be integrated exactly, various numerical integration schemes, such as the second-order Trapezoidal scheme, can be used to

Stability is critical to numerical methods used to solve time-dependent systems. In this section, we conduct Von Neumann stability analysis for the new method combined with the standard second-order central difference scheme. The following one-dimensional linear

The new method combined with the second-order central difference for solving Eq. (21) is

<sup>Δ</sup>*x*<sup>2</sup> *<sup>δ</sup>*2*un*+<sup>1</sup>

*<sup>i</sup>* +

 *u*ˆ*n v*ˆ*n* 

<sup>2</sup> ) −*β*

<sup>2</sup> ) <sup>1</sup> <sup>+</sup> <sup>2</sup>*α*sin2( *<sup>θ</sup>*

2 ) ,

<sup>Δ</sup>*x*<sup>2</sup> *<sup>δ</sup>*2*un*+<sup>1</sup>

= *M*−1*R*

*<sup>i</sup>* <sup>+</sup> *cvn*+<sup>1</sup>

*δ*2 <sup>Δ</sup>*x*<sup>2</sup> *<sup>v</sup>n*+<sup>1</sup>

*aux*| *x*1 *<sup>x</sup>*<sup>0</sup> + *cu*|

*<sup>v</sup>*<sup>0</sup> <sup>=</sup> <sup>1</sup> *a*

Solving *v*<sup>0</sup> from Eq. (19), we obtain the boundary condition

generate boundary condition for *v* in a similar way.

convection-diffusion equation is used in the analysis

*un*+<sup>1</sup> *<sup>i</sup>* <sup>−</sup> *<sup>u</sup><sup>n</sup> i* <sup>Δ</sup>*<sup>t</sup>* <sup>=</sup> <sup>1</sup>

*vn*+<sup>1</sup> *<sup>i</sup>* <sup>−</sup> *<sup>v</sup><sup>n</sup> i* <sup>Δ</sup>*<sup>t</sup>* <sup>=</sup> *<sup>c</sup>*

*M* = 

Taking the discrete Fourier transform of the above equations, we have

 *u*ˆ*n*+<sup>1</sup> *v*ˆ*n*+<sup>1</sup> 

1 + 2*α*sin2( *<sup>θ</sup>*

2*γ*sin2( *<sup>θ</sup>*

*M u*ˆ*n*+<sup>1</sup> *v*ˆ*n*+<sup>1</sup> = *R u*ˆ*n v*ˆ*n* 

**3. Stability analysis**

*auxx* + *cux* + *d* = 0, (17)

*av*<sup>1</sup> − *av*<sup>0</sup> + *cu*<sup>1</sup> − *cu*<sup>0</sup> + *d*Δ*x* = 0. (19)

(*av*<sup>1</sup> + *cu*<sup>1</sup> − *cu*<sup>0</sup> + *d*Δ*x*). (20)

*ut* = *uxx* + *cux* (21)

*<sup>i</sup>* , (22)

*<sup>i</sup>* . (23)

. (24)

, (25)

*<sup>x</sup>*<sup>0</sup> + *d*(*x*<sup>1</sup> − *x*0) = 0, (18)

coefficient and constant *f*(*u*), we have

for Solving Unsteady Convection-Diffusion Equations

or

Thus

where

Applying the operator (1 + *<sup>δ</sup>*<sup>2</sup> <sup>12</sup> ) to both sides of these two equations, we obtain

$$(1 + \frac{\delta^2}{12})\frac{u\_i^{n+1} - u\_i^n}{\Delta t} = a \frac{\delta^2 u\_i^{n+1}}{\Delta x^2} + (1 + \frac{\delta^2}{12})[c\_i v\_i + f(u\_i^{n+1})],\tag{10}$$

$$\begin{split} (1+\frac{\delta^2}{12})\frac{\upsilon\_i^{n+1}-\upsilon\_i^n}{\Delta t} &= c\_i \frac{\delta^2 u\_i^{n+1}}{\Delta x^2} + a \frac{\delta^2 \upsilon\_i^{n+1}}{\Delta x^2} \\ &+ (1+\frac{\delta^2}{12})[(c\_x)\_i + f\_u(u\_i^{n+1})] \upsilon\_i^{n+1}. \end{split} \tag{11}$$

Eqs. (10) and (11) are fourth-order accurate in space but only contains the second-order central difference operator *δ*<sup>2</sup> that is defined on a 3-point stencil. As a result, the discretized equations form a system of block tri-diagonal algebraic equations.

#### **2.2 Initial and boundary conditions**

The initial condition for *v* can be obtained by differentiating *h*(*x*), the initial condition fro *u* given in Eq. (1) with respect to *x*, i.e.

$$v(x,0) = \frac{dh(x)}{dx}.$$

The boundary conditions for *v*, or for Eq. (7), are less straight forward. In the following, we discuss three different ways to generate boundary conditions for *v* at the spatial grid point *i* = 0, assuming Dirichlet boundary conditions are given for *u*. The boundary conditions for *v* at the spatial grid point *i* = *M* and other scenarios can be dealt with in similar ways. Standard finite difference approximation: Since *v* = *ux*, we have

$$v\_0 = \frac{u\_1 - u\_0}{\Delta \mathbf{x}}\text{,}\tag{12}$$

which provides an equation to complement the equation obtained by discretizing Eq. (7) at *i* = 1. This approximation is first-order accurate in space. If necessary, higher-order one-sided finite difference schemes can be used to approximate *v* = *ux*.

Padé approximation: We can use the fourth-order Padé approximation at *i* = 1 for *v* = *ux*

$$w\_1 = \frac{\delta u\_1}{2\Delta x (1 + \frac{\delta^2}{6})} = \frac{u\_2 - u\_0}{2\Delta x (1 + \frac{\delta^2}{6})}.\tag{13}$$

Applying (1 + *<sup>δ</sup>*<sup>2</sup> <sup>6</sup> ) to both sides of Eq. (13), we have

$$(1 + \frac{\delta^2}{6})v\_1 = \frac{\mu\_2 - \mu\_0}{2\Delta x},\tag{14}$$

or

$$
\frac{1}{6}v\_0 + \frac{2}{3}v\_1 + \frac{1}{6}v\_2 = \frac{\mu\_2 - \mu\_0}{2\Delta x}.\tag{15}
$$

Solving *v*<sup>0</sup> from Eq. (15), we obtain the boundary condition

$$v\_0 = \frac{\Im(u\_2 - u\_0)}{\Delta x} - 4v\_1 - v\_2. \tag{16}$$

Direct integration of Eq. (1): This approach is particularly convenient for special cases of Eq. (1). For example, when dealing with steady state equation with a constant convection coefficient and constant *f*(*u*), we have

$$
tau\_{\rm xx} + c u\_{\rm x} + d = 0,\tag{17}$$

where *a*, *c* and *d* are constants. Integrating this equation from *x*<sup>0</sup> to *x*1, we have

$$
tau\_{\mathbf{x}}|\_{\mathbf{x}\_0}^{\mathbf{x}\_1} + c\boldsymbol{u}|\_{\mathbf{x}\_0}^{\mathbf{x}\_1} + d(\mathbf{x}\_1 - \mathbf{x}\_0) = \mathbf{0},\tag{18}$$

or

4 Will-be-set-by-IN-TECH

*δ*2*un*+<sup>1</sup> *i* <sup>Δ</sup>*x*<sup>2</sup> + (<sup>1</sup> <sup>+</sup>

> *δ*2*un*+<sup>1</sup> *i* <sup>Δ</sup>*x*<sup>2</sup> <sup>+</sup> *<sup>a</sup>*

Eqs. (10) and (11) are fourth-order accurate in space but only contains the second-order central difference operator *δ*<sup>2</sup> that is defined on a 3-point stencil. As a result, the discretized equations

The initial condition for *v* can be obtained by differentiating *h*(*x*), the initial condition fro *u*

*<sup>v</sup>*(*x*, 0) = *dh*(*x*)

The boundary conditions for *v*, or for Eq. (7), are less straight forward. In the following, we discuss three different ways to generate boundary conditions for *v* at the spatial grid point *i* = 0, assuming Dirichlet boundary conditions are given for *u*. The boundary conditions for *v*

*<sup>v</sup>*<sup>0</sup> <sup>=</sup> *<sup>u</sup>*<sup>1</sup> <sup>−</sup> *<sup>u</sup>*<sup>0</sup>

which provides an equation to complement the equation obtained by discretizing Eq. (7) at *i* = 1. This approximation is first-order accurate in space. If necessary, higher-order one-sided

Padé approximation: We can use the fourth-order Padé approximation at *i* = 1 for *v* = *ux*

6 )

<sup>6</sup> )*v*<sup>1</sup> <sup>=</sup> *<sup>u</sup>*<sup>2</sup> <sup>−</sup> *<sup>u</sup>*<sup>0</sup>

<sup>=</sup> *<sup>u</sup>*<sup>2</sup> <sup>−</sup> *<sup>u</sup>*<sup>0</sup> 2Δ*x*(1 + *<sup>δ</sup>*<sup>2</sup>

*<sup>v</sup>*<sup>2</sup> <sup>=</sup> *<sup>u</sup>*<sup>2</sup> <sup>−</sup> *<sup>u</sup>*<sup>0</sup>

6 )

at the spatial grid point *i* = *M* and other scenarios can be dealt with in similar ways.

Standard finite difference approximation: Since *v* = *ux*, we have

finite difference schemes can be used to approximate *v* = *ux*.

*<sup>v</sup>*<sup>1</sup> <sup>=</sup> *<sup>δ</sup>u*<sup>1</sup>

(1 + *δ*2

<sup>6</sup> ) to both sides of Eq. (13), we have

1 6 *v*<sup>0</sup> + 2 3 *v*<sup>1</sup> + 1 6

Solving *v*<sup>0</sup> from Eq. (15), we obtain the boundary condition

2Δ*x*(1 + *<sup>δ</sup>*<sup>2</sup>

*<sup>v</sup>*<sup>0</sup> <sup>=</sup> <sup>3</sup>(*u*<sup>2</sup> <sup>−</sup> *<sup>u</sup>*0)

*dx* .

*δ*2

+(1 +

<sup>12</sup> ) to both sides of these two equations, we obtain

*δ*2

*δ*2*vn*+<sup>1</sup> *i*

<sup>12</sup> )[(*cx*)*<sup>i</sup>* <sup>+</sup> *fu*(*un*+<sup>1</sup>

<sup>12</sup> )[*civi* <sup>+</sup> *<sup>f</sup>*(*un*+<sup>1</sup>

*<sup>i</sup>* )], (10)

<sup>Δ</sup>*x*<sup>2</sup> (11)

*<sup>i</sup>* )]*vn*+<sup>1</sup> *<sup>i</sup>* .

<sup>Δ</sup>*<sup>x</sup>* , (12)

<sup>2</sup>Δ*<sup>x</sup>* , (14)

<sup>2</sup>Δ*<sup>x</sup>* . (15)

<sup>Δ</sup>*<sup>x</sup>* <sup>−</sup> <sup>4</sup>*v*<sup>1</sup> <sup>−</sup> *<sup>v</sup>*2. (16)

. (13)

Applying the operator (1 + *<sup>δ</sup>*<sup>2</sup>

(1 + *δ*2 12 )

> (1 + *δ*2 12 )

**2.2 Initial and boundary conditions**

given in Eq. (1) with respect to *x*, i.e.

Applying (1 + *<sup>δ</sup>*<sup>2</sup>

or

*un*+<sup>1</sup> *<sup>i</sup>* <sup>−</sup> *<sup>u</sup><sup>n</sup> i* <sup>Δ</sup>*<sup>t</sup>* <sup>=</sup> *<sup>a</sup>*

*vn*+<sup>1</sup> *<sup>i</sup>* <sup>−</sup> *<sup>v</sup><sup>n</sup> i* <sup>Δ</sup>*<sup>t</sup>* <sup>=</sup> *ci*

form a system of block tri-diagonal algebraic equations.

$$av\_1 - av\_0 + cu\_1 - cu\_0 + d\Delta \mathbf{x} = 0.\tag{19}$$

Solving *v*<sup>0</sup> from Eq. (19), we obtain the boundary condition

$$v\_0 = \frac{1}{a}(av\_1 + cu\_1 - cu\_0 + d\Delta x). \tag{20}$$

This boundary condition has no truncation error since the integrations are carried out exactly. For more general cases when some of the terms in Eq. (1) cannot be integrated exactly, various numerical integration schemes, such as the second-order Trapezoidal scheme, can be used to generate boundary condition for *v* in a similar way.

#### **3. Stability analysis**

Stability is critical to numerical methods used to solve time-dependent systems. In this section, we conduct Von Neumann stability analysis for the new method combined with the standard second-order central difference scheme. The following one-dimensional linear convection-diffusion equation is used in the analysis

$$u\_t = u\_{xx} + cu\_x \tag{21}$$

The new method combined with the second-order central difference for solving Eq. (21) is

$$\frac{u\_i^{n+1} - u\_i^n}{\Delta t} = \frac{1}{\Delta x^2} \delta^2 u\_i^{n+1} + c v\_i^{n+1},\tag{22}$$

$$\frac{v\_i^{n+1} - v\_i^n}{\Delta t} = \frac{c}{\Delta x^2} \delta^2 u\_i^{n+1} + \frac{\delta^2}{\Delta x^2} v\_i^{n+1}.\tag{23}$$

Taking the discrete Fourier transform of the above equations, we have

$$\mathcal{M}\begin{bmatrix}\mathfrak{A}^{n+1}\\\mathfrak{d}^{n+1}\end{bmatrix} = \mathcal{R}\begin{bmatrix}\mathfrak{A}^{n}\\\mathfrak{d}^{n}\end{bmatrix}.\tag{24}$$

Thus

$$
\begin{bmatrix} \hat{\mathfrak{U}}^{n+1} \\ \hat{\mathfrak{O}}^{n+1} \end{bmatrix} = \boldsymbol{M}^{-1} \boldsymbol{R} \begin{bmatrix} \hat{\mathfrak{U}}^{n} \\ \hat{\mathfrak{O}}^{n} \end{bmatrix} \tag{25}
$$

where

$$M = \begin{bmatrix} 1 + 2\alpha \sin^2(\frac{\theta}{2}) & -\beta \\ 2\gamma \sin^2(\frac{\theta}{2}) & 1 + 2\alpha \sin^2(\frac{\theta}{2}) \end{bmatrix} / 2$$

Upwind Stand-cfd 4*th*-order new Error Time Error Time Error Time 9.7E-03 .0016 4.1E-04 .0011 7.4E-07 .0010 9.9E-04 .0098 9.8E-06 .0012 9.2E-09 .0012 1.0E-04 .5155 1.0E-06 .0016 1.8E-10 .0014 1.0E-05 68.05 1.0E-08 .0101 4.7E-12 .0016 1.0E-08 ∞ 1.0E-10 .4870 2.9E-13 .0022 1.0E-09 ∞ 1.5E-11 7.059 4.1E-14 .0035

<sup>87</sup> A Fourth-Order Compact Finite Difference Scheme

first-order upwind scheme, it is almost impossible to obtain an error less than 1.0*E* − 08 since

Note that the symbol ∞ in Table 1 does not really represent infinity. It just represents an excessively large number. The computation time for each test case is the average of five runs.

Several numerical examples are presented here to compare the new method with the standard

**Example 1**: In this example, we compare the accuracy of four different algorithms using the same 3-point finite difference stencil. The following steady state convection-diffusion

> *c <sup>a</sup> <sup>x</sup>* <sup>−</sup> <sup>1</sup> *e c <sup>a</sup>* − 1

Δ*x* Upwind 2*nd*-stand 2*nd*-new 4*th*-new 1/10 1.32E-01 3.45E-02 1.41E-02 7.41E-04 1/20 7.64E-02 7.90E-03 3.70E-03 4.74E-05 1/40 4.17E-02 1.90E-03 9.52E-04 2.98E-06 1/80 2.18E-02 4.79E-04 2.39E-04 1.87E-07 1/160 1.12E-02 1.20E-04 5.98E-05 1.17E-08 Table 2. Errors between the exact solution and the numerical solution of Eq.(27) with *a* = 0.1

Table 2. shows the error �*u<sup>e</sup>* <sup>−</sup> *<sup>u</sup>c*�<sup>∞</sup> between the exact solution *<sup>u</sup><sup>e</sup>* and the numerical solution *u<sup>c</sup>* calculated using the following four finite difference schemes with *a* = 0.1 and *c* = 1:

1. Upwind: First-order upwind difference for the convection term and second-order central

2. 2*nd*-stand: Second-order central difference for both the diffusion and convection terms. 3. 2*nd*-new: New method that uses second-order central difference for Eqs. (5) and (7).

*<sup>u</sup>*(*x*) = *<sup>e</sup>*

− *auxx* + *cux* = 0, *u*(0) = 0, *u*(1) = 1, (27)

. (28)

Table 1. Computational cost vs accuracy for different algorithms

where *a* and *c* are constants. The exact solution is given as

finite difference for the diffusion term.

a grid size of *h* = 1.0*E* − 08 is needed.

for Solving Unsteady Convection-Diffusion Equations

The unit for computation time is second.

**5. Numerical results**

and *c* = 1.

finite difference algorithms.

equation is used in the example:

and

$$\mathcal{R} = \begin{bmatrix} 1 \ 0 \\ 0 \ 1 \end{bmatrix}'$$

with *α* = <sup>Δ</sup>*<sup>t</sup>* <sup>Δ</sup>*x*<sup>2</sup> , *<sup>β</sup>* <sup>=</sup> *<sup>c</sup>*Δ*t*, and *<sup>γ</sup>* <sup>=</sup> *<sup>c</sup>*Δ*<sup>t</sup>* <sup>Δ</sup>*x*<sup>2</sup> .

Here *M*−1*R* is the amplification matrix at each time-step. In order for the numerical algorithm to be stable, the modulus of the eigenvalues of *M*−1*R* must be less than or equal to unity for all possible values of *θ*. For the 2 × 2 matrix, it is easy to see that it has two conjugate complex eigenvalues. If the modulus of the two eigenvalues is less than or equal to 1, then the method is unconditionally stable.

The eigenvalues of *M*−1*R* can be calculated as

$$
\omega = \frac{1 + 2a\sin^2(\frac{\theta}{2}) \pm i(c\sqrt{2a\Delta t}\sin(\frac{\theta}{2}))}{1 + 4a\sin^2(\frac{\theta}{2}) + 4a^2\sin^4(\frac{\theta}{2}) + 2c^2\Delta t a \sin^2(\frac{\theta}{2})},
$$

thus we have

$$\left|\omega\right|^2 = \frac{1}{1 + 4a\sin^2(\frac{\theta}{2}) + 4a^2\sin^4(\frac{\theta}{2}) + 2c^2\Delta t a \sin^2(\frac{\theta}{2})}$$

Obviously, |*ω*| <sup>2</sup> <sup>≤</sup> 1 since 0 <sup>≤</sup> sin2( *<sup>θ</sup>* <sup>2</sup> ) ≤ 1 and *α* > 0.

#### **4. Computational complexity**

In this section an estimation of computational cost of the new method is obtained, and computation times for three methods (upwind, standard central finite difference or standard cfd, and fourth-order new method) are compared using the following model equation

$$
\mu\_{\rm x} - \mu\_{\rm xx} = 0, \quad \mathfrak{u}(0) = 0, \quad \mathfrak{u}(1) = 1,\tag{26}
$$

with exact solution *<sup>u</sup>*(*x*) = *<sup>e</sup><sup>x</sup>*−<sup>1</sup> *<sup>e</sup>*−<sup>1</sup> .

If the above model equation is discretized over *n* subintervals with equal size, and the numerical solution at the *i*-th grid point is denoted as *ui*, the solutions can be obtained by solving the linear system *AU* = *b*. For both the upwind and standard central finite difference schemes, *A* is an *n* × *n* tridiagonal matrix and *U* = (*u*1, *u*2, ··· , *un*). For the fourth-order new method the matrix *A* is a 2*n* × 2*n* banded matrix with a bandwidth of 7, and *U* = (*u*1, *v*1, *u*2, *v*2, ··· , *un*, *vn*) contains the numerical solutions for both *u* and *v*, with *v* = *ux*.

If an efficient algorithm is used to solve the tridiagonal linear system resulted from upwind and standard central finite difference schemes,a total of 2*n* − 1 divisions, 3(*n* − 1) multiplications, and 3(*n* − 1) subtractions are needed. The computational complexity is 2(*n* − 1) + 1 + 3(*n* − 1) + 3(*n* − 1) = 8*n* − 7 = 8*n* + **O**(1) flops.

If the Gaussian elimination algorithm is used to solve the linear system resulted from the fourth-order new method, a total of 18*n* − 15 divisions, 24*n* − 22 multiplications, and 24*n* − 22 subtractions are needed. The computational complexity is 66*n* − 59 = 66*n* + **O**(1) flops.

Table 1 clearly shows that to achieve the same accuracy the fourth-order new method is much faster than both upwind and the standard second-order central finite difference schemes. For instance, to obtain an error less than 1.0*E*-05, the upwind scheme needs 68.054 seconds, and the standard second-order central finite difference scheme needs 0.0012 seconds, while the fourth-order new method needs only 0.0010 seconds to obtain an error of 7.4*E*-07. For the


Table 1. Computational cost vs accuracy for different algorithms

first-order upwind scheme, it is almost impossible to obtain an error less than 1.0*E* − 08 since a grid size of *h* = 1.0*E* − 08 is needed.

Note that the symbol ∞ in Table 1 does not really represent infinity. It just represents an excessively large number. The computation time for each test case is the average of five runs. The unit for computation time is second.

#### **5. Numerical results**

6 Will-be-set-by-IN-TECH

Here *M*−1*R* is the amplification matrix at each time-step. In order for the numerical algorithm to be stable, the modulus of the eigenvalues of *M*−1*R* must be less than or equal to unity for all possible values of *θ*. For the 2 × 2 matrix, it is easy to see that it has two conjugate complex eigenvalues. If the modulus of the two eigenvalues is less than or equal to 1, then the method

<sup>2</sup> ) ± *i*(*c*

<sup>2</sup> ) + <sup>4</sup>*α*2sin4( *<sup>θ</sup>*

<sup>2</sup> ) + <sup>4</sup>*α*2sin4( *<sup>θ</sup>*

In this section an estimation of computational cost of the new method is obtained, and computation times for three methods (upwind, standard central finite difference or standard cfd, and fourth-order new method) are compared using the following model equation

If the above model equation is discretized over *n* subintervals with equal size, and the numerical solution at the *i*-th grid point is denoted as *ui*, the solutions can be obtained by solving the linear system *AU* = *b*. For both the upwind and standard central finite difference schemes, *A* is an *n* × *n* tridiagonal matrix and *U* = (*u*1, *u*2, ··· , *un*). For the fourth-order new method the matrix *A* is a 2*n* × 2*n* banded matrix with a bandwidth of 7, and *U* = (*u*1, *v*1, *u*2, *v*2, ··· , *un*, *vn*) contains the numerical solutions for both *u* and *v*, with

If an efficient algorithm is used to solve the tridiagonal linear system resulted from upwind and standard central finite difference schemes,a total of 2*n* − 1 divisions, 3(*n* − 1) multiplications, and 3(*n* − 1) subtractions are needed. The computational complexity is

If the Gaussian elimination algorithm is used to solve the linear system resulted from the fourth-order new method, a total of 18*n* − 15 divisions, 24*n* − 22 multiplications, and 24*n* − 22 subtractions are needed. The computational complexity is 66*n* − 59 = 66*n* + **O**(1) flops. Table 1 clearly shows that to achieve the same accuracy the fourth-order new method is much faster than both upwind and the standard second-order central finite difference schemes. For instance, to obtain an error less than 1.0*E*-05, the upwind scheme needs 68.054 seconds, and the standard second-order central finite difference scheme needs 0.0012 seconds, while the fourth-order new method needs only 0.0010 seconds to obtain an error of 7.4*E*-07. For the

<sup>2</sup> ) ≤ 1 and *α* > 0.

<sup>√</sup>2*α*Δ*<sup>t</sup>* sin( *<sup>θ</sup>*

<sup>2</sup> ))

2 ) ,

2 ).

<sup>2</sup> ) + <sup>2</sup>*c*2Δ*tα*sin2( *<sup>θ</sup>*

<sup>2</sup> ) + <sup>2</sup>*c*2Δ*tα*sin2( *<sup>θ</sup>*

*ux* − *uxx* = 0, *u*(0) = 0, *u*(1) = 1, (26)

*R* = 1 0 0 1 ,

<sup>Δ</sup>*x*<sup>2</sup> .

*<sup>ω</sup>* <sup>=</sup> <sup>1</sup> <sup>+</sup> <sup>2</sup>*α*sin2( *<sup>θ</sup>*

1 + 4*α*sin2( *<sup>θ</sup>*

2(*n* − 1) + 1 + 3(*n* − 1) + 3(*n* − 1) = 8*n* − 7 = 8*n* + **O**(1) flops.

<sup>2</sup> <sup>=</sup> <sup>1</sup> 1 + 4*α*sin2( *<sup>θ</sup>*

and

with *α* = <sup>Δ</sup>*<sup>t</sup>*

thus we have

Obviously, |*ω*|

*v* = *ux*.

is unconditionally stable.

<sup>Δ</sup>*x*<sup>2</sup> , *<sup>β</sup>* <sup>=</sup> *<sup>c</sup>*Δ*t*, and *<sup>γ</sup>* <sup>=</sup> *<sup>c</sup>*Δ*<sup>t</sup>*

The eigenvalues of *M*−1*R* can be calculated as


**4. Computational complexity**

with exact solution *<sup>u</sup>*(*x*) = *<sup>e</sup><sup>x</sup>*−<sup>1</sup> *<sup>e</sup>*−<sup>1</sup> .

<sup>2</sup> <sup>≤</sup> 1 since 0 <sup>≤</sup> sin2( *<sup>θ</sup>*

Several numerical examples are presented here to compare the new method with the standard finite difference algorithms.

**Example 1**: In this example, we compare the accuracy of four different algorithms using the same 3-point finite difference stencil. The following steady state convection-diffusion equation is used in the example:

$$-\mathfrak{a}\mathfrak{u}\_{\mathfrak{X}} + \mathfrak{c}\mathfrak{u}\_{\mathfrak{X}} = \mathfrak{0}, \quad \mathfrak{u}(\mathfrak{0}) = \mathfrak{0}, \mathfrak{u}(\mathfrak{1}) = \mathfrak{1}, \tag{27}$$

where *a* and *c* are constants. The exact solution is given as

$$u(x) = \frac{e^{\frac{\varepsilon}{4}x} - 1}{e^{\frac{\varepsilon}{4}} - 1}. \tag{28}$$


Table 2. Errors between the exact solution and the numerical solution of Eq.(27) with *a* = 0.1 and *c* = 1.

Table 2. shows the error �*u<sup>e</sup>* <sup>−</sup> *<sup>u</sup>c*�<sup>∞</sup> between the exact solution *<sup>u</sup><sup>e</sup>* and the numerical solution *u<sup>c</sup>* calculated using the following four finite difference schemes with *a* = 0.1 and *c* = 1:


Fig. 2. Solution curves corresponding to *a* = 0.1, 0.05, 0.001 and *c* = 1 calculated using the

<sup>89</sup> A Fourth-Order Compact Finite Difference Scheme

Fig. 3. Solution curves corresponding to *a* = 1.0, 0.1, 0.01 and 0.001 calculated by the new

central difference scheme and the convection term is discretized by the first-order upwind scheme assuming the convection coefficient is positive. The change of sign of the convection coefficient is not taken into consideration. The solution curves correspond to *a* = 1.0 and 0.1 look reasonable. But the solution curve correspond to *a* = 0.01 is not accurate. It does not cross the solution curves corresponding to *a* = 1.0 and 0.1 at *x* = 0.5. The solution curve corresponds to *a* = 0.001 does not resemble the correct solution at all. This is not surprising since the convection coefficient changes sign but the upwind scheme used in the calculation

<sup>2</sup>*nd*−new algorithm with <sup>Δ</sup>*<sup>x</sup>* <sup>=</sup> 0.01.

for Solving Unsteady Convection-Diffusion Equations

method using central difference scheme

does not take this into consideration.

4. 4*th*-new: New method that uses fourth-order Padé approximation for Eqs. (5) and (7).

All results are obtained using a 3-point stencil. It is clear from Table 2 that the accuracies of different schemes are as expected: when the grid size Δ*x* is reduced by <sup>1</sup> <sup>2</sup> , the errors are reduced by approximately <sup>1</sup> <sup>2</sup> , ( <sup>1</sup> <sup>2</sup> )<sup>2</sup> and ( <sup>1</sup> <sup>2</sup> )<sup>4</sup> for the 1*st*-, 2*nd*- and 4*th*-order algorithms, respectively.

Fig. 1. Solution curves corresponding to *a* = 0.1, 0.05, 0.001 and *c* = 1 using the 2*nd*-stand algorithm with Δ*x* = 0.001.

Although both the 2*nd*-stand and 2*nd*-new algorithms produced comparable results for the case of *a* = 0.1 and *c* = 1 shown in Table 2, the 2*nd*-new algorithm is much more robust when the boundary layer near *x* = 1 becomes steeper as *a* decreases. Fig. 1 and Fig. 2 show the solutions calculated using the 2*nd*-stand and 2*nd*-new central difference schemes, respectively. The three curves in each figure correspond to the diffusion coefficient *a* = 0.1, 0.05, and 0.001, respectively. It is clear from Fig. 1 that the solution calculated by the standard second-order central difference scheme becomes highly oscillatory when *a* reaches 0.001, while the solution calculated by the new second-order central difference scheme shown in Fig. 2 describes the steep boundary layer near *x* = 1 very well.

**Example 2**: In this example, we compare the robustness of three different algorithms. The governing equation for this example is

$$(-\
u u\_{xx} + (x - \frac{1}{2})u\_x = 0, \quad u(0) = 0, \quad u(1) = 1,\tag{29}$$

where *a* is a constant. The convection coefficient changes sign in the middle of the domain, which makes this turning point problem difficult (Morton, 1996).

Fig. 3 shows the solution curves obtained by the new method using the second-order central difference algorithm. The four solution curves in the figure correspond to diffusion coefficient *a* = 1.0, 0.1, 0.01, 0.001, respectively, all calculated with Δ*x* = 0.001. It is clear that as *a* decreases, the boundary layers at *x* = 0 and *x* = 1 become steeper. Also note that all four solution curves pass the same point *x* = 0.5, where the convection coefficient is zero.

Fig. 4 shows the solution curves obtained by the standard finite difference schemes on the same grid with Δ*x* = 0.001. The diffusion term is discretized by the second-order 8 Will-be-set-by-IN-TECH

<sup>2</sup> , the errors

<sup>2</sup> )<sup>4</sup> for the 1*st*-, 2*nd*- and 4*th*-order algorithms,

4. 4*th*-new: New method that uses fourth-order Padé approximation for Eqs. (5) and (7). All results are obtained using a 3-point stencil. It is clear from Table 2 that the accuracies

Fig. 1. Solution curves corresponding to *a* = 0.1, 0.05, 0.001 and *c* = 1 using the 2*nd*-stand

Although both the 2*nd*-stand and 2*nd*-new algorithms produced comparable results for the case of *a* = 0.1 and *c* = 1 shown in Table 2, the 2*nd*-new algorithm is much more robust when the boundary layer near *x* = 1 becomes steeper as *a* decreases. Fig. 1 and Fig. 2 show the solutions calculated using the 2*nd*-stand and 2*nd*-new central difference schemes, respectively. The three curves in each figure correspond to the diffusion coefficient *a* = 0.1, 0.05, and 0.001, respectively. It is clear from Fig. 1 that the solution calculated by the standard second-order central difference scheme becomes highly oscillatory when *a* reaches 0.001, while the solution calculated by the new second-order central difference scheme shown in Fig. 2 describes the

**Example 2**: In this example, we compare the robustness of three different algorithms. The

where *a* is a constant. The convection coefficient changes sign in the middle of the domain,

Fig. 3 shows the solution curves obtained by the new method using the second-order central difference algorithm. The four solution curves in the figure correspond to diffusion coefficient *a* = 1.0, 0.1, 0.01, 0.001, respectively, all calculated with Δ*x* = 0.001. It is clear that as *a* decreases, the boundary layers at *x* = 0 and *x* = 1 become steeper. Also note that all four

Fig. 4 shows the solution curves obtained by the standard finite difference schemes on the same grid with Δ*x* = 0.001. The diffusion term is discretized by the second-order

solution curves pass the same point *x* = 0.5, where the convection coefficient is zero.

)*ux* = 0, *u*(0) = 0, *u*(1) = 1, (29)

of different schemes are as expected: when the grid size Δ*x* is reduced by <sup>1</sup>

<sup>2</sup> )<sup>2</sup> and ( <sup>1</sup>

<sup>2</sup> , ( <sup>1</sup>

are reduced by approximately <sup>1</sup>

algorithm with Δ*x* = 0.001.

steep boundary layer near *x* = 1 very well.

<sup>−</sup> *auxx* + (*<sup>x</sup>* <sup>−</sup> <sup>1</sup>

which makes this turning point problem difficult (Morton, 1996).

2

governing equation for this example is

respectively.

Fig. 2. Solution curves corresponding to *a* = 0.1, 0.05, 0.001 and *c* = 1 calculated using the <sup>2</sup>*nd*−new algorithm with <sup>Δ</sup>*<sup>x</sup>* <sup>=</sup> 0.01.

Fig. 3. Solution curves corresponding to *a* = 1.0, 0.1, 0.01 and 0.001 calculated by the new method using central difference scheme

central difference scheme and the convection term is discretized by the first-order upwind scheme assuming the convection coefficient is positive. The change of sign of the convection coefficient is not taken into consideration. The solution curves correspond to *a* = 1.0 and 0.1 look reasonable. But the solution curve correspond to *a* = 0.01 is not accurate. It does not cross the solution curves corresponding to *a* = 1.0 and 0.1 at *x* = 0.5. The solution curve corresponds to *a* = 0.001 does not resemble the correct solution at all. This is not surprising since the convection coefficient changes sign but the upwind scheme used in the calculation does not take this into consideration.

**Example 3**: In this example we solve a nonlinear Black-Scholes equation, which is widely used to model option price when transaction cost is considered. It is well known that in the area of financial mathematics, the calculations of both the solution and its first derivative are required. The derivative is the so-called hedge delta, which actually represents the number of shares (of stock) that one should hold or sale, in order to maximize the profit. The nonlinear Black-Scholes equation solved in this chapter is given below. More details about this model

<sup>91</sup> A Fourth-Order Compact Finite Difference Scheme

*<sup>u</sup>*(*x*, 0) = *max*(<sup>1</sup> <sup>−</sup> *<sup>e</sup>*−*x*, 0), *<sup>u</sup>*(−∞,*t*) = 0, *<sup>u</sup>*(∞,*t*) = 1,

is the volatility of the stock, and Φ is a function defined as the solution to the following ODE:

*<sup>s</sup>*Φ(*s*) <sup>−</sup> *<sup>s</sup>*

(*s*) = <sup>Φ</sup>(*s*) + <sup>1</sup> 2

Fig. 6. Solution curves for option price corresponding to various transaction cost

In the numerical solution process, the infinite domain is approximated by a finite interval of sufficient length. In financial industry, one of the widely used methods for calculating hedge delta is to first solve the Black-Scholes equation to obtain numerical solution of *u*, and then apply finite difference schemes to the numerical solution to obtain approximations to *ux*. With the new method discussed in this chapter, both the solution *u* and the hedge delta *ux* are calculated simultaneously. Fig. 6 and Fig. 7 show the option price and the hedge delta for

**Example 4**: In this example we solve a time dependent nonlinear equation, the Burgers'

*ut* + *uux* = *auxx*, *u*(0, *t*) = 0, *u*(1, *t*) = 1, (31)

rates(*a* = 0.0, 0.01, 0.02) and the pay-off curve.

various transaction cost rates, respectively.

Equation:

*<sup>a</sup>*2*E*(*uxx* <sup>+</sup> *ux*)])(*uxx* <sup>+</sup> *ux*) <sup>−</sup> *Kux*,

*σ*2 0 , *ρ* is the risk-free interest rate, *σ*<sup>0</sup>

, Φ(0) = 0, (30)

can be found in (Barles & Soner, 1998).

with initial and boundary conditions

*ut* = (1 + Φ[*e*

for Solving Unsteady Convection-Diffusion Equations

where *a* is transaction cost rate, *E* is the strike price, *K* = <sup>2</sup>*<sup>ρ</sup>*

Φ�

(*Kt*+*x*)

Fig. 4. Solution curves corresponding to *a* = 1.0, 0.1, 0.01 and 0.001 calculated by the upwind scheme without considering the change of sign of the convection coefficient.

Fig. 5. Solution curves corresponding to *a* = 1.0, 0.1, 0.01 and 0.001 calculated by the upwind scheme taking into consideration the sign change of the convection coefficient.

Fig. 5 shows the solution curves obtained by the standard finite difference algorithm. The diffusion term is discretized by the second-order central difference scheme and the convection term is discretized by the first-order upwind scheme taking into consideration the sign change of the convection coefficient. We can see the solution curves correspond to *a* = 1.0, 0.1, and 0.01 are very similar to those in Fig. 3. All three solution curves cross each other at the same point *x* = 0.5. However, the solution curve corresponding to *a* = 0.001 does not resemble the true solution at all. This is somewhat surprising but probably can be attributed to the difficulty caused by the turning point at *x* = 0.5 where *c*(*x*) = 0. The situation does not improve even when the grid spacing Δ*x* is further reduced.

10 Will-be-set-by-IN-TECH

Fig. 4. Solution curves corresponding to *a* = 1.0, 0.1, 0.01 and 0.001 calculated by the upwind

Fig. 5. Solution curves corresponding to *a* = 1.0, 0.1, 0.01 and 0.001 calculated by the upwind

Fig. 5 shows the solution curves obtained by the standard finite difference algorithm. The diffusion term is discretized by the second-order central difference scheme and the convection term is discretized by the first-order upwind scheme taking into consideration the sign change of the convection coefficient. We can see the solution curves correspond to *a* = 1.0, 0.1, and 0.01 are very similar to those in Fig. 3. All three solution curves cross each other at the same point *x* = 0.5. However, the solution curve corresponding to *a* = 0.001 does not resemble the true solution at all. This is somewhat surprising but probably can be attributed to the difficulty caused by the turning point at *x* = 0.5 where *c*(*x*) = 0. The situation does not improve even

scheme taking into consideration the sign change of the convection coefficient.

when the grid spacing Δ*x* is further reduced.

scheme without considering the change of sign of the convection coefficient.

**Example 3**: In this example we solve a nonlinear Black-Scholes equation, which is widely used to model option price when transaction cost is considered. It is well known that in the area of financial mathematics, the calculations of both the solution and its first derivative are required. The derivative is the so-called hedge delta, which actually represents the number of shares (of stock) that one should hold or sale, in order to maximize the profit. The nonlinear Black-Scholes equation solved in this chapter is given below. More details about this model can be found in (Barles & Soner, 1998).

$$\mu\_t = (1 + \Phi[e^{(Kt+x)}a^2E(\mu\_{xx} + \mu\_x)])(\mu\_{xx} + \mu\_x) - K\mu\_{xx}$$

with initial and boundary conditions

$$u(\mathbf{x},0) = \max(1 - e^{-\mathbf{x}}, 0), \quad u(-\infty, t) = 0, \quad u(\infty, t) = 1,$$

where *a* is transaction cost rate, *E* is the strike price, *K* = <sup>2</sup>*<sup>ρ</sup> σ*2 0 , *ρ* is the risk-free interest rate, *σ*<sup>0</sup> is the volatility of the stock, and Φ is a function defined as the solution to the following ODE:

$$\Phi'(s) = \frac{\Phi(s) + 1}{2\sqrt{s\Phi(s)} - s}, \quad \Phi(0) = 0,\tag{30}$$

Fig. 6. Solution curves for option price corresponding to various transaction cost rates(*a* = 0.0, 0.01, 0.02) and the pay-off curve.

In the numerical solution process, the infinite domain is approximated by a finite interval of sufficient length. In financial industry, one of the widely used methods for calculating hedge delta is to first solve the Black-Scholes equation to obtain numerical solution of *u*, and then apply finite difference schemes to the numerical solution to obtain approximations to *ux*. With the new method discussed in this chapter, both the solution *u* and the hedge delta *ux* are calculated simultaneously. Fig. 6 and Fig. 7 show the option price and the hedge delta for various transaction cost rates, respectively.

**Example 4**: In this example we solve a time dependent nonlinear equation, the Burgers' Equation:

$$
\mu\_t + \mu u\_\mathbf{x} = a u\_{\mathbf{x}\mathbf{x}} \quad \text{u}(0, t) = 0, \quad \mathbf{u}(1, t) = 1,\tag{31}
$$

shown in Fig. 10, and agree with the results reported in other works, such as (Hassanien et al.,

<sup>93</sup> A Fourth-Order Compact Finite Difference Scheme

Fig. 9. Solution curves corresponding to *a* = 0.1, 0.01, 0.001, 0.0001 using the upwind scheme

Fig. 10. Solution curves corresponding to *a* = 0.1, 0.01, 0.001, 0.0001 using the second-order

*<sup>u</sup>*(*x*, 0) = <sup>2</sup>*a<sup>π</sup>* sin(*πx*)

*<sup>u</sup>*(*x*, *<sup>t</sup>*) = <sup>2</sup>*aπe*−*π*<sup>2</sup> *at* sin(*πx*)

with *κ* > 1, for which the exact solution to Eq. (31) is known as

*κ* + cos(*πx*)

*κ* + *e*−*π*<sup>2</sup> *at* cos(*πx*)

, (32)

. (33)

2005).

for Solving Unsteady Convection-Diffusion Equations

with Δ*x* = 0.01.

new algorithm with Δ*x* = 0.01.

Next we use the initial condition

Fig. 7. Solution curves of the hedge delta *<sup>∂</sup><sup>V</sup> <sup>∂</sup><sup>S</sup>* corresponding to various transaction cost rate(*a* = 0.0, 0.01, 0.02).

where *a* is a constant. We first solve the equation with the initial condition *u*(*x*, 0) = *x*(1 − *x*), for which the analytic solution is not available. Figs. 8, 9, and 10 show the solution curves for various *a* values obtained by the standard central difference scheme (for both the convection and diffusion terms), the upwind (for convection)-central difference (for diffusion) scheme, and the new method using second-order central difference, respectively. It is clear that as *a* decreases, for instance, *a* = 0.001 , the standard central difference scheme produces oscillations in Fig. 8.

Fig. 8. Solution curves corresponding to *a* = 0.1, 0.01, 0.001, 0.0001 using the standard central difference scheme with Δ*x* = 0.01.

The new method, on the other hand, produces solutions that are oscillation free as shown in Fig. 9. These solutions are similar to those produced by the upwind-central difference scheme 12 Will-be-set-by-IN-TECH

where *a* is a constant. We first solve the equation with the initial condition *u*(*x*, 0) = *x*(1 − *x*), for which the analytic solution is not available. Figs. 8, 9, and 10 show the solution curves for various *a* values obtained by the standard central difference scheme (for both the convection and diffusion terms), the upwind (for convection)-central difference (for diffusion) scheme, and the new method using second-order central difference, respectively. It is clear that as *a* decreases, for instance, *a* = 0.001 , the standard central difference scheme produces

Fig. 8. Solution curves corresponding to *a* = 0.1, 0.01, 0.001, 0.0001 using the standard central

The new method, on the other hand, produces solutions that are oscillation free as shown in Fig. 9. These solutions are similar to those produced by the upwind-central difference scheme

*<sup>∂</sup><sup>S</sup>* corresponding to various transaction cost

Fig. 7. Solution curves of the hedge delta *<sup>∂</sup><sup>V</sup>*

rate(*a* = 0.0, 0.01, 0.02).

oscillations in Fig. 8.

difference scheme with Δ*x* = 0.01.

shown in Fig. 10, and agree with the results reported in other works, such as (Hassanien et al., 2005).

Fig. 9. Solution curves corresponding to *a* = 0.1, 0.01, 0.001, 0.0001 using the upwind scheme with Δ*x* = 0.01.

Fig. 10. Solution curves corresponding to *a* = 0.1, 0.01, 0.001, 0.0001 using the second-order new algorithm with Δ*x* = 0.01.

Next we use the initial condition

$$u(\mathbf{x},0) = \frac{2a\pi\sin(\pi\mathbf{x})}{\kappa + \cos(\pi\mathbf{x})},\tag{32}$$

with *κ* > 1, for which the exact solution to Eq. (31) is known as

$$u(\mathbf{x},t) = \frac{2a\pi e^{-\pi^2at}\sin(\pi\mathbf{x})}{\kappa + e^{-\pi^2at}\cos(\pi\mathbf{x})}.\tag{33}$$

(Gustafsson et al., 1995) should be used to first decompose the original equation into a series of one-dimensional problems. The new method discussed in this chapter can then be applied to these one-dimensional problems to calculate numerical solutions efficiently. Details will be

<sup>95</sup> A Fourth-Order Compact Finite Difference Scheme

Barles, G. & Soner, H. (1998). Option pricing with transaction costs and a nonlinear

Gu, Y.; Liao, W. & Zhu, J. (2003). An efficient high order algorithm for solving systems of 3-D

Gupta, M.; Manohar, R. & Stephenson, J. (1984). A single cell high-order scheme for the

Gustafsson, B.; Kreiss, H. & Oliger, J. (1995). *Time-Dependent problems and difference methods*,

Hassanien, I; Salaman, A. & Hosham, H. (2005). Fourth-Order finite difference method for

Hundsdorfer, W. & Verwer, J. (2003). *Numerical solution of time-dependent*

Kalita, J.; Dalal, D. & Dass, A. (2002). A class of higher order compact schemes for the

Karaa, S. & Zhang, J. (2004). High order method for solving unsteady convection-diffusion problems, *Journal of Computational Physics*, Vol. 198, no.1, 2004, pp. 1 - 9. Morton, K. (1996). *Numerical solution of convection-diffusion problems*, Chapman & Hall,

Noye, B. & Tan, H. (1989). Finite difference methods for solving the two-dimensional

Rigal, A. (1989). Numerical analysis of two-level finite difference schemes for unsteady

Rigal, A. (1990). Numerical analysis of three-time-level finite difference schemes for

Rigal, A. (1994). High order difference schemes for unsteady one-dimensional

Spotz, W. & Carey, G. (2001). Extension of high-order compact schemes to time-dependent

reaction-diffusion equations, *Journal of Computational and Applied Mathematics*, Vol.

convection-diffusion equation with variable coefficients, *International Journal for*

solving Burgers' equation, *Applied Mathematics and Computation*, Vol. 170, 2005, pp.

unsteady two-dimensional convection-diffusion equation with variable convection coefficients, *International Journal for Numerical Methods in Fluids*, Vol. 38, 2002, pp.

advection-diffusion equation, *International Journal for Numerical Methods in Fluids*,

diffusion-convection problems, *International Journal for Numerical Methods in*

unsteady diffusion-convection problems, *International Journal for Numerical Methods*

diffusion-convection problems, *Journal of Computational Physics*, Vol. 114, 1994,

problems, *Numerical Methods for Partial Differential Equations*, Vol. 17, 2001, pp. 657 -

Black-Scholes equation, *Finance Stochast*, Vol. 2, 1998, pp. 369 - 397.

*Numerical Methods in Fluids*, Vol. 4, 1984, pp. 641 - 651.

*advection-diffusion-reaction equations*, Springer, Berlin.

presented in future papers.

781 - 800.

1111 - 1131.

London.

pp. 59 - 76.

672.

Vol. 9, 1989, pp. 75 - 98.

*Engineering*, Vol. 28, 1989, pp. 1001 - 1021.

*in Engineering*, Vol. 30, 1990, pp. 307 - 330.

Seydel, R. (2002). *Tools for computational finance*, Springer, Berlin.

155, 2003, pp. 1 - 17.

John Wiley & Sons, New York.

for Solving Unsteady Convection-Diffusion Equations

**7. References**

We apply both the upwind-central scheme and the second-order new method to this example to compare the accuracy of the two algorithms. Table 3 shows the error �*u<sup>e</sup>* <sup>−</sup> *<sup>u</sup>c*�<sup>∞</sup> between the exact solution *u<sup>e</sup>* and the numerical solution *u<sup>c</sup>* calculated using the two algorithms with *κ* = 1.2, for *a* = 0.001 and *a* = 0.0001. All results are obtained using a 3-point stencil. Since the main focus of the study is to compare spatial accuracy of the two algorithms, the first-order explicit time integration is used for simplicity with Δ*t* = 0.0001 to ensure stability. It is clear from Table 3 that the accuracies of the two schemes are as expected: When the grid size Δ*x* is reduced by <sup>1</sup> <sup>2</sup> , the errors are reduced by approximately <sup>1</sup> <sup>2</sup> and ( <sup>1</sup> <sup>2</sup> )<sup>2</sup> for the upwind-central and the 2*nd*−new schemes, respectively.


Table 3. Errors between the exact solution and the numerical solution of Eq. (31) with *a* = 0.001 and *a* = 0.0001.

The numerical examples presented in this chapter show that the standard central difference scheme is second-order accurate on a 3-point stencil but produces oscillatory solutions for convection dominated problems. The upwind scheme is more robust but is only first-order accurate on a 3-point stencil. The new method discussed in this chapter appears to combine the advantages of accuracy of the standard central difference algorithm with the robustness of the upwind scheme for convection dominated equations.
