Numerical solutions of the modified Burgers’ equation by Petrov–Galerkin method

Applied Mathematics and Computation 218 (2011) 3673–3679

Contents lists available at SciVerse ScienceDirect

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

Numerical solutions of the modified Burgers’ equation by Petrov–Galerkin method Thoudam Roshan ⇑, K.S. Bhamra Department of Mathematics, Manipur University, Canchipur 795 003, Imphal, India

a r t i c l e

i n f o

a b s t r a c t The modified Burgers’ equation (MBE) is solved numerically by the Petrov–Galerkin method using a linear hat function as the trial function and a cubic B-spline function as the test function. Product approximation has been used in this method. A linear stability analysis of the scheme shows it to be unconditionally stable. The accuracy of the presented method is demonstrated by two test problems. The numerical results are found in good agreement with the exact solutions. Ó 2011 Elsevier Inc. All rights reserved.

Keywords: Cubic B-spline Petrov–Galerkin Product approximation Wave packet

1. Introduction The Burgers’ equation

ut þ uux  muxx ¼ 0;

ð1Þ

where m is a positive constant and the subscripts x and t denote space and time derivatives respectively was first introduced by Batman [1] and later treated by Burgers’ [2] as a mathematical model for turbulence. Since then the equation has found applications in fields as diverse as number theory, gas dynamics, heat conduction, elasticity, etc. The Burgers’ equation was solved numerically by various methods such as the finite difference [3], Galerkin [4,5], least squares [6] and collocation methods [7,8] etc. Indeed, the Burgers’ equation is a special case of the modified Burgers’ equation (MBE) of the form

ut þ ul ux  muxx ¼ 0;

ð2Þ

where l is a positive constant and m can be interpreted as viscosity. The MBE has the strong nonlinear aspects of the governing equation in many practical transport problems such as nonlinear waves in medium with low frequency pumping or absorption, ion reflection at quasi perpendicular shocks, turbulence transport, wave processes in thermoelastic medium, transport and dispersion of pollutants in river and sediment transport etc. The initial condition associated with Eq. (2) will be

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

a 6 x 6 b;

ð3Þ

with the boundary conditions

uða; tÞ ¼ g 1 ðtÞ and uðb; tÞ ¼ g 2 ðtÞ;

t > t0 :

ð4Þ

Various numerical methods have been proposed to solve MBE. Ramadan et al. [8] used septic B-spline collocation method for the numerical solution of the MBE with l = 1 and l = 2. Ramadan and Danaf [9] also obtained the numerical solution of MBE using quintic B-spline collocation method for the case l = 2 only. Griewank et al. [10] used non-polynomial spline-functions ⇑ Corresponding author. E-mail address: [email protected] (T. Roshan). 0096-3003/$ - see front matter Ó 2011 Elsevier Inc. All rights reserved. doi:10.1016/j.amc.2011.09.010

3674

T. Roshan, K.S. Bhamra / Applied Mathematics and Computation 218 (2011) 3673–3679

for the numerical treatment of the MBE for the cases l = 1 and l = 2. Sachdev et al. [11] obtained large-time asymptotics for periodic solutions of the MBE for the cases l = 2 and l = 3. Lattice Boltzmann model for the MBE have been studied by Duan et al. [12] for various initial conditions. In all the numerical techniques mentioned above, the small value of the viscosity m considered was up to 0.001 only. In this work, the Petrov–Galerkin method is developed for the MBE using the linear hat function as the trial function and the cubic B-spline function as the test function. Here the proposed method is shown to represent accurately the various wave packets for the cases l = 2 and 3 and it can treat small value of m up to 0.0001. The numerical results obtained by this proposed method are compared with the known solutions. 2. The Petrov–Galerkin method For convenience, the MBE (2) is rewritten as

ut þ

1

lþ1

ðulþ1 Þx  muxx ¼ 0:

ð5Þ

The space interval a 6 x 6 b is discretized with (N + 1) uniform grid points xj = a + jh, where j = 0, 1, 2, . . . , N, and the grid spacing is given by h = (b  a)/N. Let Uj(t) denote the approximation to the exact solution u(xj, t). Following the Petrov–Galerkin method used in [13], we assume that the approximate solution of Eq. (3) as

uh ðx; tÞ ¼

N X

U j ðtÞ/j ðxÞ:

ð6Þ

j¼0

The product approximation technique [14] is used for treating the nonlinear term in the following manner: lþ1

uh

N X lþ1 U j ðtÞ/j ðxÞ;

ðx; tÞ ¼

ð7Þ

j¼0

