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

Authors: Don Liu; Qin Chen; Yifan Wang
Journal: Journal of Computational Engineering
 (2013)

---

## Abstract

A system of coupled nonlinear partial differential equations with convective and dispersive terms was modified from Boussinesq-type 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.


---

## Body

## 1.  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 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 time-dependent Maxwell equation with spurious high-frequency oscillatory components of the numerical solution being suppressed by a spatial filter. Spotz and Carey extended high-order 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.


## 2. Simulation Approaches

### 2.1. Modified Boussinesq Equations

In terms of the velocityu and water depth η, Boussinesq-type equations in one horizontal dimension could be modified into below equations which are similar to [37]:
(1)∂η∂t+∂(uη)∂x=f(u,η,x,t),∂u∂t+u∂u∂x-μ3∂3u∂x2∂t+∂η∂x=g(u,η,x,t),
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 u and one for η; and boundary conditions, two for u 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 w:
(2)∂η∂t=F,∂w∂t=G,
where
(3)F=f(u,η,x,t)-∂(uη)∂x=f(u,η,x,t)-η∂u∂x-u∂η∂x,(4)w=u-μ3∂2u∂x2,G=g(u,η,x,t)-u∂u∂x-∂η∂x.
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 f and g in these equations with μ=3, the exact solutions to (1) are available. In this paper, we present the special forms of f and g, and the corresponding exact solutions as given below:
(5)f=(2x+1)e2x-2t-ex-t,g=(x2+x)e2x-2t+3ex-t,
subject to the initial conditions
(6)η(x,0)=ex,u(x,0)=xex,
and the boundary conditions
(7)η(1,t)=e1-t,u(1,t)=e1-t,∂u(2,t)∂x=3e2-t;
the exact solutions to (1) are
(8)η=ex-t,u=xex-t,1≤x≤2,t≥0.


### 2.2. Solution Ideas and Procedures

The major steps in finding numerical solution are elaborated below. These steps involveCompact 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 initialu(0) and η(0), we use Compact Scheme I to compute ∂u(0)/∂x and ∂η(0)/∂x, then by (3) and (4) to acquire F(0) and G(0); and we use Compact Scheme II to compute ∂2u(0)/∂x2 from u(0), then by (4) to obtain w(0). With F(0), G(0), η(0), and w(0) being known, by the first order Adams-Bashforth method, η(1) and w(1) could be computed. Next, we use Compact Scheme III to compute u(1). This concludes the computation in time level one.Based onu(1) and η(1), we use Compact Scheme I to compute∂u(1)/∂xand∂η(1)/∂x, then by (3) and (4) to acquire F(1) and G(1). With F(1), G(1), F(0), G(0), η(1), and w(1) being known, by the second order Adams-Bashforth method, η(2) and w(2) could be computed. Next, we use Compact Scheme III to compute u(2). This finishes the computation in time level two.Based onu(2) and η(2), repeat the same procedures and use the third order Adams-Bashforth method to compute η(3) and u(3). Continue the computations and use the third order Adams-Bashforth method from now on to compute u and η at the next time level.


### 2.3.  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≤x≤2, into N equal intervals with N+1 nodes, with the nodal index starting from 0 to the last one N.


#### 2.3.1.Compact Scheme I: Solving for ∂u/∂x,∂η/∂x with Known u,η

Denote eitheru or η with Y. The first derivative of Y relative to x at node xi is Yi′=∂Yi/∂x. The uniform grid size is h. Based on the known values of u and η at the current time level (e.g., the initial values), ∂u/∂x and ∂η/∂x could be computed accurately with compact methods. For the interior nodes, the scheme of the sixth order accuracy is
(9)12Yi-1′+36Yi′+12Yi+1′=h-1(-Yi-2-28Yi-1+28Yi+1+Yi+2)+o(h6).
For the two left end nodes, use the sixth order scheme with i=1,2, respectively:
(10)-60Yi-1′-300Yi′=h-1(197Yi-1+25Yi-300Yi+1+100Yi+2-25Yi+3+3Yi+4)+o(h6).
For the two right end nodes, when solving for ∂η/∂x, use the sixth order scheme below with i=N,N-1, respectively; when solving for ∂u/∂x, use the same scheme with i=N-1and the Neumann boundary condition:
(11)-60Yi′-300Yi+1′=h-1(-197Yi-25Yi-1+300Yi-2-100Yi-3+25Yi-4-3Yi-5)+o(h6).
Equations (9), (10), and (11) form a tridiagonal system of dimension N+1 by N+1 in terms of the nodal values of the first derivative Y′. 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].


