An implicit fourth-order compact finite difference scheme for one-dimensional Burgers’ equation

Applied Mathematics and Computation 206 (2008) 755–764

Contents lists available at ScienceDirect

Applied Mathematics and Computation journal homepage: www.elsevier.com/locate/amc

An implicit fourth-order compact finite difference scheme for one-dimensional Burgers’ equation Wenyuan Liao Department of Mathematics and Statistics, The University of Calgary, 2500 University Drive NW, Calgary, Alberta, Canada T2N 1N4

a r t i c l e

i n f o

a b s t r a c t A fourth-order compact finite difference method is proposed in this paper to solve onedimensional Burgers’ equation. The newly proposed method is based on the Hopf–Cole transformation, which transforms the original nonlinear Burgers’ equation into a linear heat equation, and transforms the Dirichlet boundary condition into the Robin boundary condition. The linear heat equation is then solved by an implicit fourth-order compact finite difference scheme. A compact fourth-order formula is also developed to approximate the Robin boundary conditions, while the initial condition for the heat equation is approximated using Simpson’s rule to maintain overall fourth-order accuracy. Numerical experiments have been conducted to demonstrate the efficiency and high-order accuracy of this method. The numerical results also show that the method is unconditionally stable, as there is no constraint on time step size. Ó 2008 Elsevier Inc. All rights reserved.

Keywords: Burgers’ equation Fourth-order Compact Finite difference Implicit Unconditionally stable

1. Introduction In this paper, we consider the one-dimensional Burgers’ equation:

ou ou o2 u þu ¼m 2; ot ox ox

x 2 ð0; 1Þ; t 2 ð0; T

ð1Þ

with initial condition

uðx; 0Þ ¼ f ðxÞ;

0
ð2Þ

and boundary conditions

uð0; tÞ ¼ aðtÞ;

uð1; tÞ ¼ bðtÞ;

t 2 ð0; T;

1 Re

ð3Þ

where m ¼ and Re is the Reynolds number characterizing the size of viscosity. Burgers’ equation (1) is an important and basic partial differential equation from fluid mechanics, and has been widely used for various applications, such as modeling of gas dynamics and traffic flow, describing wave propagation in acoustics and hydrodynamics, etc. The first attempt to analytically solve the one-dimensional Burgers’ equation was by Bateman [2], who derived the steady solution for (1), which had been used Burger in [5] to model turbulence. For a set of special initial and boundary conditions, the Burgers’ equation can be analytically solved. In the context of gas dynamics, Hopf [12] and Cole [7] have independently shown that, for any initial condition, the Burgers’ equation can be reduced to a linear homogeneous heat equation that can be solved exactly, thus the exact solution to the original Burgers’

E-mail address: [email protected] 0096-3003/$ - see front matter Ó 2008 Elsevier Inc. All rights reserved. doi:10.1016/j.amc.2008.09.037

756

W. Liao / Applied Mathematics and Computation 206 (2008) 755–764

equation can be expressed in the form of a Fourier series. Based on the previous results, Benton and Platzman [3] surveyed the exact solution of the one-dimensional Burgers’ equation. Though the exact solution is available in the form of a Fourier series, an efficient numerical method is still required, especially when the initial condition is not smooth, or is only available on discrete locations. Another issue with the Fourier series solution is its slow convergence, therefore, a considerably long series is required to obtain an accurate approximation. In the past several decades, numerical solutions to Burgers’ equation have attracted a lot of attention from both scientists and engineers and this has resulted in various finite difference, finite-element and boundary element methods for solving problems such as (1), initial condition (2) and boundary conditions (3). To name a few, Varoglu and Finn [22] developed an isoperimetrical space–time, finite-element method; Nguyen and Reynen [18] proposed a method based on the least-squares weak-formulation of the finite-element method; and Caldwell and Smith [6] developed a finite-element method that features adaptive element size at different time stages. Depending on the time integration techniques used, the finite difference methods for solving time-dependent problems fall into three categories: explicit, implicit and semi-implicit. For a comparison of various early works on the finite difference method readers are referred to [4]. Recently, several numerical finite difference schemes had been developed. For instance, Kutluay et al. [14] have developed explicit and exact-explicit finite difference methods based on the Hopf–Cole transformation. Ozis and Aslan [19] have proposed a semi-approximate approach for solving (1) with high Reynolds number, while an unconditionally stable fourth-order finite difference method has been proposed by Hassanien et al. [10]. For a more detailed survey of other numerical methods, we refer to [8,13,14,17,20]. For the sake of stability, an implicit algorithm is preferred, especially when the computation is conducted over a long time interval. However, implicit methods are usually not very efficient, since a linear or nonlinear algebraic system (for a nonlinear PDE such as Burgers’ equation) must be solved at each time step. On the other hand, explicit schemes are very efficient but have poor stability, which usually enforces a very small time step size. An alternative is to treat the nonlinear term explicitly, which leads to a semi-implicit scheme. When an implicit scheme is used to solve a linear partial differential equation, one needs to solve a linear algebraic system Au = b at each time step. If a compact scheme is used for a one-dimensional problem, the matrix A is tri-diagonal, however, if a non-compact scheme is used, the matrix A will be a band matrix with a bandwidth greater than one. Solving a tri-diagonal linear algebraic system is more efficient than solving a band linear algebraic system. There are at least two reasons why compact schemes are preferred: the linear system can be solved very efficiently and no extra points are needed near the boundary. In this paper, a high order (fourth-order in both time and space variables) compact (three-point stencil) scheme for solving Burgers’ equation (1), with initial condition (2) and boundary condition (3), is proposed. The rest of the paper is organized as follows. In Section 2, we briefly review some existing numerical methods related to the scheme in this paper. The fourthorder compact scheme is described in Section 3, followed by numerical experiments in Section 4. Finally concluding remarks are presented in Section 5. 2. Review of existing algorithms The domain [0, 1][0, T] is divided into an N  M mesh, with h ¼ N1 and Dt = T/M, respectively. xi = i  h, for i = 0, 1, 2, . . . , N is the ith node. In the temporal dimension, the uniform step size Dt is used, thus tj = j  Dt is the time level for jth step. The quantity u(xi, tj) represents the exact solution at (xi, tj) while uji represents the numerical solution at (xi, tj). A standard explicit central finite difference method for solving (1) is defined as