where

8 > < 1 þ ðx  jhÞ=h; x 2 ½xj1 ; xj  /j ðxÞ ¼ 1  ðx  jhÞ=h; x 2 ½xj ; xjþ1  > : 0; otherwise The unknown functions Uj(t), j = 0, 1, 2, . . . , N, are determined from the variational formulation

ððuh Þt ; wj Þ þ

 1  lþ1  ; wj  mððuh Þxx ; wj Þ ¼ 0; uh x lþ1

ð8Þ

where wj, j = 0, 1, 2, . . . , N, are test functions, which are taken to be the cubic B-splines given by

wj ðxÞ ¼

1

8 > ðx  xj2 Þ3 ; x 2 ½xj2 ; xj1  > > > > 3 2 2 3 > > < h þ 3h ðx  xj1 Þ þ 3hðx  xj1 Þ  3ðx  xj1 Þ ; x 2 ½xj1 ; xj  3

2

h þ 3h ðxjþ1  xÞ þ 3hðxjþ1  xÞ3  3ðxjþ1  xÞ3 ; x 2 ½xj ; xjþ1  3 h > > > > > x 2 ½xjþ1 ; xjþ2  ðxjþ2  xÞ3 ; > > : 0; otherwise

and (, ) denotes the usual inner product:

ðf ; gÞ ¼

Z

b

f ðxÞgðxÞdx:

a

Integrating by parts and using the fact that w(a) = w(b) = w0 (a) = w0 (b) = 0, Eq. (8) leads to the formulation

ððuh Þt ; wj Þ þ

1

lþ1

   lþ1 ; wj  mðuh ; ðwj Þxx Þ ¼ 0: uh x

ð9Þ

Performing the integrations on (9) will give the following system of ordinary differential equations (ODEs):

1  1 m Aþ B  2 C ¼ 0; 20 4ðm þ 1Þh h where

A ¼ U j2 þ 26U j1 þ 66U j þ 26U jþ1 þ U jþ2 ; B ¼ ðU j2 Þlþ1  10ðU j1 Þlþ1 þ 10ðU jþ1 Þlþ1 þ ðU jþ2 Þlþ1 C ¼ U j2 þ 2U j1  6U j þ 2U jþ1 þ U jþ1

ð10Þ

T. Roshan, K.S. Bhamra / Applied Mathematics and Computation 218 (2011) 3673–3679

3675

and ‘’ means the derivative with respect to time t. Now to solve the ODEs, we assume U nj to be a fully discrete approximation to the exact solution u(xj, tn), where tn = nDt and Dt is the time step size. Replacing the time derivative of the parameter U by usual difference approximation,  U ¼ ðU nþ1  U n Þ=Dt and parameter U by the Crank–Nicholson formulation, U = (Un+1 + Un)/2 the system of ODEs (10) reduces to a system of nonlinear equations

1 nþ1 Dt mDt 1 n Dt mDt þ A A  Bnþ1  2 Cnþ1 ¼ B n þ 2 Cn ; 10 4ðl þ 1Þh 10 4ðl þ 1Þh h h

ð11Þ

where j = 2, 3, . . . , N  2 and n = 1, 2, 3, . . . . The nonlinear system of Eq. (11) can be solved by Newton’s method and the desire solution of the MBE can be obtained. 3. Stability analysis To apply the Von Neumann stability for the system (11), we must first linearize this system. Assuming u in the nonlinear ~ , the linearized scheme is term ulux of the MBE as locally constant, u nþ1

nþ1

0

0

nþ1 n n 0 n 0 n 0 n aU nþ1 þ dU jþ1 þ eU nþ1 j2 þ bU j1 þ cU j jþ2 ¼ a U j2 þ b U j1 þ c U j þ d U jþ1 þ e U jþ2 ;

ð12Þ

where

~ l Dt mDt ~ l Dt mDt 1 u 1 u þ  2 ; a0 ¼  þ 2 ; 10 4h 10 4h h h l ~ Dt 2mDt ~ l Dt 2 m Dt 26 10u 26 10u 0 b¼ þ  2 ; b ¼  þ 2 ; 10 4h 10 4h h h 66 6mDt 66 6mDt 0 c¼ þ 2 ; c ¼  2 ; 10 10 h h ~ l Dt 2mDt ~ l Dt 2mDt 26 10u 26 10u 0 d¼   2 ; d ¼ þ þ 2 ; 10 4h 10 4h h h l l ~ Dt mDt ~ Dt m Dt 1 u 1 u 0   2 þ þ 2 : ;e ¼ e¼ 10 4h 10 4h h h

