A high-order compact finite difference scheme and precise integration method based on modified Hopf-Cole transformation for numerical simulation of n-dimensional Burgers' system

TL;DR

This paper introduces a high-order compact finite difference scheme combined with a precise integration method, based on a modified Hopf-Cole transformation, to accurately simulate multi-dimensional Burgers' systems.

Contribution

It develops a novel high-order numerical scheme that enhances accuracy for n-dimensional Burgers' equations using a modified Hopf-Cole transformation and advanced integration techniques.

Findings

Significantly improves computational accuracy over existing methods.

Demonstrates excellent adaptability in 2D and 3D simulations.

Achieves high accuracy in medium Reynolds number flows.

Abstract

This paper modifies a n-dimensional Hopf-Cole transformation to the n-dimensional Burgers' system. We obtain the n-dimensional heat conduction equation through the modification of the Hopf-Cole transformation. Then the fourth-order precise integration method (PIM) in combination with a spatially global sixth-order compact finite difference (CFD) scheme is presented to solve the equation with high accuracy. Moreover, coupling with the Strang splitting method, the scheme is extended to multi-dimensional (two, three-dimensional) Burgers' system. Numerical results show that the proposed method appreciably improves the computational accuracy compared with the existing numerical method.Moreover, the two-dimensional and three-dimensional examples demonstrate excellent adaptability, and the numerical simulation results also have very high accuracy in medium Reynolds numbers.

Figures

Tables

Table 1. Table 1: Comparison with the following numerical and analytical solutions of different values of x 𝑥 x and t 𝑡 t with R ​ e = 10 𝑅 𝑒 10 Re=10 for Example 1.

$x$	$t$	Ref.[58]	Ref.[42]	Ref.[45]	Ref.[44]	Proposed scheme	Analytical solution
0.25	0.4	0.30891	0.30881	0.30887	0.30889	0.308894228585555	0.308894228585318
	0.6	0.24075	0.24069	0.24070	0.24075	0.240739023291803	0.240739023291448
	0.8	0.19568	–	0.19566	0.19569	0.195675570103972	0.195675570103439
	1.0	0.16257	0.16254	0.16255	0.16258	0.162564857111346	0.162564857110671
	3.0	0.02720	0.02720	0.02721	0.02720	0.027202314473410	0.027202314472951
0.50	0.4	0.56964	0.56955	0.56956	0.56956	0.569632450695361	0.569632450693995
	0.6	0.44721	0.44714	0.44715	0.44724	0.447205521200320	0.447205521198742
	0.8	0.35924	–	0.35920	0.35927	0.359236058517410	0.359236058515669
	1.0	0.29192	0.29188	0.29188	0.29195	0.291915957127591	0.291915957125836
	3.0	0.04021	0.04021	0.04022	0.04021	0.040204924438755	0.040204924438046
0.75	0.4	0.62542	0.62540	0.62540	0.62537	0.625437893711249	0.625437893706948
	0.6	0.48721	0.48715	0.48716	0.48718	0.487214974885767	0.487214974882163
	0.8	0.37392	–	0.37389	0.37391	0.373921753212449	0.373921753209455
	1.0	0.28748	0.28744	0.28743	0.28747	0.287474405919467	0.287474405916976
	3.0	0.02977	0.02978	0.02978	0.02977	0.029772126859293	0.029772126858766

Table 2. Table 2: Comparison with the following numerical and analytical solutions of different values of x 𝑥 x and t 𝑡 t with R ​ e = 100 𝑅 𝑒 100 Re=100 for Example 1.

$x$	$t$	Ref.[58]	Ref.[42]	Ref.[45]	Ref.[44]	Proposed scheme	Analytical solution
0.25	0.4	0.34819	0.34229	0.34184	0.34191	0.341914932413026	0.341914932411983
	0.6	0.27536	0.26902	0.26891	0.26896	0.268964845317425	0.268964845316620
	0.8	0.22752	–	0.22143	0.22148	0.221481914524793	0.221481914524373
	1.0	0.19375	0.18817	0.18815	0.18820	0.188193961397110	0.188193961396738
	3.0	0.07754	0.07511	0.07510	0.07511	0.075114083887341	0.075114083887190
0.50	0.4	0.66543	0.66797	0.66060	0.66069	0.660710972121299	0.660710970851541
	0.6	0.53525	0.53211	0.52932	0.52942	0.529418263880240	0.529418263729147
	0.8	0.44526	–	0.43905	0.43914	0.439138250704212	0.374420037644682
	1.0	0.38047	0.37500	0.37436	0.37443	0.374420037662816	0.374420037644682
	3.0	0.15362	0.15018	0.15017	0.15019	0.150179005234648	0.150179005235832
0.75	0.4	0.91201	0.93680	0.91026	0.91023	0.910227039712648	0.910268136079484
	0.6	0.77132	0.77724	0.76719	0.76723	0.767241968056591	0.767243282478514
	0.8	0.65254	–	0.64745	0.64740	0.647395126988680	0.647395234822761
	1.0	0.56157	0.56157	0.55608	0.55606	0.556050682504419	0.556050704468246
	3.0	0.22874	0.22485	0.22504	0.22486	0.224811248098947	0.224811248193590

Table 3. Table 3: Numerical results of R ​ e = 5 𝑅 𝑒 5 Re=5 and R ​ e = 10 𝑅 𝑒 10 Re=10 with t = 1 𝑡 1 t=1 for Example 1.

	$N$	$11$	$21$	$41$	$81$
$Re=5$	$L_{\infty }$	1.4625E-06	4.2340E-09	2.4716E-11	3.0104E-13
	CPU(s)	0.604	0.674	0.780	1.508
	ROC	–	8.4322	7.4204	6.3340
$Re=10$	$L_{\infty }$	2.8131E-06	7.3746E-09	3.7254E-11	2.1694E-12
	CPU(s)	1.150	1.207	1.479	2.786
	ROC	–	8.5754	7.6374	4.1719

Table 4. Table 4: L ∞ subscript 𝐿 L_{\infty} , rates of convergence, and CPU time of the following numerical schemes in Matlab for Example 2.

	Ref.[11]			Proposed scheme
$h$	$L_{\infty }$	ROC	CPU(s)	$L_{\infty }$	ROC	CPU(s)
$\pi /8$	3.660E-05	-	0.1492	1.149E-05	-	0.071
$\pi /16$	2.277E-06	4.0066	0.3109	2.968E-08	8.8967	0.079
$\pi /32$	1.422E-07	4.0017	0.5419	1.329E-10	7.8030	0.081
$\pi /64$	8.882E-09	4.0004	1.0726	1.480E-12	6.4886	0.118
platform	Intel Core i7-4510U 2.60 GHz workstation			Intel Core i5-6300HQ 2.30 GHz processor

Table 5. Table 5: Comparsion of L ∞ subscript 𝐿 L_{\infty} of the different Re at t = 1 𝑡 1 t=1 for Example 3.

	$N$	$11\times 11$	$21\times 21$	$41\times 41$	$81\times 81$
$Re=100$	$L_{\infty }\left(u\right)$	5.2252E-05	5.1284E-07	2.4176E-09	1.4047E-11
	ROC	–	6.6708	7.7288	7.4272
	$L_{\infty }\left(v\right)$	1.2944E-04	1.1137E-06	5.0211E-09	4.2705E-11
	ROC	–	6.8607	7.7932	6.8761
