A Sixth Order Accuracy Solution to a System of Nonlinear Differential Equations with Coupled Compact Method

A system of coupled nonlinear partial differential equations with convective and dispersive terms was modified from Boussinesqtype equations. Through a special formulation, a system of nonlinear partial differential equations was solved alternately and explicitly in time without linearizing the nonlinearity.Coupled compact schemes of sixth order accuracy in space were developed to obtain numerical solutions.Within couple compact schemes, variables and their first and second derivatives were solved altogether. The sixth order accuracy in space is achieved with a memory-saving arrangement of state variables so that the linear system is banded instead of blocked. This facilitates solving very large systems. The efficiency, simplicity, and accuracy make this coupled compact method viable as variational and weighted residual methods. Results were compared with exact solutions which were obtained via devised forcing terms. Error analyses were carried out, and the sixth order convergence in space and second order convergence in time were demonstrated. Long time integration was also studied to show stability and error convergence rates.


Introduction
Compact difference methods, essentially the implicit versions of finite different methods, are superior to the explicit versions in achieving high order accuracy. High order compact methods directly approximate all derivatives with high accuracy without any variational operation or projection. They are very effective for complex and strongly nonlinear partial differential equations, especially when variational methods are inconvenient to implement. They have smaller stencils than explicit finite difference methods while maintaining high order accuracy.
Much development has been achieved in finite difference methods. A large number of papers related to various finite difference algorithms have been published over the past decades. Kreiss [1] was among the first who pioneered the research on compact implicit methods. Hirsh [2] developed and applied a 4th order compact scheme to Burgers' equation in fluid mechanics. Adam [3] developed 4th order compact schemes for parabolic equations on uniform grids. Hoffman [4] studied the truncation errors of the centered finite difference method on both uniform and nonuniform grids and discussed the accuracy of these schemes after applying grid transformations. Dennis and Hudson [5] studied a 4th order compact scheme to approximate Navier-Stokes type operators. Rai and Moin [6] presented several finite difference solutions for an incompressible turbulence channel flow. They also compared their solutions of high order upwind schemes with solutions obtained from spectral methods. Their work demonstrated that finite difference methods can be a good approach for the direct numerical simulation of turbulence. Lele [7] demonstrated the spectral-like resolution of the classical Padé schemes and the flexibility of compact finite difference methods. He developed compact difference methods for mid-point interpolation and filtering with high order formal accuracy. Yu et al. [8] solved an unsteady 2D Euler equation with a 6th order compact difference scheme in space and a 4th order Runge-Kutta scheme in time, which provide a good example of high order (both spatial and temporal) finite difference methods. Mahesh [9] presented high order compact schemes with both the first and second order derivatives evaluated simultaneously. Dai and Nassar [10,11] developed high order compact difference schemes 2 Journal of Computational Engineering for high order derivatives in time and space in solving a 3D heat transport problem at the microscale. Sun and Zhang [12] proposed a new high order finite difference discretizaion strategy based on the Richardson extrapolation technique and on an operator interpolation scheme for solving one-and two-dimensional convection-diffusion equations.
Compact methods are widely used in applications in engineering and physics. Ekaterinaris [13] presents fourth order, factorized, upwind compact schemes for time accurate solutions of hyperbolic equation systems. Li et al. [14] obtained solutions of Navier-Stokes equation expressed in the stream function-vorticity form, with a fourth order compact scheme for a driven cavity problem. Shang [15] solved the timedependent Maxwell equation with spurious high-frequency oscillatory components of the numerical solution being suppressed by a spatial filter. Spotz and Carey extended highorder compact schemes to time-dependent problems [16]. El-Zoheiry [17] devised a three-level iterative scheme based on the compact implicit method for solving the improved Boussinesq equation. Zhang et al. [18] derived a fourth order compact scheme with multigrid mesh refinement for convection diffusion problems. Central and upwind compact methods were analyzed by Sengupta et al. [19], and schemes with improved boundary closure and an explicit high-order upwind scheme were proposed. Shen et al. [20] demonstrated a new way of constructing high accuracy shock-capturing compact schemes. Jocksch et al. [21] developed a fifth order compact upwind-biased finite difference scheme for compressible Couette flow. Mohamad et al. [22] presented a third order compact upwind scheme for flows containing discontinuities. Anthonio and Hall [23] simulated free surface flows with a high-order compact finite difference method and a fourth order Runge-Kutta method for temporal discretization. Cruz and Chen [24] modeled porous beds with high order time-marching fourth order finite difference schemes. A dual-compact scheme, with the increased dispersive accuracy for the convective terms, was proposed by Chiu and Sheu [25]. A new family of high-order compact upwind difference schemes with good spectral resolution was introduced by Zhou et al. [26]. Sanyasiraju and Mishra [27] developed sixth order compact schemes for convection diffusion equations with constant coefficients. Dehghan and Mohebbi [28] presented unconditionally stable fourth order accuracy in both space and time solution to unsteady convection diffusion equation. Shah et al. [29] solved the Navier-Stokes equations in which the convective terms were approximated with a third order upwind compact scheme using flux differencing splitting. Xie et al. [30] presented the numerical solution of the one dimensional nonlinear Schrödinger equation. Tonelli and Petti [31] used a finite volume technique to solve the advective part of the Boussinesq equations and applied a finite difference method to discretize dispersive and source terms. The time integration was performed with the fourth order Adams-Bashforth-Moulton method.
Irregular geometries could hinder the application of compact methods. To facilitate nonuniform discretization, one way is to use a transformation so that compact methods can be implemented in the transformed space where the nominal rate of accuracy could be achieved. For example, zeros of Chebyshev polynomials could be adopted as the grid points for computation [32,33]. Gamet et al. [34] developed a 4th order compact scheme for approximating the first order derivative with a full inclusion of metrics directly on a nonuniform mesh with variable grid spacing. Using this method, a direct numerical simulation of a compressible turbulent channel flow was conducted. Vasilyev [35] presented high order finite difference simulation of turbulent flows on nonuniform meshes with good conservation properties. Hermanns and Hernández [36] developed high order finite difference methods on nonuniform grids. Usually, the actual accuracy is determined by the effect of the transformation. Another option for nonuniform grids is to use full inclusion of metrics (FIM), such as Liu et al. [33] where fourth and sixth order compact methods were used for parabolic partial differential equations encountered in geophysical modeling. Generally speaking, FIM could be difficult for deriving high accuracy schemes.
In the past two decades, Boussinesq-type equations, derived from depth-integration of conservation laws for mass and momentum, have been widely used to simulate nonlinear water waves in coastal and ocean engineering. One of the challenges of solving this class of partial differential equations is the numerical treatment of the dispersive terms described by the mixed third order derivative in space (twice) and time.
One objective of this paper is to develop 6th order coupled compact schemes for modified Boussinesq equations with analytical exact solution for comparison. These schemes demonstrate high order accuracy and fast convergence rates in both space and time. The idea in this paper could be used for solving other systems of coupled nonlinear equations.
This paper is outlined as below. In simulation approaches, the modified Boussinesq equations and their exact solutions corresponding to specially devised forcing functions are given. Then the solution ideas and procedures are presented, followed by schemes of sixth order accuracy developed for these equations. In simulation results, the solutions of all coupled variables and the convergence rates in space and time are demonstrated. In addition, the effect of long time integration is investigated to show numerical stability over time. Finally, conclusion summarises the novelty and uniqueness of the schemes developed in this paper, and future work is planned.