uj  uji1 ujþ1  uji i ¼ uji iþ1 Dt 2h

! þm

ujiþ1  2uji þ uji1 2

h

ð4Þ

with a truncation error O(Dt) + O(h2). It is well-known that this scheme is only conditionally stable and a strict constraint on time step size Dt is usually enforced. The Crank–Nicolson algorithm for solving (1) is given as

! ! " # " # jþ1 jþ1 jþ1 jþ1 j j ujþ1 þ ui1 ujiþ1  2uji þ uji1 ujþ1  uji 1 1 jþ1 uiþ1  ui1 j uiþ1  ui1 iþ1  2ui i þm þm þ : ui ui ¼ 2 2 2 2 Dt 2h 2h h h

ð5Þ

This scheme is obviously unconditionally stable and with higher order truncation error O(Dt2) + O(h2). However, it involves solving a nonlinear algebraic system which makes it inefficient in practice. To resolve this issue, we can treat part of the nonlinear convection term explicitly, resulting in the following semi-implicit scheme:

! ! " # " # jþ1 jþ1 jþ1 j j uiþ1  2ujþ1 þ ujþ1 ujiþ1  2uji þ uji1 ujþ1  uji 1 1 j uiþ1  ui1 j uiþ1  ui1 i i1 i þm þm þ ui ui ¼ 2 2 2 2 Dt 2h 2h h h

ð6Þ

or one could treat the whole convection term explicitly

! ! " # " # j j jþ1 j j uiþ1  2ujþ1 þ ujþ1 ujiþ1  2uji þ uji1 ujþ1  uji 1 1 j uiþ1  ui1 j uiþ1  ui1 i i1 i þm þm þ : ui ui ¼ 2 2 2 2 Dt 2h 2h h h

ð7Þ

W. Liao / Applied Mathematics and Computation 206 (2008) 755–764

757

Another semi-implicit scheme proposed by Ozis in [19] is given as jþ1 ujþ1  ui1 ujþ1  uji i ¼ uji iþ1 Dt 2h

! þm

jþ1 ujþ1 þ ujþ1 iþ1  2ui i1

h

2

;

ð8Þ

which has a truncation error O(Dt) + O(h2). In [14], the authors have proposed a numerical scheme based on the Hopf–Cole transformation which, instead of solving the original Burgers’ equation, first solves the linear heat equation using an explicit finite difference method, and then applies the transformation to obtain the solution of the original Burgers’ equation. The new algorithm presented in this paper is also based on the Hopf–Cole transformation but the linear heat equation is solved by an implicit, compact and fourth-order finite difference scheme. As shown later by numerical experiments, this is very efficient and accurate, and is unconditionally stable. 3. Fourth-order finite difference scheme The Burgers’ equation (1) can be transformed into a heat equation using the Hopf–Cole transformation [12,7]:

uðx; tÞ ¼ 2m

wx ðx; tÞ : wðx; tÞ

ð9Þ

It is easy to verify that w(x, t) satisfies the following heat equation:

ow o2 w ¼m 2 ; ot ox

0 < x < 1; t 2 ð0; T

ð10Þ

with initial condition

 Z wðx; 0Þ ¼ exp  0

x

 f ðsÞ ds ; 2m

