Algorithms for coupled Burgers’ equations by sharing characteristic curves within BSLM

This paper introduces a new perspective of the traditional view on the velocity of each physical particle in the coupled Burgers’ equation in the backward semi-Lagrangian method (BSLM). The proposed methods reduce the number of Cauchy problems to be solved by observing a single virtual characteristic curve with a velocity. This can drastically reduce the computational cost of determining the departure point. Then, we solve the derived system reflected by the single virtual characteristic curve. Moreover, an efficient strategy for the derived linear system of equations is provided. Four examples are tested to demonstrate the adaptability and efficiency of the proposed method. The test results show that the proposed method has third- and fourth-order accuracy in time and space, respectively. In addition, compared with the existing method of solving the problem along two particles with different velocities, we confirm that the proposed method significantly reduces computational cost while maintaining accuracy well.


Introduction
Under the effect of gravity, the sedimentation of two types of particle concentrations in fluid suspensions and colloids is described in the following viscous coupled Burgers' equations [9]: with the following initial and boundary conditions: Coefficients α k and β k (k = 1, 2) are real constants that describe the interaction of gravity and thermal fluctuation (Péclet number) and the velocities of the moving particles, respectively [6].In addition, the positive real constants, ν 1 and ν 2 , represent the viscosity coefficients related to the Reynolds numbers of u and v, respectively.The rate at which a particular reaction proceeds mainly depends on the reactant concentration.Suppose the reactant concentration is increased.In that case, the number of molecules or ions will increase, thereby increasing the number of collisions between them and the reaction rate.
Thus, the solutions to the model problem simultaneously represent particle velocity and concentration.The governing equations are of advection-diffusion type, but the nonlinear coupling terms of u and v are not a common feature of this type.
Over the last 10 years, various numerical methods have been developed to solve the viscous coupled Burgers' equations.For example, Kapoor and Joshi developed a differential quadrature method with uniform algebraic trigonometric tension B-spline [12], Başhan applied a mixed method with the finite difference and differential quadrature method based on third-order modified cubic B-spline functions [5], Hussein and Kashkool used a weak Galerkin finite element method [10], Zhang et al. introduced an improved backward substitution method [21], Abdullah et al. developed a numerical procedure based on the cubic B-spline and the Hermite formula [1], Mohanty and Sharma presented a high accuracy two-level implicit method based on cubic spline approximation [17], Bhatt and Khaliq modified exponential time-differencing Runge-Kutta method based on Padé approximation [7], Kumar and Pandit proposed a composite scheme based on finite difference and Haar wavelets [13], Jiwari and Alshomrani developed a new collocation method based on modified cubic trigonometric B-spline function [11], Mohanty et al. proposed a two-level implicit compact operator method [16].Among the abovementioned methods, there is a backward semi-Lagrangian (BSLM) to solve the model problem.It solves the system on fixed Eulerian spatial grids at every time step by tracing in the reverse time direction the trace of a particle observed from a Lagrangian viewpoint.This BSLM is called a coupled backward semi-Lagrangian method (CBSLM), which is the first attempt at using the BSLM to solve the model equations [2].The CBSLM inherits the following well-known advantages of the traditional BSLM.The BSLM guarantees unconditional stability by implicitly dealing with the particle trajectories, allowing the use of a larger time step size than that in Eulerian methods [3,4,18].It does not require domain remeshing, unlike the Lagrangian method.The CBSLM has second-order accuracy in time and fourth-order accuracy in space, and it implicitly treats the particle trajectories without any iteration process.
In the previous BSLM task, CBSLM, the virtual characteristics of π(t) and φ(t) (or particle trajectories), respectively, satisfy the following nonlinear Cauchy problems: and the model equations are reformulated by the general form of the Lagrangian formulations: where In the CBSLM [2], the particle's velocities are fixed as The CBSLM process for simultaneously solving (3) and ( 4) (or model equations) can be divided into three parts: (i) finding the departure points at the previous time steps by solving the Cauchy problems; (ii) interpolating the function values at the departure points; and (iii) solving the systems to discretize (4).Among these processes, in most BSLMs, as well as in CBSLM, interpolation processes incur the highest computational cost.In CBSLM, interpolation is required even in the process of determining the departure points.In particular, a more significant number of departure points and interpolations are inevitably required to obtain high-order accuracy in time and space.Thus, we describe a high-order (third-order) accuracy version of CBSLM, CBSLM3, and propose a faster algorithm while maintaining the accuracy of CBSLM3.Hereafter, we denote the existing CBSLM as CBSLM2 to distinguish it from CBSLM3.
Our main goal is to develop a novel methodology, without losing accuracy, than both CBSLM2 and CBSLM3 that individually consider all characteristic curves of each particle.For this, we propose a new method for choosing the velocity functions.The proposed method observes a single virtual characteristic curve with velocity s.First, we propose the choice of three versions of s and then present their algorithms for solving the corresponding Lagrangian form.Additionally, to construct high-order algorithms, the third-order backward differentiation formula and fourth-order finite difference method are applied to discretize the Lagrangian formulations.We apply the third-order multistep error correction method to solve the nonlinear Cauchy problem.Note that the typical Lagrange interpolation is used for estimating the required function values at nongrid points in all proposed methods.In particular, we improve two parts among the three parts of BSLM by (i) handling the number of Cauchy problems, and (ii) providing an efficient strategy for solving a linear system of equations.To improve these parts, we reduce the number of the Cauchy problems to be solved by sharing the characteristic curve arriving at point x as only one.We additionally design an explicit formula that converts to a pentadiagonal one for the system of equations occurring at every time integration step.This formula is used to reduce the computational cost of solving the linear system of equations.
The remainder of this paper is organized as follows.In Sect.2, we review and extend CBSM2 for solving (1) and the Cauchy problem solver with third-order accuracy, which individually determines all characteristic curves of each particle.In Sect.3, we introduce three proposed schemes for solving the model problem and provide a fast solver for derived nonlinear systems.The numerical simulations are illustrated in Sect. 4 to demonstrate the efficiency of the proposed algorithms.Finally, the conclusions are summarized in the last section.