#### 2.3.2.Compact Scheme II: Solving for ∂2u/∂x2 with Known u

Denoteu with Y, Yi′′=∂2Yi/∂x2. Based on the known values of u at the current time level (e.g., the initial time), ∂2u/∂x2 could be computed with high order compact methods. For the interior nodes, the scheme of the sixth order accuracy is
(12)8Yi-1′′+44Yi′′+8Yi+1′′=h-2(3Yi-2+48Yi-1-102Yi+48Yi+1+3Yi+2)+o(h6).
For the two left end nodes, use the following sixth order scheme with i=1,2, respectively, with the Dirichlet boundary condition Y0=e1-t at x=1 included:
(13)αYi-1′′+βYi′′=h-2(bYi-1+a0Yi+a1Yi+1+a2Yi+2+a3Yi+3+a4Yi+4+a5Yi+5)+o(h6),
where α=2823(1/3), β=32340, b=37350(19/27), a0=-75537, a1=36767(1/2), a2=4762(16/27), a3=-4620, a4=1463, and a5=(11/54)-187. For the right end node, use the fifth order mixed scheme as below with the Neumann boundary condition YN′=3e2-t at x=2 included:
(14)72YN-1′′+18YN-2′′=-30h-1YN′+h-2(127YN-216YN-1+81YN-2+8YN-3)+o(h5).
For the node next to the right end node, use the sixth order scheme below with i=N and the same coefficients as in (13):
(15)βYi-1′′+αYi′′=h-2(a5Yi-6+a4Yi-5+a3Yi-4+a2Yi-3+a1Yi-2+a0Yi-1+bYi)+o(h6).
Equations (12), (13), (14), and (15) form a quasitridiagonal system of dimension N+1 by N+1 in terms of the nodal values of the second derivative of u, 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.


#### 2.3.3.Compact Scheme III: Solving for u,∂2u/∂x2 with Known w

After the time integration with the Adams-Bashforth method,w at the next time level becomes known. Now we compute the new nodal values of u and ∂2u/∂x2 at the next time level with w, based on (4):
(16)u-∂2u∂x2=w.
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:
(17)∑i‍ai∂2ui∂x2=1h2∑k‍bkuk+o(h6).
Therefore, the resulting system is of dimension 2N+2 by 2N+2 in terms of the nodal values of ∂2u/∂x2 and u. In other words, there are two equations for each node. Denote u with Y, at the node i, ui=Yi, ∂2ui/∂x2=∂2Yi/∂x2=Yi′′. For the interior nodes, we choose (17) to be (12). Then, the system of (16) and (17) at the interior nodes becomes
(18)-Yi′′+Yi=wi,(19)8Yi-1′′+44Yi′′+8Yi+1′′=h-2(3Yi-2+48Yi-1-102Yi+48Yi+1+3Yi+2)+o(h6).
For the node next to the very end node on each side, (18) with i=1 for the left and i=N-1 for the right is implemented. The second equation at the same node is the sixth order scheme as below, with i=3 for the left one and i=N-1 for the right one:
(20)Yi-3′′-24Yi-2′′-194Yi-1′′-24Yi′′+Yi+1′′=h-2(-240Yi-2+480Yi-1-240Yi)+o(h6).
For the left end node at x=1, the Dirichlet boundary condition for Y is implemented. The second equation at the same node is (13) with i=1 and the Dirichlet boundary condition appended. For the right end node at x=2, (18) is implemented with i=N. The second equation at the same node is the fifth order mixed scheme, that is, (14) with i=N and the Neumann boundary condition YN′=3e2-t appended.We specially arrangedYi′′ and Yi alternately, that is, Y0′′, Y0, Y1′′, Y1,…, to acquire a banded linear system of dimension 2N+2 by 2N+2 in terms of the nodal values of ∂2u/∂x2 and u. 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.


