Efficient and accurate finite difference schemes for solving one-dimensional Burgers’ equation

In this paper, two efficient fourth-order compact finite difference algorithms have been developed to solve the one-dimensional Burgers' equation: u_t+u u_x=ε u_xx. The methods are based on the Hopf–Cole transformation, Richardson's extrapolation, and multilevel grids. In both methods, we first transform the original nonlinear Burgers' equation into a linear heat equation: w_t=ε w_xx using the Hopf–Cole transformation, which is given as u=−2ε (w_x/w). In the first method, the resulted heat equation is solved by the second-order accurate Crank–Nicholson algorithm while w_x is approximated by central finite difference, which is also second-order accurate. Richardson's extrapolation technique is then applied in both time and space to obtain fourth-order accuracy. In the second method, to reduce the cancellation error in approximating w_x, we derive the heat equation satisfied by w_x, which is then solved by the Crank–Nicholson algorithm. The original Dirichlet boundary condition is transformed into the Robin boundary condition, which is also approximated using second-order central finite difference. Finally, Richardson's extrapolation and multilevel grid techniques are applied in both time and space to obtain fourth-order accuracy. To study the efficiency, accuracy and robustness, we solved two numerical examples and the results are compared with those of two other higher-order methods proposed in W. Liao [An implicit fourth-order compact finite difference scheme for one-dimensional Burgers' equation, Appl. Math. Comput. 206(2) (2008), pp. 755–764] and I.A. Hassanien, A.A. Salama, and H.A. Hosham [Fourth-order finite difference method for solving Burgers' equation, Appl. Math. Comput. 170 (2005), pp. 781–800].

Keywords: Burgers' equation; Multilevel grids; Richardson's extrapolation; finite difference scheme; higher-order compact algorithm; 65M06; 65M55; 65L12; 65N06

In this paper, we consider the one-dimensional Burgers' equation:

Graph

with an initial condition

Graph

and boundary conditions

Graph

where and Re is the Reynolds number characterizing the strength of viscosity.

Burgers' equation (1) has been widely used in modelling various phenomena in science and engineering with a broad range of applications, including gas dynamics, traffic flow, and wave propagation in acoustics and hydrodynamics, etc. It also serves as a prototype in the study of nonlinear effects in fluid mechanics and in the general partial differential equation theory. The first attempt to analytically solve the one-dimensional Burgers' equation was by Bateman [2], who derived the steady-state solution for Equation (1), which had been used by Burger [5] to model turbulence.

It is well known that Burgers' equation in its original form of Equation (1) can be solved analytically only for some special initial and boundary conditions. Hopf [11] and Cole [7] have independently shown that for more general initial and boundary conditions the one-dimensional Burgers' equation can be transformed to a linear homogeneous heat equation that can be solved analytically, thus the analytical solution to the original Burgers' equation can be expressed in the form of an infinite Fourier series. However, efficient numerical methods are still required for cases in which the initial conditions are only available at discrete points, or are discontinuous or not smooth. For these cases, the Fourier series solutions may converge very slowly.

There has been intense research in the past several decades in efficient numerical methods for solving Burgers' equation in its original form (1). Various numerical algorithms have been developed to solve application problems involving Burgers' equation, including finite difference, finite-element, boundary element, and spectral methods. For instance, in [18], Ozis and Aslan developed a semi-approximate approach for solving the one-dimensional Burgers' equation with high Reynolds numbers. In [10], Hon and Mao applied the multiquadratic as a spatial approximation scheme and a low order explicit finite difference approximation for time discretization to solve Burgers' equation. The authors reported that the major numerical error is from the time integration instead of the multiquaric spatial approximation. Hassanien et al.[9] developed a two-level three-point finite difference method, which is fourth-order accurate in space and second-order accurate in time. Based on our numerical tests which are presented later in Section 4, this method can be improved to fourth order in time by applying Richardson's extrapolation. However, this method is an iterative method, thus at each time step, several iterations are needed to ensure the convergence. Nevertheless, the Hassanien method is unconditionally stable and very robust when it is used to solve Burgers' equation with large Reynolds number. In [12], Huang and Abduwali proposed a modified local Crank–Nicolson method to solve one- and two-dimensional Burgers' equation. It has been shown that the proposed method is efficient and unconditionally stable. For more details of the recent methods for solving the Burgers' equation, the readers are referred to [1][6][8][9][10][15][18][19][21].

There are two main difficulties in numerically solving Burgers' equation in its original form (1) using the finite difference method. The first difficulty is that to avoid non-physical oscillations, the convection term in Burgers' equation requires one-sided upwind finite difference schemes, which are more complicated to implement than the central schemes when the solution changes sign in the solution domain. Further, the one-sided upwind schemes typically have a lower order of accuracy compared to the central difference schemes on the same finite difference stencil. The second difficulty is that the discretization of Burgers' equation will lead to a system of nonlinear equations if an implicit time integration method is used. Since the implicit methods are usually more stable, it is often necessary to use them for difficult problems for which the explicit or semi-explicit methods are not practical. For a more detailed comparison of various early works on the finite difference method readers are referred to [3].

To overcome these difficulties, the Hopf–Cole transformation has been used by many researchers to reduce the nonlinear Burgers' equation to a linear heat equation by eliminating the nonlinear convection term. Kutluay et al.[15] have developed explicit and exact-explicit finite difference methods based on the Hopf–Cole transformation. In [14], Kadalbjoo and Awasthi reduced the nonlinear Burgers' equation to a heat equation and then solved the heat equation by the Crank–Nicholson algorithm. An efficient implicit fourth-order finite difference method based on Hopf–Cole transformation and the Padé approximation was proposed by Liao in [16]. In [23], a compact fourth-order finite difference method also based on Hopf–Cole transformation and Padé approximation of the central difference operator was developed. A detailed proof of the stability and convergence of the new method was provided as well.

In this paper, two fourth-order compact finite difference schemes based on the Hopf–Cole transformation, Richardson's extrapolation, and multilevel grids are proposed to solve the one-dimensional Burgers' equation (1)–(3). The rest of the paper is organized as follows. In Section 2, we briefly describe the Hopf–Cole transformation and derive the linear heat equation from the nonlinear Burgers' equation. The two new fourth-order methods are described in Section 3, followed by numerical experiments in Section 4. Finally concluding remarks are presented in Section 5.

Assume that u(x, t) is the solution to the one-dimensional nonlinear Burgers' equation given in Equation (1), then the following Hopf–Cole transformation [7][11]

Graph

connects u(x, t) to w(x, t), which is the solution to the following linear heat equation:

Graph

with initial condition

Graph

and boundary conditions

Graph

In [16], a fourth-order finite difference scheme based on the above Hopf–Cole transformation and the Padé approximation was presented to solve the one-dimensional Burgers' equation. The integration in Equation (6) was approximated by the fourth-order Simpson's rule and the boundary conditions (7) and (8) were approximated by the method developed in [17]. The solution w(x, t) to the heat equation (5) and the derivative wx were calculated in [16] by fourth-order compact schemes based on the Padé approximation. The resulting algorithm is fourth-order accurate, compact, and unconditionally stable. It works well for problems with low Reynolds numbers. In this paper, we follow similar procedures in [16] to approximate the initial and boundary conditions (6)–(8) and develop two new fourth-order solution algorithms for solving the one-dimensioanl Burgers' equation with large Reynolds numbers.

The two new methods are based on the Crank–Nicholson algorithm and Richardson's extrapolation in both time and space. The difference between the two methods is how wx is calculated. The first method uses the standard second-order central finite difference to calculate wx after the solution w has been calculated, while the second method solves a linear heat equation satisfied by wx.

To simplify the presentation, we divide the solution domain [0, 1]×[0, T] into a uniform M×N grid with step sizes h=1/(M−1) and Δ t=T/N. The numerical solutions of w and wx at grid point xi and time level tn are denoted as and , respectively. The second-order central finite difference operator is denoted as , which is defined as .

It is well known that the truncation error of the Crank–Nicholson algorithm is second-order in both space and time. By using Taylor expansion, we can see that the third-order terms Δ t3 and h3 are absent in the truncation error. This allows the use of the Richardson's extrapolation to achieve fourth-order accuracy in both space and time. Richardson's extrapolation method is very simple and flexible. It can be applied to either the calculated numerical solution [17] or the residuals produced by the numerical solution [4]. In the following two subsections, we introduce the two fourth-order finite difference methods, both are based on the Crank–Nicholson algorithm and Richardson's extrapolation.

Applying the well-known Crank–Nicholson algorithm to heat Equation (5), we obtain

Graph

To derive the boundary conditions for w at x=0 and x=1, we first approximate Equations (7) and (8) using second-order central finite difference:

Graph

We then solve the boundary conditions on both sides as the following

Graph

where and are the numerical solutions of w(x, t) at the two virtual points x0=−h and xM+1=1+h, respectively.

Upon solving Equation (9), we obtain that is second-order accurate in both space and time. Using the standard second-order central finite difference to approximate , we have

Graph

Substituting and into Equation (4), we obtain numerical solution to Burgers' equation at final time level T

Graph

The truncation error in algorithm (9) is second-order accurate in time and space. Moreover, it can be shown by Taylor series expansion that the third-order terms h3 and Δ t3 are absent in the truncation error. We can also verify that the boundary conditions given in Equations (12) and (13) are second-order accurate with the third-order term h3 absent in the truncation error. Similarly, the central finite difference approximation given in Equation (14) is second-order accurate with the third-order term absent in the truncation error. Based on the Lax Theorem that a consistent and stable two-step finite difference algorithm for a linear PDE is convergent with a global error that is of the same order as that of the local truncation error [13][20], we have the following error estimations for the numerical solutions of w and wx at the grid point xi and time T:

Graph

where w(xi, T) and wx(xi, T) are the exact analytic solutions, while and are the numerical solutions.

As a result, the numerical solution calculated from Equation (15) is second-order accurate in time and space with the third-order terms absent in the truncation error, which is very important when Richardson's extrapolation and multigrid techniques are applied to obtain fourth-order accuracy. We summarize the results in the following theorem:

Assume that u(xi, T) is the exact solution of Equation(1)at xiand time T, whileis the numerical solution of u(xi, T) calculated from the algorithm given in Equations(9), (12)–(15), then the numerical solution is second-order accurate in both time and space, with the following error bound

Graph

where and ẽ 2 are constants depending on T, x i and ε.

From the Hopf–Cole transformation (4) and the algorithm in Equation (15), we can derive the error bound as follows:

Graph

Using the error estimations in Equations (16) and (17), we have

Graph

where , . Assume and |wx|<ξ, we have

Graph

Obviously, when and , thus there exist δ>0, such that when h<δ and , . Therefore, when max, we have the following error estimation:

Graph

where , , , and κ= max(ζ, ξ) are constants depending on xi, T and ε.   ▪

Both w(xi, T) and are non-zero, based on the Hopf–Cole transformation given in Equation (4), thus the coefficients e1, e2, ẽ1 and ẽ2 are well defined. The feature that the third-order terms h3 and Δ t3 are absent in the truncation error of u(xi, T) allows the use of Richardson's extrapolation and multigrid techniques to obtain fourth-order accuracy in both time and space.

Denote the numerical solution obtained by using Δ t and h as . We then repeat the solution procedure using a finer mesh with step sizes (h/2) and (Δ t/2), and denote the numerical solution on the finer mesh as . Applying Richardson's extrapolation to the resulting numerical solutions, we have

Graph

which is fourth-order accurate in both space and time, with the leading error term given as . Note that the linear combination in Equation (21) is normally performed on the coarser grid with step sizes h and Δ t. If the extrapolated solution is desired on the finer grid, the solution can be first interpolated to the finer mesh with second-order accuracy and then combined with using Equation (21) on the finer mesh. It is worth pointing out that the third-order terms h3 and Δ t3 should also be absent in the interpolation error.

Although the first new method is fourth-order accurate, we noted that it involves the central finite difference (14) when wx is calculated. When a central finite difference is directly used to calculate the first derivative, , there could be significant cancellation errors and loss of significant digits if values of wi+1 and wi−1 are close to each other. For example, if and both have eight significant digits, then will only have four significant digits. To reduce the cancellation error in calculating wx, we use a different way to approximate wx in this method. Setting v(x, t)=wx, and differentiating both sides of the heat equation (5) with respect to x, we derive the following heat equation satisfied by v(x, t):

Graph

We first solve the two heat equations (5) and (22) by the Crank–Nicholson algorithm:

Graph

The initial condition for (22) is obtained from (4) as

Graph

where w(xi, 0) is the initial condition for w at xi, which is available through (6). The boundary condition for v can be derived from the boundary conditions given in Equations (7) and (8). For instance, at x=0, we have

Graph

while at the right boundary x=1, we have

Graph

Since the values and used in the calculation of boundary conditions (26) and (27) are from algorithm (23), there will be truncation errors of order O(h2), with the third-order term absent. This, however, does not change the order of the truncation error for all grid points in any norm. As a result, the calculated solutions for w and v are second-order accurate in space and time, and again the third-order terms h3 and Δ t3 are absent in the global errors.

The solution u to the original Burgers' equation is then obtained by substituting the calculated w and wx into Equation (4). Repeating this process on two grids, one with Δ t and h, and one with Δ t/2 and h/2, we obtain and , which are second-order accurate in time and space with no Δ t3 and h3 terms in the truncation error. Richardson's extrapolation given in Equation (21) can then be used to obtain u that is fourth-order accurate in space and time.

Although the extrapolations in space and time using two grids require as much as five times the amount of computation as solving the equations on the coarse grid alone, the resulting higher-order accuracy allows the use of much larger steps in space and time. Thus, the overall computational cost is reduced significantly.

Since the Crank–Nicholson algorithm is unconditionally stable, the two new methods are also unconditionally stable since they are based on the Crank–Nicholson algorithm. The numerical results in the next section also show that there is no restriction on the time step.

Two numerical examples are presented to verify the order of convergence of the new methods, and to compare accuracy and efficiency of the new methods discussed in the previous section, denoted as new method 1 and new method 2, with those of the fourth-order method proposed in [16], denoted as the Padé method, and the fourth-order method proposed in [9], denoted as the Hassanien method. For clarity of presentation, we give below a brief description of the Padé method and the Hassanien method before the two numerical examples.

Padé method: The Padé method uses the Hopf–Cole transformation to turn Equation (1) into a linear heat equation (5) and then achieves fourth-order accuracy in space by using the following Padé approximation to the second-order derivative:

Graph

where . Replacing in Equation (9) by Equation (28) and ignore the truncation error term, we have

Graph

Applying the operator on both sides of Equation (29) and moving all unknown terms to the left-hand side of the equation, we have

Graph

Equation (30) involves only the standard 3-point finite difference stencil , but is fourth-order accurate in space and second-order accurate in time. Richardson's extrapolation is then used to improve the results to fourth-order accurate in time.

Hassanien method: The Hassanien method solves the nonlinear Burgers equation (1) directly. The nonlinear convective term is linearized by lagging the multiplier u one iteration behind and discretizing fully implicitly. The algorithm can be written as

Graph

where m represents iterations within each time step to deal with the nonlinear convection term, a1, a0, and a−1 are coefficients depending on , , and , and b1, b0, and b−1 are coefficients depending on , , and . All these coefficients are calculated by maximizing the order of truncation errors in approximating Equation (1) by Equation (31). In each iteration, Equation (31) forms a system of tridiagonal equations that can be solved efficiently by the Thomas algorithm.

Note that the Hassanien method is fourth-order accurate in space but only second-order accurate in time, therefore the comparison of accuracy of this method with other methods in this section is made only in spatial dimension. However, when the overall efficiencies of the four methods are compared, we apply Richardson's extrapolation in time to the original Hassanien method to make it fourth-order accurate in time and space. The revised method is still denoted as Hassanien method.

The first example is used to demonstrate that while all four methods are fourth-order accurate, the new method 2 appears to be more accurate than the Padé method, which in turn appears to be more accurate than the new method 1 and Hassanien method. The second example is used to show that the efficiency and robustness of the new methods when they are used to solve Burgers' equation with high Reynolds numbers.

We first solve Burgers' equation (1) for which the exact solution is available [22] as

Graph

where σ>1 is a parameter. The initial and boundary conditions are taken from this exact solution.

We first verify that the two new methods are fourth-order accurate in both space and time. As shown in Table 1, the maximum errors for both methods are reduced by a factor of about 24=16 when h and Δ t are reduced by a factor of 2, which apparently indicates the fourth-order accuracy in space and time.

Table 1. Fourth-order convergence of the two new methods.

We then compare the two new methods with the other two methods in terms of accuracy in the spatial dimension. To reduce the effect from temporal dimension, the time step size Δ t is set extremely small () so that discretization error in time is negligible. We initially set h=1/10, then reduced it by a factor of 2 in each refinement. The results in Table 2 confirm that all of the four methods are fourth-order accurate in the spatial dimension, since the errors are reduced by a factor of about 16 in each refinement. It is clear that the new method 2 that approximates wx by solving the heat equation (22) is the most accurate while the Hassanien method is the least accurate, though all of them are fourth-order accurate. When the same step size h is used, the maximal error of the new method 2 is about 1/5 of the maximum error of the Padé method, about 1/25 of the maximum error of the new method 1, and about 1/30 of the maximum error of the Hassanien method.

Table 2. Comparison of numerical errors in space of the four fourth-order methods for Example 1.

It is well known that the error of a fourth-order method can be expressed as , where the constant c can be estimated as . For the four methods in Table 2, the constant c is , , and , respectively. Obviously, the smaller the constant is, the more accurate the method is.

To compare the overall efficiency of these methods, we adjust h and Δ t for each method, solve this example with such that the maximal error for each method is around , and then record the CPU time (in seconds) for each method. The results are displayed in Table 3, which shows that the Padé method proposed in [16] is the most efficient method, which is about 30 times faster than the least efficient method, the Hassanien method. This is reasonable, as we can see that the Hassanien method is an iterative method, and in each step, the coefficient matrix needs to be evaluated. Between these two methods are the two new methods, with the new method 1 being about 50% faster than the new method 2.

Table 3. Comparison of efficiency of the four fourth-order methods when used to solve Example 1 with ε=0.1, σ=2.

In this example, we solve Burgers' equation (1) with the following boundary and initial conditions:

Graph

Note that the exact solution to this problem is available as the following infinite series:

Graph

where the coefficients are defined as

Graph

To verify the accuracy of the numerical methods, the infinite series solution (34) needs to be evaluated to high accuracy. In this example, the number of terms N is chosen such that the coefficient cN is less than .

It is well known that one of the difficulties in solving Burgers' equation with large Reynolds numbers (small ε) is that the solution may become almost discontinuous with a sharp front after some time, even if the initial data are smooth. A robust and accurate numerical algorithm should be able to capture the sharp front and produces numerical solution that exhibits the correct physical behaviour of the exact solution.

We plot the solutions at final time Tf=1.0 for and 0.0001 by the two new methods, the Padé method, the Hassanien method, and the analytical formula (34) in Figure 1. It is clear that when ε is large (the Reynolds number is small), all four numerical solutions are very accurate, smooth, and hardly distinguishable from the exact solutions. There is no sharp front in the solutions. When , it is clear that a sharp front develops around x=0.9 in all four numerical solutions and the exact solution, as can be seen from Figure 1(b). All four numerical algorithms are still accurate and capture the sharp front very well.

Graph: Figure 1. Comparison of exact solution with the numerical solutions from the four numerical methods for Example 2, with h=Δ t=0.01 and various ε.

When ε is reduced to 0.001, the Hassanien method begins to produce non-physical oscillations and fails to capture the sharp front of the solution, while the solution by Padé method blows up (not plotted in the figure). However, it seems that the Hassanien method is pretty accurate before the sharp front. It suggests that this method is suitable for solving problem without sharp front. Clearly, both the new method 1 and new method 2 are still able to capture the sharp front reasonably well. It seems that the new method 1 is more accurate and robust for problems with smaller ε (large Reynolds number).

If we further reduce ε to 0.0001, the solution by the Padé method blows up (not plotted in the figure). The solution from new method 2, although not accurate, remains oscillation-free and captures the sharp front at the location where a sharp front is expected. One interesting observation is that the Hassanien method, which is very accurate before the sharp front, failed to capture the sharp front, instead it produces a smooth solution where a sharp front is expected. Finally, the new method 1 still produces a oscillation-free and reasonably good numerical solution to the problem.

It is well known that grid refinement is an effective way to resolve the oscillation of numerical method. In the following test, we fix , and solve this example using various methods with different h and Δ t. The results are displayed in Figure 2. It shows that when the grid is very fine, all numerical methods are accurate and can capture the sharp front. As can be seen from Figure 2(a), all numerical solutions are almost identical to the exact solution. When h and Δ t increase to 0.002, we can see in Figure 2(b), that all numerical methods are still accurate, though it is clear that the new method 2 is the least accurate method. If we further increase h and Δ t to 0.0025, we noted in Figure 2(c), that both the new method 2 and the Hassanien method produced solutions that are visibly different from the exact solution, however all four numerical methods still capture the sharp front as expected. When , the Padé method begins to produce oscillation in the solution, while the other three methods still produce reasonable solutions, as can be seen from Figure 2(d). Note that the solutions by Padé method are not plotted in part (e) and (f) of Figure 2, as they blow up. When h and Δ t are increased to 0.01, there is noticeable difference among the solutions of the three numerical methods. As we can see from Figure 2(e), the Hassanien method produces oscillation around the location where a sharp front is expected, the new method 2, is not accurate although it captures the sharp front at the right location. It seems that the new method 1 is the most robust method among the four methods, which captures the sharp front and produces reasonably accurate numerical solutions. Finally, we tested these numerical methods with a very coarse grid, with and the results are plotted in Figure 2(f). It clearly shows that the new method 1 is the most robust method, which still produces numerical solution that is reasonably close to the exact solution, although there is a noticeable difference. On the other hand, the solutions by Hassanien method and new method 2 are not acceptable at all.

Graph: Figure 2. Comparison of exact solution with the numerical solutions from the four numerical methods for Example 2, with ε=0.001 and various stepsizes (h=Δ t).

To further study and compare the behaviours of the two new methods and the Hassanien method near the right boundary, we plot the numerical solutions for Figures 2(c) and (d) near the grid point x=0.99 in Figure 3. The detailed pictures shows that the first new method is more accurate than both the second new method and the Hassanien method, although all of the three methods are very accurate.

Graph: Figure 3. Comparison of exact solution with the numerical solutions from the four numerical methods for Example 2, with ε=0.001 near the sharp front (h=Δ t).

To further compare the four methods, the L2 errors of the numerical solutions for Example 2 are calculated and displayed in Table 4.

Table 4. L2 errors of the numerical solutions for Example 2 with ε=0.001 by the four methods using various step sizes.

In summary, it seems that when ε is small (large Reynolds number), fine grid is preferred, which of course will make the methods very time-consuming, and requires much more computing resource. When both accuracy and robustness are considered, it seems that the new method 1 is the best method among all of the four numerical methods.

In this example, we solve Burgers' equation (1) with the following boundary and initial conditions:

Graph

Note that the exact solution is available as

Graph

where the coefficients are defined as

Graph

To further compare these methods in accuracy, we solve this example for and present the results in Table 5. It is clear that all numerical methods are accurate. It seems that the new method 1 is the most accurate one, while the Hassanien method is the least accurate one. The Padé method and the new method 2 are somewhere between them Here the exact solution is obtained by evaluating the series given in Equation (38) for which the integer N is chosen such that the coefficient cN is less than 1.0e−15.

Table 5. Comparison of numerical errors(in maximum norm) of the four fourth-order methods for Example 3.

Two fourth-order compact finite difference methods extended from the fourth-order Padé method in [16] are developed in this paper. In the first new method, the heat equation (5) is solved by the standard Crank–Nicholson algorithm, and wx is approximated by the standard second-order central finite difference using the calculated numerical solution w. Richardson's extrapolation is then applied to solutions from two grids to obtain fourth-order accuracy. In the second new method, both w and wx are calculated by solving Equations (5) and (22) using the standard Crank–Nicholson algorithm. Richardson's extrapolation is then used to obtain fourth-order accuracy. Numerical examples are used to verify the fourth-order accuracy of the two new methods. The fourth-order Padé method developed in [16] and Hassanien method developed in [9] are implemented to solve these numerical examples and the efficiency, accuracy, and robustness of the four methods are compared.

The numerical results in Examples 1 and 3 have shown that when the Reynolds number is small and the solution has no sharp front, the Padé method and the two new methods are more accurate than the Hassanien method. Since the Padé method is computationally more efficient than the other methods to achieve a given accuracy, it should be the preferred method for problems with a small Reynolds number. Example 2 shows that the new method 1 and the Hassanien method are more robust than the Padé method and the new method 2 when the Reynolds number is large. Since the new method 1 is more computationally efficient than the Hassanien method, it is suggested that the new method 1 should be used for problems with a large Reynolds number.

The numerical tests show that both new methods are fourth-order accurate in time and space, however it seems that the parameter ε plays an important role in the truncation error. In the future, we plan to further investigate the accuracy of these methods, and extend these higher-order methods to multidimensional problems.

The work of the first author is supported by the Natural Sciences and Engineering Research Council of Canada (NSERC). The authors would like to thank the anonymous referees for their efforts and constructive comments on the revision of the manuscript.