06x61

ð11Þ

and Robin boundary conditions

2mwx ð0; tÞ þ aðtÞwð0; tÞ ¼ 0;

t > 0;

ð12Þ

2mwx ð1; tÞ þ bðtÞwð1; tÞ ¼ 0;

t > 0:

ð13Þ

In the following sections, we will derive the compact fourth-order scheme for solving (10), the fourth-order approximations for initial condition (11) and boundary conditions (12) and (13), respectively. 3.1. Fourth-order scheme for solving (10) The well-known Crank–Nicolson algorithm had been widely used to solve (10) because of its unconditional stability and second order accuracy. Applying the Crank–Nicolson algorithm to (10), we obtain

winþ1  wni m ¼ 2 ðd2x winþ1 þ d2x wni Þ; Dt 2h

ð14Þ

where d2x is the second-order central finite difference operator, which is defined as: d2x wi ¼ wiþ1  2wi þ wi1 . The Crank–Nicolson algorithm is compact and unconditionally stable but is only second-order accurate in both the temporal and spatial dimensions [9]. One way to obtain a fourth-order approximation is to use Padé approximation, which is imple1 2 2 dx Þdx , thus the following fourth-order scheme is mented by replacing the central finite difference operator d2x with ð1  12 obtained:

        1 2 2 1 2 2 1r 1 ¼ 1 þ r 1  dx dx wnþ1 d d wni ; i 12 12 x x

ð15Þ

1 2 2 where r ¼ 2m hD2t. However the scheme requires a five-point stencil. To maintain a compact three-point stencil, ð1  12 dx Þdx is 2 approximated by 1þd1x d2 , which is fourth-order accurate and eventually leads to the following fourth-order and compact x 12 scheme:

1r

d2x 1 2 1 þ 12 dx

! wnþ1 i

¼

1þr

d2x 1 2 1 þ 12 dx

! wni :

ð16Þ

2

1 2 In real computations, the operator 1þd1x d2 is not implemented, so we multiply 1 þ 12 dx to both sides of (16) and obtain the fol12 x

lowing compact and fourth-order scheme:

    1 2 1 2 2 n 1þ ¼ 1 þ þ rd dx  rd2x wnþ1 d x wi ; i 12 12 x

i ¼ 1; 2; . . . ; N:

ð17Þ

758

W. Liao / Applied Mathematics and Computation 206 (2008) 755–764

3.2. Fourth-order approximation of the initial condition be expressed explicitly. Thus, the fourth-order Simpson’s integration rule is In general, the anti-derivative of f2ðxÞ m cannot Rx used here to calculate wðxi ; 0Þ ¼ expf 0 i f2ðsÞ m dsg, for i = 1, 2, 3, . . . , N. w(x0, 0) = 1. When i = 1, the following fourth-order approximation is used:

  h wðx1 ; 0Þ ¼ exp  ðf ðx0 Þ þ 4f ðx1 Þ þ f ðx1 ÞÞ ; 2 12m

ð18Þ

x0 þx1 2

where x1 ¼ is the middle point of [x0, x1]. 2 For i = 2, 3, . . . , N, based on the following integration decomposition:

Z 0

xi

f ðsÞ ds ¼ 2m

Z

xi2

f ðsÞ ds þ 2m

0

Z

xi

xi2

f ðsÞ ds 2m

ð19Þ

the initial condition is approximated as

  h wðxi ; 0Þ ¼ wðxi2 ; 0Þ þ exp  ðf ðxi2 Þ þ 4f ðxi1 Þ þ f ðxi ÞÞ : 6m

ð20Þ

3.3. Compact fourth-order approximation for boundary conditions As can be seen, the heat equation is associated with Robin boundary conditions (12) and (13), so a compact fourth-order approximation for the boundary condition is necessary to maintain an overall fourth-order accuracy. We will use the left boundary to demonstrate the fourth-order approximation. At x = 0, 2mwx(0, t) + a(t)w(0, t) = 0. We first 2 x w0 þ Oðh Þ. Here use a second order central finite difference operator to approximate the first order derivative, i.e. wx ð0; tÞ  D2h Dx is a finite difference operator defined as Dxwi = wi+1  wi1. Therefore, the second order approximation for the left boundary condition is given as

2mDx wðx0 ; tÞ þ aðtÞwðx0 ; tÞ ¼ 0: 2h

ð21Þ

The second order approximation defined above can be improved to fourth-order by substituting Dx with 1þD1xd2 , which leads to 6 x

mðwðx1 ; tÞ  wðx1 ; tÞÞ   h 1 þ 16 d2x

þ aðtÞwðx0 ; tÞ ¼ 0:

ð22Þ

Then, apply the method introduced in [16], and solve for w(x1, t), we obtain the following boundary condition:



wðx1 ; tÞ ¼

m þ aðtÞ

h



6

wðx1 ; tÞ þ 23 aðtÞwðx0 ; tÞ m  aðtÞ

h

ð23Þ

:

6

Similarly, we can derive the approximation for the boundary condition at x = 1 as



wðxNþ1 ; tÞ ¼



m  bðtÞ

h

6

wðxN1 ; tÞ  23 bðtÞwðxN ; tÞ m þ bðtÞ

h

ð24Þ

:

6

where x1 = h and xN+1 = 1 + h are two virtual grid points from outside of the left and right sides respectively. w1 and wN+1 can be used to complete the linear system (17). The two functions a(t) and b(t) are Dirichlet boundary conditions of u(x, t), i.e, u(0, t) = a(t) and u(1, t) = b(t), respectively. For the numerical examples in this paper, we assume that both functions are zero, so the boundary conditions we actually obtain are w(x1, t) = w(x1, t) and w(xN+1, t) = w(xN1, t). Combining the algorithm (17), initial conditions (18) and (20), and the boundary conditions (23) and (24), the problem can be solved efficiently with fourth-order accuracy in space. 3.4. Richardson’s extrapolation It is well-known that the truncation error of the Crank–Nicolson algorithm in the temporal dimension has a form of c1Dt2 + c2Dt4, so it is natural to apply the Richardson’s extrapolation here to improve accuracy in the temporal dimension. Assume wDt is the numerical solution obtained by using Dt, while wDt is the solution obtained using time step size D2t, the new solution of Richardson’s extrapolation is defined as w ¼ order accurate in the temporal dimension.

4wDt wDt 2

3

2

, which has an error of the form  14 c2 Dt4 and thus is fourth-

3.5. Compact fourth-order approximation for wx x ðx;tÞ Note that, uðx; tÞ ¼ 2m wwðx;tÞ , to reconstruct u(x, t), accurate numerical solutions for both w(x, t) and wx(x, t) are required. Here, we first solve the heat equation (10) with fourth-order accuracy, then solve for wx(x, t) by the compact fourth-order

759

W. Liao / Applied Mathematics and Computation 206 (2008) 755–764

method (26). wx(x, t) must be solved with a fourth-order method, otherwise the overall order of accuracy for solving u(x, t) will be reduced. As u(0, t) and u(1, t) are already specified, we need only to approximate wx(xi, t) for i = 1, 2, . . . , N  1. Again, we apply the second order central finite difference operator Dx on w(xi, t), then replace it with 1þD1xd2 , which leads to the fol6 x lowing fourth-order accurate scheme:

wx ðxi ; tÞ ¼

Dx   wðxi ; tÞ; 2h 1 þ 16 d2x

i ¼ 1; 2; . . . ; N  1:

ð25Þ

We multiply 1 þ 16 d2x to both sides of (25), resulting in the following linear system:

  1 Dx 1 þ d2x wx ðxi ; tÞ ¼ wðxi ; tÞ; 6 2h

i ¼ 1; 2; . . . ; N  1;

ð26Þ

which can be solved efficiently with fourth-order accuracy, and the solution wx will be used to reconstruct the solution u(x, t). 4. Numerical results 4.1. Example 1 We first solve the Burgers’ equation (1) and boundary conditions

uð0; tÞ ¼ uð1; tÞ ¼ 0;

t>0

ð27Þ

with the following initial condition:

uðx; 0Þ ¼ 2m

p sinðpxÞ ; x 2 ð0; 1Þ; r þ cosðpxÞ

ð28Þ

where r > 1 is a parameter. Here the initial condition is chosen such that w(x, 0) can be derived and expressed explicitly as

wðx; 0Þ ¼ r þ cosðpxÞ;

x 2 ½0; 1

ð29Þ

and the closed form of the exact solution is available as [23] 2

uðx; tÞ ¼

2mpep mt sinðpxÞ ; r þ ep2 mt cosðpxÞ

0 < x < 1:

ð30Þ

With the closed form of the exact solution available, we can calculated the numerical error and convergence order of the present method. In Table 1, the maximal error between the numerical solution by the higher order method without Richardson’s extrapolation and the exact solution is calculated for r = 2, m = 0.1, T = 1.0. To demonstrate the order of accuracy in the spatial dimension, Dt is set as 0.000001 such that the effect of the truncation error from time discretization can be ignored. Clearly, when h is reduced by a factor of 2, the error is reduced by a factor of 16, which confirms the conclusion that the algorithm is fourth-order accurate in the spatial dimension. To show that the method is only second-order accurate in the temporal dimension without Richardson’s extrapolation, we fix h = 0.001 and reduce Dt by a factor of 2. The maximal errors with various Dt at T = 1.0 for m = 0.1, r = 2 are displayed in Table 2, which confirms the conclusion claimed above.