Preliminary results
In this section, we introduce CBSLM3, which is an extended version of CBSLM2 with a third-order temporal accuracy for solving the model equations.To do that, we first define time sequences t n := t 0 + nh, n = 0, 1, 2, . . ., N , where h := (Tt 0 )/N is the time step size.Furthermore, the uniformly divided spatial grid sequence is defined as x j = a + j x, x = (ba)/J.The time variable evaluated function u(t n , x) is denoted by superscript u n (x).

BSLM with BDF
Now, we provide the BSLM algorithm with third-order temporal accuracy to obtain the solutions at time t n+1 .The first step of BSLM is evaluating x = x j at t = t n+1 for the Lagrangian formulations (4) with (5).As mentioned in the Introduction, CBSLM chooses the fixed virtual velocities: Then, along with the characteristic curves π(t) and φ(t) satisfying nonlinear Cauchy problems (3) with initial values π j (t n+1 ) = x j and φ j (t n+1 ) = x j , the Lagrangian formulations (4) are given by To obtain the solutions u and v at (t, x) = (t n+1 , x j ), we employ the third-order backward differentiation formula for the total derivatives du/dt and dv/dt in (6).Then, we obtain the following system of nonlinear boundary value problems: where To discretize the nonlinear boundary value problems (7), we use the fourth-order finite difference for first and second partial derivatives with weight matrices D 1 and D 2 and vectors d 2,u and d 2,v given in the appendix.Then, we obtain the semidiscretized matrix system where Here, the notation diag(w) denotes a diagonal matrix, the diagonal entries of which are those of vector w and vector d 1,w is defined in the appendix.The hat notation on vector ŵ means a vector without both ends of column vector w.
Note that the departure points π j (t n-l ) and φ j (t n-l ), l = 0, 1, 2, do not usually coincide with the spatial grid point; consequently, the function values at these points u n-k (π j (t n-l )) and v n-k (φ j (t n-l )) require an appropriate interpolation.Herein, we use the typical Lagrange interpolation with sixth-order accuracy.The approximate vector, R n+1 κ , of r n+1 κ is defined as follows: where L is the interpolation operator using typical Lagrange polynomials of the sixth degree and is defined in the appendix.After applying the Lagrange interpolation to approximate r n+1 κ and dropping the truncation error, we achieved a fully-discretized system for the approximation U n+1 and V n+1 of u n+1 and v n+1 , respectively, namely with Vn+1 ex = 3 Vn -3 Vn-1 + Vn-2 .The final process for implementing this system involves the determination of the locations of departure points.This process is essential not only for CBSLM3 but also for the proposed methods.