a¼

Substituting the Fourier mode

U nj ¼ nn expðijhÞ;

h ¼ hk and i ¼

pffiffiffiffiffiffiffi 1;

ð13Þ

where k is the mode number and h the space step size into the linearized system, we obtained the amplification factor, g as

g¼

nnþ1 X þ iY ¼ ; X 1  iY nn

ð14Þ

where

    1 m Dt 26 2mDt 66 6mDt cos h þ þ 2 cos 2h þ 2 þ 2  2 ; X¼2 10 10 10 h  h  h   1 mDt 26 2mDt 66 6mDt cos h þ  2 cos 2h þ 2  2 þ 2 X1 ¼ 2 10 10 10 h h h and

Y¼

~ l Dt u ðsin 2h þ 10 sin hÞ: 2h

As h, Dt and m are all positive, therefore

X1  X ¼

4mDt h

2

ð3  cos 2h  2 cos hÞ P 0 i:e: X 1 P X:

So,

pffiffiffiffiffiffi jgj ¼ g g ¼

sffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffi X2 þ Y 2 X 21 þ Y 2

6 1:

Hence, the proposed method is unconditionally stable in the linear sense.

3676

T. Roshan, K.S. Bhamra / Applied Mathematics and Computation 218 (2011) 3673–3679 Table 1 The L2 and L1 error norm with l = 2, m = 0.01, h = 0.005 and Dt = 0.01. Time

L2  103

L1  103

Time

L2  103

L1  103

2 3 4 5 6

0.375515 0.344421 0.316752 0.294779 0.27492

0.81766 0.712453 0.608055 0.52857 0.467517

7 8 9 10 11

0.255093 0.234735 0.214154 0.193914 0.174518

0.418355 0.376199 0.338116 0.302855 0.270127

Table 2 The L2 and L1 error norm with l = 2, m = 0.005, h = 0.005 and Dt = 0.01. Time

L2  103

L1  103

Time

L2  103

L1  103

2 3 4 5 6

0.223323 0.205476 0.189344 0.176734 0.166348

0.580808 0.506182 0.432122 0.375678 0.332598

7 8 9 10 11

0.157487 0.149747 0.142853 0.13659 0.130782

0.298876 0.271735 0.249436 0.230739 0.214837

4. Numerical results We now obtain the numerical solutions of the MBE for two standard problems. In the first example, the case l = 2 for the MBE is consider and in this case the numerical solution are compare with the analytic solution and available numerical results obtained in[8,9]. In the second example l = 3 is consider and we show that our results agree with the asymptotic solution obtained in [11]. Our numerical results are also compared with that of the results obtained in [12]. In both the cases we used L2 and L1-error norms for testing the accuracy of the scheme. 4.1. Example 1 The analytic solution of the MBE, for l = 2, which was obtained in [15] take the form:

uðx; tÞ ¼

1þ

ðx=tÞ pffiffi  ; t =c0 expðx2 =4mtÞ

ð15Þ

where c0 is constant and 0 < c0 < 1. In this case we take c0 = 0.5. The initial condition u(x,1) is given from equation (15), while the boundary conditions (4) are g1(t) = 0, g2(t) = u(1, t) with u(1, t) given by (15). For comparison with the relevant known results in [8,9] the values m = 0.01, 0.005 and 0.001 are used. The obtained results are summarized in Tables 1–3. From these tables, we observed that as the value of m decreases the values of the error norms are also decrease. The results obtained in Tables 1–3 are compare with the results in [8,9] and it is tabulated in Table 4. Form this table we observed that our method is superior than the methods found in [8,9]. Figs. 1–3 illustrate the behaviors of the numerical solutions for m = 0.01, 0.005 and 0.001, where h = 0.005,Dt = 0.01 at some different times t = 1, 4, 7 and 10 respectively. The top curve is at t = 1 and the bottom curve is at t = 10. The Figures show that as the time increases the curve of the numerical solution decays. Also note that the decay gets fast as the viscosity m decreases. We also consider the case, where m = 0.0001. This case was not considered in [8,9]. Here, we take h = 0.005, Dt = 0.01. The simulations are done up to t = 11. Table 5 represents the error norms for the given method for this case m = 0.0001. 4.2. Example 2 In this example, we solved the MBE, forl = 3. Following [11], the initial and boundary conditions are used as

Table 3 The L2 and L1 error norm with l = 2, m = 0.001, h = 0.005 and Dt = 0.01. Time

L2  103

L1  103

Time

L2  103

L1  103

2 3 4 5 6

