Burgers’ equation is a nonlinear partial differential equations (PDEs) which was first introduced by BatemanBATEMAN , and was later treated as the turbulence of the mathematical model book:968450 ; BURGERS1948171 . Burgers’ equation is an especially important PDEs in fluid mechanics, which combines the characteristics of the first order wave equation and heat conduction equation. Burgers’ equation is used as a tool to describe the interaction between convection and diffusion. Over the decades, Burgers’ equation has a large variety of applications in the modeling of water in dynamic soil water, surface disturbances electromagnetic waves, density waves, statistics of flow problems, mixing and turbulent diffusion, cosmology and seismology Zabusky1965 ; Dubansky2017 , etc. Hopf Hopf1950 and Cole Cole1951 showed independently that for any given initial conditions the Burgers’ equation can be reduced to a linear homogeneous heat equation that can be solved analytically and the analytical solution of the original Burgers’ equation can be expressed in the form of Fourier series. Even though the analytical solution is available in the form of the Fourier series, accurate and efficient numerical schemes are still required to solve the Burgers’ equation which consists of multi-dimensional system or the complex initial condition. In such situations, the Fourier series solutions for the practical applications are very limited, which converges slowly or diverges in many cases. The analytical or numerical solutions are essential for the corresponding Burgers’ equations. Apart from the limited number of these problems, most of them do not have exact analytical solutions, so it is imperative to get a satisfactory solution of Burgers’ equation. Here, we first analyze one-dimensional coupled Burgers’ equation. Owing to the nonlinear convection term and viscous term, the coupled Burgers’ equation can be studied as a simple example of the Navier-Stokes equation.
The one-dimensional coupled nonlinear Burgers’ equation Kumar2014
subject to the initial conditions:
and boundary conditions
where is kinematic viscosity parameters of the fluid, which correspond to an inverse of Reynolds number (if , then ) . are real constants and are arbitrary constants. and are given smooth functions.
It is pervasively acknowledged that the nonlinear coupled Burgers’ equation (1) does not have precise analytic solutions. Researchers are interested in using various numerical techniques to study the properties of Burgers’ equation, because it has wide applicability in various fields of science and engineering. Due to the existence of nonlinear terms and viscosity parameters, numerical approximation of nonlinear coupled Burgers’ equation is a challenging task. In Burgers’ equation, discontinuities may appear in finite time, even if the initial condition is smooth. They give rise to the phenomenon of shock waves which have important applications in physics BREZIS199876 ; Frisch2001 . Recently, some contributions related to the time-dependent coupled viscous Burgers’ equation have been published, which analyze the theoretical and numerical aspects. Several numerical experiments on non-coupled and coupled Burgers’ equation were run to compare the accuracy of the proposed schemes with other existing methods Bhatt2016 ; Gao2016 ; Khater2008 ; Mittal2011 ; Jaradat2018 . For the sake of clarity, a brief description of the comparison method is provided below.
Bhatt et al. Bhatt2016 proposed A-stable and L-stable Fourth-order exponential time difference Runge-Kutta schemes in combination with a global fourth-order CFD scheme for the numerical solution of the coupled Burgers’ equations.
In Ref. Gao2016 , the analytical solutions of two-dimensional and three-dimensional Burgers’ equations are derived. For multi-dimensional problems, these solutions can describe the shock wave phenomenon in large Reynolds numbers (), which is very useful for testing analytical solutions of numerical methods.
In Ref. Khater2008
the author developed a Chebyshev spectral collocation method for solving approximate solutions of nonlinear PDEs. Using Chebyshev spectral collocation method, this problem is reduced to a set of ordinary differential equations (ODEs), and then solved with Runge-Kutta fourth-order method.
In Ref. Mittal2011 the author uses cubic B-spline function to construct a collocation method for numerical simulation of coupled Burgers’ equation. The time derivative term is discretized by the conventional Crank-Nicolson (C-N) scheme, while the space derivative term is discretized by the cubic B-spline method. The results obtained by finite difference cubic B-spline show that the accuracy of the solution decreases with the increase of time due to the time truncation error of the time derivative term.
Jaradat et al. Jaradat2018 establish new two-mode coupled Burgers’ equations which are introduced for the first time. The authors find the necessary conditions in which the multiple kink and multiple singular kink solutions exist and present the two-front solutions.
In this paper, we mainly discuss the numerical scheme of n-dimensional Burgers’ system. The n-dimensional Burgers’ system includes non-coupled and coupled problems, and multi-dimensional problems are for coupled Burgers’ equations.
where is the fluid velocity fields, is dimension of space, is the kinematic viscosity of the fluid and are real constants. denotes the Laplace operator, and is the Hamilton gradient operator.
The system of Eq. (4) is Burgers’ equations which includes non-coupled() and coupled() problems. When and , the system of Eq.(4) respectively becomes one-dimensional, two-dimensional and three-dimensional Burgers’ equation
Chen et al. Chen2016 show an n-dimensional Hopf-Cole transformation between the n-dimensional Burgers’ system and an n-dimensional heat equation under an irrotational condition. Motivated by this idea, the purpose of this paper is to intend to extend the Hopf-Cole transformation to linearize the n-dimensional Burgers’ equation (4); After obtaining the n-dimensional heat conduction equation, the CFD scheme with high precision and high efficiency is used to solve it.
Currently, there are many numerical methods for heat conduction equation Andrea2001 ; Marsden1994 , such as finite difference method (FDM), finite element method (FEM), finite volume method (FVM) and spectrum method, etc. The traditional FDM shows great limitations in accuracy. An important measure to improve the accuracy of the traditional FDM is to refine mesh, which in turn will increase the amount of storage and calculating speed, especially in high-dimensional cases. Therefore, it is of great theoretical significance and practical value to construct a scheme with high accuracy and excellent stability in time and space.
The CFD scheme is one of the most studied FDM at present. Experience proves that the compact scheme is much more accurate than the corresponding explicit scheme of the same order Li2008 . Over the past three decades, the methods for developing high-order CFD scheme have made great progress. Dennis et. al. proposed the fourth-order CFD scheme for convection-diffusion problems Dennis1989 , this scheme can get more accurate results with a thicker grid. Lele Lele1992 developed CFD scheme with pseudo spectral resolution on the basis of summarizing the previous work and proposed a linear sixth-order central CFD scheme, which can achieve the accuracy of the spectral method. Subsequently, many scholars constructed different schemes of CFD scheme and solved many types of partial differential equations Zhao2013 ; Lai2007 ; Nihei2003 ; Sutmann2007 , such as integro-differential equations, three-dimensional Poisson equations, the shallow water equations, and the Helmholtz equations, they all achieved better numerical results. Sengupta et. al. developed a class of upwind compact difference schemes, and such schemes could be applied to different fields Wang2005 . In that same year, Kumar Kumar2009 discussed a high-order compact difference scheme for singularly perturbed reaction diffusion problems on a new Shish Kin mesh. Sen Mehra2017 ; Sen2016 discussed the fourth-order exact compact difference scheme for mixed derivative parabolic problems with variable coefficients.
The CFD scheme is a widely used method for spatial discretization of heat conduction equations to obtain the ODEs, and then other methods of time discretization are used for discretizing the ODEs, such as Euler method, multistep methods and Runge-Kutta method. The exact solution of heat conduction equation contains the calculation of exponential matrix. How to accurately calculate the exponential matrices is an essential problem in solving PDEs. Moler et al. Moler2003 summarized nineteen schemes for calculating the exponential matrices. These nineteen schemes are aimed at different practical problems, and their numerical solutions also have corresponding advantages and disadvantages. In 1994, Zhong Wan-Xie2004
proposed the precise integration method (PIM) of exponential matrices to solve the initial value problem of linear ODEs. PIM is an approximated method to calculate the exponential matrices, which contains Taylor approximation and Padé approximation. The PIM avoided the computer error caused by fine division and improved the numerical solution of exponential matrices by the accuracy of the computation.
Alternating Direction Implicit (ADI) method is a classical numerical scheme for solving multi-dimensional heat conduction equation. ADI, such as Peacemen-Rachford scheme, D’Yakonov scheme and Douglas scheme, are only the second-order accuracy schemes Zhang2016 ; Karaa2004 ; Journal2013 . ADI often fail to meet the accuracy requirements of practical problems. Strang splitting method (SSM) is a numerical method for solving differential equations that are decomposable into a sum of differential operators, which is to solve multi-dimensional PDEs by reducing their dimensionality to a sum of one-dimensional problems Gilbert2013 . This is a scheme of operator splitting method. If the differential operators of the SSM commute, then it will lead to no loss of accuracy. Therefore, the proposed schemes will extend to multi-dimensional heat conduction equation through SSM.
The remainder of the paper is arranged as follows. The n-dimensional Hopf-Cole transformation between the n-dimensional Burgers’ system and n-dimensional heat conduction equation are presented in Section 2; Moreover, we give the modification of the Hopf-Cole transformation. The high-order exponential time differencing PIM in combination with a spatially global sixth-order CFD scheme for solving n-dimensional heat condution equations are presented in Section 3. In Section 4, the Strang splitting method is described and the proposed schemes are extended to multi-dimensional problems. In Section 5, numerical examples are carried out to test the accuracy and adaptability of the proposed schemes. The conclusions are drawn in Section 6.
2 n-dimensional Hopf-Cloe transformation
The purpose of n-dimensional Hopf-Cole transformation is to convert Eq. (4) into n-dimensional heat equation
by the n-dimensional Hopf-Cole transformation
where , .
When and , the system of Eq. (8) respectively becomes one-dimensional, two-dimensional and three-dimensional heat equations
The initial and boundary conditions are
2.1 the modification of Hopf-Cole transformation
Based on this method, we intend to extend Hopf-Cole transformation to n-dimensional Burgers’ system. Assuming that the n-dimensional heat conduction equation has the irrotational condition
where are the basis of n-dimensional Euclidean space.
Applying Hopf-Cole transformation, Eq. (19) will becomes the n-dimensional heat conduction equations. Introduce for the system of Eq. (4), we obtain Eq. (8). It is especially noted that disappears in Hopf-Cole transformation.
With the Development of Hopf-Cole transformation in past decades, Kadalbajoo et al. Kadalbajoo2006 proposed that the C-N scheme for Eq. (5) based on Hopf-Cole transformation. They discretized the space twice with C-N scheme and central difference. Due to twice spatial dispersion in of Eq. (9), the numerical solution results in loss of accuracy. In 2015, Mukundan et al. Mukundan2015 present numerical techniques for Burgers’ equation, which use backward difference and central difference for . The accuracy of these numerical schemes will decline because of the two discretizations of . We have improved Hopf-Cole transformation, which will only be one dispersed in space. Firstly, Eq. (18) can be written as
Substituting into Eq. (20)
Then the solution of Eq. (21) can be obtained utilizing high precision numerical schemes such as CFD scheme. Thus, Eq. (21) will get after a spatial discretization for n-dimensional Burgers’ equations. In this way, the improvement of Hopf-Cole transformation avoids the truncation error of twice spatial difference and can obtain the higher precision the first derivative of .
2.2 The simplification of initial value problem
For some initial value problems, Fourier series solutions of Hopf-Cole transformation will converge very slowly, which dramatically increases the complexity of the calculation. In this subsection, we simplify the initial value condition of n-dimensional(one-dimensional, two-dimensional, three-dimensional) Burgers’ equation.
2.2.1 one-dimensional modification
It is widely noted that the analytical solution of one-dimensional heat conduction equation can be written in the standard form of the Fourier series
where is Fourier coefficient.
The initial conditions of one-dimensional heat conduction equation are extracted from the Eq. (32)
with boundary conditions
by one-dimensional Hopf-Cole transformation
The challenge of the initial value problem is to calculate the coefficient of the Eq. (29), which is difficult for that the single integral consisting of the exponential and trigonometric function. To simplify one-dimensional problem, the main work is to convert the calculations of into more efficient kind. can be written as
where is the modified Bessel function of the first kind and of order . In 2016, Gao et al. Gao2016 give numerical modification of analytical solution for two and three dimensional Burgers’ equation. Their modification is similar to our simplification, but Gao et al. did not provide one-dimensional case. Therefore, the following two-dimensional and three-dimensional improvements refer to the ideas put forward by Gao et al.
2.2.2 two-dimensional modification
For the two-dimensional Burgers’ equation (6), we select the modified Bessel function of the first kind and of order to replace Fourier coefficient under the initial condition. Considering the following initial and boundary conditions Gao2016 ; Zhang2009 ; Zhang2010 ; Siraj-ul-Islam2012
where the space domain is , and the time domain is .
It is widely noted that the analytical solution of two-dimensional heat conduction equation can be written in the standard form of the Fourier series
where is Fourier coefficient.
The initial conditions of two-dimensional heat conduction equation are extracted from the Eq. (32)
and the boundary conditions
by two-dimensional Hopf-Cole transformation
Applying the Fourier transformation to Eq. (33), we will obtain Fourier coefficient
The challenge of the initial value problem is to calculate the coefficient of the Eq. (38), which is difficult for that the double integral consisting of the exponential and trigonometric function. To simplify two-dimensional problem, the main work is to convert the calculations of into more efficient kind. can be written as
2.2.3 three-dimensional modification
Because the three-dimensional Burgers’ equation is too complex, its improvement is quite different from the one-dimensional and two-dimensional ones. Considering the following initial and boundary conditions Gao2016
where the space domain is ,and the time domain is .
It is widely noted that the analytical solution of three-dimensional heat conduction equation can be written in the standard form of the Fourier series
where is Fourier coefficient.
The initial conditions of three-dimensional heat conduction equation are extracted from the Eq. (43)
and the boundary conditions
By three-dimensional Hopf-Cole transformation
Applying the Fourier transformation to Eq. (44), we will obtain Fourier coefficient
The challenge of the initial value problem is to calculate the coefficient of the Eq. (49), which is difficult for that the triple integral consisting of the exponential and trigonometric function. To simplify three-dimensional problem, the main job is to convert the calculations of into more efficient kind. can be written as
where is the generalized hypergeometric series. Ref. Perger1993 defines the generalized hypergeometric series
in which is the Pochhammer symbol and is defined as:
The coefficient in Eq. (51) is defined by the following equation:
3 high-order numerical scheme
In this section, to solve the n-dimensional heat conduction equation obtained by Hopf-Cole transformation, we will present the sixth-order CFD scheme and the precise integration method (PIM). For simplicity, we consider one-dimensional heat conduction equation (10) with mesh size , in which , where is spatial step size. We firstly apply the sixth-order CFD scheme to discretization in space. If and represent the second derivative of at , then an approximation of the second derivatives at interior nodes may be expressed as
In order to make those near-boundary points have the same order accuracy as interior nodes, they should be obtained by matching Taylor series expansions to the order of at boundary points and , hence we get the following formulae Li2008
Therefore the sixth-order compact finite difference approximation of second derivatives is given by
3.1 Precise integration method
After the spatial discretization, the governing PDEs become the following ODEs
Giving as the temporal step size, then integrating Eq. (65) directly, the following recurrence formula is obtained
where is an exponential matrix.
The present work will focus on how to compute the exponential matrix very precisely. Moler et al. Moler2003 had discussed nineteen dubious ways to compute the exponential matrix, they pointed out that the problem of calculating exponential matrix had not been fully solved. In this paper, we apply the PIM to calculate the exponential matrix, which was proposed by Zhong et al. Wan-Xie2004
. The PIM is a algorithm of high precision for computing exponential matrix, which avoids the computer truncation error caused by the fine division and improves the numerical accuracy of the exponential matrix. In short, PIM is a series of matrix or vector multiplication calculations. Therefore, the main problem is how to calculate the exponential matrix. The precise computation of exponential matrix has two core contents Zhong2001 :
(1)The incremental part of the exponential matrix is calculated separately, rather than as a whole.
(2)The addition theorem of exponent is achieved by algorithm.
Using the addition theorem for the exponential matrix , the following equation is obtained
where is a relatively large positive integer. In order to ensure computational accuracy, Ref. Wan-Xie2004 suggested , is defined as bisection order.
3.1.1 Taylor approximation methods
for eigenvalues ofclose to 0. To overcome this problem and other numerical issues associated with it, many researchers have proposed different methods. This study focuses on the new technique and develops CFD scheme based on Taylor approximation of in order to alleviate computational difficulties associated with them. The Taylor expansion formula to the exponential matrix is defined as
where is an extremely short time, so the truncation error in time can be ignored, the fourth-order Taylor expansion can have high precision. Hence
Because is very small, it is enough to expand only the the first five terms of the series. The exponential matrix departs from the unit matrix to a very small extent. Hence it should be distinguished as
In order to obtain exponential matrix , we need to use algorithm for the matrix .
3.1.2 algorithm of the exponential matrix
. One of the core contents of PIM is the identity matrixcannot be directly added to the incremental matrix for Eq. (70). Because is a miniature matrix. If they add up directly, becomes the mantissa of in the process of computer operation. Thus, will become an appended part and its precision will seriously drop during the round-off operation in computer arithmetic. As a matter of fact, is an incremental part, which must be calculated and stored separately. Therefore, we will apply algorithm to calculate .
For computing the matrix , Eq. (67) should be factored as
Because has the following equation relation of factorization
Thus, Eq. (71) can be written as
The factorization (72) should be iterated times for . Then, no longer has a small value after such an iteration circulated times according to the following computer cycle language
At the end of the cycles, the computer stores . At this point, can be directly added to the identity matrix to obtain the exponential matrix
Therefore, the Taylor approximation of the PIM can be combined with algorithm to calculate the exponential matrix to obtain a high precision numerical solution . In the same way, we can use the CFD based on PIM (CFD-PIM) to obtain the solution of the first derivative of the one-dimensional heat conduction equation. According to one-dimensional Hopf-Cole transformation (26), and are submitted into Eq. (26) to get the solution of one-dimensional Burgers’ equation.
3.2 Stability analysis
To study the stability of our scheme, we only consider the periodic boundary condition for simplicity.
In Eq. (64), if is the eigenvalue of matrix , then is the eigenvalue of exponential matrix
with the same corresponding eigenvector. To prensent that CFD-PIM scheme is unconditionally stability, we need to prove the spectral radius of matrix is less than 1. To this end, the following two lemmas are needed.
Lemma 1. If is an eigenvalue of matrix with its corresponding eigenvector , then the eigenvalue is real number and .
Proof. By the definitions of eigenvalue and eigenvector, we may write , implying that Zhou2018 . This gives
Here, for periodic boundary condition the matrix and are as follows
Obviously, the matrix , and are really symmetrical, so the eigenvalue is real number. Meanwhile, for arbitrary , the right-hand side of the Eq. (76) is
Using the inequality , we obtain
and the left-hand side of the Eq. (76) is