Modified Boussinesq Equations.
In terms of the velocity and water depth , Boussinesq-type equations in one horizontal dimension could be modified into below equations which are similar to [37]: where is a dimensionless parameter measuring the wave dispersion at a given still-water depth. The above equations are closed with two initial conditions, one for and one for ; and boundary conditions, two for and one for . This is a system of nonlinearly coupled partial differential equations with strong convection and dispersivity. To solve this system of equations, by the method of lines [38], (1) could be formulated into a system of time-dependent equations in terms of and a new variable : Now we could take care of the system of equations in time and space separately. To march (2) in time, a variety of methods are available [32,38,39]. Since implicit methods are used in space and the system is strongly nonlinear, we select explicit methods for time evolution. In this paper, up to third order Adams-Bashforth methods are used to integrate in time. During the initial time steps, lower order Adams-Bashforth methods are used instead. In order to compare results and investigate the rate of convergence, exact solutions to (1) are needed. With some specific forms of and in these equations with = 3, the exact solutions to (1) are available. In this paper, we present the special forms of and , and the corresponding exact solutions as given below: subject to the initial conditions and the boundary conditions the exact solutions to (1) are

Solution Ideas and Procedures.
The major steps in finding numerical solution are elaborated below. These steps involve Compact Schemes I, II, and III, which are described in the next subsection. We use superscripts to denote the time levels of variables. Based on the initial (0) and (0) , we use Compact Scheme I to compute (0) / and (0) / , then by (3) and (4) to acquire (0) and (0) ; and we use Compact Scheme II to compute 2 (0) / 2 from (0) , then by (4) to obtain (0) . With (0) , (0) , (0) , and (0) being known, by the first order Adams-Bashforth method, (1) and (1) could be computed. Next, we use Compact Scheme III to compute (1) . This concludes the computation in time level one.

Coupled Compact Schemes of High
Accuracy. Now we focus on the spatial discretization and the specific high order compact schemes mentioned earlier. We partition the domain of interest, 1 ≤ ≤ 2, into equal intervals with +1 nodes, with the nodal index starting from 0 to the last one .

