Fourth-order compact schemes for the numerical simulation of coupled Burgers’ equation

Accepted Manuscript Fourth-order compact schemes for the numerical simulation of coupled Burgers’ equation H.P. Bhatt, A.Q.M. Khaliq PII: DOI: Reference:

S0010-4655(15)00418-X http://dx.doi.org/10.1016/j.cpc.2015.11.007 COMPHY 5808

To appear in:

Computer Physics Communications

Received date: 4 May 2015 Revised date: 11 November 2015 Accepted date: 16 November 2015 Please cite this article as: H.P. Bhatt, A.Q.M. Khaliq, Fourth-order compact schemes for the numerical simulation of coupled Burgers’ equation, Computer Physics Communications (2015), http://dx.doi.org/10.1016/j.cpc.2015.11.007 This is a PDF file of an unedited manuscript that has been accepted for publication. As a service to our customers we are providing this early version of the manuscript. The manuscript will undergo copyediting, typesetting, and review of the resulting proof before it is published in its final form. Please note that during the production process errors may be discovered which could affect the content, and all legal disclaimers that apply to the journal pertain.

*Manuscript

1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60 61 62 63 64 65

Fourth-order compact schemes for the numerical simulation of coupled Burgers’ equationI H. P. Bhatt∗, A. Q. M. Khaliq Department of Mathematical Sciences and Center for Computational Science, Middle Tennessee State University, Murfreesboro, TN 37132, USA

Abstract This paper introduces two new modified fourth-order exponential time differencing Runge-Kutta (ETDRK) schemes in combination with a global fourth-order compact finite difference scheme (in space) for direct integration of nonlinear coupled viscous Burgers’ equations in their original form with out using any transformations or linearization techniques. One scheme is a modification of the Cox and Matthews ETDRK4 scheme based on (1, 3)−Pad´e approximation and other is a modification of Krogstad’s ETDRK4-B scheme based on (2, 2)−Pad´e approximation. Efficient versions of the proposed schemes are obtained by using a partial fraction splitting technique of rational functions. The stability properties of the proposed schemes is studied by plotting the stability regions, which provide an explanation of their behavior for dispersive and dissipative problems. The order of convergence of the schemes is examined empirically and found that the modification of ETDRK4 converges with the expected rate even if the initial data are nonsmooth. On the other hand, modification of ETDRK4-B suffers with order reduction if the initial data are nonsmooth. Several numerical experiments are carried out in order to demonstrate the performance and adaptability of the proposed schemes. The numerical results indicate that the proposed schemes provide better accuracy than other schemes available in the literature. Moreover, the results show that the modification of ETDRK4 is reliable and yields more accurate results than modification of ETDRK4-B, while solving problems with nonsmooth data or with high Reynolds number. Keywords: Compact scheme; Exponential time differencing scheme; Coupled viscous Burgers’ equation; Partial fraction splitting technique

Pad´e approximation;

1. Introduction Burgers’ equation is one of the most common nonlinear time dependent partial differential equation (PDE) in fluid mechanics, which consists both of nonlinear propagation and diffusive effects. The equation has been widely used as a simplified model for several physical phenomena, such as Work done while the second author was on Sabbatical leave at the Department of Mathematics and Statistics, King Fahd University of Petroleum and Minerals, Saudi Arabia. ∗ Corresponding author. E-mail address: [email protected](H. P. Bhatt), [email protected] (A. Q. M. Khaliq). I

Preprint submitted to Computer Physics Communications

November 11, 2015

1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60 61 62 63 64 65

modeling of traffic flow and gas dynamics, shock waves, shallow water waves, describing wave propagation in acoustics, hydrodynamics, and most commonly to test various numerical schemes. In addition, it is one of the few nonlinear PDEs that can be solved analytically for a restricted set of arbitrary initial conditions. Hopf [1] and Cole [2] showed independently that for any given initial conditions the Burgers’ equation can be reduced to a linear homogeneous heat equation that can be solved analytically and the analytical solution of original Burgers’ equation can be expressed in the form of Fourier series. Even though the exact solution is available in the form of the Fourier series, accurate and efficient numerical schemes are still required to solve the Burgers’ equation which consists of non smooth data or the initial condition is only available on discrete locations. In such situations, the Fourier series solutions may converge very slowly [3]. Numerical solution of Burgers’ equation has been carried out with various numerical techniques during the past several decades due to equation’s wide range of applicability in various field of science and engineering and is still an active field of research. To name a few, in [4] the authors constructed an explicit and exact-explicit finite difference schemes based on the Hopf-Cole transformation. Fourthorder finite difference methods have been proposed in [3, 5]. A hybrid numerical scheme based on Euler implicit method, quasilinearization and uniform Haar wavelets has been introduced in [6]. The coupled viscous Burgers’ equation was derived by Esipov [7] and is a simple model of sedimentation or evolution of scaled volume concentrations of two kinds of particles in fluid suspensions or colloids under the effect of gravity [7]. In recent years, numerical solution to coupled Burgers’ equation has attracted a lot of attention from both scientists and engineers and this has resulted in a number of numerical algorithms. For instance, cubic B-spline collocation scheme [8], Chebyshev spectral collocation method [9], fully implicit finite difference scheme [10], Fourier Pseudospectral method [11], differential quadrature method [12], compact operator method [13], mesh free interpolation method [14], the variational iteration method [15], composite numerical scheme [16], ChebyshevLegendre Pseudo-Spectral method [17], Haar wavelet-based scheme [18], finite element method [19], and the general distributed approximating functional (DAF) method [20] have been developed for the numerical solution of coupled Burgers’ equations. To the best knowledge of authors, no numerical scheme of order four in time and four in space for the coupled Burgers’ equation given in (4) has been developed so far. Hence, the main goal of this paper is to develop A−stable and L−stable Fourth-order ETD schemes in combination with a global fourth-order compact finite difference scheme for the numerical solution of the following Burgers’ and coupled Burgers’ equation in their original form with out using any transformations or linearization techniques: • Burgers’ Equation:

∂u ∂ 2u ∂u + ηu = δ 2 , x ∈ Ω = [a, b], t ∈ [0, T ], ∂t ∂x ∂x with the initial condition: u(x, 0) = g1 (x),

x ∈ Ω = [a, b],

(1)

(2)

and the boundary conditions: u(a, t) = f1 (t), u(b, t) = f2 (t), 2

t ∈ [0, T ],

(3)

1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60 61 62 63 64 65

• Coupled Burgers’ Equation [7]: ∂u ∂ 2u ∂(uv) ∂u − δ 2 + ηu +α = 0, x ∈ Ω = [a, b], t ∈ [0, T ], ∂t ∂x ∂x ∂x ∂v ∂ 2v ∂(uv) ∂v − µ 2 + ξv +β = 0, x ∈ Ω = [a, b], t ∈ [0, T ], ∂t ∂x ∂x ∂x

(4)

subject to the initial conditions: u(x, 0) = g1 (x), v(x, 0) = g2 (x),

x ∈ Ω = [a, b], x ∈ Ω = [a, b],

(5)

and boundary conditions: u(a, t) = f1 (t), u(b, t) = f2 (t), v(a, t) = f3 (t), v(b, t) = f4 (t),

t ∈ [0, T ], t ∈ [0, T ],

(6)

where δ, µ are positive viscosity parameters which correspond to an inverse of Reynolds number Re ( if δ = µ, then δ = µ = 1/Re ), η, ξ are real constants and α, β are arbitrary constants depending upon the system parameters such as the P´eclet number, f1 (t), f2 (t), f3 (x), f4 (x), g1 (x), g2 (x) are given functions. Several numerical experiments on Burgers’ and coupled Burgers’ equations were run in order to compare the accuracy of the proposed schemes with other existing methods [8, 9, 11–13, 19, 20]. The numerical results show that the proposed schemes provide better solutions in most of the cases than the methods discussed below. For clarity a brief description of compared methods are provided below. Authors in [8] presented a method based on Crank-Nicolson and cubic B-spline functions scheme for the numerical solution of the coupled Burgers’ equation. In [9] the authors developed a Chebyshev spectral collocation method for the approximate solution of nonlinear PDEs. Utilizing the Chebyshev spectral collocation method, the problem is reduced to a system of ODEs and is then solved by Runge-Kutta method of order four. Rashid et al. [11] proposed an approximate solution of coupled Burgers’ equation with a set of initial and periodic boundary conditions through Fourier pseudo-spectral method in space and classical fourth-order Runge-Kutta method in time. In [12] the authors analyzed a coupled Burgers’ equation with a polynomial differential quadrature method (PDQM). The PDQM transformed the coupled Burgers’ equation in to system of nonlinear ODEs and is then solved by fourth-order Runge-Kutta method. Mohanty et al. [13] proposed a new two-level implicit compact operator method of order two in time and four in space for the solution of coupled Burgers’ equations. In their method Mohanty et al. did not use any transformation or linearization techniques to handle nonlinearity. In [19], Arminjon and Beauchamp introduced a finite element method for the 2D coupled Burgers’ equation and concluded that the method is efficient compared to the other methods, namely the method of lines and Runge-Kutta-type method. Wei et al. [20] proposed distributed approximating functional (DAF) scheme for spatial discritization and Taylor series expansion for time discritization to solve one and two-dimensional Burgers’ 3

1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60 61 62 63 64 65