## 2.1. Modified Boussinesq Equations

In terms of the velocityu and water depth η, Boussinesq-type equations in one horizontal dimension could be modified into below equations which are similar to [37]:
(1)∂η∂t+∂(uη)∂x=f(u,η,x,t),∂u∂t+u∂u∂x-μ3∂3u∂x2∂t+∂η∂x=g(u,η,x,t),
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 u and one for η; and boundary conditions, two for u 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 w:
(2)∂η∂t=F,∂w∂t=G,
where
(3)F=f(u,η,x,t)-∂(uη)∂x=f(u,η,x,t)-η∂u∂x-u∂η∂x,(4)w=u-μ3∂2u∂x2,G=g(u,η,x,t)-u∂u∂x-∂η∂x.
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 f and g in these equations with μ=3, the exact solutions to (1) are available. In this paper, we present the special forms of f and g, and the corresponding exact solutions as given below:
(5)f=(2x+1)e2x-2t-ex-t,g=(x2+x)e2x-2t+3ex-t,
subject to the initial conditions
(6)η(x,0)=ex,u(x,0)=xex,
and the boundary conditions
(7)η(1,t)=e1-t,u(1,t)=e1-t,∂u(2,t)∂x=3e2-t;
the exact solutions to (1) are
(8)η=ex-t,u=xex-t,1≤x≤2,t≥0.


## 2.2. Solution Ideas and Procedures

The major steps in finding numerical solution are elaborated below. These steps involveCompact 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 initialu(0) and η(0), we use Compact Scheme I to compute ∂u(0)/∂x and ∂η(0)/∂x, then by (3) and (4) to acquire F(0) and G(0); and we use Compact Scheme II to compute ∂2u(0)/∂x2 from u(0), then by (4) to obtain w(0). With F(0), G(0), η(0), and w(0) being known, by the first order Adams-Bashforth method, η(1) and w(1) could be computed. Next, we use Compact Scheme III to compute u(1). This concludes the computation in time level one.Based onu(1) and η(1), we use Compact Scheme I to compute∂u(1)/∂xand∂η(1)/∂x, then by (3) and (4) to acquire F(1) and G(1). With F(1), G(1), F(0), G(0), η(1), and w(1) being known, by the second order Adams-Bashforth method, η(2) and w(2) could be computed. Next, we use Compact Scheme III to compute u(2). This finishes the computation in time level two.Based onu(2) and η(2), repeat the same procedures and use the third order Adams-Bashforth method to compute η(3) and u(3). Continue the computations and use the third order Adams-Bashforth method from now on to compute u and η at the next time level.


## 2.3.  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≤x≤2, into N equal intervals with N+1 nodes, with the nodal index starting from 0 to the last one N.


### 2.3.1.Compact Scheme I: Solving for ∂u/∂x,∂η/∂x with Known u,η

Denote eitheru or η with Y. The first derivative of Y relative to x at node xi is Yi′=∂Yi/∂x. The uniform grid size is h. Based on the known values of u and η at the current time level (e.g., the initial values), ∂u/∂x and ∂η/∂x could be computed accurately with compact methods. For the interior nodes, the scheme of the sixth order accuracy is
(9)12Yi-1′+36Yi′+12Yi+1′=h-1(-Yi-2-28Yi-1+28Yi+1+Yi+2)+o(h6).
For the two left end nodes, use the sixth order scheme with i=1,2, respectively:
(10)-60Yi-1′-300Yi′=h-1(197Yi-1+25Yi-300Yi+1+100Yi+2-25Yi+3+3Yi+4)+o(h6).
For the two right end nodes, when solving for ∂η/∂x, use the sixth order scheme below with i=N,N-1, respectively; when solving for ∂u/∂x, use the same scheme with i=N-1and the Neumann boundary condition:
(11)-60Yi′-300Yi+1′=h-1(-197Yi-25Yi-1+300Yi-2-100Yi-3+25Yi-4-3Yi-5)+o(h6).
Equations (9), (10), and (11) form a tridiagonal system of dimension N+1 by N+1 in terms of the nodal values of the first derivative Y′. 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].


### 2.3.2.Compact Scheme II: Solving for ∂2u/∂x2 with Known u