A Cauchy problem solver finding departure points
In this section, we introduce a strategy for solving the Cauchy problems to complete the discretized system (12).Hereafter, we use the notation θ j to represent a characteristic curve π j or φ j .Recall that the position θ j (t n-l ) := θ (t n+1 , x j ; t n-l ) (l = 0, 1, 2) of the particle arriving at x j (j = 1, . . ., J -1) on time level t = t n+1 is referred to as the departure point.The trajectory of the departure point θ j (t) satisfies the Cauchy problem of Herein, the particle velocity s is combined with the solutions, u and v of the model problem, which is strong nonlinear.To resolve this nonlinearity, we employ the multistep error correction method introduced in [2].The mechanism of the error correction method starts with defining a linear local approximation where S n (x j ) is an approximation for s(t n , x j ).The local approximation leads a perturbation ψ(t) := θ j (t)y(t) to be corrected, which satisfies the following asymptotic initial value problem: with the initial value ψ(t n+1 ) = 0. Here, S n is a column vector defined as S n := [S n (x 0 ), S n (x 1 ), . . ., S n (x J )] T and L x is introduced in the appendix as the derivative operator of L.
To estimate numerical solutions for the asymptotic initial value problem ( 16), the Radau IIA-type Runge-Kutta method with third-order accuracy is expressed as presented below The above is applied to the problem, inducing the following linear system: where Finally, the definition of the perturbation ψ allows us to obtain the approximation θ n-l j of θ j (t n-l ) by Remark 1 To find the approximation of departure points π n j and φ n j in (11), two linear systems of the form (17) are solved for π and φ, respectively, in CBSLM3.CBSLM3 requires 2 × J Jacobian calculations and 6 × J function evaluations at each time step to construct these systems.The proposed schemes, introduced in the next section, require exactly half of the interpolation process of CBSLM3 to construct the system by sharing a single virtual characteristic curve.In addition, the computational cost to solve this system is also reduced by half.

Sharing characteristic curve schemes for solving the model problem
In the introduction, we have described the general form of the Lagrangian formulations (4) with ( 5), for coupled Burgers' equations.The Lagrangian formulations depend on how the velocity functions s 1 and s 2 are set.This section presents a sharing characteristic curve semi-Lagrangian (SC-SL) scheme by choosing the velocity functions.The proposed method observes a single virtual characteristic curve θ j (t) with velocity s.Here, velocity s is chosen by considering the simplicity of the resulting Lagrangian formulations, (4) with (5).That is, the velocity s is chosen to reduce the computational cost by selecting a single virtual characteristic curve and to minimize the number of terms to be approximated in (5).We propose the choice of three versions of s in the following subsection and then present their algorithms for solving the Lagrangian formulations accordingly.

Velocity type 1: SC-SL 1
The first choice of the velocity s(t, θ j (t)) along the single characteristic curve θ j (t) for each grid x j is defined by Based on this choice, the Lagrangian formulations can be defined as with where γ 12 := α 1β 2 and γ 21 := α 2β 1 .After evaluating t = t n+1 in (21), the third-order backward differentiation formula is applied to total derivatives in (21).From the fact that θ j (t n+1 ) = x j , we can obtain the following system of nonlinear boundary value problems: To discretize ( 22), the fourth-order finite difference weight matrices D 1 and D 2 , defined in the appendix, are used for first and second partial derivatives.Then, the semidiscretized matrix system for (22) is obtained by where To approximate r n+1 κ (x j ), we use the Lagrange interpolation and drop the truncation errors in (23).Then, the fully-discretized system is obtained by with