Table 1 Maximal error ku  uexactk1 between the exact solution (30) and the numerical solution by the fourth-order finite-difference scheme at T = 1 for r = 2, m = 0.1 and Dt = 0.000001 h E(h) EðhÞ r ¼ Eðh=2Þ r log2

1/10 1.5763e05 – –

1/20 9.8231e07 16.0473 4.0042586

1/40 6.1273e08 16.0318 4.0028645

1/80 3.8303e09 15.9969 3.9997205

1/160 2.3954e10 15.9901 3.9991071

Table 2 Maximal error ku  uexactk1 between the exact solution (30) and the numerical solution by the fourth-order finite-difference scheme without Richardson’s extrapolation at T = 1 for m = 0.1, r = 2 and h = 0.001

Dt E(Dt) tÞ r ¼ EðEðDDt=2Þ r log 2

1/16 3.9298e05 – –

1/32 9.8217e06 4.001177 2.000424

1/64 2.4552e06 4.000293 2.000106

1/128 6.1380e07 4.00007 2.00003

1/256 1.5345e07 4.000017 2.000006

760

W. Liao / Applied Mathematics and Computation 206 (2008) 755–764

To show that the method has an overall fourth-order convergence rate with Richardson’s extrapolation, we initially set h = 1/20 and Dt = 1/16, then reduced them both by a factor of 2 each time. The results displayed in Table 3 confirm that the new method is fourth-order accurate in both the temporal and spatial dimensions. 4.2. Example 2 In this example, we solve the same Burgers’ equation (1) and boundary condition (27), but with a different initial condition

uðx; 0Þ ¼ sinðpxÞ;

x 2 ½0; 1:

ð31Þ

This example was also solved by other numerical methods and the detailed numerical results can be found in [1,11]. For the initial condition defined above, the initial condition w(x, 0) can be expressed explicitly as

  1  cosðpxÞ wðx; 0Þ ¼ exp  ; 2pm

x 2 ½0; 1;

ð32Þ

so the high-order approximation described in Section 3.2 is not necessary. Using the Hopf–Cole transformation, the exact Fourier solution to the problem is available as

uðx; tÞ ¼ 2pm

P1 2 2 n¼1 c n expðn p mtÞn sinðnpxÞ P ; 2 p2 mtÞ cosðnpxÞ c0 þ 1 c expðn n¼1 n

ð33Þ

where the coefficients are defined as

  1  cosðpxÞ dx; exp  2pm 0   Z 1 1  cosðpxÞ cn ¼ 2 cosðnpxÞdx; exp  2pm 0

c0 ¼

Z

1

ð34Þ ðn ¼ 1; 2; 3; . . .Þ:

ð35Þ

To calculate the error of the numerical solution, the exact Fourier series solution (33) needs to be evaluated. In real computation, a large number of terms should be evaluated to ensure high accuracy. In this example, a number N is chosen such that the coefficient cN is less than 1.0e15. To compare the present fourth-order scheme with the explicit finite difference scheme and the exact-explicit finite difference scheme proposed in [14], the same problem is solved at T = 0.1 for m = 1 and Dt = 0.00001. Clearly, the fourth-order algorithm is much more accurate. As shown in Table 4, at the nine grid points: x = 0.1, 0.2, . . . , 0.9, the numerical solution by the fourth-order method with h = 1/10 is more accurate than the numerical solutions using the other two algorithms (explicit finite difference and exact-explicit finite difference methods) with h = 1/80. The numerical solutions calculated by the explicit finite difference scheme and exact-explicit finite difference scheme can be found in Tables 1 and 2 of [14]. Table 5 shows that the new fourth-order scheme works well and results in accurate numerical solution with large Dt for which the explicit finite difference scheme is unstable. Table 3 Maximal error ku  uexactk1 between the exact solution (30) and the numerical solution by fourth-order finite-difference scheme with Richardson’s extrapolation at T = 1 for m = 0.1, r = 2 while h and Dt are reduced at the same rate simultaneously

E Eðh;DtÞ r ¼ Eðh=2; Dt=2Þ r log 2

h = 1/20 Dt = 1/16

h = 1/40 Dt = 1/32

h = 1/80 Dt = 1/64

h = 1/160 D t = 1/128

h = 1/320 Dt = 1/256

2.6121e07 – –

1.6179e08 16.1448427 4.0130015

1.0105e09 16.011112 4.0010016

6.3190e11 15.991599 3.99924

3.8831e12 16.2730177 4.0244099

Table 4 Comparison of the solution by the fourth-order finite-difference scheme and the exact solution provided by (33) at T = 0.1 for m = 1 and Dt = 0.00001 x

0.1 0.2 0.3 0.4 0.5 0.6 0.7 0.8 0.9

Numerical solution

Exact solution