equations. The remainder of the paper is organized as follows. In section 2 compact fourth-order schemes are described to approximate uxx and ux . In addition, in the same section, two modified fourth-order ETD schemes are derived. The stability of the proposed schemes are discussed in section 3. In section 4, numerical experiments are performed to test the accuracy and adaptability of the proposed schemes. The conclusions are presented in section 5, which briefly summarizes the numerical outcomes. 2. Fourth-order numerical Schemes This paper concerned with a uniform grid, thus we descritize the solution domain Ω × [0, T ] = {(x, t)/ a ≤ x ≤ b, 0 ≤ t ≤ T } into grids described by the set of nodes {(xi , tj )}, in which xi = a + ih, i = 0, 1, · · · , N and tj = jk, j = 0, 1, · · · , M = Tk , where h and k are spatial and temporal step sizes, respectively. 2.1. Compact finite difference formulae There are many methods used to generate compact finite difference formulae to approximate first and second-order spatial derivatives. Readers are referred to [21, 22] for more details on how to generate compact finite difference formulae. In this study, the spatial derivatives in Burgers’ equation (1) are approximated with the formulae introduced by Zhao in [23]. Below the formulae in [23] are stated without change. The standard compact finite difference formula for first derivatives of u(x, t) at interior points is: u0i−1 + 4u0i + u0i+1 = where ui = u(xi ) and u0i = i = 1 use: 4u01

+

u02

1 = h

du(xi ) . dx



3 (ui+1 − ui−1 ), i = 2, 3, · · · , N − 2, h

(7)

The truncation error for formula (7) is O(h4 ). At boundary, when

 11 4 1 − u0 − 4u1 + 6u2 − u3 + u4 , 12 3 4

and when i = N − 1, the formula is:   1 4 11 1 0 0 uN −2 + 4uN −1 = − uN −4 + uN −3 − 6uN −2 + 4uN −1 + uN . h 4 3 12

(8)

(9)

The truncation errors in the formula (8) and (9) are O(h4 ). Writing (7)-(9) into matrix form as: A1 U 0 = M1 U + H1 ,

(10)

4

1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60 61 62 63 64 65

where  4 1  1 4  . A1 =  0 . . .  .. 0 ···

 ··· 0 ..  1 . 0 ... ...  , 0   1 4 1 0 1 4 (N −1)×(N −1) 0

  −u0  0   11   ..  , H1 =   12h  .   0  uN (N −1)×1

 4  1 − 3 2 − 49 12 0 ··· 0  ..   −1 0 1 0 0 .   .  . . ...  0 −1 0 1 0   3  . , M1 =  0 . . 0 0 −1 0 1  h  ..  .. . . .. .. ..  .  . . . . .    0 0 −1 0 1  4 1 −2 43 (N −1)×(N −1) 0 0 0 − 12 9 

u01 0 .. .



      0 U = ,     0  u0N −1 (N −1)×1



u1 0 .. .



      U = .     0  uN −1 (N −1)×1

Hence the fourth-order compact finite-difference approximation of ux is given by: U 0 = A−1 1 (M1 U + H1 ).

(11)

The standard compact finite difference formula for second derivatives of u(x, t) at interior points is: u00i−1 + 10u00i + u00i+1 =

12 (ui−1 − 2ui + ui+1 ), i = 2, 3, · · · , N − 2, h2

(12)

2

u(xi ) where u00i = d dx . The truncation error for formula (12) is O(h4 ). 2 At boundary, when i = 1, apply:

14u001 − 5u002 + 4u003 − u004 =

12 (u0 − 2u1 + u2 ), h2

(13)

and when i = N − 1, − u00N −4 + 4u00N −3 − 5u00N −2 + 14u00N −1 =

12 (uN −2 − 2uN −1 + uN ). h2

(14)

The truncation error in both of the formulae is also O(h4 ). Writing (12)-(14) in matrix form yields: A2 U 00 = M2 U + H2 ,

(15)

5

1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60 61 62 63 64 65

where   −14 −5 4 −1 0 · · · 0    ..  −2 1 0 ··· 0  1  10 1 0 0 .   ..   ...    1 −2 1 . 0 1 10 1 0 0  0   12    ... ... ... .. .. ..  A2 =  , M = ,  2 0 0  . . .   h2    .  ..   .. 1 −2 1   . 0 1 10 1 0    0 ··· 0 1 −2 (N −1)×(N −1)  0 0 1 10 1  0 · · · 0 −1 4 −5 14 (N −1)×(N −1) 

 u0 0  12   ..  , H2 = 2  .  h   0 uN (N −1)×1



u001 0 .. .



      00 U = .     0  u00N −1 (N −1)×1

Hence the fourth-order compact finite difference approximation of uxx is given by: U 00 = A−1 2 (M2 U + H2 ).

(16)

2.2. Exponential time differencing procedure This section includes brief discussion about ETD procedure, which motivated the present study. The Burgers’ equation (1) is written as: ∂u + Au = F(u, t), ∂t

(17)

where A and F represent linear and nonlinear operators of Burgers’ equation, respectively. Approximating Au with fourth-order compact approximation (16) and F(u, t) with fourth-order compact approximation (11), the following system of ODEs is obtained: ∂u + Au = F (u, t), u(x, 0) = g1 (x), ∂t

(18)

η −1 −1 2 where A = −δA−1 2 M2 and F (u, t) = − 2 A1 (M1 u + H1 ) + δA2 H2 . Let k = tn+1 − tn be the time step size, then using a variation of constant formula, the following recurrence formula is obtained: Z 1 −kA un+1 = e un + k e−kA(1−τ ) F (u(tn + τ k), tn + τ k)dτ, (19) 0

where un = u(tn ). The expression (19) is an exact solution of system (18) and its approximation leads to different ETD schemes. This paper considers two popular ETD schemes of Runge-Kutta type,

6

1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60 61 62 63 64 65

namely ETDRK4 [24] and ETDRK4-B [25] and propose the new modified forms for their general applicability. The Cox and Matthews ETDRK4 scheme is given as:

un+1 = ϕ0 (kA)un +kϕ1 (kA)Fn +kϕ2 (kA) −3Fn + 2Fna + 2Fnb − Fnc +4kϕ3 (kA) Fn − Fna − Fnb + Fnc , (20) where k k Fn = F (un , tn ), Fna = F (an , tn + ), Fnb = F (bn , tn + ) and Fnc = F (cn , tn + k), 2 2 k an = ϕ0 (kA/2)un + ϕ1 (kA/2)Fn , 2 k bn = ϕ0 (kA/2)un + ϕ1 (kA/2)Fna , 2
k cn = ϕ0 (kA/2)an + ϕ1 (kA/2) 2Fnb − Fn , 2 The ETDRK4-B [25] scheme, which has been shown to possess an even smaller error compared (20) is given as:

un+1 = ϕ0 (kA)un +kϕ1 (kA)Fn +kϕ2 (kA) −3Fn + 2Fna + 2Fnb − Fnc +4kϕ3 (kA) Fn − Fna − Fnb + Fnc , (21) where k k Fn = F (un , tn ), Fna = F (an , tn + ), Fnb = F (bn , tn + ) and Fnc = F (cn , tn + k), 2 2 k an = ϕ0 (kA/2)un + ϕ1 (kA/2)Fn , 2 k bn = ϕ0 (kA/2)un + ϕ1 (kA/2)Fn + kϕ2 (kA/2) (Fna − Fn ) , 2
cn = ϕ0 (kA)un + kϕ1 (kA)Fn + 2kϕ2 (kA) Fnb − Fn . Note that the schemes (20) and (21) are same but they contain different definitions of an , bn and cn . In addition, both of the schemes (20) and (21) contain matrix functions: ! µ−1 j X (−kA) , µ = 1, 2, 3. (22) ϕ0 (kA) = e−kA , ϕµ (kA) = (−kA)−µ e−kA − j! j=0 The accurate computation of (22) is a challenging problem in numerical analysis and is discussed in the literature (see [26, 27]). One major issue is the cancellation error arising during the direct computation of (22) for eigenvalues of A close to zero. To overcome this and other numerical issues associate with it, many researchers have proposed different methods. One well-known methods is the explicit computation of ϕµ (kA) proposed by Kassam and Trefethen [26]. This method includes a contour integration using contour that must encircle the spectrum 7

1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60 61 62 63 64 65

of the discretization matrix A, but that makes their approach unrealistic for the problems with large spectral radii [28]. Generally, for increasingly finer discretization, the spectrum of the discretization matrix is not easily known and it is typically unbounded. This is the primary limitation of the Kassam and Trefethen approach because the contour varies from problem to problem with dependence in spatial mesh. This limitation makes their approach problem dependent. A more serious and challenging issue in Kassam and Trefethen approach is computing the matrix exponential, when a matrix is large. Kassam-Trefethen used Matlab’s expm function, which has O(n3 ) complexity, to compute the matrix exponential. Computation of a matrix exponential with an expm function will make their scheme computationally inefficient and unrealistic for larger problems in multiple dimensions. Another well-known methods is the computation of the matrix vector product ϕµ (kA)u without the explicit computation of ϕµ (kA). There are three main techniques to the evaluation: the rational approximation of ϕµ (kA) [28–30], the Krylov subspace method [31, 32], and polynomial approximation of ϕµ (kA) [33, 34]. This study focuses on the first technique and develops a modified version of ETDRK4 and ETDRK4-B schemes based on Pad´e approximation of ϕµ (kA) in order to alleviate ˜ r,s (z) for (r, s)-Pad´e computational difficulties associated with them. The notations Rr,s (z) and R − z2 −z approximation are utilized to approximate e and e , respectively. The (r + s)th order rational Pad´e approximation to e−z is defined as: Rr,s (z) =

Pr,s (z) , Qr,s (z)

where

Pr,s (z) =

r X j=0

s

X (s + r − j)! (s + r − j)!r! (−z)j and Qr,s (z) = (z)j . (s + r)!j!(r − j)! (s + r)!j!(s − j)! j=0