0.066071 0.061793 0.0573984 0.053741 0.0506308

0.261856 0.228837 0.195767 0.170233 0.150895

7 8 9 10 11

0.0479362 0.0455715 0.0434759 0.0416037 0.0399192

0.135618 0.123286 0.113173 0.104701 0.097458

3677

T. Roshan, K.S. Bhamra / Applied Mathematics and Computation 218 (2011) 3673–3679 Table 4 Comparisons of results with l = 2, h = 0.005 and D t = 0.01. L2  103 t=2

L1  103 t=2

L2  103 t=6

L1  103 t=6

L2  103 t = 10

L1  103 t = 10

m = 0.01 Ours [9] [8]

0.37552 0.52308 0.79043

0.81766 1.21698 1.70309

0.27492 0.49023 0.57672

0.46752 0.72249 0.76105

0.19391 0.64007 0.80026

0.23074 1.28124 1.80329

m = 0.005 Ours [9]

0.22332 0.25786

0.58081 0.72264

0.16635 0.22569

0.33260 0.43082

0.13659 0.18735

0.23074 0.30006

m = 0.001 Ours [9] [8]

0.06607 0.06703 0.18355

0.26186 0.27967 0.81862

0.05063 0.06046 0.08142

0.15090 0.17176 0.21348

0.04160 0.05010 0.05512

0.10470 0.12129 0.13943

m-values & methods

Fig. 1. Solution behavior of the MBE with l = 2, m = 0.01, h = 0.005, Dt = 0.01.

Fig. 2. Solution behavior of the MBE with l = 2, m = 0.005, h = 0.005, Dt = 0.01.

uðx; 0Þ ¼ A sin

px ; d

uð0; tÞ ¼ uðd; tÞ ¼ 0

in which d = p and A = 1, then the MBE has an asymptotic solution[11] of the form

uðx; tÞ ¼ ekt f0 ðx; tÞ þ e4kt f1 ðx; tÞ þ e7kt f2 ðx; tÞ þ    ;

ð15Þ

3678

T. Roshan, K.S. Bhamra / Applied Mathematics and Computation 218 (2011) 3673–3679

Fig. 3. Solution behavior of the MBE with l = 2, m = 0.001, h = 0.005, Dt = 0.01.

Table 5 The L2 and L1 error norm with l = 2, m = 0.0001, h = 0.005 and Dt = 0.01 Time

L2  104

L1  104

Time

L2  104

L1  104

2 3 4 5 6

0.112962 0.112783 0.10709 0.100867 0.095058

0.818695 0.743028 0.640708 0.557252 0.498194

7 8 9 10 11

0.089841 0.0852 0.0810706 0.0773826 0.0740724

0.448047 0.405932 0.373655 0.344082 0.321101

Fig. 4. Solution behavior of the MBE with l = 3, m = 0.005, h = 0.005, Dt = 0.01.

Table 6 Comparisons of results with l = 3, m = 0.005, h = 0.005 and Dt = 0.01. Time, t

L2  103(Ours)

L2  103 [12]

L1  103(Ours)

L1  103 [12]

150 200 250 300 350 400 450

6.128090 2.229470 0.914486 0.416609 0.234546 0.165391 0.131525

3.227 0.9912 0.5031 0.5939 0.6940 0.7567 0.7990

6.845380 2.045900 0.839436 0.399913 0.221936 0.144269 0.105834

5.172 1.671 1.400 1.452 1.488 1.513 1.531

T. Roshan, K.S. Bhamra / Applied Mathematics and Computation 218 (2011) 3673–3679

3679

where

    px A4 p 2px A4 d 4px ; f 1 ðx; tÞ   1 sin þ 1 sin ; d d 96mp d 4d       px 3px 5p x 7p x þ g 4 ðtÞ sin þ g 5 ðtÞ sin þ g 6 ðtÞ sin ; f2 ðx; tÞ ¼ g 3 ðtÞ sin d d d d " # " # 2 2 2 2 d d D1 d d D2 g 3 ðtÞ ¼  D1 t þ E1 þ D2 t þ E2  ; g 4 ðtÞ ¼ 2 2 6mp 6m p 2mp 2mp2 " # 2 2 2 d d D3 d E4 mp2 ; g D t þ E  ðtÞ ¼ ; k ¼ ; g 5 ðtÞ ¼ 3 3 6 2 18mp2 18mp2 42mp2 d

f0 ðx; tÞ ¼ A1 sin