h = 1/10

h = 1/20

h = 1/40

h = 1/80

0.1095377613 0.2097910776 0.2918941470 0.3479201929 0.3715721789 0.3590391476 0.3098983995 0.2278119285 0.1206835728

0.1095381295 0.2097920869 0.2918962199 0.3479236879 0.3715771535 0.3590451860 0.3099045949 0.2278170693 0.1206864989

0.1095381499 0.2097921451 0.2918963428 0.3479238985 0.3715774561 0.3590455555 0.3099049754 0.2278173856 0.1206866791

0.10953815119 0.20979214867 0.29189635032 0.34792391150 0.37157747490 0.35904557846 0.30990499905 0.22781740532 0.12068669034

0.1095381512705 0.2097921489100 0.2918963508255 0.3479239123656 0.3715774761468 0.3590455799850 0.3099050006311 0.2278174066274 0.1206866910894

761

W. Liao / Applied Mathematics and Computation 206 (2008) 755–764 Table 5 Comparison of the solution by the fourth-order finite difference scheme and the exact solution provided by (33) at T = 1 for m = 0.01 x

Numerical solution

0.1 0.2 0.3 0.4 0.5 0.6 0.7 0.8 0.9 0.925 0.95 0.975 1.00

Exact solution

1 Dx ¼ 20 1 Dt ¼ 16

1 Dx ¼ 40 1 Dt ¼ 32

1 Dx ¼ 80 1 Dt ¼ 64

1 Dx ¼ 160 1 Dt ¼ 128

0.0754108599 0.1506563858 0.2256097221 0.3002190934 0.3745366036 0.4487357325 0.5231375441 0.5983301760 0.6751424000 N/A 0.6881108787 N/A 0.0

0.075383758 0.150645871 0.225662289 0.300303279 0.374427793 0.447873253 0.520438673 0.591859310 0.660803353 0.672942844 0.659340057 0.511045515 0.0

0.0753820250 0.1506452228 0.2256656626 0.3003084997 0.3744205283 0.4478196667 0.5202789585 0.5914997849 0.6600672865 0.6719821691 0.6578703069 0.5089060830 0.0

0.0753819162 0.1506451824 0.2256658746 0.3003088255 0.3744200684 0.4478163167 0.5202690473 0.5914776670 0.6600223774 0.6719234253 0.6577797390 0.5087732703 0.0

0.07538190894 0.15064517967 0.22566588871 0.30030884722 0.37442003764 0.44781609332 0.52026838747 0.59147619875 0.66001946026 0.67191967836 0.65777406985 0.50876500576 0.0

It is well-known that one of the difficulties in solving the Burgers’ equation is that shock of the solution may occur after some time, even if the initial data is smooth. When the characteristic curves of the Burgers’ equation (1) cross, a shock of the solution occurs. If the initial condition u(x, 0) is given, the time when the shock forms can be determined exactly as 1 , as mentioned in [15]. From this formula, we can see that, if Ts is positive, the solution will break and form T s ¼ Minx ðu x ðx;0ÞÞ a shock at some point. On the other hand, if Ts is negative, the solution will not break, and no shock forms. A robust and accurate numerical algorithm should be able to capture the shock and the numerical solution should exhibit the correct physical behavior. From the last three rows in Table 5, we can see that a shock of the solution occurs. Numerical solutions and exact solutions at four points (x = 0.925, 0.95, 0.975, 1.0) are inserted in the table, so the shock is more visible from the data. It is 1 (the two worth pointing out that the numerical solutions at x = 0.925 and x = 0.975 are not available for the column Dx ¼ 20 blanks are marked as N/A). To obtain a better view of the shock, we plot the numerical solutions for m = 0.01 together with the exact solutions at different times in Fig. 1. As we can see from Table 5, the numerical solutions are very accurate, so it is

1

1 Numer.

0.6 0.4 0.2 0

Exact

0.8 u at T = 0.6

0.8

u at T = 0.2

Numer.

Exact

0.6 0.4 0.2

0

0.2

0.4

0.6

0.8

0

1

0

0.2

0.4

x 0.7

0.8

1

0.6

0.8

1

0.4 Numer.

Numer.

0.6

Exact

Exact

0.3 u at T = 2

0.5

u at T = 1

0.6

x

0.4 0.3 0.2

0.2

0.1

0.1 0

0

0.2

0.4

0.6

x

0.8

1

0

0

0.2

0.4

x

Fig. 1. Numerical and exact solutions of example 2 for m = 0.01, h = 0.0025 and Dt = 0.0001.

762

W. Liao / Applied Mathematics and Computation 206 (2008) 755–764