Compact Scheme I: Solving for
Denote either or with . The first derivative of relative to at node is = / . The uniform grid size is ℎ. Based on the known values of and at the current time level (e.g., the initial values), / and / could be computed accurately with compact methods. For the interior nodes, the scheme of the sixth order accuracy is For the two left end nodes, use the sixth order scheme with = 1, 2, respectively: For the two right end nodes, when solving for / , use the sixth order scheme below with = , − 1, respectively; when solving for / , use the same scheme with = − 1 and the Neumann boundary condition: Equations (9), (10), and (11) form a tridiagonal system of dimension + 1 by + 1 in terms of the nodal values of the first derivative . The coefficient matrix of the linear system is stored as three vectors in order to save storage. The linear system is solved with the Thomas Algorithm [40].

Compact Scheme II: Solving for 2 / 2 with Known .
Denote with , = 2 / 2 . Based on the known values of at the current time level (e.g., the initial time), 2 / 2 could be computed with high order compact methods. For the interior nodes, the scheme of the sixth order accuracy is For the two left end nodes, use the following sixth order scheme with = 1,2, respectively, with the Dirichlet boundary condition 0 = 1− at = 1 included: For the node next to the right end node, use the sixth order scheme below with = and the same coefficients as in (13): Equations (12), (13), (14), and (15) form a quasitridiagonal system of dimension + 1 by + 1 in terms of the nodal values of the second derivative of , since only one nonzero element in the last row lies outside the tridiagonal lines. The coefficient matrix of the linear system could still be stored as three vectors. The linear system is solved by the Thomas Algorithm with a minor modification.

Compact Scheme
Since two variables are solved for at the same time, another equation is needed here. A compact scheme of the sixth order accuracy as below could be the supplementary equation: Therefore, the resulting system is of dimension 2 + 2 by 2 +2 in terms of the nodal values of 2 / 2 and . In other words, there are two equations for each node. Denote with , at the node , = , 2 / 2 = 2 / 2 = . For the interior nodes, we choose (17) to be (12). Then, the system of (16) and (17) at the interior nodes becomes For the node next to the very end node on each side, (18) with = 1 for the left and = −1 for the right is implemented. The second equation at the same node is the sixth order scheme as below, with = 3 for the left one and = − 1 for the right one: For the left end node at = 1, the Dirichlet boundary condition for is implemented. The second equation at the same node is (13) with = 1 and the Dirichlet boundary condition appended. For the right end node at = 2, (18) is implemented with = . The second equation at the same node is the fifth order mixed scheme, that is, (14) with = and the Neumann boundary condition = 3 2− appended. We specially arranged and alternately, that is, 0 , 0 , 1 , 1 ,. . ., to acquire a banded linear system of dimension 2 + 2 by 2 + 2 in terms of the nodal values of 2 / 2 and . Without such deliberate alternating arrangement, the coefficient matrix would be blocked with sparsity. This system with a small bandwidth could be efficiently solved via the banded Gaussian elimination with scaled partial pivoting.

Approximate and Exact Solutions.
We compare our approximate (numerical) solutions with exact solutions for and its derivatives in Figure 1(a), and its derivative in Figure 1(b), at the dimensionless time = 0.1. The time step size is = 0.0001 and the space step size is = 0.05. No visual difference could be seen from these figures between the approximate and exact solutions. To investigate errors, we need to check convergence rates in both space and time.    drop to the machine zero in double precision. The data points demonstrate the sixth order convergence rate for both / and / .

Compact Scheme II.
To see the performance of Compact Scheme II, we use the known exact solution to compute the Euclidean norm of the relative errors for 2 / 2 at all nodes is. In the log-log scale, Figure 2(b) shows the 2 norm of the relative error versus the number of nodes used in space. Before the spatial step reduces to 0.01, the 2 norm of the relative errors of 2 / 2 drops to the machine zero. This figure clearly demonstrates the sixth order convergence rate.

Compact Scheme III. The performance of Compact
Scheme III is demonstrated below. Taking advantage of the known exact solution, we compute the Euclidean norms of the relative errors for both and 2 / 2 at all nodes. Figure 3(a) plots the 2 norms of the relative errors versus the number of nodes used in discretization. The linear system generated from the Compact Scheme III is no longer tridiagonal, but it is still banded. Because a fifth order scheme is used at the right end node = in order to adopt the Neumann boundary condition, the convergence rates for and 2 / 2 are slightly different but about sixth order.

Time
Step = 10 −4 . Now we examine the convergence rates for all three compact schemes together over 1,000 time steps after the initial conditions are imposed. During this period of time, Compact Schemes I, II, and III are used one after another to solve these modified Boussinesq equations. Figure 3(b) presents the overall spatial convergence rates (indicated in the legend, resp.) for all variables after 1,000 time steps, which was set to be = 10 −4 . The lowest convergence rate of 4.84 is observed for 2 / 2 because the variable carries a second derivative and the boundary condition is not explicitly given in terms of 2 / 2 . The rest of the variables have the sixth order convergence rate.