Denoteu with Y, Yi′′=∂2Yi/∂x2. Based on the known values of u at the current time level (e.g., the initial time), ∂2u/∂x2 could be computed with high order compact methods. For the interior nodes, the scheme of the sixth order accuracy is
(12)8Yi-1′′+44Yi′′+8Yi+1′′=h-2(3Yi-2+48Yi-1-102Yi+48Yi+1+3Yi+2)+o(h6).
For the two left end nodes, use the following sixth order scheme with i=1,2, respectively, with the Dirichlet boundary condition Y0=e1-t at x=1 included:
(13)αYi-1′′+βYi′′=h-2(bYi-1+a0Yi+a1Yi+1+a2Yi+2+a3Yi+3+a4Yi+4+a5Yi+5)+o(h6),
where α=2823(1/3), β=32340, b=37350(19/27), a0=-75537, a1=36767(1/2), a2=4762(16/27), a3=-4620, a4=1463, and a5=(11/54)-187. For the right end node, use the fifth order mixed scheme as below with the Neumann boundary condition YN′=3e2-t at x=2 included:
(14)72YN-1′′+18YN-2′′=-30h-1YN′+h-2(127YN-216YN-1+81YN-2+8YN-3)+o(h5).
For the node next to the right end node, use the sixth order scheme below with i=N and the same coefficients as in (13):
(15)βYi-1′′+αYi′′=h-2(a5Yi-6+a4Yi-5+a3Yi-4+a2Yi-3+a1Yi-2+a0Yi-1+bYi)+o(h6).
Equations (12), (13), (14), and (15) form a quasitridiagonal system of dimension N+1 by N+1 in terms of the nodal values of the second derivative of u, 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.


### 2.3.3.Compact Scheme III: Solving for u,∂2u/∂x2 with Known w

After the time integration with the Adams-Bashforth method,w at the next time level becomes known. Now we compute the new nodal values of u and ∂2u/∂x2 at the next time level with w, based on (4):
(16)u-∂2u∂x2=w.
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:
(17)∑i‍ai∂2ui∂x2=1h2∑k‍bkuk+o(h6).
Therefore, the resulting system is of dimension 2N+2 by 2N+2 in terms of the nodal values of ∂2u/∂x2 and u. In other words, there are two equations for each node. Denote u with Y, at the node i, ui=Yi, ∂2ui/∂x2=∂2Yi/∂x2=Yi′′. For the interior nodes, we choose (17) to be (12). Then, the system of (16) and (17) at the interior nodes becomes
(18)-Yi′′+Yi=wi,(19)8Yi-1′′+44Yi′′+8Yi+1′′=h-2(3Yi-2+48Yi-1-102Yi+48Yi+1+3Yi+2)+o(h6).
For the node next to the very end node on each side, (18) with i=1 for the left and i=N-1 for the right is implemented. The second equation at the same node is the sixth order scheme as below, with i=3 for the left one and i=N-1 for the right one:
(20)Yi-3′′-24Yi-2′′-194Yi-1′′-24Yi′′+Yi+1′′=h-2(-240Yi-2+480Yi-1-240Yi)+o(h6).
For the left end node at x=1, the Dirichlet boundary condition for Y is implemented. The second equation at the same node is (13) with i=1 and the Dirichlet boundary condition appended. For the right end node at x=2, (18) is implemented with i=N. The second equation at the same node is the fifth order mixed scheme, that is, (14) with i=N and the Neumann boundary condition YN′=3e2-t appended.We specially arrangedYi′′ and Yi alternately, that is, Y0′′, Y0, Y1′′, Y1,…, to acquire a banded linear system of dimension 2N+2 by 2N+2 in terms of the nodal values of ∂2u/∂x2 and u. 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.


## 2.3.1.Compact Scheme I: Solving for ∂u/∂x,∂η/∂x with Known u,η