hard to tell the difference between the exact solutions and the numerical solution from the figure. Since the initial condition is given as u(x, 0) = sin(px), from the formula, we have T s ¼ p1  0:32. We can see from Fig. 1 that when t = 0.2 < Ts, the numerical solution and exact solution are smooth and no shock forms, however when t > Ts, both the numerical solution and exact solution form shocks at the location around x = 0.9. It is clear that the fourth-order method presented in this paper precisely captures the shock of the solution. 4.3. Example 3 For this example, we solve the Burgers’ equation (1) with boundary condition (27) but with the following initial condition:

uðx; 0Þ ¼ 4xð1  xÞ;

0 < x < 1:

ð36Þ

The initial condition for w(x, 0) can be expressed as

  3x2 þ 2x3 ; wðx; 0Þ ¼ exp  3m

0 < x < 1:

ð37Þ

The exact solution to this example is given by (33), but in this case the Fourier coefficients are

  3x2  2x3 dx; exp  3m 0   Z 1 3x2  2x3 cn ¼ 2 cosðnpxÞdx; exp  3m 0

c0 ¼

Z

1

ð38Þ ðn ¼ 1; 2; 3; . . .Þ:

ð39Þ

The solution obtained using the fourth-order finite difference scheme is compared with the solution obtained using the explicit finite difference scheme and the solution obtained using exact-explicit finite difference scheme from [14]. The results are displayed in Table 6. Note that the exact solution is available as a Fourier series, so the exact solution in the table is actually obtained by evaluating the Fourier series (33) with Fourier coefficients (38) and (39). In the real computation, the same h(0.0125) is used for all three schemes but time step size Dt is different. For the fourth-order method, Dt = 0.01 is used, while for the other two algorithms, Dt = 0.001 is used. Here different Dt0 s are chosen for the following three reasons. Firstly, with Dt = 0.01, the two algorithms in [14] are not stable, while the fourth-order method still produces an accurate numerical solution, so the fourth-order method is more robust than the other two methods. Secondly, we want to show that, even with larger Dt, the present fourth-order method produces numerical solutions that are more accurate than the solutions obtained using the other two algorithms. Thirdly, since Dx = 0.0125, and the new method is fourth-order accurate in both time and space, a smaller Dt (<0.01) does not make much difference in improving the overall accuracy, because the dominant error is from the discretization in space. To save computational time, a larger Dt is chosen. The numerical solutions for different times at different locations are displayed in Table 6, which shows that the fourthorder scheme is very accurate and robust.

Table 6 Comparison of solutions by exact-explicit finite-difference presented in [13], the solution by the fourth-order finite difference scheme and the exact solution 1 provided by (33) (with Fourier coefficients (38) and (39)) at different times for m = 0.01, with Dx ¼ 80 x

0.25

0.50

0.75

Tend

0.4 0.6 0.8 1.0 3.0 0.4 0.6 0.8 1.0 3.0 0.4 0.6 0.8 1.0 3.0

Numerical solution

Exact solution

Explicit

Exact-explicit

Fourth-order

0.36296 0.28217 0.23043 0.19463 0.07611 0.69591 0.55351 0.45625 0.38705 0.15220 0.95925 0.80197 0.67267 0.57501 0.22796

0.36185 0.28193 0.23046 0.19474 0.07617 0.67851 0.54508 0.45176 0.38446 0.15215 0.91169 0.77402 0.65617 0.56478 0.22746

0.36225863703815 0.28203607105989 0.23045089854492 0.19469030365258 0.07605681071165 0.68368800021020 0.54832361620351 0.45371620836465 0.38567653580000 0.15202591034708 0.92040036375162 0.78303771863201 0.66275532819945 0.56933797105110 0.22751146310186

0.36225937607348 0.28203659151197 0.23045114915252 0.19469040825968 0.07613409779556 0.68367860034256 0.54831636831404 0.45371356235593 0.38567577349503 0.15217998215782 0.92050204057761 0.78299398529618 0.66272038276880 0.56931867464330 0.22774304791098

763

W. Liao / Applied Mathematics and Computation 206 (2008) 755–764

0.7 T = 0.6 T=1 T=2 T=4

0.6 0.5

u

0.4 0.3 0.2 0.1 0

0

0.1

0.2

0.3

0.4

0.5

0.6

0.7

0.8

0.9

1

x Fig. 2. Numerical solutions of example four at different times.

4.4. Example 4 For this example, we solve the Burgers’ equation (1) with boundary condition (27) and the initial condition is defined as 2

uðx; 0Þ ¼

ex 3 þ cosðxÞ

ð40Þ