Time
Step = 10 −5 and = 10 −6 . Similar to Figures  3(b) and 4(a) with = 10 −5 and Figure 4(b) with = 10 −6 present the spatial convergence rates (indicated in the legend) involving Compact Schemes I, II, and III for all variables after 1,000 time steps. The lowest convergence rate about 4.9 is expected for 2 / 2 for the same reason as before. For the rest variables, the sixth order convergence rates are clearly demonstrated. From Figures 3(b), 4(a), and 4(b), we observe that as the time step size reduces, high order convergence rates are maintained for higher spatial resolution (larger in figures). This is because a smaller time step size produces less error during time integration. However, as the spatial resolution ( ) increases, the rounding error increases accordingly, and eventually overtakes the truncation error. Therefore, we do not see the errors keep decreasing as increases.

Convergence Rate in Time Integration.
After we examined the spatial convergence rates, now we demonstrate the convergence rate in time integration. To this end, we fix the spatial step size to = 0.02 ( = 50) and vary the size of time step to plot the 2 norms of the relative errors in all variables. Figure 5  the size of time step for all variables. A second order time convergence is observed when is less than 10 −5 . Since the spatial step size is still large, at = 10 −5 , the time convergence rate is impaired, especially for . After we reduce to 0.0125, the second order time convergence rate is restored upto = 10 −5 for all variables. When is further diminished to 10 −6 , the time convergence rate becomes less sensitive to again, but the 2 norm of the relative error of is still better than the case with a larger spatial step = 0.02.

The Effect of Long Time Integration.
To examine the effect of long time integration, we fix the spatial resolution to = 50, that is, = 0.02, and use = 10 −5 to plot the overall convergence rates for all variables with Compact Schemes I, II, III implemented one after another during computations. Since all three schemes are of sixth order accuracy, the magnitude of discretization error is (0.02) 6 = 6.4 × 10 −11 , the time step size has to be small enough so that the time integration error will be of the same order as in space. Figure 6(a) shows the evolution of convergence rates for and over 168,000 steps in time integration with = 10 −5 . It is observed that during this long time integration, the 2 norms of the relative errors in and are bounded below 10 −9 , although slightly growing over time. Figure 6(b) shows the evolution of convergence rates for the derivative of and in the same run. The 2 norms of the relative errors in , , and are below 10 −8 . The 2 norm of the relative error in is almost one order higher than those in and . The reason is that, although all derivatives of and are treated as the state variables just like and themselves and solved in Compact Schemes I, II, III simultaneously, no Dirichlet boundary condition for is specified and is solved after the time integration; therefore, the relative error in is superposed by the time integration error. Next we reduce the time step size to = 10 −6 . Figure 7(a) shows the evolution of convergence rates for and after a time integration period over 220,000 steps. It is observed that during this time period, the 2 norms of the relative errors in and are well bounded below 10 −10 , one order lower than in Figure 6(a). Figure 7(b) presents the evolution of convergence rates for the derivatives of and in the same run. The 2 norms of the relative errors in and are below 10 −10 while the relative error in is stable around 10 −9 . There is no indication of error growth since the time step size is small enough to stabilize the growth of the time integration error.

Conclusion
In this paper, coupled compact methods are used to solve nonlinear differential equations with third order mixed derivatives in space and time. We have developed sixth order compact schemes for numerical solutions to modified Boussinesq-type of equations, a coupled system of nonlinear differential equations with strong convective, and dispersive terms. The novelty of this work is summarized as below. First, the state variables include the original unknowns, , , and their first and second derivatives. However, all state variables are not solved at once, because solving , , , , and simultaneously in one linear system requires to triple the size of the linear system for solving and together. To reduce memory usage for potentially very large systems, we presolve the first and the second derivatives of in Journal of Computational Engineering tridiagonal systems. After the time integration, we solve for and . In other words, , , and their derivatives are solved in two consecutive time levels. Second, the coupled partial differential equations are solved in space and time alternately. This is especially effective for strong nonlinear terms. No explicit linearization is needed. Third, for the nonlinear partial differential equations with mixed third order derivative, the implicit differencing is applied in space and explicit differencing is used in time integration. This unique treatment avoids the linearization of any nonlinear terms. This could be very useful for strongly nonlinear and time-dependent partial differential equations.
With the aid of convergence analyses in space and time, by fixing either the time step or the spatial step, the expected accuracy of these schemes is demonstrated. Long time integrations were performed to show the stability of these schemes.
In the subsequent work, we investigate the effectiveness of the proposed schemes for solutions with oscillatory signs, such as the velocity of water waves. Upwind treatment will be adopted for strong hyperbolic systems.