Denote eitheru or η with Y. The first derivative of Y relative to x at node xi is Yi′=∂Yi/∂x. The uniform grid size is h. Based on the known values of u and η at the current time level (e.g., the initial values), ∂u/∂x and ∂η/∂x could be computed accurately with compact methods. For the interior nodes, the scheme of the sixth order accuracy is
(9)12Yi-1′+36Yi′+12Yi+1′=h-1(-Yi-2-28Yi-1+28Yi+1+Yi+2)+o(h6).
For the two left end nodes, use the sixth order scheme with i=1,2, respectively:
(10)-60Yi-1′-300Yi′=h-1(197Yi-1+25Yi-300Yi+1+100Yi+2-25Yi+3+3Yi+4)+o(h6).
For the two right end nodes, when solving for ∂η/∂x, use the sixth order scheme below with i=N,N-1, respectively; when solving for ∂u/∂x, use the same scheme with i=N-1and the Neumann boundary condition:
(11)-60Yi′-300Yi+1′=h-1(-197Yi-25Yi-1+300Yi-2-100Yi-3+25Yi-4-3Yi-5)+o(h6).
Equations (9), (10), and (11) form a tridiagonal system of dimension N+1 by N+1 in terms of the nodal values of the first derivative Y′. 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].


## 2.3.2.Compact Scheme II: Solving for ∂2u/∂x2 with Known u

Denoteu with Y, Yi′′=∂2Yi/∂x2. Based on the known values of u at the current time level (e.g., the initial time), ∂2u/∂x2 could be computed with high order compact methods. For the interior nodes, the scheme of the sixth order accuracy is
(12)8Yi-1′′+44Yi′′+8Yi+1′′=h-2(3Yi-2+48Yi-1-102Yi+48Yi+1+3Yi+2)+o(h6).
For the two left end nodes, use the following sixth order scheme with i=1,2, respectively, with the Dirichlet boundary condition Y0=e1-t at x=1 included:
(13)αYi-1′′+βYi′′=h-2(bYi-1+a0Yi+a1Yi+1+a2Yi+2+a3Yi+3+a4Yi+4+a5Yi+5)+o(h6),
where α=2823(1/3), β=32340, b=37350(19/27), a0=-75537, a1=36767(1/2), a2=4762(16/27), a3=-4620, a4=1463, and a5=(11/54)-187. For the right end node, use the fifth order mixed scheme as below with the Neumann boundary condition YN′=3e2-t at x=2 included:
(14)72YN-1′′+18YN-2′′=-30h-1YN′+h-2(127YN-216YN-1+81YN-2+8YN-3)+o(h5).
For the node next to the right end node, use the sixth order scheme below with i=N and the same coefficients as in (13):
(15)βYi-1′′+αYi′′=h-2(a5Yi-6+a4Yi-5+a3Yi-4+a2Yi-3+a1Yi-2+a0Yi-1+bYi)+o(h6).
Equations (12), (13), (14), and (15) form a quasitridiagonal system of dimension N+1 by N+1 in terms of the nodal values of the second derivative of u, 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.


## 2.3.3.Compact Scheme III: Solving for u,∂2u/∂x2 with Known w

After the time integration with the Adams-Bashforth method,w at the next time level becomes known. Now we compute the new nodal values of u and ∂2u/∂x2 at the next time level with w, based on (4):
(16)u-∂2u∂x2=w.
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:
(17)∑i‍ai∂2ui∂x2=1h2∑k‍bkuk+o(h6).
Therefore, the resulting system is of dimension 2N+2 by 2N+2 in terms of the nodal values of ∂2u/∂x2 and u. In other words, there are two equations for each node. Denote u with Y, at the node i, ui=Yi, ∂2ui/∂x2=∂2Yi/∂x2=Yi′′. For the interior nodes, we choose (17) to be (12). Then, the system of (16) and (17) at the interior nodes becomes
(18)-Yi′′+Yi=wi,(19)8Yi-1′′+44Yi′′+8Yi+1′′=h-2(3Yi-2+48Yi-1-102Yi+48Yi+1+3Yi+2)+o(h6).
For the node next to the very end node on each side, (18) with i=1 for the left and i=N-1 for the right is implemented. The second equation at the same node is the sixth order scheme as below, with i=3 for the left one and i=N-1 for the right one:
(20)Yi-3′′-24Yi-2′′-194Yi-1′′-24Yi′′+Yi+1′′=h-2(-240Yi-2+480Yi-1-240Yi)+o(h6).
For the left end node at x=1, the Dirichlet boundary condition for Y is implemented. The second equation at the same node is (13) with i=1 and the Dirichlet boundary condition appended. For the right end node at x=2, (18) is implemented with i=N. The second equation at the same node is the fifth order mixed scheme, that is, (14) with i=N and the Neumann boundary condition YN′=3e2-t appended.We specially arrangedYi′′ and Yi alternately, that is, Y0′′, Y0, Y1′′, Y1,…, to acquire a banded linear system of dimension 2N+2 by 2N+2 in terms of the nodal values of ∂2u/∂x2 and u. 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.