Definition 2.2.1. A rational approximation Rr,s (z) of ez is said to be A−acceptable, if |Rr,s (z)| < 1, whenever <(z) < 0 and L−acceptable if, in addition, |Rr,s (z)| → 0 as <(z) → −∞. The rational Pad´e approximation Rr,s (z) to ez is : • A−acceptable if r = s. • L−acceptable if r = s − 1 or s − 2. The Fig. 1 demonstrates the behavior of the (1, 3)-Pad´e approximation given by R1,3 (z) = (1− z4 ) (12−6z+z 2 ) −z and (2, 2)-Pad´e approximation given by R2,2 (z) = (12+6z+z . As can be seen 2 ) of e 3 z2 z3

(1+ 4 z+

4

+ 24 )

from the Fig. 1, the (1, 3)-Pad´e approximation, possessing growth factors, which tend to zero asymptotically and satisfies the maximum modulus theorem [35] for L-acceptability given in definition 2.2.1. On the other hand, (2, 2)-Pad´e does not converge to zero for increasing z and was unable to satisfy maximum modulus theorem [35] for L-acceptability.

8

1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60 61 62 63 64 65

(a)

(b)

Fig. 1 Behavior of the functions R1,3 (z), and R2,2 (z) for z ∈ [0, 25] × [−10, 10]. 2.3. A modified ETDRK4-B scheme This subsection present a modified version of ETDRK4-B scheme by utilizing Fourth-order (2, 2)Pad´e approximation to e−kA . The modified scheme avoids direct computation of the matrix exponential and higher powers of matrix inverse. An advantage we found in utilizing (2, 2)-Pad´e approximation is that the factors A−1 and A−3 cancel out. To modify the Krogstad’s scheme (21), a (2, 2)-Pad´e approximation R2,2 (kA) = (12I − 6kA + 2 2 k A )(12I + 6kA + k 2 A2 )−1 is utilized into (21) to approximate matrix exponential functions, which yields:

un+1 = R2,2 (kA)un + P1 (kA)Fn + P2 (kA) −3Fn + 2Fna + 2Fnb − Fnc + P3 (kA) Fn − Fna − Fnb + Fnc , (23) where P1 (kA) = 12k(12I + 6kA + k 2 A2 )−1 , P2 (kA) = k(6I + kA)(12I + 6kA + k 2 A2 )−1 , P3 (kA) = 2k(4I + kA)(12I + 6kA + k 2 A2 )−1 .

In addition: ˜ 2,2 (kA)un + P˜1 (kA)Fn , an = R ˜ 2,2 (kA)un + P˜1 (kA)Fn + P˜2 (kA) (Fna − Fn ) , bn = R
cn = R2,2 (kA)un + P1 (kA)Fn + 2P2 (kA) Fnb − Fn , 9

1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60 61 62 63 64 65

with ˜ 2,2 (kA) = (48I − 12kA + k 2 A2 )(48I + 12kA + k 2 A2 )−1 , R P˜1 (kA) = 24k(48I + 12kA + k 2 A2 )−1 , P˜2 (kA) = 2k(12I + kA)(48I + 12kA + k 2 A2 )−1 . 2.3.1. The partial fraction form of modified ETDRK4-B scheme The scheme (23) introduced in subsection 2.3 consists of high order matrix polynomials to invert, which would be computationally burdensome and numerically unstable, if the matrices have high condition numbers. Moreover, round off error in computing the power of the matrices can produce bad ˜ 2,2 (kA) will not be computed approximations [36]. In order to handle this difficulty, R2,2 (kA) and R directly. Instead,the problem of stably computing the inverse of matrix polynomials inherent to Pad´e approximation is handled by a partial fraction decomposition approach as suggested by [29, 37]. This decomposition does reduce the computational complexity to just two LU decompositions over the entire time interval (provided the space step h, and time step k are held constant). In addition, this decomposition approach is important in alleviating ill-conditioning problems because only implicit Euler type solvers are required. To compute un+1 utilize: 2

R2,2 (kA) = (−1) +

q1 X j=1

  q1 +q2 X wj wj +2 < , kA − cj I kA − cj I j=1+q 1

and corresponding {Pi (z)}3i=1 takes the form: Pi (kA) = k

q1 X j=1

  q1 +q2 X wij wij + 2k < , kA − cj I kA − cj I j=1+q

i = 1, 2, 3,

1

where R2,2 as well as of Pi have q1 real and 2q2 complex poles {cj } with q1 + 2q2 = 2 and wj and wij are corresponding weights respectively. To evaluate an , bn , and cn employ: ˜ 2,2 (kA) = (−1)2 + R

q1 X j=1

  q1 +q2 X w˜j w˜j +2 < , kA − c˜j I kA − c˜j I j=1+q 1

and corresponding P˜i (z) as: P˜i (kA) = k

q1 X j=1

q1 +q2 X Ω˜j + 2k < kA − c˜j I j=1+q 1

Ω˜j kA − c˜j I

!

,

i = 1, 2,

˜ 2,2 as well as of P˜i have q1 real and 2q2 complex poles {˜ where R cj } with q1 + 2q2 = 2 and w˜j and Ω˜j are corresponding weights respectively. 10

1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60 61 62 63 64 65

2.3.2. Modified ETDRK4-B algorithm To implement the scheme (23) in a computationally efficient manner, a description of the algorithm is presented by implementing a partial fraction splitting technique. For i = 1, · · · , q1 + q2 , where q1 = 0, and q2 = 1.0 : Step 1. Solve (kA − c˜i I)Ra = w˜i un + k Ω˜i Fn , i

for Rai and then compute an as: an = un +

q1 X

q1 +q2

Rai + 2

i=1+q1

i=1

Step 2. Solve

X

<(Rai ).

  (kA − c˜i I)Rbi = w˜i un + k Ω˜i − Ω˜i+1 Fn + k Ω˜i+1 Fna ,

for Rb and then compute bn as: bn = un +

q1 X

q1 +q2

Rbi + 2

i=1

Step 3. Solve

X

i=1+q1

<(Rbi ).

(kA − ci I)Rci = wi un + k (w1i − 2w2i ) Fn + 2kw2i Fnb ,

for Rci and then compute cn as: cn = un +

q1 X

q1 +q2

Rci + 2

i=1

Step 4. Solve

X

i=1+q1

<(Rci ).


(kA − ci I)Rui = wi un + k (w1i − 3w2i + w3i ) Fn + k (2w2i − w3i ) Fna + Fnb − k (w2i − w3i ) Fnc ,

for Rui and then compute un+1 as: un+1 = un +

q1 X i=1

q1 +q2

Rui + 2

X

i=1+q1

<(Rui ).

In order to implement this Fourth-order scheme in its partial fraction form, poles and correspond˜ 1,3 (z), and {P˜i (z)}2 , which are as follows: ing weights were computed for R1,3 (z), {Pi (z)}3i=1 , R i=1 c1 = −3.0 + i1.7320508075688772935, w11 = −i3.4641016151377545871, w31 = 1.0 − i0.57735026918962576452,

w˜1 = −12.0 − i20.784609690826527522, ˜ 2 = 1.0 − i1.7320508075688772935. Ω

w1 = −6.0 − i10.39230484541326376, w21 = 0.5 − i0.8660254037844386467, c˜1 = −6.0 + i3.4641016151377545871, ˜ 1 = −i3.4641016151377545870, Ω

From this point the algorithm presented in subsection 2.3.2 is referenced with Krogstad-P22. 11

1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60 61 62 63 64 65