for which an explicit expression for w(x, 0) is not available, thus, a higher-order approximation based on Simpson’s rule, such as that given in Section 3.2 is used here. The solutions for different times are plotted in Fig. 2. 5. Conclusion An efficient fourth-order numerical algorithm based on the Hopf–Cole transformation, the Padé approximation and Richardson’s extrapolation has been developed in this paper. The nonlinear Burgers’ equation was transformed into a homogeneous heat equation using the Hopf–Cole transformation. A compact fourth-order finite difference scheme was then developed to solve the heat equation. To maintain overall fourth-order accuracy, some novel higher-order compact approximation schemes for boundary and initial conditions are developed as well. Finally, based on the Hopf–Cole transformation, a compact fourth-order finite difference scheme is developed to reconstruct the solution to u(x, t) from the numerical solution of the heat equation. Numerical experiments have shown that this new scheme is robust, unconditionally stable based on the strict von Neumann condition for stability [21], and very accurate (fourth-order in both temporal and spatial dimensions). Acknowledgement This work is supported by the Natural Sciences and Engineering Research Council of Canada (NSERC). The author would like to thank the anonymous referees for their time, effort, and extensive comments on the revision of the manuscript. References [1] [2] [3] [4] [5] [6] [7] [8] [9] [10]

E.N. Aksan, A. Ozdes, A numerical solution of Burgers’s equation, Appl. Math. Comput. 156 (2004) 395–402. H. Bateman, Some recent researches on the motion of fluids, Monthly Weather Rev. 43 (1915) 163–170. E. Benton, G.W. Platzman, A table of solutions of the one-dimensional Burgers’ equations, Quart. Appl. Math. 30 (1972) 195–212. S. Bringen, A. Saati, Comparison of several finite-difference methods, J. Aircraft 27 (1990) 90–99. J.M. Burger, A Mathematical Model Illustrating the Theory of Turbulence, Adv. in Appl. Mech., vol. I, Academic Press, New York, 1948. pp. 171–199. J. Caldwell, P. Smith, Solution of Burgers’ equation with a large Reynolds number, Appl. Math. Model. 6 (1982) 381–385. J.D. Cole, On a quasilinear parabolic equations occurring in aerodynamics, Quart. Appl. Math. 9 (1951) 225–236. D.J. Evans, A.R. Abdullah, The group explicit method for the solution of Burgers’ equation, Quart. Appl. Math. 30 (1984) 239–253. B. Gustafson, H. Kress, J. Oliger, Time Dependent Problems and Difference Methods, John Wiley & Sons, New York, 1995. I.A. Hassanien, A.A. Salama, H.A. Hosham, Fourth-order finite difference method for solving Burgers’ equation, Appl. Math. Comput. 170 (2005) 781– 800. [11] Y.C. Hon, X.Z. Mao, An efficient numerical scheme for Burgers’ equation, Appl. Math. Comput. 95 (1998) 37–50. [12] E. Hopf, The partial differential equation ut + uux = luxx, Commun. Pure Appl. Math. 3 (1950) 201–230. [13] S. Kutluay, A. Rsen, A linearized numerical scheme for Burgers-like equations, Appl. Math. Comput. 156 (2004) 295–305.

764

W. Liao / Applied Mathematics and Computation 206 (2008) 755–764

[14] S. Kutluay, A.R. Bahadir, A. Ozdes, Numerical solution of one-dimensional Burgers equation: explicit and exact-explicit finite difference methods, J. Comput. Appl. Math. 103 (1999) 251–261. [15] R. Leveque, Numerical Methods for Conservation Laws, Birkhauser, 1992. [16] W. Liao, J. Zhu, A.Q.M. Khaliq, A fourth-order compact algorithm for nonlinear reaction–diffusion equations with Neumann boundary conditions, Numer. Methods Partial Differen. Equat. 22 (2006) 600–616. [17] R.C. Mitta, P. Signnal, Numerical solution of Burgers’ equation, Commun. Numer. Methods Eng. 9 (1993) 397–406. [18] W. Nguyen, J. Reynen, A space–time finite element approach to Burgers’ equation, in: C. Taylor et al. (Eds.), Proceedings of the International Conference on Numerical Methods for Non-linear Problems, Universal Poltecnica de Barcelona, Spain, vol. 2, Pineridge Press, Swansea, UK, 1984. [19] T. Ozis, Y. Aslan, The semi-approximate approach for solving Burgers’ equation with high Reynolds number, Appl. Math. Comput. 163 (2005) 131–145. [20] T. Ozis, A. Ozdes, A direct variational methods applied to Burgers’ equation, J. Comput. Appl. Math. 71 (1996) 163–175. [21] G.D. Smith, Numerical solution of partial differential equation: finite difference Methods, third ed., Clarendon Press, Oxford, 1987. [22] E. Varoglu, W.D. Finn, Space–time finite elements incorporating characteristics for the Burgers’ equation, Int. J. Numer. Methods Eng. 16 (1980) 171– 184. [23] W.L. Wood, An exact solution for Burgers’ equation, Commun. Numer. Methods Eng. 22 (7) (2006) 797–798.