## 3.  Simulation Results

### 3.1. Approximate and Exact Solutions

We compare our approximate (numerical) solutions with exact solutions foruand its derivatives in Figure 1(a), η and its derivative in Figure 1(b), at the dimensionless time T=0.1. The time step size is dt=0.0001 and the space step size is dx=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.Comparison of the approximate and exact solutions foru (a) and η (b).
(a)
(b)


### 3.2. Convergence Rates of Separate Compact Schemes in Space

#### 3.2.1.Compact Scheme I

The performance ofCompact Scheme I is examined here. Since the exact solutions are available, we compute the Euclidean (L2) norms of the relative errors for the first derivative of u and η at all nodes. Figure 2(a) presents, in the log-log scale, the L2 norm of the relative errors versus the number of nodes used in discretization. When the spatial resolution goes down to 0.01, that is, N=100, the Euclidean norms of the relative errors drop to the machine zero in double precision. The data points demonstrate the sixth order convergence rate for both ∂u/∂x and ∂η/dx.Convergence rates forCompact Scheme I (a) and Compact Scheme II (b) in space.
(a)
(b)


#### 3.2.2.Compact Scheme II

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


#### 3.2.3.Compact Scheme III

The performance ofCompact Scheme III is demonstrated below. Taking advantage of the known exact solution, we compute the Euclidean norms of the relative errors for both u and ∂2u/∂x2 at all nodes. Figure 3(a) plots the L2 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 i=N in order to adopt the Neumann boundary condition, the convergence rates for u and ∂2u/∂x2 are slightly different but about sixth order.(a) Convergence rates forCompact Scheme III. (b) Convergence rates for combined schemes in space with dt=10-4.
(a)
(b)


### 3.3. Convergence Rates of Combined Schemes in Space

#### 3.3.1. Time Stepdt=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 dt=10-4. The lowest convergence rate of 4.84 is observed for ∂2u/∂x2 because the variable u carries a second derivative and the boundary condition is not explicitly given in terms of ∂2u/∂x2. The rest of the variables have the sixth order convergence rate.


#### 3.3.2. Time Stepdt=10-5 and dt=10-6

Similar to Figures3(b) and 4(a) with dt=10-5 and Figure 4(b) with dt=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 ∂2u/∂x2 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 N in figures). This is because a smaller time step size produces less error during time integration. However, as the spatial resolution (N) increases, the rounding error increases accordingly, and eventually overtakes the truncation error. Therefore, we do not see the errors keep decreasing as N increases.Convergence rates for combined schemes in space with (a)dt=10-5 and (b) dt=10-6.
(a)
(b)


### 3.4.  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 todx=0.02 (N=50) and vary the size of time step to plot the L2 norms of the relative errors in all variables. Figure 5(a) presents the rate of convergence versus the size of time step for all variables. A second order time convergence is observed when dt is less than 10-5. Since the spatial step size is still large, at dt=10-5, the time convergence rate is impaired, especially for uxx. After we reducedxto 0.0125, the second order time convergence rate is restored upto dt=10-5 for all variables. When dt is further diminished to 10-6, the time convergence rate becomes less sensitive to dt again, but the L2 norm of the relative error of uxx is still better than the case with a larger spatial step dx=0.02.Convergence rates for combined schemes in time:L2 norm of the relative errors versus the size of time step at fixed spatial resolution: (a) dx=0.02 and (b) dx=0.0125.
(a)
(b)


### 3.5.  The Effect of Long Time Integration

