A High-Order Weakly L -Stable Time Integration Scheme with an Application to Burgers’ Equation

: In this paper, we propose a 7th order weakly L -stable time integration scheme. In the process of derivation of the scheme, we use explicit backward Taylor’s polynomial approximation of sixth-order and Hermite interpolation polynomial approximation of ﬁfth order. We apply this formula in the vector form in order to solve Burger’s equation, which is a simpliﬁed form of Navier-Stokes equation. The literature survey reveals that several methods fail to capture the solutions in the presence of inconsistency and for small values of viscosity, e.g., 10 − 3 , whereas the present scheme produces highly accurate results. To check the effectiveness of the scheme, we examine it over six test problems and generate several tables and ﬁgures. All of the calculations are executed with the help of Mathematica 11.3. The stability and convergence of the scheme are also discussed.


Local Trunction Error
By applying Taylor's series expansion, we have then it follows that where t n (h) = O(h 8 ).
Thus, the scheme (29) is seventh order convergent. Accordingly, the order of the proposed method is reduced by one, but it is weakly stable, which is a great advantage.
From Figure 1, it can be seen that Ψ(s) ≮ 1 and, hence, our scheme is not A-stable. Since Ψ(s) → 0 as s → ∞, the scheme is weakly L-stable (see [35]).

Stability Region
To find the boundary of the stability region, we apply the boundary locus method (p. 64, Chapter 7, Ref. [38]). It can be easily seen that, outside of the region (see Figure 2), it is unconditionally stable.

The Numerical Scheme
Here, we consider solution space with uniform nodes expressed as Σ T i,j = {(x i , t j ) : i = 0(1)N, j = 0(1)M}. For that, we partition the interval [α 0 , α 1 ] into N (a positive integer) equal sub intervals with the spatial point x i = i∆x, i = 0(1)N, where ∆x is the spatial step.
Additionally, dividing the interval [0, T] into M equal subintervals with the temporal point t j = jτ, j = 0(1)M, where τ = T/M and M is a positive integer. Now, define ψ i (t) = ψ(x i , t) and consider equation (5) and compute the solution ψ(x i , t) for a given t and for x i on [α 0 , α 1 ]. We use (8)- (9) to deduce the following formula for computing the ω(x i , t j ) Now we approximate second order spatial derivative by fourth order central finite difference formula which is given by and convert the linearized Burgers' equation into an initial value problem in vector form. Now, we apply the above finite difference discretization on (5) with the Neumann boundary conditions we get the following equation where and Let ρ = ν d τ/24∆x 2 , then applying the time integration formula on the initial value problem (45), we obtain Now, we use the above defined Ψ j+1/6 , Ψ j+2/6 , Ψ j+3/6 , Ψ j+4/6 , Ψ j+5/6 in Equation (47) and deduce our final formula used to compute numerical solutions This method is of order O(∆x 4 ) + O(τ 7 ). By using (47), we can compute Ψ j+1 and, hence, w ij is computed at different x i 's for a given time level t j . The physical properties of the solutions are discussed later in the form of figures and tables.

Stability Analysis
Equation (47) can be written as where .
The matrix P is similar to a symmetric matrix.

Proof. Let us introduce a diagonal matrix
i.e., D is similar to a symmetric matrixD. Now, we will show that P is similar to symmetric matrix. Let but matricesL −1 1 andL 2 are symmetric and commute and therefore P is similar to a symmetric matrix P and therefore all the eigenvalues of the matrix P are real.

Lemma 2.
All of the eigenvalues of the matrix D are non-negative.

Numerical Experiment
To confirm the effectiveness of the proposed weakly L-stable scheme, we apply it to some examples and compute the numerical solutions and depict them in tables and figures. To check the accuracy of the proposed scheme, we also analyze the following type of errors (i) Mean root square error norm (L 2 ) and (ii) Maximum error norm (L ∞ ) where ω num. is numerical solution by present method and ω exact is the exact solution and j indicates solution at jth grid point.