2.4. Modified ETDRK4 algorithm 2 2 )(I + 34 kA + k 4A + In order to modify ETDRK4, a (1, 3)-Pad´e approximation R1,3 (kA) = (I − kA 4 k3 A3 −1 ) is utilized into (21) to approximate matrix exponential functions. Following the same pro24 cedure mentioned in the previous subsections, the following algorithm results: For i = 1, · · · , q1 + q2 where q1 = 1.0, and q2 = 1.0 : Step 1. Solve (kA − c˜i I)Rai = w˜i un + k Ω˜i F (un , tn ), for Rai and then compute an as: an =

q1 X

q1 +q2

X

Rai + 2

i=1

i=1+q1

<(Rai ).

Step 2. Solve k (kA − c˜i I)Rbi = w˜i un + k Ω˜i F (an , tn + ), 2 for Rb and then compute bn as: bn =

q1 X

q1 +q2

Rbi + 2

i=1

X

i=1+q1

<(Rbi ).

Step 3. Solve k (kA − c˜i I)Rci = w˜i an + k Ω˜i (2F (bn , tn + ) − F (un , tn )), 2 for Rci and then compute cn as: cn =

q1 X

q1 +q2

Rci + 2

i=1

X

i=1+q1

<(Rci ).

Step 4. Solve
(kA − c1 I)Rui = wi un + k (w1i − 3w2i + w3i ) Fn + k (2w2i − w3i ) Fna + Fnb − k (w2i − w3i ) Fnc ,

for Rui and then compute un+1 as: un+1 =

q1 X i=1

q1 +q2

Rui + 2

X

i=1+q1

<(Rui ). 12

1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60 61 62 63 64 65

In order to implement this Fourth-order scheme in its partial fraction form, poles and correspond˜ 1,3 (z), and P˜ (z), which are as follows: ing weights were computed for R1,3 (z), {Pi (z)}3i=1 , R c1 = −2.6258168189584667160, c2 = −1.6870915905207666420 − i2.5087317549248805108, w1 = 5.5407990186788211678, w2 = −2.7703995093394105839 − i0.1591864442851235025, w11 = 2.1101239731096647932, w12 = −0.55506198655483239663 + i0.73103036863338010984, w21 = 0.80360669406735231150, w22 = 0.098196652966323844250 + i0.28728808194688171542, w31 = 1.2241626122055290024, w32 = 0.38791869389723549879 + i0.10430280315974445268, c˜1 = −5.2516336379169334320, c˜2 = −3.3741831810415332840 − i5.0174635098497610217, w˜1 = 11.081598037357642334, w˜2 = −5.5407990186788211672 − i0.31837288857024700490, ˜ 1 = 2.1101239731096647932, Ω ˜ 2 = −0.55506198655483239660 + i0.73103036863338010983. Ω From this point the algorithm presented in subsection 2.4 is referenced with an ETDRK4-P13.

3. Stability Analysis In this section the linear stability of the schemes ETDRK4-P13 and Krogstad-P22 were investigated by utilizing an approach suggested and discussed in [24, 38, 39] (and references therein) for the non-linear autonomous ODE: ut = −cu + F (u), (24) with F (u) as a nonlinear part, we assume that there exist a fixed point u0 such that −cu0 +F (u0 ) = 0. Linearizing about this fixed point, we thus obtain: ut = −cu + λu,

(25)

0

where u is perturbation of u0 , and λ = F (u0 ). If (24) represents a system of ODEs, then λ is a diagonal or a block diagonal matrix containing the eigenvalues of F . To keep the fixed point stable, we need that Re(λ − c) < 0, for all λ (see[24]). This approach only provides an indication to how stable a numerical scheme is, since in general one cannot linearize both terms simultaneously [25]. In general, the parameters c and λ may both be complex-valued. The stability region of the schemes ETDRK4-P13 and Krogstad-P22 is four dimensional and therefore difficult to represent [24]. The two-dimensional stability region is obtained, if both c and λ are purely imaginary or purely real [40], or if λ is complex and c is fixed and real [38]. The application of the schemes ETDRK4-P13 and Krogstad-P22 to the linearized Eq. (25) leads to a recurrence relation involving un , and un+1 . By letting r = uun+1 , x = λk, and y = −ck, we come n up with the following amplification factor: r(x, y) = c0 + c1 x + c2 x2 + c3 x3 + c4 x4 .

(26)

13

1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60 61 62 63 64 65

where for Krogstad-P22 1 1 1 1 5 c0 = 1 + y + y 2 + y 3 + y 4 + y + O(y 6 ), 2 6 24 144 1 1 7 4 5 5 c1 = 1 + y + y 2 + y 3 + y + y + O(y 6 ), 2 6 192 2304 1 1 1 23 3 11 4 37 5 c2 = + y + y 2 + y + y − y + O(y 6 ), 2 2 4 288 768 27648 1 1 47 2 77 3 5 4 515 5 c3 = + y + y + y + y − y + O(y 6 ), 6 6 576 3456 9216 165888 1 1 35 2 1 169 4 457 5 c4 = + y+ y − y3 − y − y + O(y 6 ). 24 32 3456 13824 82944 331776 and for ETDRK4-P13 1 1 1 1 c0 = 1 + y + y 2 + y 3 + y 4 + y 5 + O(y 6 ), 2 6 24 96 1 2 1 3 3 4 53 5 c1 = 1 + y + y + y + y + y + O(y 6 ), 2 6 64 3072 1 2 53 3 311 4 3323 5 1 1 y + y + y + O(y 6 ), c2 = + y + y + 2 2 4 576 9216 221184 1 1 25 2 13 3 193 4 401 5 c3 = + y + y + y + y + y + O(y 6 ), 6 6 288 384 13824 55296 1 1 5 2 13 3 131 4 223 5 c4 = + y+ y + y + y + y + O(y 6 ). 24 32 384 2304 36864 98304 The boundaries of the stability regions of the ETDRK4-P13 and Krogstad-P22 schemes are obtained by substituting r = eiθ , θ ∈ [0, 2π] into the Eq. (26) and solving for x, but unfortunately we do not know the explicit expression for |r(x, y)| = 1. We will only be able to plot it. In this paper various examples of stability regions of interest are demonstrated. • y ∈ <− : This study concentrates on the case where λ is complex and c is fixed and real. The analysis begins by selecting several real negative values of y and looking for a region of stability in the complex x plane where |r(x, y)| = 1. In the Fig. 2, the corresponding families of stability regions of ETDRK4-P13 and Krogstad-P22 in the complex x plane are plotted. According to Beylkin et al.[38] for method to be useful, it is important that stability regions grow as y → −∞. Clearly, as shown in Fig. 2, the stability regions for the schemes ETDRK4P13 and Krogstad-P22 grow larger as y → −∞. These regions give an indication of the stability of the schemes. The red curve corresponds to the case y = 0, where the stability region of the schemes ETDRK4-P13 and Krogstad-P22 coincides with that of the corresponding fourth-order Runge-Kutta (RK4) scheme. • Re(y) = 0 : This is the case where λ is complex and c is purely imaginary. Regions for values of y in the interval [−2πi, 2πi] are plotted in Fig. 3. As seen from Fig. 3, both of the schemes have almost similar intersection and include an interval of the imaginary axis. This argument gives 14

1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60 61 62 63 64 65

an indication of stability of the scheme ETDRK4-P13 and Krogstad-P22 to solve dispersive PDEs.

(a) ETDRK4-P13

(b) Krogstad-P22

Fig. 2 Stability regions for different y ∈ <− .

(a) ETDRK4-P13

(b) Krogstad-P22

Fig. 3 Stability regions for different y ∈ [−2πi, 2πi], and their intersection. Actually, for y < 0, |y| >> 1 the boundaries of the ETDRK4-P13 and Krogstad-P22 tend to ellipses, whose parameters have been fitted numerically with the following result:  2 Im(x) 2 (Re(x)) + = y2. (27) 0.7 15

1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60 61 62 63 64 65

For y = −25, in Fig. 4 the ellipse and experimental boundaries of the proposed schemes are drawn. It is clear from the Fig. 4 that the boundary of ETDRK4-P13 fits to ellipse better than Krogstad-P22. According to (27), when Re(x) = 0 the intersection with the imaginary axis increase as |y|. Since the boundary of stability of ETDRK4-P13 tends to ellipse closer than Krogstad-P22 and at the same time boundary of stability grows faster than x, the ETDRK4-P13 scheme should perform better than Krogstad-P22 to solve the dissipative PDES like (for example, Burgers’ equation or Kuramoto-Sivashinsky equation).

(a) ETDRK4-P13

(b) Krogstad-P22

Fig. 4 Ellipse and experimental boundary for y=-25. 4. Numerical experiments and discussions In this section several test examples are considered in order to illustrate the adaptability of the proposed schemes and to compare their accuracy with those which are already available in literature for solving Burgers’ and coupled Burgers’ equations. The accuracy of the schemes is measured in terms of maximum error norm L∞ , which is defined as: L∞ = max |ei |, i

where ei = ui − Ui , ui , and Ui are the ith exact and numerical solutions, respectively. The empirical order of spatial and temporal convergence of the schemes, when the exact solution of the problem(s) is/are available are calculated by the formulae: log10 (ku − uhi k/ku − uhi+1 k) log10 (hi /hi+1 )

and

log10 (ku − uki k/ku − uki+1 k) , log10 (ki /ki+1 )

where k · k represents L∞ norm, u is the exact solution, uhi and uki are the numerical solutions with spatial step size hi and temporal step size ki , respectively. 16

1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60 61 62 63 64 65

On the other hand, when the exact solution of the problem(s) is/are unavailable, the following formula is used to measure the temporal convergence rates of the schemes: log10 (Ek /E k ) 2 , log10 (2) where Ek = kuk − u2k k and E k = ku k − uk k are discrete maximum norm errors at k and k2 . All the 2 2 numerical experiments are conducted in MATLAB 7.12 platforms based on an Intel Core i7-4510U 2.60 GHz workstation. Example 1. Consider Burgers’ equation (1) over a domain Ω = [0, 2π] with periodic boundary conditions and following low regularity initial condition:  sin(x), 0 ≤ x < π, u(x, 0) = 0, π ≤ x < 2π. Limited work has been available in literature dealing numerically with Burgers’ equation involving periodic boundary conditions in contrast to Burgers’ equation with nonperiodic boundary conditions. This example is considered as a benchmark problem in order to show the adaptability of the proposed schemes to solve Burgers’ equation with periodic boundary conditions and to show how L−stable and A−stable schemes behave with problem involving low regularity initial condition. Three sets of numerical experiments on this example were conducted. In the first sets of experiment, the order of accuracy in temporal direction of the proposed schemes was checked by running the simulation until time T = 1.0 with the parameters δ = 0.1 and η = 1.0. To reduce the effect of spatial dimension, spatial step size length h = π/304 was selected small enough so that discretization error in space is negligible compared to the temporal error. Initially k = 0.2 was set and repeatedly halved it at each time. The computed results are captured in Table 1 and from the results, it is clear that L−stable scheme ETDRK4-P13 is able to maintain the expected fourth-order accuracy in time. On the other hand, A−stable scheme Krogstad-P22 faces order reduction, while dealing with the problem involving low regularity initial condition. Table 1 The maximum error, time rates of convergence, and CPU time of ETDRK4-P13 and Krogstad-P22 for Example 1 with δ = 0.1, η = 1.0, and h = π/304 at T = 1.0.

ETDRK4-P13

Krogstad-P22

k Ek Ratio Order CPU(s) k Ek Ratio Order CPU(s)

0.1 3.116E-04 4.5431 0.1 2.819E-03 2.0915

17

0.1/2 2.250E-05 13.8467 3.7915 8.1900 0.1/2 1.534E-03 1.8383 0.8784 4.4543

0.1/4 1.975E-06 14.8513 3.8925 15.6797 0.1/4 5.270E-04 2.9104 1.5412 8.5095

0.1/8 1.000E-07 15.1503 3.9213 29.7746 0.1/8 2.112E-05 24.9514 4.6410 16.6095

1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60 61 62 63 64 65

In the second sets of experiment, the simulation was run until T = 5.0 by ETDRK4-P13 and Krogstad-P22 with h = π/100, k = 0.1, δ = 1.0, and η = 1.0. The simulation profile of u(x, t) at T = 5.0 is captured in Fig. 5 (a) and (b). It is clear from the Fig. 5 (b) that the numerical solution contains unwanted oscillations at x = π, when Krogstad-P22 is used, which are killed of by ETDRK4-P13 scheme due to its ability to damp spurious oscillation cased by high frequency components in the solution.

(a) ETDRK4-P13

(b) Krogstad-P22

Fig. 5 (color online) Numerical solution of Example 1 via (a) ETDRK4-P13, (b) Krogstad-P22 with h = π/100, k = 0.1, δ = 1.0, η = 1.0 at T = 5.0.

One of the well known difficulties in solving Burgers’ equation with larger Reynolds numbers (small δ) is that as time advances, the solution curves steepen and develops a shock like discontinuity. So, in the third sets of experiment, simulation was run via the proposed schemes until final time T = 5.0 with δ = 0.1, 0.05, 0.01 and 0.005 to observe how they behave with mild to larger Reynolds numbers. The wave propagation scenario of u(x, t) via ETDRK4-P13 is displayed in Fig. 6 (a)-(d). It is clear from Fig. 6 (a) and (b) that there is no sharp front in the solutions for δ = 0.1 and δ = 0.05, which becomes sharper for δ = 0.01 and δ = 0.005 and eventually leads to a discontinuity for smaller δ. As can be seen from Fig. 6 (c) and (d), the scheme ETDRK4-P13 is able to capture the sharp front very well. Meanwhile, for δ = 0.1 and δ = 0.05, the Krogstad-P22 scheme also able to produce solution with no sharp front but, for δ = 0.01 and δ = 0.005 Krogstad-P22 scheme blows up.

18

1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60 61 62 63 64 65

(a)

(b)

(c)

(d)

Fig. 6 (color online) Numerical solution of Example 1 via ETDRK4-P13 with (a) h = π/40, k = 0.1, δ = 0.1, η = 1.0, (b) h = π/40, k = 0.1, δ = 0.05, η = 1.0, (c) h = π/100, k = 0.01, δ = 0.01, η = 1.0, (d) h = π/200, k = 0.005, δ = 0.005, η = 1.0. The final time T = 5.0 and the time interval between the curves in (a)-(b) is 0.5 and (c)-(d) is 0.25.

Example 2. In what follows, the Burgers’ equation given in Eq. (1) is solved for different values of δ and η over a domain Ω = [0, 1], with the homogeneous Dirichlet boundary conditions and following discontinuous initial conditions:   0, 0 ≤ x < 21 , 1, 12 ≤ x ≤ 43 , u(x, 0) =  0 34 < x ≤ 1. As can be seen from the Table 2, the scheme ETDRK4-P13 is much more accurate than the 19

1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60 61 62 63 64 65

scheme Krogstad-P22 and not suffer with order reduction, while solving problem with discontinuous initial condition. In the Fig. 7 (a) and (b), the space-time evolution of the component u(x, t) obtained by ETDRK4-P13 and Krogstad-P22, respectively is captured. From the Fig. 7 (b), it can be seen that the solution obtained by Krogstad-P22 oscillate very severely at the points of discontinuity, which was not case for ETDRK4-P13. The ETDRK4-P13 produces very smooth solution, which is depicted in Fig. 7 (a). Table 2 The maximum error, time rates of convergence, and CPU time of ETDRK4-P13 and Krogstad-P22 for Example 2 with δ = 1.0, η = 0.001, and h = 0.01 at T = 1.0.

ETDRK4-P13

Krogstad-P22

k Ek Ratio Order CPU(s) k Ek Ratio Order CPU(s)

0.05 3.0264E-07 0.2699 0.05 1.4840E-01 0.1762

(a) ETDRK4-P13

0.05/2 2.2626E-08 13.3759 3.7416 0.51182 0.05/2 1.3301E-01 1.1157 0.1581 0.2959

0.05/4 1.5651E-09 14.4575 3.8537 0.8909 0.05/4 1.4453E-01 0.9202 -0.1199 0.5802

0.05/8 1.0466E-10 15.9536 3.9024 1.7406 0.05/8 4.6879E-02 3.0831 1.6243 1.2185

(b) Krogstad-P22

Fig. 7 (color online) Numerical solution of Example 2 via (a) ETDRK4-P13, (b) Krogstad-P22 with h = 0.001, k = 0.001, δ = 1.0, η = 1.0. at T = 0.05.

Example 3. In this example, the coupled Burgers’ equation given in Eq.(4) is considered in the region −π ≤ x ≤ π with an initial and boundary conditions extracted from the following exact solution given by Kaya [42] for δ = µ = 1.0, η = ξ = −2.0, α = β = 1.0: 20

1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60 61 62 63 64 65

u(x, t) = v(x, t) = e−t sin(x). To investigate the order of accuracy of the proposed schemes for coupled Burgers’ equation, a numerical test on Example 3 was performed. In the computation, a simulation was run up to T = 1.0 by initially setting h = π/8 and k = 0.1 then reduced both of them by a factor of 2 in each refinement. The maximum error, rates of convergence and CPU times are listed in Table 3. As can be seen from Table 3 that the computed convergence rates of the proposed schemes apparently demonstrate the fourth-order accuracy in both time and space. Meanwhile, Table 3 shows that both the schemes have same accuracy but the scheme Krogstad-P22 is computationally cheaper than scheme ETDRK4-P13. Table 3 The maximum error, rates of convergence, and CPU time of ETDRK4-P13 and Krogstad-P22 for Example 3 at T = 1.0. ETDRK4-P13 h π/8 π/16 π/32 π/64 π/128

k 0.1 0.1/2 0.1/4 0.1/8 0.1/16

L∞ (u) 3.660E-05 2.277E-06 1.422E-07 8.882E-09 5.552E-10

Ratio 16.0739 16.0187 16.0048 15.9973

Krogstad-P22 order 4.0066 4.0017 4.0004 3.9998

CPU(s) 0.1492 0.3109 0.5419 1.0726 3.4372

L∞ (u) 3.673E-05 2.285E-06 1.426E-07 8.913E-09 5.571E-10

Ratio 16.0731 16.0183 16.0046 15.9991

order 4.0666 4.0017 4.0004 3.9999

CPU(s) 0.1170 0.1704 0.2980 0.6868 2.1609

To compare the accuracy of the proposed schemes with other existing schemes in terms of maximum errors at various time levels from T = 0.5 to 3.0, the maximum errors of the components u(x, t) and v(x, t) were computed with exactly same values of N and k considered in [12, 13] via the proposed schemes and compared the results with results provided in [12, 13]. The computed results were tabulated in Table 4. From the Table 4, it can be seen that the proposed schemes provide more accurate solutions than the solutions obtained by the schemes in [12, 13]. Table 4 Comparison of accuracy at different times with k = 0.01 for Example 3. ETDRK4-P13 N = 31

T = 0.5

T = 1.0

Results reported in [12]

Krogstad-P22

T = 2.0

T = 3.0

T = 0.5

T = 1.0

T = 2.0

T = 3.0

T = 0.5

T = 1.0

T = 2.0

T = 3.0

L∞ (u)

2.422E-06 2.938E-06

2.162E-06

1.193E-06

1.517E-04

1.841E-04

1.352E-04

7.460E-05

2.422E-06

2.938E-06

2.162E-06

1.193E-06

L∞ (v)

2.422E-06

2.162E-06

1.193E-06

1.517E-04

1.841E-04

1.352E-04

7.460E-05

2.422E-06

2.938E-06

2.162E-06

1.193E-06

T = 2.0

T = 3.0

T = 0.5

T = 1.0

T = 2.0

T = 3.0

T = 0.5

T = 1.0

T = 2.0

T = 3.0

2.938E-06

ETDRK4-P13 N = 101 T = 0.5

T = 1.0

Results reported in [13]

Krogstad-P22

L∞ (u)

1.969E-08 2.389E-08

1.757E-08

9.698E-08

2.507E-06

3.042E-06

2.2380E-06

1.235E-06

1.969E-08

2.389E-08

1.757E-08

9.698E-08

L∞ (v)

1.969E-08

1.757E-08

9.698E-08

2.507E-06

3.042E-06

2.2380E-06

1.235E-06

1.969E-08

2.389E-08

1.757E-08

9.698E-08

2.389E-08

In the Table 5, the maximum errors computed by the proposed schemes with different values of N and fixed value of k = 0.0001 are compared with results given in [11]. As can be seen from the Table 21

1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60 61 62 63 64 65

5, the scheme introduced in [11] provided better accuracy than our proposed schemes for fewer mesh points. For the refined mesh size the proposed schemes provided far better accuracy than scheme in [11]. For the scheme introduced in [11], the maximum errors stay at the same order of magnitude while refining the mesh size. On the other hand, maximum errors getting smaller and smaller with refining mesh size for the proposed schemes. Over all the proposed schemes provided more accurate solutions in spatial direction than scheme in [11]. Table 5 Comparison of errors for u(x, t) and v(x, t) with k = 0.0001 for Example 3 at T = 1.0.

ETDRK4-P13

N L∞ (u) L∞ (v)

4 2.934E-02 2.934E-02

8 1.000E-03 1.000E-03

16 4.726E-05 4.726E-05

32 2.588E-06 2.588E-06

64 1.517E-07 1.517E-07

128 9.184E-09 9.184E-09

Results reported in [11]

N L∞ (u) L∞ (v)

4 1.542E-03 1.542E-03

8 1.165E-05 1.165E-05

16 1.164E-05 1.164E-05

32 1.163E-05 1.163E-05

64 1.160E-05 1.160E-05

128 1.159E-05 1.159E-05

Krogstad-P22

N L∞ (u) L∞ (v)

4 2.934E-02 2.934E-02

8 1.000E-03 1.000E-03

16 4.726E-05 4.726E-05

32 2.588E-06 2.588E-06

64 1.517E-07 1.517E-07

128 9.184E-09 9.184E-09

In the Fig. 8 (a) the evolution profile of u(x, t) until final time T = 3.0 was illustrated. In the Fig 8 (b), the comparisons of analytical solutions with numerical solutions of u(x, t) obtained via ETDRK4-P13 at different time T = 0.0, T = 1.0, T = 2.0, and T = 3.0 were visualized. From the Fig. 8 (b) it is observed that the numerical solutions show great agreement with the analytical solutions. The evolution profile clearly describes the wave propagation, where it can be seen that the amplitude of the wave reduces with time and wavelength stays the same.

(a)

(b)

Fig. 8 (color online) (a) Evolution profile of u(x, t) via ETDRK4-P13 until final time T = 3.0, (b) A comparison between numerical and exact solutions at different time levels, the solid line represents an exact solution. The parameters are taken as h = π/40 and k = 0.1. 22

1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60 61 62 63 64 65

Example 4. For this example, the coupled Burgers’ equation given in Eq. (4) is considered with the exact solution given in [43] for δ = µ = 1.0 and η = ξ = 2.0:   2α − 1 u(x, t) = a0 − 2A tanh(A(x − 2At)), 4αβ − 1     2β − 1 2α − 1 v(x, t) = a0 − 2A tanh(A(x − 2At)), 2α − 1 4αβ − 1
and a0 , α, β are arbitrary constants. The boundary and initial conditions are where A = a0 4αβ−1 4α−2 taken from the exact solution. In the Example 4, the computations have been carried out on the fixed domain [−10, 10] with k = 0.01, a0 = 0.05 and different values of α, β and N as in [9, 11, 12] for comparisons. Tables 6 and 7 depict the comparisons of maximum errors of u(x, t) and v(x, t) at different time levels with those already available in the literature. As can be seen from tables, the maximum errors given by the proposed schemes are of the higher order of magnitude than Rashid [11] and are of the same order of magnitude than Khater [9] and Mittal [12]. In some cases the proposed schemes in this paper show better results especially, when α = 0.1 and β = 0.3. Table 6 Comparison of errors at different time for u(x, t) with k = 0.01, δ = µ = 1.0, η = ξ = 2.0 and a = 0.05 for Example 4.

β

Results reported in [9] L∞ (u)

Results reported in [11] L∞ (u)

Results reported in [12] L∞ (u)

ETDRK4-P13 L∞ (u)

Krogstad-P22 L∞ (u)

T = 0.5 -

T = 0.5 T = 1.0 9.619E-04 1.153E-03 4.310E-04 1.268E-03

T = 0.5 -

T = 0.5 T = 1.0 4.183E-05 8.235E-05 4.589E-05 9.179E-05

T = 0.5 T = 1.0 4.183E-05 8.235E-05 4.589E-05 9.179E-05

4.166E-05 8.274E-05 4.588E-05 9.169E-05

4.166E-05 8.274E-05 4.588E-05 9.169E-05

N

α

16

0.1 0.3 0.3 0.03

20

0.1 0.3 4.38E-05 8.66E-05 0.3 0.03 4.58E-05 9.16E-05

T = 1.0 -

-

-

T = 1.0 -

4.173E-05 8.275E-05 4.585E-05 9.167E-05

Table 7 Comparison of errors at different time for v(x, t) with k = 0.01, δ = µ = 1.0, η = ξ = 2.0 and a = 0.05 for Example 4.

β

Results reported in [9] L∞ (v)

Results reported in [11] L∞ (v)

Results reported in [12] L∞ (v)

ETDRK4-P13 L∞ (v)

Krogstad-P22 L∞ (v)

T = 0.5 -

T = 0.5 T = 1.0 3.332E-04 1.162E-03 1.148E-03 1.639E-03

T = 0.5 -

T = 0.5 T = 1.0 2.141E-05 4.188E-05 1.807E-04 3.612E-04

T = 0.5 T = 1.0 2.141E-05 4.188E-05 1.807E-04 3.612E-04

2.173E-05 4.155E-05 1.809E-04 3.616E-04

2.173E-05 4.155E-05 1.809E-04 3.616E-04

N

α

16

0.1 0.3 0.3 0.03

20

0.1 0.3 4.99E-05 9.92E-05 0.3 0.03 1.81E-04 3.62E-04

T = 1.0 -

-

-

T = 1.0 -

5.418E-05 1.074E-04 2.826E-05 5.673E-05

Further, the evolution profiles of u(x, t) and v(x, t) are depicted in Fig. 9 (a) and (c). In addition, the comparison between the analytical and numerical solutions are captured in Fig. 9 (b) and (d) 23

1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60 61 62 63 64 65

for N = 100, k = 0.01, α = 1.0, β = 2.0, a0 = 0.05, and T = 1.0. It is clear from the Fig. 9 (b) and (d) that there is a good agreement between numerical and analytical solutions.

(a)

(b)

(c)

(d)

Fig. 9 (color online) (a) Evolution profile of u(x, t) via ETDRK4-P13 for Example 4, (b) A comparison between

numerical and analytical solution of u(x, t) (c) Evolution profile of v(x, t) via ETDRK4-P13 for Example 4, (d) A comparison between numerical and analytical solution of v(x, t). The parameters are taken as N = 100, k = 0.01, and final time T = 1.0. Notice the different scales on the vertical axes.

Example 5. The coupled Burgers’ equation given in (4) is considered with the following initial conditions:  sin(2πx), 0 ≤ x ≤ 12 , u(x, 0) = 1 0, ≤ x ≤ 1. 2  0, 0 ≤ x ≤ 12 , v(x, 0) = −sin(2πx), 21 ≤ x ≤ 1, 24

1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60 61 62 63 64 65

and absorbing boundary conditions. In order to show that L−stable methods are always preferable among the A−stable methods for the problem like Example 5, which consists of low regularity initial condition, maximum errors and time rates of convergence were computed in the domain [0, 1] with δ = µ = 1.0, η = ξ = 20, α = β = 1.0 via ETDRK4-P13 and Krogstad-P22 for a comparison. A simulation was run until T = 0.1 and recorded the maximum errors, rates of convergence, and CPU times in Table 8. In the simulation, spatial step size was fixed at h = 0.005 and reduced time step k by the factor of 2 at each time. Comparing the results in Table 8, it can be seen that ETDRK4-P13 performs better than Krogstad-P22 in terms of accuracy and was able to maintain fourth-order accuracy in time. The accuracy of the Krogstad-P22 scheme is badly affected by high frequency components in the solutions, as a consequence it suffers with order reduction. The time rates of convergence of the proposed schemes are visualized in Fig. 10 with log-log scale graphs. Table 8 The maximum error, rates of convergence, and CPU time of ETDRK4-P13 and Krogstad-P22 for Example 5 with δ = µ = 1.0, η = ξ = 20, α = β = 1.0, and h = 0.005 at T = 0.1.

ETDRK4-P13

Krogstad-P22

k Ek Ratio Order CPU(s) k Ek Ratio Order CPU(s)

0.025 1.9598E-02 0.2797 0.025 1.7715E-01 0.2179

0.025/2 1.1363E-03 17.2474 4.1083 0.4909 0.025/2 2.3929E-02 4.6971 2.2318 0.3559

(a) Time rate of convergence of ETDRK4-P13

0.025/4 1.0461E-04 10.8633 3.4414 0.8681 0.025/4 1.0368E-02 2.0844 1.0597 0.5906

0.025/8 6.7428E-06 15.5128 3.9554 1.5861 0.025/8 5.8418E-03 1.7789 0.8310 1.0717

0.025/16 3.5468E-07 19.0111 4.2488 2.9965 0.025/16 2.7046E-03 2.1596 1.1108 2.0044

(b) Time rate of convergence of Krogstad-P22

Fig. 10 (color online) (a) The log-log graph of L∞ error VS. time step k of ETDRK4-P13, (b) The log-log graph of L∞ error VS. time step k of Krogstad-P22. 25

1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60 61 62 63 64 65

In the Table 9, the maximum values of u(x, t) and v(x, t) were computed at different time levels for α = β = 10.0, N = 51, and k = 0.01 via the ETDRK4-P13 and were compared the results with those already available in the literature [8, 13]. From the Table 9, it is clear that the maximum values and their positions at different time levels obtained with ETDRK4-P13 are in good agreement with results reported in [8, 13]. The scheme in [13] are only second-order in time so that, the authors used extremely small temporal step (k = 0.00001) to get the same results as we obtained with 100 times larger (k = 0.01) temporal step. This results justify that the higher-order schemes allow us to use large discretization step in comparison with second-order scheme. The computed results are demonstrated graphically in Fig. 11 (a) and (b). The dots on solution curves of Fig. 11 (a) and (b) represent a position of the maximum values of the components u(x, t) and v(x, t), respectively for α = β = 10.0. Table 9 Maximum values of u(x, t), and v(x, t) at different time with α = β = 10.0, η = ξ = 2.0, and N = 51.

ETDRK4-P13

k=0.01

Results reported in [8]

k=0.01

Results reported in [13]

k=0.00001

T 0.1 0.2 0.3 0.4 T 0.1 0.2 0.3 0.4 T 0.1 0.2 0.3 0.4

Max(u) 0.1445 0.0524 0.0193 0.0072 Max(u) 0.1446 0.0524 0.0193 0.0072 Max(u) 0.1445 0.0524 0.0193 0.0072

26

At point(x) 0.58 0.54 0.52 0.50 At point(x) 0.58 0.54 0.52 0.50 At point(x) 0.58 0.54 0.52 0.50

Max(v) 0.1432 0.0470 0.0173 0.0064 Max(v) 0.1431 0.0471 0.0172 0.0064 Max(v) 0.1432 0.0470 0.0173 0.0064

At point (x) 0.66 0.56 0.52 0.50 At point (x) 0.66 0.56 0.52 0.50 At point (x) 0.66 0.56 0.52 0.50

1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60 61 62 63 64 65

(a)

(b)

Fig. 11 (color online) (a) The solution profile of u(x, t) of Example 5 via ETDRK4-P13 at different time levels for α = β = 10.0, η = ξ = 2.0, k = 0.01, and N = 51, (b) The solution profile of v(x, t) of Example 5 via ETDRK4-P13 at different time levels for α = β = 10.0, η = ξ = 2.0, k = 0.01, and N = 51. To examine the performance of L−stable and A−stable schemes for solving Example 5, the solution profile of the component u(x, t) with h = 0.01, k = 0.02, η = ξ = 2.0, and α = β = 10.0 was simulated until different time levels and capture the space-time evolution profile in the Fig. 12 (a) and (b). As can be seen from the Fig 12 (b), the Krogstad-P22 scheme suffers with spurious oscillation close to x = 0.5 due to its inablity to damp high frequency components in the solution. On the other hand, scheme ETDRK4-P13 does not suffer with this phenomena and provides smooth solution.

(a) ETDRK4-P13

(b) Krogstad-P22

Fig. 12 Evolution profile of u(x, t) of Example 5 via (a) ETDRK4-P13, (b) Krogstad-P22, at different time. The parameters are takes as h = 0.01, k = 0.02, η = ξ = 2.0, and α = β = 10.0.

27

1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60 61 62 63 64 65

In order to see how the schemes ETDRK4-P13 and Krogstad-P22 approximate the solution of coupled Burgers’ equation with high Reynolds numbers, a space-time evolution profile of the component u(x, t) with δ = µ = 0.01 and 0.005 was simulated until the final time T = 1.0 and captured the wave propagation scenario in Fig. 13. Form the Fig. 13 (a) and (b), it can be seen that, as the the time goes on, the maximum amplitude of u(x, t) keeps decreasing, leading to steeper the solution profile for larger Reynolds numbers and the scheme ETDRK4-P13 is able to capture the sharp front of the solution reasonably well. On the other hand, for δ = µ = 0.01 and 0.005, scheme Krogstad-P22 blows up hence, plot of the wave propagation scenario of the component u(x, t) simulated via Krogstad-P22 is not provided here.

(a)

(b)

Fig. 13 (color online) Numerical solution of Example 5 via ETDRK4-P13 with (a) N = 101, k = 0.001, δ = µ = 0.01, (b) N = 161, k = 0.001, δ = µ = 0.005. The final time T = 1.0 and the time interval between the curves in (a) and (b) is 0.05. Example 6 ( 2D Coupled Burgers’ Equation ). In order to show the adaptability of the proposed schemes to solve Burgers’ equation in higher spatial dimensions, we considered the following 2D coupled Burgers’ equation defined over the square domain Ω = [0, 1]2 :  ∂ 2u ∂ 2u  ∂u ∂u ∂u + +u +v =δ , ∂t ∂x ∂y ∂x2 ∂y 2  ∂ 2v ∂ 2v  ∂v ∂v ∂v + +u +v =µ , ∂t ∂x ∂y ∂x2 ∂y 2

(28)

subject to the initial conditions: u(x, y, 0) = sin(πx) sin(πy), v(x, y, 0) = [sin(πx) + sin(2πx)] [sin(πy) + sin(2πy)],

28

(29)

1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60 61 62 63 64 65

and boundary conditions: u(0, y, t) = u(1, y, t) = 0, v(0, y, t) = v(1, y, t) = 0, t ∈ [0, T ], u(x, 0, t) = u(x, 1, t) = 0, v(x, 0, t) = v(x, 1, t) = 0, t ∈ [0, T ],

(30)

1 where δ = µ = Re . No closed form solution has been reported for this problem. In order to decide whether the proposed scheme ETDRK4-P13 provides reliable and accurate results for 2D coupled Burgers equation, we compared the results computed by ETDRK4-P13 using central difference approximation in space with results provided in [19]. Computations are carried out with exactly same values of the parameters Re and N those considered in [19]. We selected time step length k large in comparison to values considered in [19] in order to show the strength of higher-order time stepping scheme over low order accurate time stepping schemes. Table 10 reported the comparison of the solutions obtained at final time T = 0.01 by using ETDRK4-P13 with results given in [19] using finite element and method of lines. It is seen from the Table 10 that the present results computed for large time step are in excellent agreement with results given in [19] for refined time steps. Moreover, we have computed solutions at different time with large Reynolds number Re = 100 and compared the results with those given in [20]. The results computed with ETDRK4P13 and results given in [20] are stored in Table 11. As we can see from the Table 11, the results obtained with ETDRK4-P13 scheme are in good agreement with results given in [20].

Table 10 Comparison of the results obtained by using the scheme ETDRK4-P13, method of lines and the finite element method [19] at T = 0.01 and Re = 1.0 for Example 6. Points

ETDRK4-P13

Method of Lines [19]

N = 10 N = 20 N = 40 k = 1/200 N = 1/200 N = 1/200 u (0.1, (0.2, (0.4, (0.7, (0.9,

0.1) 0.8) 0.4) 0.1) 0.9)

0.07277 0.28888 0.72315 0.20139 0.07956

0.07258 0.28846 0.72206 0.20113 0.07949

0.1) 0.8) 0.4) 0.1) 0.9)