To examine the effect of long time integration, we fix the spatial resolution toN=50, that is, dx=0.02, and use dt=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 u and η over 168,000 steps in time integration with dt=10-5. It is observed that during this long time integration, the L2 norms of the relative errors in u and η are bounded below 10-9, although slightly growing over time. Figure 6(b) shows the evolution of convergence rates for the derivative of u and η in the same run. The L2 norms of the relative errors in ux, uxx, and ηx are below 10-8. The L2 norm of the relative error in uxx is almost one order higher than those in ux and ηx. The reason is that, although all derivatives of u and η are treated as the state variables just like u and η themselves and solved in Compact Schemes I, II, III simultaneously, no Dirichlet boundary condition for uxx is specified and uxx is solved after the time integration; therefore, the relative error in uxx is superposed by the time integration error.Long time integration over 168,000 time steps foru and η (a) and for the derivatives of u and η (b), the time step size being 10-5.
(a)
(b)Next we reduce the time step size todt=10-6. Figure 7(a) shows the evolution of convergence rates for u and η after a time integration period over 220,000 steps. It is observed that during this time period, the L2 norms of the relative errors in u 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 u and η in the same run. The L2 norms of the relative errors in ux and ηx are below 10-10 while the relative error in uxx 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.Long time integration over 220,000 time steps foru and η (a) and for the derivatives of u and η (b), the time step size being 10-6.
(a)
(b)


## 3.1. Approximate and Exact Solutions

We compare our approximate (numerical) solutions with exact solutions foruand its derivatives in Figure 1(a), η and its derivative in Figure 1(b), at the dimensionless time T=0.1. The time step size is dt=0.0001 and the space step size is dx=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.Comparison of the approximate and exact solutions foru (a) and η (b).
(a)
(b)


## 3.2. Convergence Rates of Separate Compact Schemes in Space

### 3.2.1.Compact Scheme I

The performance ofCompact Scheme I is examined here. Since the exact solutions are available, we compute the Euclidean (L2) norms of the relative errors for the first derivative of u and η at all nodes. Figure 2(a) presents, in the log-log scale, the L2 norm of the relative errors versus the number of nodes used in discretization. When the spatial resolution goes down to 0.01, that is, N=100, the Euclidean norms of the relative errors drop to the machine zero in double precision. The data points demonstrate the sixth order convergence rate for both ∂u/∂x and ∂η/dx.Convergence rates forCompact Scheme I (a) and Compact Scheme II (b) in space.
(a)
(b)


### 3.2.2.Compact Scheme II

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


### 3.2.3.Compact Scheme III

The performance ofCompact Scheme III is demonstrated below. Taking advantage of the known exact solution, we compute the Euclidean norms of the relative errors for both u and ∂2u/∂x2 at all nodes. Figure 3(a) plots the L2 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 i=N in order to adopt the Neumann boundary condition, the convergence rates for u and ∂2u/∂x2 are slightly different but about sixth order.(a) Convergence rates forCompact Scheme III. (b) Convergence rates for combined schemes in space with dt=10-4.
(a)
(b)


## 3.2.1.Compact Scheme I

The performance ofCompact Scheme I is examined here. Since the exact solutions are available, we compute the Euclidean (L2) norms of the relative errors for the first derivative of u and η at all nodes. Figure 2(a) presents, in the log-log scale, the L2 norm of the relative errors versus the number of nodes used in discretization. When the spatial resolution goes down to 0.01, that is, N=100, the Euclidean norms of the relative errors drop to the machine zero in double precision. The data points demonstrate the sixth order convergence rate for both ∂u/∂x and ∂η/dx.Convergence rates forCompact Scheme I (a) and Compact Scheme II (b) in space.
(a)
(b)


## 3.2.2.Compact Scheme II

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


## 3.2.3.Compact Scheme III

The performance ofCompact Scheme III is demonstrated below. Taking advantage of the known exact solution, we compute the Euclidean norms of the relative errors for both u and ∂2u/∂x2 at all nodes. Figure 3(a) plots the L2 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 i=N in order to adopt the Neumann boundary condition, the convergence rates for u and ∂2u/∂x2 are slightly different but about sixth order.(a) Convergence rates forCompact Scheme III. (b) Convergence rates for combined schemes in space with dt=10-4.
(a)
(b)


## 3.3. Convergence Rates of Combined Schemes in Space

### 3.3.1. Time Stepdt=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 dt=10-4. The lowest convergence rate of 4.84 is observed for ∂2u/∂x2 because the variable u carries a second derivative and the boundary condition is not explicitly given in terms of ∂2u/∂x2. The rest of the variables have the sixth order convergence rate.