Example 1
Consider the Equation (1) with Dirichlet BCs and IC where ν d is the coefficient of viscosity. Using the transformation Equation (1) is transformed into with IC and BCs The solution of the above initial boundary value problem is defined by Equation (11), where In this example, we depict several tables and figures and list numerical solutions and exact solutions in order to exhibit the correctness of the scheme. Additionally, we have computed the L 2 and L ∞ error defined by Equations (61) and (62), respectively. In Table 1, we take ν d = 2 and time step τ = 0.0001 with ∆x = 0.0125. It can be noticed that computed solutions are very close to analytical solutions. It can be seen that, at a specified location, the solution decreases when time passes out. Additionally, it can be seen that, at a fixed moment, the numerical solution first increases and then decreases with changing location from 0 to 1. Table 2 comprise the computed solutions and analytical solutions for ν d = 0.2, ∆x = 0.0125 and with the time step τ = 0.0001.
In Table 2, we can see that the results provided by the present scheme are better than the results in [39]. Additionally, our results are actually better than those of [40]. Table 3 represents that the current scheme gives satisfactory results for viscous coefficient ν d = 0.01 with ∆x = 0.0125 and τ = 0.01 at different times T. Figure 3 illustrates the accuracy of the present scheme at different times for ν d = 0.2 and we are not able to distinguish between the analytical solutions and computed solutions. It is known that the Fourier series solution fails to converge for ν d < 0.01 due to the slow convergence rate of infinite series (11). In Figure 4, the analytical solution shows high oscillation, while the computed solutions follow physical behavior. Figure 5 represents the physical behavior of the computed solutions for the different small values of ν d . In Figure 6, we illustrate the physical behavior of the computed results for a small value of ν d in the three-dimensional mode.

Example 2
Consider Equation (1) with Dirichlet BCs and IC where ν d is the viscous coefficient. Using the transformation we see that the Equation (11) represents the analytical solution where In Table 4, we depict the numerical results and compare them with exact solutions for τ = 0.0001, ν d = 2 with ∆x = 0.0125. L ∞ and L 2 -error indicate that the difference between the analytical solution and numerical solution is very less. For the comparison purpose, in Table 5, we take ν d = 0.2, τ = 0.0001 and ∆x = 0.0125 and notice that the results by the present method are slightly more close to the exact solutions than that given in [39] and [40]. It can be seen that, at a specified location, the solution decreases when time passes out. Additionally, it can be seen that, at a fixed moment, the numerical solution first increases and then decreases with changing location from 0 to 1. In Table 6, we take ∆x = 0.0125, τ = 0.01, and ν d = 0.01. It can be noticed that solutions that are produced by the present method and analytical solutions are very close to each other.
In Figure 7, we can see that the exact solution starts oscillating between x = 0.8 to x = 1 for a small value of ν d = 0.001 due to slow rate of convergence of infinite series, but the results obtained by this method follow the parabolic profile. Figure 8 demonstrates the accuracy of the method for ν d = 0.1 and it can be seen that the analytical solution and numerical solution are almost the same at different times throughout the domain. Figure 9 shows that the results that are projected by the present scheme follow the nature of the solutions for different small values of ν d . Figure 10 illustrates the physical nature of the solution in three dimensions.

Example 3
Here, we take the exact solution (shock-like) [41] of Equation (1) ω(x, t) = x t where t 0 = e 1 4ν d having BCs and IC In Table 7, for the comparison, we take time step τ = 0.01, spatial step ∆x = 0.0005 and small value of viscous coefficient ν d = 0.002. The numerical solutions at the different discrete points are compared with the existing results presented in [41] and also with the exact solutions. It can be noticed that the results by the present technique are slightly more close to exact solutions than that of given in [41]. For this example, discrete L ∞ and L 2 -error norms are also given and compared with the error that is given in [41]. It can be seen that the error produced by the method in [41] is very high when compared to the error provided by the current scheme. Figure 11 shows the nature of the solutions by present method for small value of ν d = 0.001.
The physical behavior of the obtained results by current scheme for ν d = 0.001, τ = 0.01 and ∆x = 0.025 is exhibited in the Figure 12 left. The absolute errors are presented in Figure 12 right and it is clear that the absolute errors are ≤10 −3 for different times. Hence, the numerical results that are obtained by the present method are acceptable.
Equation (11) represents the analytical solution of the above problem, where This example shows inconsistent IC and BCs at the boundary point x = 1. Figure 13 shows high oscillation near the boundary point x = 1 by the CN method, while the present method gives accurate and stable numerical solutions throughout the domain.
Equation (11) represents the analytical solution of the above problem, where x cos lπx dx.
This example shows inconsistent IC and BCs at both the boundary points x = 0 and x = 1. From Figure 14, it is clear that the CN method produces high oscillation near both boundary points, while the present method gives accurate and stable numerical solutions throughout the domain.

Conclusions
In the present paper, we have used explicit backward Taylor's series approximation formula of order six and Hermite interpolation polynomial of order five to derive 7th order time integration formula, which is weakly L-stable. The present method is tested over some problems and observed that the approximated results are quite satisfactory and also provides good results when compared to the existing results. It is also observed that the analytical and computed results are very close to each other for small values of viscosity. The strength of this method is that it is easy to apply and takes very little time for computation. It will be interesting to see whether we can give general n − 1th order convergent weakly L-stable Newton-Cotes formulae by using nth order convergent Newton-Cotes formula. Though the order of the weakly L-stable method is reduced by 1, it is very fruitful for a certain class of nonlinear initial boundary value problems.