A41 p A4 d A3 B1 p 9A3 B1 p 5A3 B1 p ; B2 ¼ 1 ; D1 ¼ 1 ; D2 ¼  1 ; D3 ¼ 1 ; 96mp 4d 4d 8d 8d 3 3 3 3 A B2 p 9A B2 p 15A1 B2 p 7A B2 p ; E2 ¼ 1 ; E3 ¼  ; E4 ¼ 1 E1 ¼  1 8d 8d 8d 8d

B1 ¼ 

and A1 is an arbitrary constant, called the old age constant; this represents the memory of the initial conditions [11] and it converge to 0.365366. For the numerical test we choose m = 0.005, h = 0.005 and Dt = 0.01. We conduct a comparison between our numerical scheme and the asymptotic as well as the numerical solutions found in [11,12]. Fig. 4 shows the numerical results and the asymptotic solutions at different time stages. Here the doted line represents the asymptotic solution and the solid line represents our numerical solution. At time t = 150, there is little difference between the asymptotic and numerical solution. From the graph we observed that our numerical result is in good agreement with the asymptotic solution as time increases. L2 and L1-error norms are used to measured the computational efficiency and accuracy of the numerical scheme. Our numerical results and results obtained in [12] at some different time levels are tabulated in Table 6. From the table it is clear that as time increases, the error norms decrease. This shows that our numerical scheme is quite stable. 5. Conclusion The Petrov–Galerkin method using the linear hat function and cubic B-spline function as trial and test functions respectively has been successfully implemented to study the numerical solutions of the MBE. It has been shown that our scheme is unconditionally stable. According to the two test problems, the calculations with the analytic/asymptotic solution shows the propose method is capable of solving MBE accurately. The further comparison of the presented numerical results with methods mentioned in this paper shows that our method is better than the methods cited in[8,9,12]. Lastly, we conclude that our scheme is capable of solving MBE for small values of viscosity. Acknowledgements The authors thank the referee for giving valuable comments and suggestions. References [1] H. Batman, Some recent researches on the motion of fluids, Monthly Weather Rev. 43 (1915) 163–170. [2] J.M. Burger, A Mathematical Model Illustrating the Theory of Turbulence, in: Adv. Appl. Mech., Academic Press, New York, 1948, pp. 177–199. [3] S. Kutluay, A.R. Bahader, A. Ozdes, Numerical solution of one dimensional Burgers’ equation explicit and exact-explicit finite difference method, J. Comput. Appl. Math. 109 (1948) 251–261. [4] A. Dogan, A Galerkin finite element approach to Burgers’ equation, Appl. Math. Comput. 157 (2004) 331–346. [5] T. Özis, E.N. Aksan, A. Özdes, A finite element approach for solution of Burgers’ equation, Appl. Math. Comput. 139 (2003) 417–428. _ Dag˘, Numerical solution of Burgers’ by the least-squares quadratic B-spline finite element method, J. Comput. Appl. Math. 167 [6] S. Kutluay, A. Aksan, I. (2004) 21–33. _ [7] I. Dag˘, D. Irk, B. Saka, Numerical solution of the Burgers’ equation using cubic B-splines, Appl. Math. Comput. 163 (2005) 199–211. [8] M.A. Ramadan, T.S. El-Danaf, Abd. Alalel El, A numerical solution of the Burgers’ equation using septic B-splines, Chaos Solitons & Fractals 26 (2005) 1249–1258. [9] M.A. Ramadan, T.S. El Danaf, Numerical treatment for the modified Burgers’ equation, Math. Comput. Simul. 70 (2005) 90–98. [10] A. Griewank, S. Talaat, El-Danaf, Efficient accurate numerical treatment of the modified Burgers’ equation, Appl. Anal. 88 (1) (2009) 75–87. [11] P.L. Sachdev, Ch. Srinivasa Rao, B.O. Enflo, Large-time asymptotics for periodic solutions of the modified Burgers’ equation, Stud. Appl. Math. 114 (2005) 307–323. [12] Y. Duan, R. Liu, Y. Jiang, Lattice Boltzmann model for the modified Burgers’ equation, Appl. Math. Comput. 202 (2008) 489–497. [13] J.M. Sanz-Serna, I. Christie, Petrov–Galerkin methods for nonlinear dispersive waves, J. Comput. Phys. 39 (1981) 94–103. [14] I. Christie, D.F. Griffths, A.R. Mitchell, J.M. Sanz-Serna, Product approximation for nonlinear problems in finite element methods, IMA J. Numer. Anal. 1 (1981) 253–266. [15] S.L. Harris, Sonic shocks governed by the modified Burgers’ equation, EJAM 6 (1996) 75–107.