### 3.3.2. Time Stepdt=10-5 and dt=10-6

Similar to Figures3(b) and 4(a) with dt=10-5 and Figure 4(b) with dt=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 ∂2u/∂x2 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 N in figures). This is because a smaller time step size produces less error during time integration. However, as the spatial resolution (N) increases, the rounding error increases accordingly, and eventually overtakes the truncation error. Therefore, we do not see the errors keep decreasing as N increases.Convergence rates for combined schemes in space with (a)dt=10-5 and (b) dt=10-6.
(a)
(b)


## 3.3.1. Time Stepdt=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 dt=10-4. The lowest convergence rate of 4.84 is observed for ∂2u/∂x2 because the variable u carries a second derivative and the boundary condition is not explicitly given in terms of ∂2u/∂x2. The rest of the variables have the sixth order convergence rate.


## 3.3.2. Time Stepdt=10-5 and dt=10-6

Similar to Figures3(b) and 4(a) with dt=10-5 and Figure 4(b) with dt=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 ∂2u/∂x2 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 N in figures). This is because a smaller time step size produces less error during time integration. However, as the spatial resolution (N) increases, the rounding error increases accordingly, and eventually overtakes the truncation error. Therefore, we do not see the errors keep decreasing as N increases.Convergence rates for combined schemes in space with (a)dt=10-5 and (b) dt=10-6.
(a)
(b)


## 3.4.  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 todx=0.02 (N=50) and vary the size of time step to plot the L2 norms of the relative errors in all variables. Figure 5(a) presents the rate of convergence versus the size of time step for all variables. A second order time convergence is observed when dt is less than 10-5. Since the spatial step size is still large, at dt=10-5, the time convergence rate is impaired, especially for uxx. After we reducedxto 0.0125, the second order time convergence rate is restored upto dt=10-5 for all variables. When dt is further diminished to 10-6, the time convergence rate becomes less sensitive to dt again, but the L2 norm of the relative error of uxx is still better than the case with a larger spatial step dx=0.02.Convergence rates for combined schemes in time:L2 norm of the relative errors versus the size of time step at fixed spatial resolution: (a) dx=0.02 and (b) dx=0.0125.
(a)
(b)


## 3.5.  The Effect of Long Time Integration

To examine the effect of long time integration, we fix the spatial resolution toN=50, that is, dx=0.02, and use dt=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 u and η over 168,000 steps in time integration with dt=10-5. It is observed that during this long time integration, the L2 norms of the relative errors in u and η are bounded below 10-9, although slightly growing over time. Figure 6(b) shows the evolution of convergence rates for the derivative of u and η in the same run. The L2 norms of the relative errors in ux, uxx, and ηx are below 10-8. The L2 norm of the relative error in uxx is almost one order higher than those in ux and ηx. The reason is that, although all derivatives of u and η are treated as the state variables just like u and η themselves and solved in Compact Schemes I, II, III simultaneously, no Dirichlet boundary condition for uxx is specified and uxx is solved after the time integration; therefore, the relative error in uxx is superposed by the time integration error.Long time integration over 168,000 time steps foru and η (a) and for the derivatives of u and η (b), the time step size being 10-5.
(a)
(b)Next we reduce the time step size todt=10-6. Figure 7(a) shows the evolution of convergence rates for u and η after a time integration period over 220,000 steps. It is observed that during this time period, the L2 norms of the relative errors in u 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 u and η in the same run. The L2 norms of the relative errors in ux and ηx are below 10-10 while the relative error in uxx 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.Long time integration over 220,000 time steps foru and η (a) and for the derivatives of u and η (b), the time step size being 10-6.
(a)
(b)


## 4.  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,u, η, and their first and second derivatives. However, all state variables are not solved at once, because solving u, ux, uxx, η, and ηx simultaneously in one linear system requires to triple the size of the linear system for solving u and η together. To reduce memory usage for potentially very large systems, we presolve the first and the second derivatives of u in tridiagonal systems. After the time integration, we solve for u and uxx. In other words, u, η, 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.



---

*Source: 432192-2013-12-04.xml*