Velocity type 2: SC-SL 2
For each x j , the second choice of the velocity s(t, Then, the Lagrangian formulations are given by with Similar to the process of deriving (24) from ( 21), the following fully-discretized system can be obtained from (26): where (28)

Velocity type 3: SC-SL 3
The last choice of the velocity s(t, θ j (t)) is Then, the Lagrangian formulations are given by with Similar to the process of deriving ( 24) from ( 21) in the previous subsection, we can obtain the following fully-discretized system from (30): where Remark 2 We set the virtual particle velocity s according to the amplitudes of the initial values of u and v.For example, if the amplitude of u is larger, then we choose SC-SL 3 .In the opposite case, we can define the velocity s(t, Then, Ûn+1 , Vn+1 , and the coefficients of the system (31) are naturally exchanged in the reversed order.

Solver of the nearly pentadiagonal system
In this subsection, we discuss a strategy for solving the nonlinear system of equations ( 9), ( 23), (31), and (27).The A 1 (•) and A 2 (•) matrices of these systems have different elements but the same size and shape.Representing those systems, we define the following system: where the coefficient matrix and loading vector are given as The matrix A is named a nearly pentadiagonal matrix for convenience.Now, we design an explicit formula for converting the nearly pentadiagonal system (33) to a fully pentadiagonal system, allowing it to be solved using an advanced Thomas algorithm [8].First, we reduce the bandwidth of system (33) by 2. The formula is based on applying a part of Gaussian elimination to p 4 , p 5 , q 4 , and q 5 .Then, we can easily reduce the nearly pentadiagonal system (33) to a pentadiagonal system as follows: where Here, the quantities pk , pk , qk , and qk , k = 1, 2, 3, 4, are defined as (37)

Numerical results and discussions
To verify the accuracy and computational costs of the proposed schemes, we test four coupled Burgers' systems.In particular, we will show the adaptability of the proposed method by dealing with cases where the two solutions u and v are the same and different.For this, we record the maximum error E ∞ (w) given by where w(t N , x j ) and W N j represent the analytical and numerical solutions at (t, x) = (t N , x j ), respectively.To assess the efficiency of the numerical scheme, we count the computational time in seconds for solving the model problem over t ∈ [0, t N ], denoted by CPU.

Example 1
The first example to be solved is the coupled Burgers' equations given by The above problem admits the following analytical traveling wave solutions [15]: where κ 0 is a real constant that determines the initial slope and the advection velocity of the solution.The initial and Dirichlet boundary conditions can be recovered from the above analytical solutions.First, to observe the adaptability of the proposed method to various initial conditions and different sharpness of the traveling wave, we plot analytical and numerical solutions of u during t ≤ 5.0 for different values of κ 0 in Fig. 1.The numerical solutions are obtained by SC-SL 3 .For three cases κ 0 = 0.5, 1.0, and 2.0, the resolutions and spatial domains (h, x, ) are employed as (0.1, 0.1, [-20, 20]), (0.05, 0.05, [-10, 30]), and (0.005, 0.005, [0, 40]), respectively.The figure shows that the wave has a faster moving speed and a sharper gradient as κ 0 is larger.It can be observed that the numerical solutions obtained by the proposed SC-SL 3 matches the analytical solutions well.
To examine the temporal and spatial convergence rates of the proposed methods, we first measure the maximum errors E ∞ (u) and E ∞ (v).The errors are computed with κ 0 = 0.5 at t = 1.0 on = [-20, 20].The temporal convergence rates are obtained by varying time step size h = 2 -(k+3) (k = 1, 2, 3, 4) with a fixed spatial grid size x = 1/60.The spatial convergence rates are also obtained by varying x = 1/2 k (k = 1, 2, 3, 4) with a fixed h = 0.005.The results are illustrated in Fig. 2 for u, where the black dotted lines denote the reference slope for intuitive comparison confirms.Similar results can be obtained for v; our results for v have been omitted.In the figures, all proposed methods have third-order temporal and fourth-order spatial convergence rates.However, SC-SL 1 has the lowest accuracy than other methods, whereas the other methods have similar accuracy.To investigate the accuracy of the CBSLM3 and SC-SLs, we compare the errors of numerical solutions obtained by the proposed methods and other methods [12,14] including CBSLM2 [2].In Table 1, the errors are listed at different time levels for two κ 0 = 0.1 and 0.5 on = [-20, 20], where temporal step size and spatial resolution are employed as (h, J) = (1/1000, 320).Here, (h, J) = (1/1000, 320) is set the same as those used in the works of Kapoor et al. [12] and Lai et al. [14], and the results are used as those recorded in the literature [12,14].In this table, the CBSLM3 and SC-SLs are about three and four digits more accurate than those obtained by Kapoor [12] and Lai [14].SC-SLs (s = 2, 3) are 2 and 10 times more accurate for κ 0 = 0.1 and 0.5, respectively, compared to those in CBLMS2.In addition, the accuracy of CBLSM3 is maintained, despite the use of a single characteristic curve.However, SC-SL 1 shows slightly insufficient accuracy, similar to the result of Fig. 2.

Example 2
For the second example, we consider the model equations on (-π, π) with the analytical solutions [1,2,5,15,20,21]  The initial and Dirichlet boundary conditions are imposed by the above analytical solutions.The solution in this example has the property of gradual diffusion over time.In Fig. 3, we plot the behaviors of analytical and numerical solutions for u using SC-SL 3 at different time levels from t = 0.0 to t = 3.0 with an interval of 0.5.According to this figure, the solution to this problem begins with a sine function and decays over time, and the SC-SL 3 works well in this case.Moreover, the numerical solution U matches u well at different time levels.
To investigate the numerical accuracy of the proposed SC-SLs, we compare the errors obtained by SC-SLs and other methods introduced in [2,5,16,21].The methods in [16,21] were compared only for the time levels at which the results were valid.The numerical results are obtained with spatial resolutions J = 100 when h = 0.01 and J = 130 when h = 0.001, except J = 190 for Başhan [5] with h = 0.01.Tables 2-3 report the results at different time levels t ≤ 10.0.In particular, CPUs for CBSLM3, SC-SLs, and CBSLM2 [2] are also reported in Tables 2-3.The tables indicate that SC-SL 2 and SC-SL 3 achieved at least one-to two-digit higher accuracy than the other methods [2,5,16,21] in the entire time intervals.Moreover, SC-SL 2 and SC-SL 3 can maintain the accuracy of CBSLM3 at a lower computational cost.In contrast, the proposed SC-SL 1 shows poor accuracy compared to CBSLM3, with the results being similar to those of Example 1.As shown Tables 2-3, as time t increases, the error differences between SC-SL 1 and other SC-SLs (s = 2, 3) decrease because this is a special situation in which the solution scale converges to 0. According to the results of Examples 1 and 2, the results of SC-SL 1 will be excluded from the comparison.We notice that the main observation point of the results obtained using the proposed method is the accuracy of the numerical results.Therefore, an important concern is whether the accuracy of the existing method, which does not share even if the virtual characteristic curve is shared, can be well maintained.Hence, we will examine the adaptability and numerical performance of the proposed method even for solutions with different behaviors in the following two examples.

Example 3
As the third example, the coupled Burgers' equation on the spatial domain is as follows: where coefficients, α k and β k are arbitrary real constants.In this experiment, we take the analytical solution as follows [19]: where and κ 0 is a nonzero constant.The initial and Dirichlet boundary conditions are imposed by the above analytical solutions.
To verify the effect of parameters and the adaptability of the proposed methods, we plot the behaviors of the numerical solutions in Figs.4-5.The parameters considered herein are wavelet speed, viscosity, the difference in magnitude between the two solutions, and nonlinearity (coupling).The numerical experiments are tested according to viscosity ν = 1.0, 0.1, and 0.001.We use the parameters (α 1 , α 2 , β 1 , β 2 ) = (2.0,2.0, 0.2, 0.8) and employ the spatial domains = [-10, 15], [-10, 10], and [-10, 10], considering the behaviors of the solutions.Fig. 4 depicts the numerical solution obtained by SC-SL 3 and compares it with the exact solution.The figure illustrates that the smaller the viscosity, the sharper the gradient of the solution.The proposed SC-SL 3 is observed to derive a fairly accurate solution even when ν becomes small.In addition, we observe numerical solutions for two  5, when κ 0 is positive, the wave propagation direction is to the left, and when κ 0 is negative, the wave propagation direction is to the right.Therefore, the sign of κ 0 determines the advection direction of the solution.Finally, the ratio of β 2 and β 1 represents the difference between the magnitudes of u and v.

Example 4
To investigate the proposed SC-SL 3 , we consider a model with shock initial conditions in our last example: This model has an initial condition as and the homogeneous boundary condition.Now, we examine the adaptability of the SC-SL 3 when solutions u and v move in opposite directions.In this case, the solutions u and v have strong nonlinearity with each other.In Figs.6-8, we plot the behaviors of the numerical solutions using CBSLM3 and SC-SL 3 at different time levels from t = 0.0 to t = 1.0 with an interval of 0.2.Here, we note that the numerical solutions of u and v move in opposite directions, i.e., to the right and left, respectively.In addition, the CPUs for each method are described in the legends of the figures.In Fig. 4, we use the parameters (α 1 , α 2 , β 1 , β 2 , ν) = (2, -2, 0.2, -0.8, 0.01) and resolution (h, x) = (1/150, 200).According to Fig. 6, the numerical solution of SC-SL 3 has the same performance as the numerical solution of CBLSM3 with less computational cost compared.Fig. 7 is the result of the simulation conducted to verify the performance for the quite small viscosity ν = 0.0005.The parameters and resolutions are set by (α 1 , α 2 , β 1 , β 2 ) = (2, -2, 0.2, -0.8) and (h, x) = (1/500, 1/800), respectively.As can be observed in Fig. 7, the numerical solutions of SC-SL 3 have a good agreement with those of CBSLM3, even if the solutions have a sharp gradient.Thus, the efficiency of SC-SL 3 can also be observed in this simulation.
Finally, to simulate the case of a significant difference in the scales of u and v, we use the following parameters and resolutions (α 1 , α 2 , β 1 , β 2 , ν) = (2, -2, 2, 8, 0.01) and (h, x) = (1/600, 1/600), respectively.In Fig. 8, the result of SC-SL 3 exhibits stable solution behavior that follows that of CBSLM3 well.Furthermore, SC-SL 3 requires fewer CPUs than CBSLM3 to obtain similar results in terms of accuracy as in the previous test.

Conclusion
This paper introduced three fast BSLMs for solving the viscous coupled Burgers' equations.The proposed methods simultaneously handle two Cauchy problems by observing a single virtual characteristic curve with a velocity.This can dramatically reduce the interpolation process and the computational cost of determining the departure point.We used four examples to assess the adaptability and efficiency of the proposed methods.In particular, we proved that SC-SL 3 displays a lower cost consumption than CBSLM3 while maintaining a similar accuracy to CBSLM3 for various types of solutions.This finding is expected to provide a reference about solving various other coupled equations with a new perspective that deviates from the traditional view of the velocity for each physical particle in the Burgers' equation coupled with the BSLM.
Consequently, the matrix form of the finite difference method for first-order spatial derivatives is given by where Next, the finite difference method for second order spatial derivative of w(x) at interior points x i (i = 2, . . ., M -2) is given by w xx (x i ) = 1 12 x 2 -w(x i-2 ) + 16w(x i-1 ) -30w(x i ) + 16w(x i+1 ) -w(x i+2 ) + O x 4 .
The finite difference method for a second-order spatial derivative also has the following matrix form: where .

Figure 2
Figure 2The temporal (left) and spatial (right) convergence rates of u for Example 1 with κ 0 = 0.5 at t = 1.0

Figure 3
Figure 3 Comparison of the numerical and analytic solutions obtained by SC-SL 3

Table 1
Error comparison of u at different time level with parameters h = 1/1000, J = 320 and

Table 2
Comparison of E ∞ (u, t) for Example 2 with (h, J) at different times t

Table 3
Comparison of E ∞ (u, t) for Example 2 with (h, J) at different times t