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

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. Thismethod 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 twoand 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 oneand 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). 4


Will-be-set-by-IN-TECH
For simplicity, we will use the following one-dimensional equation to describe the new method developed in this chapter: 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): where i = 0, 1, ··· , M,a n dn = 0, 1, ··· , N are indices for spatial and temporal grid points, respectively, u n i is the numerical approximation to the exact solution u(x i , t n ), and the central difference operator δ 2 is defined as δ 2 u n i = u n i−1 − 2u n i + u n i+1 . Applying the operator (1 + δ 2 12 ) to both sides of Eq. (2), we have Eq.
(3) is fourth-order accurate in space but contains only the second-order central difference operator δ 2 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 where the central difference operator δ is defined as δu n i = u n i+1 − u n i−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 u x from Eq. (1) and solves v = u x 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 82 Computational Simulations and Applications www.intechopen.com of convection-diffusion equations. It is unconditionally stable, and is particularly suitable for problems that require the solution of both u and u x . For example, when solving the Black-Scholes option pricing model (Seydel, 2002) both the solution V and its derivative V S are desired. The solution V is the price of an option and its derivative V S 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.

The new method
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.

Description of the new method
Setting v = u x in Eq.
(1), we have Differentiating both sides of Eq. (1) with respect to x leads to which, considering v = u x , can be written as Eq. (5) and (7) now form a system of equations for u and v. They only involve diffusion term u xx and v xx , 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 Applying the operator (1 + δ 2 12 ) to both sides of these two equations, we obtain Eqs. (10) and (11) are fourth-order accurate in space but only contains the second-order central difference operator δ 2 that is defined on a 3-point stencil. As a result, the discretized equations form a system of block tri-diagonal algebraic equations.

Initial and boundary conditions
The initial condition for v can be obtained by differentiating h(x), the initial condition fro u giveninEq. ( 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:S in c ev = u x ,wehave 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 = u x . Padé approximation: We can use the fourth-order Padé approximation at Applying (1 + δ 2 6 ) to both sides of Eq. (13), we have or Solving v 0 from Eq. (15), we obtain the boundary condition where a, c and d are constants. Integrating this equation from x 0 to x 1 ,wehave or Solving v 0 from Eq. (19), we obtain the boundary condition 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.

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 The new method combined with the second-order central difference for solving Eq. (21) is Taking the discrete Fourier transform of the above equations, we have Thus where , with α = ∆t ∆x 2 , β = c∆t,andγ = c∆t ∆x 2 . 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 θ.Fort h e2× 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 , thus we have

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 with exact solution u(x)= e x −1 e−1 . 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 u i , 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 , ··· , u n ).F o r t h e fourth-order new method the matrix A is a 2n × 2n banded matrix with a bandwidth of 7, and U =( u 1 , v 1 , u 2 , v 2 , ··· , u n , v n ) contains the numerical solutions for both u and v,w i t h v = u x . If an efficient algorithm is used to solve the tridiagonal linear system resulted from upwind and standard central finite difference schemes,a total of 2n − 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)=8n − 7 = 8n + 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 18n − 15 divisions, 24n − 22 multiplications, and 24n − 22 subtractions are needed. The computational complexity is 66n − 59 = 66n + 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.0E-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. first-order upwind scheme, it is almost impossible to obtain an error less than 1.0E − 08 since agridsizeofh = 1.0E − 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.

Numerical results
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: where a and c are constants. The exact solution is given as T able2.showstheerror u e − u c ∞ between the exact solution u e and the numerical solution u c 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 finite difference for the diffusion term.
2. 2 nd -stand: Second-order central difference for both the diffusion and convection terms.
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 1 2 ,t h ee r r o r s are reduced by approximately 1 2 , ( 1 2 ) 2 and ( 1 2 ) 4 for the 1 st -, 2 nd -a n d4 th -order algorithms, respectively. 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 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 = 0a n dx = 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      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.

www.intechopen.com
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).
where a is transaction cost rate, E is the strike price, K = 2ρ σ 2 0 , ρ is the risk-free interest rate, σ 0 is the volatility of the stock, and Φ is a function defined as the solution to the following ODE: 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 u x . With the new method discussed in this chapter, both the solution u and the hedge delta u x 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:  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. 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 shown in Fig. 10, and agree with the results reported in other works, such as (Hassanien et al., 2005).  Next we use the initial condition with κ > 1, for which the exact solution to Eq. (31) is known as u(x, t)= 2aπe −π 2 at sin(πx) κ + e −π 2 at cos(πx) .

A Fourth-Order Compact Finite Difference Scheme for Solving Unsteady Convection-Diffusion Equations
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 e − u c ∞ between the exact solution u e and the numerical solution u c 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 1 2 , the errors are reduced by approximately 1 2 and ( 1 2 ) 2 for the upwind-central and the 2 nd −new schemes, respectively.
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.

Conclusions
The method discussed in this chapter eliminates the convection term in Eq. (1) and makes it feasible to use central difference schemes to solve convection-diffusion equations accurately. The new method, combined with the central difference schemes, can achieve better accuracy than the upwind schemes on the same finite difference stencil, and is shown in the examples presented here to be as robust as the upwind schemes for convection dominated problems. It can also be easily combined with the Padé approximation to achieve fourth-order accuracy in space on a 3-point finite difference stencil. The new method does incur a modest increase in computational complexity. Instead of just solving Eq. (1), the new method requires solving Eqs. (5) and (7). With a 3-point stencil, the standard upwind-central difference schemes will generate a system of tri-diagonal algebraic equations, while the new method discussed in this chapter will lead to a system of block tri-diagonal algebraic equations with 2 × 2 blocks. This increased complexity, however, can be compensated by the use of fewer grid points with the increased order of accuracy of the new method. Furthermore, for problems that require calculations of both the solutions and their derivatives, the new method eliminates the need to calculate the derivatives after solving Eq. (1). The discussions of the new method in this chapter are based on one-dimensional problems. For higher-dimensional problems, a straightforward application of this method will lead to systems of four equations for two-dimensional problems and seven equations for three-dimensional problems. For better computational efficiency, operator splitting 94 Computational Simulations and Applications www.intechopen.com