$Re=1000$	$L_{\infty }\left(u\right)$	1.8558E-06	4.3507E-08	4.0918E-10	2.0677E-12
	ROC	–	5.4146	6.7323	7.6286
	$L_{\infty }\left(v\right)$	2.6879E-06	5.8267E-08	5.4288E-10	2.7184E-12
	ROC	–	5.5277	6.7459	7.6423

Table 6. Table 6: Comparison with the proposed scheme and the analytical solutions of different coordinate positions ( x , y ) 𝑥 𝑦 (x,y) with R ​ e = 100 𝑅 𝑒 100 Re=100 , at t = 0.25 , 0.5 𝑡 0.25 0.5 t=0.25,0.5 for Example 4.

	$t=0.25$		$t=0.5$
$(x,y)$	Analytical solution[12]	Proposed scheme	Analytical solution[12]	Proposed scheme
$(0.25,0.25)$	0.3935490117704355	0.393549011771103	0.2911828920816955	0.291182892082807
$(0.50,0.25)$	0.6861822403861596	0.686182243447071	0.5605467081571704	0.560546708625940
$(0.75,0.25)$	0.3935490117704355	0.393524523103096	0.2911828920816955	0.291180266677443
$(0.25,0.50)$	0.2619506158413131	0.261950614753384	0.2619829318741635	0.261982931825796
$(0.50,0.50)$	3.433140255652116E-70	0	6.648346088881589E-70	0
$(0.75,0.50)$	-0.2619506158413131	-0.261950614753467	-0.2619829318741635	-0.261982931825921
$(0.25,0.75)$	-0.3935490117704355	-0.393513909760747	-0.2911828920816955	-0.291179204683133
$(0.50,0.75)$	-0.6861822403861596	-0.686182243439029	-0.5605467081571704	-0.560546708617027
$(0.75,0.75)$	-0.3935490117704355	-0.393549011774452	-0.2911828920816955	-0.291182892085496

Table 7. Table 7: Comparison with the proposed scheme and the exact solutions of different coordinate positions of ( x , y ) 𝑥 𝑦 (x,y) with R ​ e = 100 𝑅 𝑒 100 Re=100 and t = 0.75 , 1.0 𝑡 0.75 1.0 t=0.75,1.0 for Example 4.

	$t=0.75$		$t=1.0$
$(x,y)$	Analytical solution[12]	Proposed scheme	Analytical solution[12]	Proposed scheme
$(0.25,0.25)$	0.2273774661403168	0.227377466141544	0.1858035619888798	0.185803561989915
$(0.50,0.25)$	0.4479960634879613	0.447996063594949	0.3687873026249988	0.368787302662550
$(0.75,0.25)$	0.2273774661403168	0.227377282040238	0.1858035619888798	0.185803538891363
$(0.25,0.50)$	0.2180447995061038	0.218044799502132	0.1818603555820653	0.181860355581896
$(0.50,0.50)$	2.301061067978292E-70	0	7.783304094972026E-71	0
$(0.75,0.50)$	-0.2180447995061038	-0.218044799502212	-0.1818603555820653	-0.181860355581938
$(0.25,0.75)$	-0.2273774661403168	-0.227377196569641	-0.1858035619888798	-0.185803526285714
$(0.50,0.75)$	-0.4479960634879613	-0.447996063584584	-0.3687873026249988	-0.368787302652006
$(0.75,0.75)$	-0.2273774661403168	-0.227377466144151	-0.1858035619888798	-0.185803561992506

Table 8. Table 8: Numerical results of κ = 1 , 2 , 5 , 10 𝜅 1 2 5 10 \kappa=1,2,5,10 with t = 1 𝑡 1 t=1 and R ​ e = 100 𝑅 𝑒 100 Re=100 for Example 5.

$\kappa$	$N$	$11\times 11$	$21\times 21$	$41\times 41$	$81\times 81$
1	$L_{\infty }$	1.1483E-07	2.7964E-10	8.8463E-13	7.2616E-15
	ROC	–	8.6817	8.3051	7.6316
2	$L_{\infty }$	5.7416E-08	1.3982E-10	4.4232E-13	3.6308E-15
	ROC	–	8.6817	8.3051	7.6316
5	$L_{\infty }$	2.2966E-08	5.5929E-11	1.7693E-13	1.4520E-15
	ROC	–	8.6817	8.3051	7.6316
10	$L_{\infty }$	1.1483E-08	2.7964E-11	8.8463E-14	7.2598E-16
	ROC	–	8.6817	8.3051	7.6316

Table 9. Table 9: The solutions of the three-dimensional heat conduction equation with R ​ e = 10 𝑅 𝑒 10 Re=10 at t = 0.1 𝑡 0.1 t=0.1 for Example 6.(The spatial step sizes of the FEM are h x = h y = h z = 10 − 2 subscript ℎ 𝑥 subscript ℎ 𝑦 subscript ℎ 𝑧 superscript 10 2 h_{x}=h_{y}=h_{z}=10^{-2} )

	$\varphi (x,y,z,t)$
$(x,y,z)$	Analytical solution[12]	FEM[12]	Proposed scheme
$(0.25,0.00,0.00)$	0.5121139094268445	0.5121574025971580	0.512112211115770
$(0.25,0.25,0.00)$	0.4021754605524777	0.4021832308799626	0.402175090853434
$(0.25,0.25,0.25)$	0.3332558679440735	0.3332624499504631	0.333255848934007
$(0.25,0.50,0.00)$	0.2221335878372145	0.2221195126283573	0.222134351589450
$(0.25,0.50,0.25)$	0.2145599507332575	0.2145578143852998	0.214560278902764
$(0.25,0.75,0.00)$	0.1232224673074507	0.1232115116009095	0.123222177969065
$(0.50,0.00,0.00)$	0.2350307713064107	0.2349996700538997	0.235032545451199
$(0.50,0.25,0.00)$	0.2221335878372145	0.2221195160032975	0.222134351729473
$(0.50,0.25,0.25)$	0.2145599507332575	0.2145578324085150	0.214560278977476
$(0.75,0.00,0.00)$	0.09757645937131027	0.09756832939249252	0.0975749310522441
$(0.75,0.25,0.00)$	0.1232224673074507	0.1232114495876140	0.123222177969033

Equations