0.43857 -0.13200 1.66510 0.06138 0.01459

0.43302 -0.12386 1.65573 0.06571 0.01372

N = 10 N = 20 N = 40 N = 1/2000 N = 1/1500 N = 1/2000

u 0.07253 0.28836 0.72178 0.20106 0.07947

0.07277 0.28887 0.72315 0.20139 0.07956

v (0.1, (0.2, (0.4, (0.7, (0.9,

Finite Element Method [19]

N = 10 N = 20 N = 40 N = 1/200 N = 1/2000 N = 1/4000

0.07257 0.28846 0.72205 0.20112 0.07948

u 0.07253 0.28836 0.72178 0.20106 0.07947

0.07279 0.28867 0.72370 0.20157 0.07951

v 0.43162 -0.12182 1.65336 0.06681 0.01349

0.43857 -0.13200 1.66509 0.06137 0.01459

0.43302 -0.12387 1.65571 0.06571 0.01372

29

0.07257 0.28842 0.72210 0.20113 0.07947

0.07252 0.28835 0.72179 0.20107 0.07946

v 0.43173 -0.12184 1.65335 0.06679 0.01349

0.44130 -0.13172 1.66212 0.06306 0.01459

0.44336 -0.12366 1.65499 0.06621 0.01367

0.43178 -0.12180 1.65316 0.06692 0.01349

1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60 61 62 63 64 65

Table 11 Comparison of the results computed via the scheme ETDRK4-P13 and results reported in [20] at different times and Re = 100.0 for Example 6.

Points

(0.1, 0.1) (0.2, 0.8) (0.4, 0.4) (0.7, 0.1) (0.8, 0.8) (0.9, 0.9)

u v u v u v u v u v u v

ETDRK4-P13 N = 80, k = 1/500

Results reported in [20] N = 80, k = 1/500

t = 0.5

t = 1.0

t = 0.5

t = 1.0

0.01510 0.12169 0.15863 0.99333 0.12819 0.70023 0.13362 0.09995 0.56408 1.18625 0.27905 0.21695

0.00727 0.05546 0.08073 0.58211 0.07043 0.36899 0.06823 0.07444 0.29573 0.69686 0.36783 0.75440

0.01509 0.12162 0.15839 0.98654 0.12822 0.70021 0.13353 0.09999 0.56378 1.18512 0.28128 0.22196

0.00726 0.05542 0.08076 0.58111 0.07045 0.36900 0.06816 0.07446 0.29571 0.69677 0.36648 0.75237

To demonstrate the strength of ETDRK4-P13 scheme further, the simulation was run until final time T = 1.0 with large Reynolds number, namely Re = 200. The solution profiles of u at different time levels corresponding to Re = 200 are depicted in Fig. 14 (a)-(d). Like as the case in Example 1, as the time goes on the solution profiles of u getting steeper and develop sharp wave front, which was captured by the proposed scheme ETDRK4-P13 reasonably well.

30

1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60 61 62 63 64 65

(a) t=0.0

(b) t=0.25

(c) t=0.5

(d) t=1.0

Fig. 14 (color online) Evolution profiles of u via ETDRK4-P13 for Example 6 with Re = 200 and (a) h = 0.01, k = 0.01 (b) h = 0.01, k = 0.01 (c) h = 0.005, k = 0.005 (d) h = 0.005, k = 0.005 at different time levels.

In order to do an empirical convergence analysis of the scheme ETDRK4-P13, we ran a simulation until final time T = 0.5. In the simulation, the spatial step size length h = 0.01 was selected small enough so that the discritization error in space is negligible. Moreover, initially k = 0.01 was set and repeatedly reduced it by the factor of 2 each time of simulation. The maximum errors, rates of convergence, and CPU times for the scheme ETDRK4-P13 are reported in Table 12. As can be seen from the Table 12, the proposed scheme is able to maintain the expected fourth-order accuracy in time.

31

1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60 61 62 63 64 65

Table 12 The maximum error, time rates of convergence, and CPU time of ETDRK4-P13 for Example 6 with Re = 100 and h = 0.01 at T = 0.5.

ETDRK4-P13

k Ek Ratio Order CPU(s)

0.01/2 2.03E-03 33.1815

0.01/4 1.36E-04 14.8882 3.8961 62.0093

0.01/8 9.03E-06 15.1181 3.9182 119.2379

0.01/16 5.66E-07 15.9534 3.9958 229.1333

Acknowledgment The authors are grateful to the anonymous referees and the editor for their valuable comments and suggestions, which have improved the quality of the paper. 5. Conclusions In this article two new modified fourth-order ETDRK schemes in combination with global fourthorder compact finite difference scheme have been developed for the numerical solution of nonlinear coupled viscous Burgers’ equations without using any transformations or linearization techniques. Stability of the proposed schemes was examined by plotting their stability regions, which provided an explanation of their behavior for dispersive and dissipative problems. The performance and adaptability of the proposed schemes have been tested with considered test problems. The numerical results showed that the proposed schemes provide better accuracy in comparison with other existing schemes [8, 9, 11–13, 19, 20] available in literature. Meanwhile, the successful implementation of the scheme ETDRK4-P13 to 2D problem showed its ability to solve problems posed in higher spatial dimensions. For the test problems with nonsmooth initial data and high Reynolds number, the computed results showed that the ETDRK4-P13 is more beneficial than Krogstad-P22 because of its quality to damp spurious oscillations caused by high frequency components in the solution. For the test problems with smooth initial data and mild Reynolds number, the numerical results demonstrated that the Krogstad-P22 is more desirable than ETDRK4-P13 because it is much more efficient and provides same accuracy as ETDRK4-P13. As the authors in [44] successfully implemented linear implicit explicit Runge-Kutta scheme to nonlinear degenerated convection-diffusion, it is believed that the proposed scheme ETDRK4-P13 with high order accurate weighted essentially non-oscillatory (WENO) schemes will also applicable for solving such PDEs. In the future, we plan to implement ETDRK4-P13 to degenerated convection-diffusion PDEs with modification in fourth-order compact scheme and approximating convective term with WENO scheme. References [1] E. Hopf, The partial differential equation ut + uux = νuxx , Commun. Pure Appl. Math. 3 (1950) 201-230. [2] J. D. Cole, On a quasilinear parabolic equations occuring in aerodynamics, Quart. Appl. Math. 9 (1951) 225-236. 32

1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60 61 62 63 64 65

[3] W. Liao, An implicit fourth-order compact finite difference scheme for one-dimensional Burgers’ equation, Appl. Math. Comput. 206 (2008) 755-764. [4] 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 (2007) 251-256. [5] W. Liao, J. Zhu, Efficient and accurate finite difference schemes for solving one-dimensional Burgers’ equation, Int. J. Comput. Math. 88 (2011) 2575-2590. [6] R. Jiwari, A hybrid numerical scheme for the numerical solution of the Burgers equation, Comput. Phys. Comm. 188 (2015) 59-67. [7] S. E. Esipov, Coupled Burgers’ equations: a model of polydispersive sedimentation, Phys. Rev. E 52 (1995) 3711-3718. [8] R. C. Mittal, G. Arora, Numerical solution of the coupled viscous Burgers’ equation, Commun. Nonlinear Sci. Numer. Simul. 16 (2011) 1304-1313. [9] H. Khater, R. S. Temash, M. M. Hassan, An Chebyshev spectral collocation method for solving Burgers’-type equations, J. Comput. Appl. Math. 222 (2008) 333-350. [10] V. K. Sirvastava, M. K. Awasthi, M. Tamsir, A fully implicit finite difference solution to one dimensional coupled Burgers’ nonlinear equation, Int. J. Math. Sci. 7 (2013) 23-28. [11] A. Rashid, A. I. Md. Ismail, A Fourier Pseudospectral method for solving coupled viscous Burgers’ equations, Comput. Methods Appl. Math. 9 (2009) 412-420. [12] R. C. Mittal, R. Jiwari, A differential quadrature method for numerical solutions of Burgers’type equations, Int. J. Numer. Methods Heat Fluid Fow 22 (2012) 880-895. [13] R. K. Mohanty, W. Dai, F. Han, Compact operator method of accuracy two in time and four in space for the numerical solution of coupled viscous Burgers’ equations, Appl. Math. Comput. 256 (2015) 381-393. [14] Siraj-ul-Islam, S. Haq, M. Uddin, A mesh free interpolation method for the numerical solution of the coupled nonlinear partial differential equations, Eng. Anal. Boundary Elements 33 (2009) 399-409. [15] M. A. Abdou, A. A. Soliman, Variational iteration method for solving Burgers’ and coupled Burgers’ equations, J. Comput. Appl. Math. 181 (2005) 245-251. [16] M. Kumar, S. Pandit, A composite numerical scheme for the numerical simulation of coupled Burgers’ equation, Comput. Phys. Comm. 185 (2014) 809-817. [17] A. Rashid, M. Abbas, A. I. Md. Ismail, A. A. Majid, Numerical solution of the coupled viscous Burgers equations by Chebyshev-Legendre Pseudo-Spectral method, Appl. Math. Comput. 245 (2014) 372-381. 33

1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60 61 62 63 64 65

[18] R. C. Mittal, H. Kaur, V. Mishra, Haar wavelet-based numerical investigation of coupled viscous Burgers’ equation, Int. J. Comput. Math. 92 (2015) 1643-1659. [19] P. Arminjon and C. Beauchamp, Numerical solution of Burgers equations in two-space dimensions, Comput. Meth. Appl. Mech. Eng. 19 (1979) 351-365. [20] G. W. Wei, D. S. Zhang, D. J. Kouri, D. K. Hoffman, Distributed approximation functional approach to Burgers’ equation in one and two space dimensions, Comput. Phys. Commun. 111 (1998), 93-109. [21] B. Fornberg, Calculation of weights in finite difference formulas, SIAM Rev. 40 (3) (1998) 685691. [22] J. Zhao1, W. Dai, T. Niu, Fourth-order compact schemes of a heat conduction problem with Neumann boundary conditions, Numer. Methods Partial Differ. Eqs. 23 (2007) 949-959. [23] J. Zhao, Highly accurate compact mixed methods for two point boundary value problems, Appl. Math. Comput. 188 (2007) 1402-1418. [24] S. M. Cox, and P. C. Matthews, Exponential time differencing for stiff systems, J. Comp. Phys. 176 (2002) 430-455. [25] S. Krogstad, Generalized integrating factor methods for stiff PDEs, J. Comput. Phys. 203 (2005) 72-88. [26] A. K. Kassam, L. N. Trefethen, Fourth-order time stepping for stiff PDEs, SIAM J. of Sci. Comput. 26 (4) (2005) 1214-1233. [27] N. J. Higham, Functions of matrices: Theory and computation, SIAM, Philadelphia (2008). [28] A. Q. M. Khaliq, J. M. Vaquero, B. A. Wade, M. Yousuf, Smoothing schemes for reactiondiffusion systems with nonsmooth data, J. Comput. Appl. Math. 223 (2009) 374-386. [29] H. P. Bhatt, A. Q. M. Khaliq, Higher order exponential time differencing scheme for system of coupled nonlinear Schr¨odinger equations, Appl. Math. Comput. 228(2014) 271-291. [30] H. P. Bhatt, A. Q. M. Khaliq, The locally extrapolated exponential time differencing LOD scheme for multidimensional reactiondiffusion systems, J. Comput. Appl. Math. 285 (2015) 256278. [31] M. Hochburck, C. Lubich, On Krylov subspace approximations to the matrix exponential operator, SIAM J. Numer. Anal. 34 (1997) 1911-1925. [32] H. Tal-Ezer, On restart and error estimation for Krylov approximation of w = f (A)v, SIAM J. Sci. Comput. 29 (2007) 2426-2441. [33] M. Caliari, A. Ostermann, Implementation of exponential Rosenbrock type integrals, App. Numer. Math. 59 (2009) 568-581. 34

1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60 61 62 63 64 65

[34] Y. A. Suhov, A spectral method for the time evolution in parabolic problems, J. Sci. Comput. 29 (2006) 201-217. [35] S. P. Norsett, A. Wolfbrandt, Attainable order of rational approximations to the exponential function with only real poles, BIT 17 (1977) 200-208. [36] C. Moler, C. Van Loan, Nineteen dubious ways to compute the exponential of a matrix, twentyfive years later, SIAM Rev. 45 (2003) 3-49. [37] A. Q. M. Khaliq, E. H. Twizell, D. A. Voss, On the parallel algorithms for semidiscretized parabolic partial differential equations based on subdiagonal Pad´e approximations, Num. Meth. Par. Diff. Eqs. 9 (1993), 107-116. [38] G. Beylkin, J. M. Keiser, L. Vozovoi, A new class of time discretization schemes for the solution of nonlinear PDEs, J. Comp. Phys. 147 (1998) 362-387. [39] Q. Du, W. Zhu, Analysis and applications of the exponential time differencing schemes and their contour integration modifications, BIT Numer. Math. 45 (2005) 307-328. [40] B. Fornberg, T. A. Driscoll, A fast spectral for nonlinear wave equations with linear dispresion, J. Comp. Phys. 155 (1999) 456-467. [41] W. L. Wood, An exact solution for Burgers’ equation, Commun. Numer. Methods Eng. 22 (7) (2006) 797-798. [42] D. Kaya, An explicit solution of coupled viscous Burgers’ equation by the decomposition method, IJMMS 27 (11) (2001) 675-680. [43] A. A. Soliman, The modified extended tanh-function method for solving Burgers-type equations, Physica A 361 (2006) 394-404. [44] S. Boscarino, R. Burger, P. Mulet, G. Russo and L.M. Villada, Linearly implicit IMEX RungeKutta methods for a class of degenerate convection-diffusion problems, SIAM J. Sci. Comput. 37 (2015) B305-B331.

35