$$\left\{\begin{array}[]{l}\frac{{\partial u_{1}}}{{\partial t}}=\omega_{1}\frac{{{\partial^{2}}u_{1}}}{{\partial{x^{2}}}}-\kappa_{1}u_{1}\frac{{\partial u_{1}}}{{\partial x}}-\delta_{1}\frac{{\partial({u_{1}}{u_{2}})}}{{\partial x}}\\ \frac{{\partial u_{2}}}{{\partial t}}=\omega_{2}\frac{{{\partial^{2}}u_{2}}}{{\partial{x^{2}}}}-\kappa_{2}u_{2}\frac{{\partial u_{2}}}{{\partial x}}-\delta_{2}\frac{{\partial({u_{1}}{u_{2}})}}{{\partial x}}\end{array}\right.,{\kern 1.0pt}{\kern 1.0pt}{\kern 1.0pt}{\kern 1.0pt}{x\in\Omega=[a,b]},{\kern 1.0pt}{\kern 1.0pt}t\in[0,T]$$

$$\begin{array}[]{l}u_{1}(x,0)={g_{1}}(t){\kern 1.0pt}{\kern 1.0pt},{\kern 1.0pt}{\kern 1.0pt}{\kern 1.0pt}{x\in\Omega=[a,b]}\\ u_{2}(x,0)={g_{2}}(t){\kern 1.0pt}{\kern 1.0pt},{\kern 1.0pt}{\kern 1.0pt}{\kern 1.0pt}{x\in\Omega=[a,b]}\end{array}$$

$$\begin{array}[]{l}u_{1}(a,t)={f_{1}}(x){\kern 1.0pt}{\kern 1.0pt},{\kern 1.0pt}{\kern 1.0pt}u_{1}(a,t)={f_{2}}(x){\kern 1.0pt}{\kern 1.0pt},{\kern 1.0pt}{\kern 1.0pt}{\kern 1.0pt}{x\in\Omega=[a,b]}{\kern 1.0pt}{\kern 1.0pt}{\kern 1.0pt},{\kern 1.0pt}{\kern 1.0pt}{\kern 1.0pt}{\kern 1.0pt}{\kern 1.0pt}t\in[0,T]\\ u_{2}(b,t)={f_{3}}(x){\kern 1.0pt}{\kern 1.0pt},{\kern 1.0pt}{\kern 1.0pt}{\kern 1.0pt}u_{2}(b,t)={f_{4}}(x){\kern 1.0pt}{\kern 1.0pt},{\kern 1.0pt}{\kern 1.0pt}{\kern 1.0pt}{\kern 1.0pt}{x\in\Omega=[a,b]}{\kern 1.0pt}{\kern 1.0pt},{\kern 1.0pt}{\kern 1.0pt}{\kern 1.0pt}{\kern 1.0pt}{\kern 1.0pt}t\in[0,T]\end{array}$$

$${\mathbf{U}}_t+\kappa (\mathbf{U}\cdot \nabla )\mathbf{U}=\omega \Delta \mathbf{U}$$

$$u_t+u{u}_x=\omega u_{xx}$$

$$\begin{array}[]{l}{u_{t}}+u{u_{x}}+v{u_{y}}=\omega({u_{xx}}+{u_{yy}})\\ {v_{t}}+u{v_{x}}+v{v_{y}}=\omega({v_{xx}}+{v_{yy}})\end{array}$$

$$\begin{array}[]{l}{u_{t}}+u{u_{x}}+v{u_{y}}+w{u_{z}}=\omega({u_{xx}}+{u_{yy}}+{u_{zz}})\\ {v_{t}}+u{v_{x}}+v{v_{y}}+w{v_{z}}=\omega({v_{xx}}+{v_{yy}}+{v_{zz}})\\ {w_{t}}+u{w_{x}}+v{w_{y}}+w{w_{z}}=\omega({w_{xx}}+{w_{yy}}+{w_{zz}})\end{array}$$

$${\varphi }_t-\omega △\varphi =0$$

$$u_i=-2\frac{\omega }{\kappa }{\partial }_{x_i}\ln \varphi =-2\frac{\omega }{\kappa }\frac{{\varphi }_{x_i}}{\varphi }$$

$${\varphi }_t-\omega {\varphi }_{xx}=0$$

$${\varphi }_t-\omega ({\varphi }_{xx}+{\varphi }_{yy})=0$$

$${\varphi }_t-\omega ({\varphi }_{xx}+{\varphi }_{yy}+{\varphi }_{zz})=0$$

$$\varphi (x_i,0)=\exp ({\int }_0^{{\mathbf{x}}_{\mathbf{i}}}-\frac{u_i(\xi ,0)}{2\omega }d\xi ),x_i\in \Omega =[a,b]$$

$${\varphi }_{x_i}(a,t)={\varphi }_{x_i}(b,t)=0,t\in [0,T]$$

$$\nabla \times \mathbf{U}=\sum _{i,j=1}^n(\frac{\partial u_j}{\partial x_i}-\frac{\partial u_i}{\partial x_j})e_i\wedge e_j=0$$

$$\frac{\partial u_i}{\partial t}+\kappa \sum _{j=1}^n{u}_j\frac{\partial u_i}{\partial x_j}=\omega \Delta u_i,i=1,2,\dots ,n$$

$$\frac{\partial u_i}{\partial x_j}=\frac{\partial u_j}{\partial x_i},i,j=1,2,\dots ,n(i\ne j)$$

$$\frac{{\partial }^2\phi }{\partial x_j\partial t}+\kappa \sum _{j=1}^n\frac{\partial \phi }{\partial x_j}\frac{{\partial }^2\phi }{\partial x_i\partial x_j}-\omega \frac{\partial \Delta \phi }{\partial x_j}=0,i=1,2,\dots ,n$$

$$\frac{\partial \phi }{\partial t}+\frac{\kappa }{2}\sum _{j=1}^n{(\frac{\partial \phi }{\partial x_j})}^2-\omega \Delta \phi =0$$

$$\frac{\partial \phi }{\partial t}+\frac{\kappa }{2}\sum _{j=1}^n{\left(\frac{\partial \phi }{\partial x_j}\right)}^2-\omega \sum _{j=1}^n\frac{{\partial }^2\phi }{\partial x_j^2}=0$$

$$-\frac{2\omega }{\kappa }\frac{{\varphi }_t}{\varphi }+\frac{\kappa }{2}{\left(\frac{-2\omega }{\kappa }\frac{1}{\varphi }\sum _{j=1}^n\frac{\partial \varphi }{\partial x_j}\right)}^2+\frac{2{\omega }^2}{\kappa }\frac{1}{{\varphi }^2}\left(\sum _{j=1}^n\left(\frac{{\partial }^2\varphi }{\partial x_j^2}\right)\varphi -\sum _{j=1}^n{\left(\frac{\partial \varphi }{\partial x_j}\right)}^2\right)=0$$

$$-2\omega \frac{{\varphi }_t}{\varphi }+2{\omega }^2\frac{1}{{\varphi }^2}\sum _{j=1}^n{\left(\frac{\partial \varphi }{\partial x_j}\right)}^2-2{\omega }^2\frac{1}{\varphi }\sum _{j=1}^n\left(\frac{{\partial }^2\varphi }{\partial x_j^2}\right)-2{\omega }^2\frac{1}{{\varphi }^2}\sum _{j=1}^n{\left(\frac{\partial \varphi }{\partial x_j}\right)}^2=0$$

$${\varphi }_t+\omega \sum _{j=1}^n\frac{{\partial }^2\varphi }{\partial x_j^2}=0\Rightarrow {\varphi }_t+\omega \Delta \varphi =0$$

$$\frac{\partial {\phi }_x}{\partial t}+\frac{\kappa }{2}\sum _{j=1}^N{(\frac{\partial {\phi }_x}{\partial x_j})}^2-\omega \Delta {\phi }_x=0$$

$$\frac{\partial {\varphi }_x}{\partial t}-\omega △{\varphi }_x=0$$

$$\begin{array}[]{l}u(x,0)=u_{0}(x)=\sin\pi x{\kern 1.0pt}{\kern 1.0pt}{\kern 1.0pt},{\kern 1.0pt}{\kern 1.0pt}{\kern 1.0pt}x\in[0,{\rm{1}}]\\ u(0,t)=u(1,t)=0{\kern 1.0pt}{\kern 1.0pt},{\kern 1.0pt}{\kern 1.0pt}t>0\end{array}$$

$$\varphi (x,t)=\sum _{\alpha =0}^{\infty }{C}_{\alpha }\exp (-{\alpha }^2{\pi }^2\omega t)\cos (\alpha \pi x)$$

$$\varphi (x,0)=\sum _{\alpha =0}^{\infty }{C}_{\alpha }\cos (\alpha \pi x),x\in \Omega =[0,1]$$

$${\varphi }_x(0,t)={\varphi }_x(1,t)=0,t\in [0,T]$$

$$u=-2\omega \frac{{\varphi }_x}{\varphi }$$

Full text

††journal: Applied Mathematics and Computation

1 Introduction

Burgers’ equation is a nonlinear partial differential equation (PDE) which was first introduced by Bateman [1], and was later treated as the turbulence of the mathematical model [2, 3]. 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 [4, 5], etc. Hopf [6] and Cole [7] 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 old 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 a 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 [8]

[TABLE]

subject to the initial conditions:

[TABLE]

and boundary conditions

[TABLE]

where $\omega_{1},\omega_{2}$ is kinematic viscosity parameters of the fluid, which correspond to an inverse of Reynolds number $Re$ (if $\omega_{1}=\omega_{2}$, then $\omega_{1}=\omega_{2}=\frac{1}{Re}$ ) . $\kappa_{1},\kappa_{2}$ are real constants and $\delta_{1},\delta_{2}$ are arbitrary constants. $f(x)$ and $g(t)$ 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 [9, 10]. 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 [11, 12, 13, 14, 15, 16]. For the sake of clarity, a brief description of the comparison method is provided below.

Bhatt et al. [11] 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. [12], 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 ($Re\geq 100$), which can be used as a reference for testing numerical methods.

In Ref. [13], the authors develop 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. [14], the author uses a 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 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. [15] establish new two-mode coupled Burgers’ equations which are introduced. The authors find the necessary conditions in which the multiple kinks and multiple singular kink solutions exist and present the two-front solutions.

Jiwari et al. [16, 17] developed a differential quadrature method to solve time dependent Burgers’ equation. And results are accurately produced by such two numerical schemes provided by the authors. These two schemes are also found quite easy to implement.

In Ref. [18], the authors proposed an algorithm based on exponential modified cubic-B-spline differential quadrature method for Burgers’ equation. With some modifications, such method is flexible enough to solve model equations in multi-dimensional problems including mechanical, physical or biophysical effects.

In Ref. [19], finite element analysis and approximation of Burgers’-Fisher equation with non-smooth initial data is presented, which provides results with high accuracy and efficiency.

In this paper, we mainly discuss the numerical scheme of n-dimensional Burgers’ system. The n-dimensional ($n\geq 1$) Burgers’ system includes low-dimensional ($n=1$) Burgers’ equations and multi-dimensional ($n\geq 2$) Burgers’ equations.

The n-dimensional Burgers’ system [20, 21, 22]

[TABLE]

where $\bm{U}=(u_{1},u_{2},...,u_{n})^{T}$ is the fluid velocity fields, $n$ is dimension of space, $\omega$ is the kinematic viscosity of the fluid and $\kappa$ are real constants. $\Delta=\frac{{{\partial^{2}}}}{{\partial x_{1}^{2}}}+\frac{{{\partial^{2}}}}{{\partial x_{2}^{2}}}+\cdots+\frac{{{\partial^{2}}}}{{\partial x_{n}^{2}}}$ denotes the Laplace operator, and $\nabla={(\frac{\partial}{{\partial{x_{1}}}},\frac{\partial}{{\partial{x_{2}}}},\cdots,\frac{\partial}{{\partial{x_{n}}}})^{T}}$ is the Hamilton gradient operator.

The system of Eq. (4) is Burgers’ equations which includes non-coupled and coupled problems. When $n=1,2,3$ and $\kappa=1$, the system of Eq.(4) respectively becomes one-dimensional, two-dimensional and three-dimensional Burgers’ equation

[TABLE]

[TABLE]

[TABLE]

Chen et al. [21] 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 [23, 24], 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 [25]. 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 [26], this scheme can get more accurate results with a thicker grid. Lele [27] 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 [28, 29, 30, 31], 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 [32]. In that same year, Kumar [33] discussed a high-order compact difference scheme for singularly perturbed reaction diffusion problems on a new Shish Kin mesh. Sen [34, 35] 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. [36] 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 [37] 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 [38, 39, 40]. 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 [41]. 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 The n-dimensional Hopf-Cole transformation

The purpose of n-dimensional Hopf-Cole transformation is to convert Eq. (4) into the n-dimensional heat equation

[TABLE]

by the n-dimensional Hopf-Cole transformation

[TABLE]

where $i=1,2,...,n$. Note: $x_{1}=x,x_{2}=y,x_{3}=z$; $u_{1}=u,u_{2}=v,u_{3}=w$.

When $n=1,2,3$ and $\kappa=1$, the system of Eq. (8) respectively becomes one-dimensional, two-dimensional and three-dimensional heat equations

[TABLE]

[TABLE]

[TABLE]

The initial and boundary conditions are

[TABLE]

[TABLE]

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

[TABLE]

where ${e_{i}}{\kern 1.0pt}{\kern 1.0pt},{\kern 1.0pt}{\kern 1.0pt}{e_{j}}$ are the basis of n-dimensional Euclidean space.

To facilitate readers to understand the derivation process, Eqs. (4) and (15) can be written as the following scalar forms

[TABLE]

[TABLE]

Let $u_{i}={\frac{{\partial\varphi}}{{\partial{x_{i}}}}}$, substituting in Eq. (16), we obtain

[TABLE]

then Eq. (18) can be written as

[TABLE]

Applying Hopf-Cole transformation, Eq. (19) will becomes the n-dimensional heat conduction equations. The detailed derivation process is as follows:

(1) Eq. (19) can be written as

[TABLE]

where $\sum\limits_{j=1}^{n}{\frac{{{\partial^{2}}}}{{\partial x_{j}^{2}}}}=\Delta$.

(2) Introduce $\varphi=-2\frac{\omega}{\kappa}\ln{\bm{\phi}}$ for the system of Eq. (20)

[TABLE]

(3) The two sides of Eq.(21) are multiplied by $\kappa$ and then simplified.

[TABLE]

It is especially noted that $\kappa$ disappears in Eq. (22).

(4) And further simplify to obtain

[TABLE]

2.1 The modification of Hopf-Cole transformation

With the Development of Hopf-Cole transformation in the past decades, Kadalbajoo et al. [42] proposed the C-N scheme based on the Hopf-Cole transformation for Eq. (5) . They discretized the space twice with C-N scheme and central difference. Due to the twice spatial dispersions of ${\phi_{x}}$ of Eq. (9), the numerical solution results were in loss of accuracy. In 2015, Mukundan et al. [43] presented numerical techniques for Burgers’ equation, which use backward difference and central difference for ${\phi_{x}}$. The accuracy of these numerical schemes will decline because of the twice discretizations of ${\phi_{x}}$. We have improved Hopf-Cole transformation, which will only be dispersed once in space. Hopf-Cole transformation is used again, but the object to be solved this time is the first derivative ${\phi_{x}}$ of the heat conduction equation. Firstly, Eq. (18) can be written as

[TABLE]

Substituting $\varphi=-2\frac{\omega}{\kappa}\ln{{{\phi_{x}}}}$ into Eq. (24)

[TABLE]

Initial and boundary conditions of Eq. (25) can be obtained from Eqs. (13) and (14).

Then the solution of Eq. (25) can be obtained by utilizing high precision numerical schemes such as CFD scheme. Thus, Eq. (25) will get ${\phi_{x}}$ after a spatial discretization for n-dimensional Burgers’ equations. In this way, the modification of n-dimensional Hopf-Cole transformation avoids the truncation error of twice spatial difference and can obtain the first derivative of ${\phi}({x},t)$ with higher precision.

The modification of n-dimensional Hopf-Cole transformation design in this section lies in two points:

(1) The modification of Hopf-Cole transformation is more general and suitable for Burgers’ system where $\kappa$ is a variable;

(2) Hopf-Cole transformation is used twice to solve the first derivative ${\phi_{x}}$ and solution ${\phi}$ of the heat conduction equation, thus avoiding the second truncation error.

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. In 2016, Gao et al. [12] gave 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.1 One-dimensional modification

Researchers have proposed the one-dimensional Burgers’ equation with the following initial and boundary condition [42, 44, 45, 46, 47]

[TABLE]

It is widely noted that the analytical solution of the one-dimensional heat conduction equation can be written in the standard form of the Fourier series

[TABLE]

where ${C_{\alpha}}$ is Fourier coefficient.

The initial conditions of the one-dimensional heat conduction equation are extracted from the Eq. (27)

[TABLE]

with boundary conditions

[TABLE]

by the one-dimensional Hopf-Cole transformation

[TABLE]

Applying the Fourier transformation to Eq. (28), we will obtain Fourier coefficient ${C_{\alpha}}$

[TABLE]

where

[TABLE]

[TABLE]

The challenge of the initial value problem is to calculate the coefficient $B_{\alpha}$ of the Eq. (33), which is difficult for the single integral consisting of the exponential and trigonometric function. To simplify the one-dimensional problem, the main work is to convert the calculations of $B_{\alpha}$ into more efficient kind. $B_{\alpha}$ can be written as

[TABLE]

where $I_{(n)}(1/2\omega\pi)$ is the modified Bessel function of the first kind and of order $n$.

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 $n$ to replace Fourier coefficient under the initial condition. Considering the following initial and boundary conditions [12, 48, 49, 50]

[TABLE]

where the space domain is $(x,{\kern 1.0pt}{\kern 1.0pt}y)\in\Omega=[0,1]\times[0,1]$, and the time domain is $t>0$.

It is widely noted that the analytical solution of the two-dimensional heat conduction equation can be written in the standard form of the Fourier series

[TABLE]

where ${C_{\alpha\beta}}$ is Fourier coefficient.

The initial conditions of the two-dimensional heat conduction equation are extracted from the Eq. (36)

[TABLE]

and the boundary conditions

[TABLE]

by the two-dimensional Hopf-Cole transformation

[TABLE]

Applying the Fourier transformation to Eq. (37), we will obtain Fourier coefficient ${C_{\alpha\beta}}$

[TABLE]

where

[TABLE]

[TABLE]

where

[TABLE]

The challenge of the initial value problem is to calculate the coefficient $B_{\alpha\beta}$ of the Eq. (42), which is difficult for the double integral consisting of the exponential and trigonometric function. To simplify the two-dimensional problem, the main work is to convert the calculations of $B_{\alpha\beta}$ into more efficient kind. $B_{\alpha\beta}$ can be written as

[TABLE]

where $I_{(n)}(\textstyle\frac{1}{4\omega\pi})$ is the modified Bessel function of the first kind and of order $n$. Gao [12] proved that Eq. (42) and Eq. (44) are equal, the scheme will get more precise and efficient analytical solution through the modification.

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 [12]

[TABLE]

[TABLE]

where the space domain is $(x,y,z)\in\Omega=[0,1]\times[0,1]\times[0,1]$,and the time domain is $t>0$.

It is widely noted that the analytical solution of the three-dimensional heat conduction equation can be written in the standard form of the Fourier series

[TABLE]

where ${C_{\alpha\beta\gamma}}$ is Fourier coefficient.

The initial conditions of the three-dimensional heat conduction equation are extracted from the Eq. (47)

[TABLE]

and the boundary conditions

[TABLE]

by the three-dimensional Hopf-Cole transformation

[TABLE]

Applying the Fourier transformation to Eq. (48), we will obtain Fourier coefficient ${C_{\alpha\beta\gamma}}$

[TABLE]

where

[TABLE]

[TABLE]

where

[TABLE]

The challenge of the initial value problem is to calculate the coefficient $B_{\alpha\beta\gamma}$ of the Eq. (53), which is difficult for that the triple integral consisting of the exponential and trigonometric function. To simplify the three-dimensional problem, the main job is to convert the calculations of $B_{\alpha\beta\gamma}$ into more efficient kind. $B_{\alpha\beta\gamma}$ can be written as

[TABLE]

where $G={}_{3}{F_{4}}({\textstyle{{\beta+1}\over 2}},{\textstyle{\beta\over 2}}+1,{\textstyle{\beta\over 2}}+j;{\textstyle{{\beta+\alpha}\over 2}}+1,{\textstyle{{\beta-\alpha}\over 2}}+1,m+1,{\textstyle{{\beta-\alpha}\over 2}}+j+1;{({\textstyle{1\over{4\omega\pi}}})^{2}})$ is the generalized hypergeometric series. Ref. [51] defines the generalized hypergeometric series

[TABLE]

in which $(a)_{k}$ is the Pochhammer symbol and is defined as:

[TABLE]

The coefficient $\mu_{j}$ in Eq. (55) is defined by the following equation:

[TABLE]

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 ${h=x_{i+1}-x_{i}}$, in which $x_{i}=ih,i=1,2,...,N$, where $h$ is spatial step size. We firstly apply the sixth-order CFD scheme to discretization in space. If $\phi=\phi(x_{i})$ and $\phi{}^{\prime\prime}_{i}$ represent the second derivative of $\phi(x)$ at $x_{i}$, then an approximation of the second derivatives at interior nodes may be expressed as

[TABLE]

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 $O(h^{6})$ at boundary points $1,2,N-1$ and $N$, hence we get the following formulae [25]

[TABLE]

[TABLE]

[TABLE]

[TABLE]

Writing Eqs. (59,60,61,62,63) in matrix form as

[TABLE]

where

[TABLE]

[TABLE]

[TABLE]

Therefore the sixth-order compact finite difference approximation of second derivatives ${\bm{\Phi^{\prime\prime}}}$ is given by

[TABLE]

where ${\bm{H}}={\bm{A^{-1}B}}$.

3.1 Precise integration method

After the spatial discretization, the governing PDEs become the following ODEs

[TABLE]

Giving $\tau{\rm{=}}{t_{k+1}}-{t_{k}}$ as the temporal step size, then integrating Eq. (69) directly, the following recurrence formula is obtained

[TABLE]

where ${\bm{T}}\left(\tau\right)={e^{{\bm{H}}\tau}}$ is an exponential matrix.

The present work will focus on how to compute the exponential matrix ${\bm{T}}$ very precisely. Moler et al. [36] 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. [37]. 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 ${e^{{\bm{H}}\tau}}$. The precise computation of exponential matrix has two core contents [52]:

(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 $2^{n}$ algorithm.

Using the addition theorem for the exponential matrix ${e^{{\bm{H}}\tau}}$, the following equation is obtained

[TABLE]

where $m$ is a relatively large positive integer and $\Delta t=\frac{\tau}{m}$. Thus $\Delta t$ is an extremely short time. In order to ensure computational accuracy, Ref. [37] suggested $m={2^{n}},n=20,m=1048576$ , $n$ is defined as bisection order.

3.1.1 Taylor approximation methods

The accurate computation of ${e^{{\bm{H}}\tau}}$ is a challenging problem for the numerical analysis [11, 36]. The major issue is the cancellation error arising during the direct computation of ${e^{{\bm{H}}\tau}}$ for eigenvalues of ${e^{{\bm{H}}}}$ close 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 ${e^{{\bm{H}}\tau}}$ in order to alleviate computational difficulties associated with them. The Taylor expansion formula to the exponential matrix ${e^{{\bm{H}}\Delta t}}$ is defined as

[TABLE]

where $\Delta t$ is extremely short, so the truncation error in time can be ignored, the fourth-order Taylor expansion can have high precision. Hence

[TABLE]

Because $\tau$ is very small, it is enough to expand only the the first five terms of the series. The exponential matrix ${\bm{T}}\left({\Delta t}\right)$ departs from the unit matrix ${\bm{I}}$ to a very small extent. Hence it should be distinguished as

[TABLE]

In order to obtain exponential matrix ${\bm{T}}(\tau)$, we need to use $2^{n}$ algorithm for the matrix ${\bm{T}}({\Delta t})$.

3.1.2 $2^{n}$ algorithm of the exponential matrix

PIM has the problem of complete loss of precision in the exponential additional theorem [52, 53, 54]. One of the core contents of PIM is the identity matrix ${\bm{I}}$ cannot be directly added to the incremental matrix ${\bm{T_{a}}}$ for Eq. (74). Because ${\bm{T_{a}}}$ is a miniature matrix. If they add up directly, ${\bm{T_{a}}}$ becomes the mantissa of ${\bm{I+T_{a}}}$ in the process of computer operation. Thus, ${\bm{T_{a}}}$ will become an appended part and its precision will seriously drop during the round-off operation in computer arithmetic. As a matter of fact, ${\bm{T_{a}}}$ is an incremental part, which must be calculated and stored separately. Therefore, we will apply $2^{n}$ algorithm to calculate ${\bm{T_{a}}}$.

For computing the matrix ${\bm{T}}(\tau)={e^{{\bm{H}}\tau}}$, Eq. (71) should be factored as

[TABLE]

Because ${\bm{T}}(\alpha)$ has the following equation relation of factorization

[TABLE]

Thus, Eq. (75) can be written as

[TABLE]

The factorization (76) should be iterated $n$ times for ${\bm{T}}(\tau)$. Then, ${\bm{T_{a}}}$ no longer has a small value after such an iteration circulated $n$ times according to the following computer cycle language

[TABLE]

At the end of the $n$ cycles, the computer stores ${\bm{T_{a}}}$. At this point, ${\bm{T_{a}}}$ can be directly added to the identity matrix ${\bm{I}}$ to obtain the exponential matrix ${\bm{T}}(\tau)$

[TABLE]

Therefore, the Taylor approximation of the PIM can be combined with $2^{n}$ algorithm to calculate the exponential matrix ${\bm{T}}(t)$ to obtain a high precision numerical solution $\bm{\Phi}(x,t)$. In the same way, we can use the CFD based on PIM (CFD-PIM) to obtain the solution $\bm{\Phi}_{x}(x,t)$ of the first derivative of the one-dimensional heat conduction equation. According to one-dimensional Hopf-Cole transformation (30), $\bm{\Phi}(x,t)$ and $\bm{\Phi}_{x}(x,t)$ are submitted into Eq. (30) to get the solution $u(x,t)$ of one-dimensional Burgers’ equation.

3.2 Stability analysis

3.2.1 Stability of Periodic boundary condition

To study the stability of our scheme, we only consider the periodic boundary condition for simplicity.

In Eq. (68), if ${\lambda_{i}}(i=1,2,...,N-1)$ is the eigenvalue of matrix ${\bm{H}}$, then ${e^{{\lambda_{i}}\tau}}$ is the eigenvalue of exponential matrix ${e^{{\bm{H}}\tau}}$ with the same corresponding eigenvector $\bm{x}={({x_{1}},{x_{2}},\ldots,{x_{N-1}})}$. To prensent that CFD-PIM scheme is unconditionally stability, we need to prove the spectral radius of matrix ${e^{{\bm{H}}\tau}}$ is less than 1. To this end, the following two lemmas are needed.

Lemma 1. If ${\lambda_{i}}$ is an eigenvalue of matrix ${\bm{H}}={{\bm{A}}^{-1}}{\bm{B}}$ with its corresponding eigenvector $\bm{x}$, then the eigenvalue ${\lambda_{i}}$ is real number and ${\lambda_{i}}\leq 0$ .

Proof. By the definitions of eigenvalue and eigenvector, we may write , implying that [55]. This gives

[TABLE]

Here, for periodic boundary condition the matrix ${\bm{A}}$ and ${\bm{B}}$ are as follows

[TABLE]

[TABLE]

Obviously, the matrix ,${\bm{A}}$ and ${\bm{B}}$ are really symmetrical, so the eigenvalue ${\lambda_{i}}$ is real number. Meanwhile, for arbitrary ${\bm{x}}\neq{\bm{0}}$, the right-hand side of the Eq. (80) is

[TABLE]

Using the inequality $2xy<{x^{2}}+{y^{2}}$, we obtain

[TABLE]

and the left-hand side of the Eq. (80) is

[TABLE]

Similarly, using the inequality $2xy<{x^{2}}+{y^{2}}$, we obtain

[TABLE]

From Eqs. (84) and (86), we can see that the right-hand side of the Eq. (80) is ${{\bm{x}}^{T}}{\bm{Ax}}>0$ and the left-hand side of the Eq. (80) is ${{\bm{x}}^{T}}{\bm{Bx}}<0$. Thus, ${\lambda_{i}}\leq 0$.

Lemma 2. Let ${\bm{W}}$ be an arbitrary square matrix. Then for any operator matrix norm $\left\|\cdot\right\|$, we obtain ${\lambda_{i}}(\bm{W})\leq\left\|\bm{W}\right\|$, where ${\lambda_{i}}(\bm{W})$ is the spectral radius of matrix ${\bm{W}}$ [56].

Theorem 1. The Taylor approximation of CFDS-PIM is unconditionally stability.

Proof. the Taylor series approximation of ${e^{{\lambda_{i}}\tau}}$ is defined as

[TABLE]

Using Lemma 1, we obtain ${\lambda_{i}}\tau\leq 0$, thus ${e^{{\lambda_{i}}\tau}}\leq 1$. For fourth-order Taylor series approximation of the PIM, ${e^{{\lambda_{i}}\tau}}<1$. We use Lemma 2 to get the spectral radius of matrix ${e^{{\bm{H}}\tau}}$ is less than 1 in fourth-order Taylor approximation of the PIM. Thus, the Taylor approximation of CFD-PIM scheme is unconditionally stability.

3.2.2 Amplification symbol

Definition 1. The rational approximation $R_{r,s}(z)$ to the exponential $e^{-z}$ is called A-acceptable when $\left|{{R_{r,s}}(-z)}\right|<1$ holds for all $-z$ with negative real part. The approximation is called L-acceptable when it is A-acceptable and it also satisfies $\left|{{R_{r,s}}(-z)}\right|\to 0$ as $\Re(-z)\to-\infty$.

In Fig. 1, we compare the behavior of $e^{-z}$ and Taylor approximation. It can be observed from the traces that Taylor approximation are A-acceptable.

Fig. 2 illustrates the traces of $e^{-z}$ and Taylor approximation for the different complex planes. Since the results of the functions are complex, we plot their real parts. It can be seen from the plots that Taylor approximation conforms to A-acceptable Definition 1.

3.2.3 Stability region

The stability of the CFDS-PIM can be observed from the plots of their stability region [11, 57].The linear ordinary differential equation (69) can be rewrite as

[TABLE]

We assume that a fixed point $u_{0}$ satisfying $cu_{0}=0$ exists, and $u$ is the perturbation of $u_{0}$. If $Re(c)<0$, then we can say the fixed point $u_{0}$ is stable. We denote $x=c\tau$, with $\tau$ being a single time step, and then apply Taylor approximation to Eq. (88). The amplification factors $r(x)$ of Taylor approximations can be calculated in the following way:

[TABLE]

Notice that we assumed $r(x)<1$ to obtain the stability region. Suppose that $x$ is complex. As can be seen in Fig. 3, the stability regions of proposed scheme is plotted. The axes of Fig. 3 are real and imaginary parts of $z$. It can be observed from Fig. 3 that the stability region of Taylor approximation are in good agreement with exponential approximation $e^{z}$.

4 The n-dimensional numerical method

Strang splitting method (SSM) is a numerical method for solving differential equations that are decompose multi-dimensional problems into a sum of differential operators. This method is named after Gilbert Strang. It is used to speed up the calculation for problems involving operators on very different time scales, and to solve the multi-dimensional PEDs by reducing them to a sum of one-dimensional problems.

4.1 Extensions to two-dimensional case

We consider the two-dimensional heat conduction equation (11). As a precursor to Strang splitting, Eq. (11) can be written as

[TABLE]

where ${\bm{H_{x}}}$ and ${\bm{H_{y}}}$ are difference operator in the $x$ and $y$ direction. The right side of Eq. (90) is already split into a sum $a+b$ of relatively simple expressions. Due to one of the properties of difference operator is the distributive law of multiplication, we obtain the following equations

[TABLE]

For Eq. (91), the analytical solution to the associated initial value problem would be

[TABLE]

This section focuses on how to calculate the exponential matrix $e^{\bm{({H_{x}}+{H_{y}})}t}$, and the calculation of $e^{\bm{({H_{x}}+{H_{y}})}t}$ is too complicated. Thus, we convert it into calculating the product of $e^{\bm{{H_{x}}}t}$ and $e^{\bm{{H_{y}}}t}$, but $e^{\bm{{H_{x}}}t}$ and $e^{\bm{{H_{y}}}t}$ must satisfy the commutativity of the addition theorem

[TABLE]

Besides, the exponentials of ${\bm{H_{x}}}$ and ${\bm{H_{y}}}$ are related to that of ${\bm{H_{x}+H_{y}}}$ by the Trotter product formula

[TABLE]

Gottleib et al. [36] suggested that the Trotter result can be used to approximated $e^{\bm{H}}$ by splitting $\bm{H}$ into $\bm{H_{x}+H_{y}}$, because $m=2^{20}$ is already very large that was proposed. Thus, we use the following approximation

[TABLE]

This approach to calculate ${e^{\bm{H}}}$ is of potential interest when the exponentials of ${\bm{H_{x}}}$ and ${\bm{H_{y}}}$ can be accurately and efficiently computed. If ${\bm{H_{x}}}$ and ${\bm{H_{y}}}$ commute, we rewrite Eq. (92) as follows

[TABLE]

Thus, the two-dimensional heat conduction equation becomes two one-dimensional problems. For each one-dimensional problem, it can be solved by the CFD-PIM scheme which introduced in Sec. 3.1.

4.2 Extensions to three-dimensional case

For the three-dimensional heat conduction equation, we can also use CFD-PIM based on the SSM(CFD-PIM-SSM) to decompose it into the sum of differential operators of three one-dimensional problems. The CFD-PIM scheme can be extended to three-dimensional case (12).

As a precursor to Strang splitting, we rewrite Eq. (12) as follows

[TABLE]

where ${\bm{H_{x}}}$, ${\bm{H_{y}}}$ and ${\bm{H_{z}}}$ are difference operators in the $x$-direction, $y$-direction, and $z$-direction, respectively. The right side of Eq. (97) is already split, which become a sum $a+b+c$. We obtain the following equations

[TABLE]

For Eq. (91), the analytical solution to the associated initial value problem would be

[TABLE]

If ${\bm{H_{x}}}$, ${\bm{H_{y}}}$ and ${\bm{H_{z}}}$ commute for Eq. (99), the analytical solution to the associated initial value problem would be

[TABLE]

Because we apply SSM to three-dimensional case, we obtain a sum of difference operator of three one-dimensional parabolic problem, the scheme has the same accuracy as one-dimensional cases. For each one-dimensional problem, it can be solved by the CFD-PIM scheme which introduced in Sec. 3.1.

5 Numerical Result

In this section, we give the six numerical examples to validate the adaptability of the proposed schemes and compare their accuracy with those which are already available in the literature for solving n-dimensional Burgers’ system. The accuracy of the schemes is measured in terms of $L_{2}$ errors, $L_{\infty}$ errors, computing time and the rate of convergence of the scheme. In our Tables, CPU(s) is computing time. The rate of convergence(ROC) of proposed schemes is defined as

[TABLE]

where ${L_{\infty}}(2h)$ and ${L_{\infty}}(h)$ are discrete maximum absolute errors at $2h$ and $h$. All the numerical experiments are conducted on MATLAB R2016a platforms based on an Intel Core i5-6300HQ 2.30 GHz processor.

Example 1. To verify the effectiveness of the modified algorithm, we test the one-dimensional Burgers’ equation (5) proposed in Sec. 2.2.1, over a domain $[0,1]$, with the initial and boundary conditions (26), according to the one-dimensional Hopf-Cole transformation (30), the analytical (Fourier series) solution of Eq. (5) is

[TABLE]

The numerical solutions are reported in Tables 1,2,3 and Fig. 4 for the different values of $t$ with $Re=100$ and $\tau{\rm{=}}1\times{10^{-5}}$. The results of the numerical solution are compared with those of Refs. [58, 44, 45, 59], and the numerical results of the proposed scheme are better than their results. It is evident that the proposed scheme has high accuracy and efficiency than other numerical schemes. Fig. 4 exhibit that as the increase of time the numerical solution of partial regions becomes steeper and steeper, and the decreasing rate of the approximate solution increases. This physical phenomenon validates the fact that the numerical solution is capable of describing the shock wave. For the numerical solutions at different time demonstrated in Fig. 4 are fantastically analogical as depicted in the figures given in Refs. [44, 45, 59, 60, 47]. According to the values of ROC in Table 3 ( ROC$>4$ ), the proposed scheme can be verified as a high order scheme

Example 2. To validate the order of convergence of the proposed scheme, a numerical experiment of coupled Burgers’ equation (1) was carried out in Example 2 with the region $x\in\Omega=[-\pi,\pi]$ with initial conditions

[TABLE]

extracted from the following exact solution given by Refs. [11, 61] for $\omega_{1}=\omega_{2}=1.0$, $\kappa_{1}=\kappa_{2}=-2.0$, $\delta_{1}=\delta_{2}=1.0$:

[TABLE]

The boundary conditions are extracted from the analytical solution(In this example, $\tau=4\times 10^{-4}$). The spatiotemporal evolution of the numerical solution is shown in the left part of Fig. 5. To intuitively observe the physical phenomena of the example, The left figure of Fig. 5 is depicted to visually compare the analytical solutions with the numerical solutions of $u(x,t)$ at different time $t$. It is observed that the numerical results show great agreement with the analytical solutions. The numerical results reflect the motion characteristics of wave propagation: the amplitude of the wave decreases with time while the wavelength remains unchanged. In addition, the errors of $L_{\infty}$, rates of convergence and CPU running time are listed in Table 4. As can be observed in Table 4, the accuracy and efficiency of our numerical scheme are much higher than that of Ref. [11] under the same spatial step size. It is observed that the error of $L_{\infty}$ becomes smaller as the mesh size is refined. The computer operating environment of the algorithm is worse than that of Ref. [11], and the computation time and convergence order of our numerical scheme are much better than that of Ref. [11], which shows that the proposed scheme has excellent adaptability. The proposed scheme presents more accurate and high-efficient solutions in spatial direction than the scheme in Refs .[11, 62].

Example 3. To test the applicability of the proposed scheme for two-dimensional problems, we consider the two-dimensional Burgers’ equations (6) over a square domain $[0,1]\times[0,1]$, with the initial conditions

[TABLE]

for which the analytical solutions [63] are

[TABLE]

[TABLE]

The boundary conditions are extracted from the analytical solution. This example uses the CFD scheme with sixth-order accuracy in space and the PIM with fourth-order accuracy in time. The numerical simulation results are shown in Table 5 and Figs. 6,7. To prove that the method is sixth-order accuracy in space, the time step $\tau$ is fixed to $5\times 10^{-4}$, thus the time truncation error can be ignored. To simplify the demonstration, we use the same grid size in the $x$ and $y$ directions. As can be seen from Table 5, when $h$ is reduced by a factor of $2$, the maximal errors for both $u(x,y,t)$ and $v(x,y,t)$ are reduced by a factor of $6$, which indicates that the method is sixth-order accurate in space. When $h$ is reduced by a factor of $2$, the $L_{\infty}$ of $u(x,y,t)$ and $v(x,y,t)$ are both reduced by a factor of $2^{6}$, which indicates that the method is sixth-order accurate in space. It is observed that the numerical solutions present excellent agreement with the analytical solutions. In large Reynolds numbers, the ROC of the proposed scheme has little influence, which means the scheme still is a high order.

Example 4. In order to verify the effectiveness of the improvement proposed above, we consider the system of the two-dimensional Burgers’ equations (6) proposed in Sec. 2.2.2, over a square domain $[0,1]\times[0,1]$, with the initial and boundary conditions (35) and the analytical solution

[TABLE]

The numerical and analytical solutions of two-dimensional examples are present in Tables 6,7 with $Re=100$ and $\tau=5\times 10^{-5}$. The numerical solutions of the different time are presented in Figs. 8,9 with $N\times N=81\times 81$. From the Tables 6,7 and Figs. 10 of numerical simulation results, it can be seen that the proposed numerical scheme has high accuracy under the condition of large Reynolds number ($Re=100,200$). It is observed that the numerical solutions show great agreement with the analytical solutions. The Figs. 8,9 exhibit that the numerical solution of partial regions becomes steeper and steeper as the increase of time. This physical phenomenon validates the fact that the numerical solution is capable of describing shock wave. Moreover, note that Tables 6,7 and Figs. 8,9 indicate the property (the boundary condition (38)) of the solution of the two-dimensional Burgers’ equation. The physical phenomena depicted in the Figs. 8,9 are analogical to those in Refs. [12, 48, 49, 50].

Example 5. In order to test the applicability of the proposed scheme for multi-dimensional problems, we consider the three-dimensional Burgers’ equations (7) over a domain $[-1,1]\times[-1,1]\times[-1,1]$, with the initial conditions

[TABLE]

The analytical solution for this problem is given by

[TABLE]

The boundary conditions are extracted from the analytical solution. It is observed that the numerical solutions present great agreement with the analytical solutions. The slices of the four-dimensional images are used for observing three-dimensional Burgers’ equations, which are depicted in Figs. 11,12,13. Figs. 14,15 present the numerical solutions and the errors between the exact solutions and the numerical solutions at $t=1$ with $z=0.25$ and $Re=100,1000$. The numerical results of this example show that the CFD-PIM-SMM scheme base on Hopf-Cole transformation is a numerical method with high precision and high efficiency for solving n-dimensional Burgers’ system. The comparison is done with solutions obtained by the numerical solutions and the exact solutions for the three-dimensional Burgers’ equations, which is presented in on left-hand side Figs. 14,15. All numerical figures of this example show that the proposed scheme has excellent accuracy and efficiency. Figs. 14,15 exhibit that the numerical solution of partial regions becomes steeper and steeper as $Re=100\longrightarrow 1000$. This physical phenomenon validates the fact that the numerical solution is capable of describing shock wave. Table 8 shows the numerical results of different $\kappa$ for Burgers’ equation. Table 8 mainly shows two phenomena: 1. With the increase of $\kappa$, the corresponding error will also increase; 2. Table 8 transformation does not cause ROC transformation; Phenomenon 2 is consistent with n-dimensional Hopf-Cole transformation (Eqs. (20) to (23)) which has no relationship between $\kappa$ and ROC.

Example 6. We consider the system of the three-dimensional Burgers’ equation (6) in Sec.2.2.3, over a square domain $[0,1]\times[0,1]\times[0,1]$, with the initial and boundary conditions (45,46), and the analytical solution

[TABLE]

In order to verify the validity and practicability of the proposed scheme of three-dimensional Burgers’ equation, the numerical solution is compared with the analytical solution. The numerical solution of the three-dimensional Burgers’ equation (12) is given by the three-dimensional Hopf-Cole transformation (50). The correctness of the numerical solution of the heat conduction equation ensures the correctness of the numerical solution of the Burgers’ equation. The higher the accuracy of numerical solution of the heat conduction equation, the higher the accuracy of the numerical solution of the Burgers’ equation will be. The numerical and analytical solutions of the three-dimensional heat conduction equation are presented in Table 9 with $Re=10$ and $\tau=5\times 10^{-5}$. From the Table 9 of numerical simulation results, it is observable that the numerical solutions present excellent consistency with the analytical solutions. The slices of the four-dimensional images are used for observing three-dimensional Burgers’ equation, which are depicted in Figs. 17,18,19,20,21. The numerical solutions of the different time are presented in Figs. 20,21,22 with $N\times N\times N=81\times 81\times 81$. It can be observed from Figs. 20,21,22 when $Re=10$, the numerical solution of the Burgers’ equation will generate shock wave. This physical phenomenon validates the fact that the numerical solution is capable of describing shock wave. The physical phenomena depicted in the Figs. 20,21,22 are similar to those in Ref.[12]. We make a three-dimensional Burgers’ equation with numerical solutions and errors in Fig. 23. It is observed that the numerical solution is very consistent with the analytical solution.

6 Conclusion

This paper proposed the modified Hopf-Cole transformation for n-dimensional Burgers’ system. After obtaining the n-dimensional heat conduction equations, we present a new high-order exponential time differencing precise integration method schemes in combination with the sixth-order compact finite difference scheme, which have been developed for the numerical solutions of n-dimensional heat conduction equations. The proposed scheme is tested on six examples, and the results are entirely satisfactory in comparison with the analytical solutions. The findings can be summarized as follows

For the examples of the ordinary initial and boundary conditions, the solutions can describe shock wave phenomena for large Reynolds numbers ($Re\geq 1000$), which is characterized by high precision and high efficiency. 2. 2.

For the complicated initial and boundary conditions problems, the proposed scheme obtains excellent accuracy and high efficiency in comparison to the other numerical techniques available in the Refs. [5, 11, 42, 44, 45, 46, 47, 59, 64, 58, 61]. 3. 3.

After the modification of Hopf-Cole transformation and optimization of computer programming, the present scheme has commendable adaptability and high efficiency in the calculation examples. 4. 4.

It is found that the proposed results are in good agreement with the analytical solutions for two-dimensional and three-dimensional problems. This scheme includes linear problems and nonlinear problems, which can be easily extended to solve model equations of high-dimensional problems.

Code

We are happy to share our research results with others. We have uploaded some example codes (Example 1, 4, and 5) on Github. This is Github’s link: https://github.com/LzEfreet/CFD-PIM.

Acknowledgement

This work is financially supported by the National Natural Science Foundation of China(No. 11826208).