Numerical solution of a coupled Burgers' equation via barycentric interpolation collocation method

1.   Introduction

Nonlinear partial differential equations are widely encountered in various scientific and engineering fields, including mathematical biology, plasma physics, chemical physics, and fluid dynamics. Burgers' equation, a prominent partial differential equation, is extensively utilized in modeling viscous fluid flow, shock waves, shallow water waves, and traffic flow, etc. ^[1,2]. It also serves a benchmark in testing various numerical algorithms. The resolution of this equation is crucial for the development of mathematical models and achieving a comprehensive understanding of physical phenomena. Numerous mathematicians and researchers are exploring a variety of methods, including not only numerical techniques ^[3,4], such as the differential quadrature method ^[5] and the finite-difference method ^[6], but also analytical approaches ^[7,8,9,10], such as the decomposition method ^[11] and the homotopy perturbation method^[12], to deepen their understanding of the dynamics and behavior of the model. For example, Kumar and Pandit ^[4] introduced a composite numerical format for solving the one-dimensional ($1$D) coupled Burgers' equation by finite difference methods and Haar wavelets. Soliman ^[7] applied the variational iteration method to derive accurate solutions for both the Burgers' equation and the coupled Burgers' equation, while Kaya ^[11] obtained exact solutions for the coupled Burgers' equation through the Adomian decomposition method. The coupled Burgers' equation is often used to simulate various coupled transport phenomena. However, due to its nonlinear nature, analytical solution for the coupled Burgers' equation can not be obtained easily. As a result, numerical solutions for this type of equation remain a vital method for investigating its dynamic behavior.

In this paper, we examine a two-dimensional ($2$D) coupled nonlinear Burgers' equation^[13,14]

where $q_{t} = {\partial q}/{\partial t}$, $q_{x} = {\partial q}/{\partial x}$, $q_{y} = {\partial q}/{\partial y}$, $q_{xx} = {\partial^{2} q}/{\partial x^{2}}$, and $q_{yy} = {\partial^{2} q}/{\partial y^{2}}$ ($q = u, v$). Furthermore, the variables $u$ and $v$ denote the velocity components that require determination, and $r = 1/R$, where $R$ is the Reynolds number. The initial conditions of Eq (1.1) is denoted as

where $D = [a, b]\times [c, d]$, $f_{0}, g_{0}$ are given functions. In addition, the boundary conditions of Eq (1.1) are denoted as

where $t > 0$, $\partial D$ denotes the boundary of region $D$, and $f, g$ are all known functions. This coupled viscous Burgers' equation offers an intuitive and straightforward model for investigating the sedimentation or scaling volume concentration evolution of two types of particles in fluid suspensions or colloids under the influence of gravitational forces.

The analysis of a coupled Burgers' equation is instrumental in elucidating the convection–diffusion phenomenon, which is influenced by the viscosity coefficients associated with the nonlinear convection and diffusion terms. Hopf ^[15] and Cole ^[16] have demonstrated that the $1$D Burgers' equation can be effectively addressed by transforming it into a linear homogeneous heat equation, with solutions represented through Fourier series. In ^[17], the exact solution of the $2$D coupled Burgers' equation could generally be obtained only under specific initial and boundary conditions. Furthermore, the inclusion of nonlinear convection terms and viscous parameters complicates the numerical analysis of the viscous Burgers' equation, making it a consistently challenging task. In the past few years, several significant numerical approaches have been developed for the coupled Burgers' equations. For instance, Chai and Ouyang ^[18] explored methods for reducing shock oscillations, including the streamline upwind/Petrov-Galerkin method and other related techniques. The authors of ^[19] introduced a weak Galerkin finite element method for solving the $2$D coupled Burgers' equation. A novel fourth-order compact difference scheme is proposed in ^[20], along with its application to solving the mixed-type time-fractional Burgers' equation. The work in ^[21] presented a new numerical method for addressing the coupled Burgers' equation, utilizing the collocation finite element method with cubic trigonometric and quintic B-splines as approximation functions. The performance and accuracy of this method are demonstrated through the calculation of the error norms $L_{2}$ and $L_{\infty}$. Another study ^[22] examined the application of the RBF-QR method in obtaining numerical solutions for the $1$D coupled viscous Burgers' equation. Lastly, Park et al. ^[23] developed an algorithm designed to reduce the computational cost associated with the backward semi-Lagrangian method (BSLM) for solving nonlinear convection-diffusion equations.

In recent years, the barycentric interpolation method has gained prominence due to its robust node adaptability and excellent numerical stability for both special and equidistant node distributions. Consequently, a substantial body of literature has employed this method to address various numerical challenges. Notable applications include Volterra integro-differential equations ^[24], nonlinear high-dimensional Fredholm integral equations ^[25], the heat conduction equation ^[26], plane elasticity problems ^[27], and semi-infinite domain problems ^[28], and so on. Furthermore, Li and Wang ^[29] and Wang and Li ^[30] have successfully applied the barycentric interpolation method to a wide range of practical engineering issues, including dynamic diffusion problems, vibration problems, and beam structural mechanics, as documented in their monographs.

The high computational efficiency, numerical stability, flexibility, and ease of implementation of the barycentric interpolation collocation method have been demonstrated. When appropriate nodes, such as Chebyshev nodes, are selected, it produces highly accurate results. Some researchers have extended the method to address $1$D Burgers' equations and modified Burgers' equations ^[30,31]. However, the numerical analysis and solutions for $2$D coupled Burgers' equations remain relatively limited. Therefore, in this article, we aim to apply the barycentric interpolation collocation method to solve $2$D coupled Burgers' equations, with the goal of further exploring its effectiveness and feasibility in numerically solving such coupled equations.

The primary objective of this paper was to develop a novel differential matrix scheme for the $2$D coupled Burgers' equation. A direct linearization approach is employed to address the nonlinear terms, and a barycentric interpolation method is introduced to resolve the differentiation matrix formulation for this coupled Burgers' equation. The computational results are compared with analytical solutions across various values of the Reynolds number $R$ in selected examples. The validity and feasibility of this method are demonstrated through several illustrative cases.

The following sections are organized as follows: Section 2 presents the application of the barycentric interpolation method to develop the numerical scheme for the coupled Burgers' equation. The validity and feasibility of this method are demonstrated through several illustrative examples in Section 3. Finally, Section 4 offers conclusions and outlines potential directions for future research.

2.   The numerical scheme for the $2$D coupled Burgers' equations

In this section, we utilize the barycentric interpolation method to aid in the resolution of the differentiation matrix formulation for the $2$D coupled Burgers' equations.

2.1.   Linearization for the $2$D coupled Burgers' equation

In the following, the direct linearization method is used to transform the $2$D coupled Burgers' equation into a linear equation. Namely, the initial value of nonlinear term $uu_{x}+vu_{y}$ and $uv_{x}+vv_{y}$ is changed to $u^{(0)}u_{x}+v^{(0)}u_{y}$ and $u^{(0)}v_{x}+v^{(0)}v_{y}$. In this case, we have

By substituting $u^{(0)}$ and $v^{(0)}$ into Eq (2.1), we can calculate $u$ and $v$, denoting them as $u^{(1)}$ and $v^{(1)}$. Thus, by substituting $u ^{(1)}$ and $v^{(1)}$ as the new $u^{(0)}$ and $v^{(0)}$ in Eq (2.1), we can compute new values for $u$ and $v$, which we denote as $u^{(2)}$ and $v^{(2)}$. Similarly, by substituting $u^{(s-1)}$ and $v^{(s-1)}$ in Eq (2.1), we can calculate new values for $u$ and $v$, denoted as $u ^{(s)}$ and $v ^{(s)}$, respectively. Thus, given the iterative representation $u \rightarrow u^{(s)}$, $u^{(0)} \rightarrow u^{(s-1)}$, $v \rightarrow v^{(s)}$, and $v^{(0)} \rightarrow v^{(s-1)}$ ($s \in N^{+}$ represents the number of iterations), we can construct a direct linearization scheme for Eq (2.1)

Here, $u^{(s-1)}$ and $v^{(s-1)}$ denote the values of $u$ and $v$ after the $s-1$ iterations, while $\frac{\partial (u^{(s)})}{\partial t}$ and $\frac{\partial (v^{(s)})}{\partial t}$ represent the partial derivatives of $u$ and $v$ with respect to $t$ after $s$ iterations, respectively. Similarly, $\frac{\partial (u^{(s)})}{\partial x}$ and $\frac{\partial (v^{(s)})}{\partial x}$ represent the partial derivatives of $u$ and $v$ with respect to $x$ after $s$ iterations, $\frac{\partial (u^{(s)})}{\partial y}$ and $\frac{\partial (v^{(s)})}{\partial y}$ represent the partial derivatives of $u$ and $v$ with respect to $y$ after $s$ iterations, $\frac{\partial^{2} (u^{(s)})}{\partial x^{2}}$ and $\frac{\partial^{2} (v^{(s)})}{\partial x^{2}}$ represent the second-order partial derivatives of $u$ and $v$ with respect to $x$ after $s$ iterations, and $\frac{\partial^{2} (u^{(s)})}{\partial y^{2}}$, $\frac{\partial^{2} (v^{(s)})}{\partial y^{2}}$ represent the second-order partial derivatives of $u$ and $v$ with respect to $y$ after $s$ iterations, respectively.

2.2.   Barycentric interpolation method

A barycentric interpolation function can be employed to deduce the approximate function of $u(x)$ for a set of $m$ distinct interpolation nodes $x_{i}$, $i = 1, 2, \cdots, m$, i.e.,

and

is the interpolation basis function with the interpolation weight ^[29,30]

and $J_{i} = \{s\in I_{x}:i-d^{*}\leq s\leq i\}$, $0\leq d^{*}\leq m$, $I_{x} = \{0, 1, 2, \cdots, m\}$, so the barycentric rational interpolation is constructed. In the following discussion, the variable $d^{*}$ in the $x, y$, and $t$ directions is denoted as $d^{*}_{x}, d^{*}_{y}$ and $d^{*}_{t}$, respectively. For the other case

we get the barycentric Lagrange interpolation. By observing the weights, we can notice their direct association with dispersed nodes. The weights provided are applicable to any node distribution. If commonly used node distributions such as equispaced nodes or Chebyshev nodes are selected, the computation of the barycentric weights becomes simpler. The equidistant nodes is

and its weight function is $\omega_{i} = (-1)^{m-i-1}C^{i-1}_{m-1}$, and $C^{i-1}_{m-1} = \frac{(m-1)!}{(m-i)!(i-1)!}$. The Chebyshev nodes is

and we then find that its weight function is

It is also worth noting that the Chebyshev nodes in the interval of $-1$ to $1$ can be transformed to any interval from $a$ to $b$ through a linear transformation $\tilde{x} = \frac{b-a}{2}x+\frac{b+a}{2}$.

Solving differential equations requires the derivatives of Eq (2.3). The $p$-th derivative at $x_{i}$ is expressed as follows:

where the first-order derivative value and $p$-th-order derivative value of $l_{i}(x)$ at the $h$-th interpolation node are ^[29,30]

Going a step further, if we substitute the discrete interpolation node $x_{h}$, $h = 1, 2, \cdots, m$ into Eq (2.3), we can obtain the discrete form as follows:

Furthermore, the matrix forms of Eq (2.7) can be expressed as

where $\boldsymbol{u} = [u_{1}, u_{2}, \cdots, u_{m}]^T$, $\boldsymbol{u}^{(p)} = [u^{(p)}_{1}, u^{(p)}_{2}, \cdots, u^{(p)}_{m}]^T$ and

In addition, $\boldsymbol{C}^{(0)} = \boldsymbol{I}_{m}$, and $\boldsymbol{I}_{m}$ is an $m$-th-order identity matrix. The detailed derivation process of Eq (2.8) can be seen in ^[29].

2.3.   Differential matrices for the coupled Burgers' equation

In the following, we illustrate the process of establishing the differential matrix for the $2$D coupled Burgers' equation. Let us consider the functions $u(x, y, t)$ and $v(x, y, t)$ defined in a spatial-temporal domain $[a, b]\times [c, d] \times (0, T]$. The spatial-temporal domain is divided by $m \times n \times l$ interpolation nodes, and the divisions of interpolation nodes in different directions are represented as

These interpolation nodes $(x_{i}, y_{j}, t_{k})$ across the entire computing area $[a, b]\times [c, d] \times (0, T]$ can be uniform or non-uniform collocation points.

Making use of the barycentric interpolation function, we obtained the approximation function for $u(x, y, t)$ and $v(x, y, t)$, which are expressed as

where the symbols $u_{ijk}$ and $v_{ijk}$ represent the function values $u(x_{i}, y_{j}, t_{k})$ and $v(x_{i}, y_{j}, t_{k})$, respectively. Additionally, the symbols $l_{qx_{i}}(x)$, $l_{qy_{j}}(y)$, and $l_{qt_{k}}(t) \; (q = u, v)$ stand for the barycentric interpolation basis functions in the directions of $x$, $y$, and $t$, respectively, and the representations of $l_{qx_{i}}(x)$, $l_{qy_{j}}(y)$, and $l_{qt_{k}}(t) \; (q = u, v)$ are as follows:

where $\omega_{x_{i}}, \omega_{y_{j}}$, and $\omega_{t_{k}}$ are the corresponding weight functions.

Furthermore, we can derive the $(p_{1}+p_{2}+p_{3})$-order partial derivative values of $u(x, y, t)$ and $v(x, y, t)$ at the interpolation nodes $\{ (x_{i}, y_{j}, t_{k}), \, \, i = 1, 2, \cdots, m; \, j = 1, 2, \cdots, n; \, k = 1, 2, \cdots, l\}$ and present them as

Thus, Eqs (2.12) and (2.13) can be expressed in matrix form as

where $\otimes$ denotes the Kronecher product of the matrix; $\boldsymbol{C}^{(p_{1}00)}_{q}$, $\boldsymbol{C}^{(0p_{2}0)}_{q}$, and $\boldsymbol{C}^{(00p_{3})}_{q} \; (q = u, v)$ are differential matrices corresponding to $x, y$, and $t$, respectively. Their expressions are

We thus obtain the matrix $\boldsymbol{C}^{(p_{1}00)}_{q}\otimes\boldsymbol{C}^{(0p_{2}0)}_{q}\otimes\boldsymbol{C}^{(00p_{3})}_{q}$ of the order $(m \times n \times l)\times (m \times n \times l)$. In addition, for $x$, $\boldsymbol{C}^{(000)}_{q} = \boldsymbol{I}_{m}$; for $y$, $\boldsymbol{C}^{(000)}_{q} = \boldsymbol{I}_{n}$; for $t$, $\boldsymbol{C}^{(000)}_{q} = \boldsymbol{I}_{l}$, and $\boldsymbol{I}_{m}, \boldsymbol{I}_{n}$, and $\boldsymbol{I}_{l}$ are all identity matrices.

Going a step further, Eq (2.2) can be written in the format of a barycentric interpolation differential matrix

where, $\boldsymbol{D}^{(001)}_{q} = \boldsymbol{I}_{m}\otimes\boldsymbol{I}_{n}\otimes\boldsymbol{C}^{(001)}_{q}$, $\boldsymbol{D}^{(100)}_{q} = \boldsymbol{C}^{(100)}_{q}\otimes\boldsymbol{I}_{n}\otimes\boldsymbol{I}_{l}$, $\boldsymbol{D}^{(010)}_{q} = \boldsymbol{I}_{m}\otimes\boldsymbol{C}^{(010)}_{q}\otimes\boldsymbol{I}_{l}$, $\boldsymbol{D}^{(200)}_{q} = \boldsymbol{C}^{(200)}_{q}\otimes\boldsymbol{I}_{n}\otimes\boldsymbol{I}_{l}$, $\boldsymbol{D}^{(020)}_{q} = \boldsymbol{I}_{m}\otimes\boldsymbol{C}^{(020)}_{q}\otimes\boldsymbol{I}_{l}$, $q = u, v$, and $\boldsymbol{0}$ is a zero column vector. In this case, Eq (2.19) can be expressed as a block matrix expression, i.e.,

where,

and $\boldsymbol{O}$ is a zero matrix.

Typically, the replacement method and the addition method are frequently utilized for the implementation of the initial-boundary conditions, as detailed in ^[30]. In this article, we employ the replacement method to tackle the initial and boundary conditions. After applying the boundary conditions to Eq (2.20), numerical computations can begin by setting an initial condition of $s = 1$ on the left side of Eq (2.20). As a result, the numerical solutions $\boldsymbol{u}^{(s)}$ and $\boldsymbol{v}^{(s)}$ can be computed iteratively.

3.   Numerical experiments and discussion

In this section, we present a series of computational experiments designed to demonstrate the accuracy and effectiveness of the proposed methodology.

Within this context, we define the subsequent error standards $E^{p}_{q\infty}$, $E^{p}_{qr}$, $(q = u, v), (p = L, R)$ to illustrate the accuracy of the proposed methodology, i.e., the error norms

Among these, $\boldsymbol{q}^{c}$ and $\boldsymbol{q}^{e}$ denotes the numerical and analytical solutions, respectively, and $\boldsymbol{q} = \boldsymbol{u, v}$. Specifically, the parameter $p$ is defined as $L$ in relation to barycentric Lagrange interpolation and as $R$ in relation to barycentric rational interpolation. In all subsequent tables, we denote the absolute error under equidistant nodes and Chebyshev nodes as $E^{pd}_{q\infty}$ and $E^{pc}_{q\infty}$, respectively, where $q = u, v$ and $p = L, R$. Similarly, the relative error under equidistant nodes and Chebyshev nodes is represented as $E^{pd}_{qr}$ and $E^{pc}_{qr}$, respectively, with the same definitions for $q$ and $p$.

Example 1. Eq (1.1) with exact solutions

$(x, y)\in [0, 1]\times [0, 1]$, $t\geq 0$ is under consideration in this example, and its initial–boundary conditions can be established on the basis of the exact solution.

In Example 1, the numerical solutions obtained from the Chebyshev and equidistant node calculations are illustrated in Figures 1 and 2. Furthermore, we present the error distribution calculated by barycentric interpolation, with the parameters set to $m = n = l = 18$, $R = 1000$, $t = 1$, and $d^{*}_{x} = d^{*}_{y} = d^{*}_{t} = 6$ in Figures 3 and 4. The variations in the absolute and relative errors associated with equidistant nodes and Chebyshev nodes, determined using different quantities of nodes, are summarized in Tables 1-4. Additionally, Table 5 offers a comprehensive overview of the errors calculated for Chebyshev nodes under varying values of $R$ with $m = n = l = 18$, $d^{*}_{x} = d^{*}_{y} = d^{*}_{t} = 6$, and $t = 1$. This example illustrates that our algorithm maintains a high level of computational accuracy, even for larger values of $R$.

Additionally, to quantitatively assess the convergence order ($C_o$) of the numerical scheme, the formula $C_o = {\log\frac{E_{i}}{E_{i+1}}}/{\log\frac{h_{i}}{h_{i+1}}}$ is applied to analyze the data presented in Table 3, where $E_{i}$ represents the numerical error for the $i$-th mesh resolution, and $h_{i}$ denotes the maximum stride length of the partition corresponding to the variables in the $i$-th mesh. The results indicate that the calculation of this method is convergent. The CPU times are also given in Table 4, which represent the computational efficiency of this method, and the results indicate that the computational cost of this method is lower than that of other methods.

Example 2. We consider Eq (1.1) with exact conditions

$(x, y)\in [0, 1]\times [0, 1]$ and $t\geq 0$. Furthermore, the initial and boundary conditions are delineated in accordance with the exact solutions.

The error distribution associated with the barycentric Lagrange collocation method, calculated using Chebyshev nodes with the parameters $m = n = l = 18$ and varying values of $R$, is illustrated in Figures 5–7 for Example 2. The numerical solutions obtained for $R = 10, 20, 50,200$ demonstrate a high degree of accuracy, making it challenging to differentiate between the exact solutions and the numerical approximations in the figures, particularly in Figures 5 and 6. In this context, the proposed algorithm exhibits effective performance for smaller values of $R$, as indicated by the numerical results. However, its effectiveness appears to diminish for larger values of $R$. It is plausible that despite the initial data being smooth, a solution shock may develop as time progresses. In the event of a solution shock, it is possible to adjust the number of computational nodes to manage the shock and enhance the accuracy of the algorithm. To further illustrate the accuracy, the numerical and exact solutions for $R = 500$ at $t = 1$, with $(x, y)\in [0, 0.5]\times [0, 0.5]$, are presented in Figure 8.

Example 3. Here, the $2$D Burgers' equation ^[14]

is considered with the exact condition

Similarly, the initial condition and boundary conditions are all specified by the given exact solution. We can bring in a complementary function $u(x, y, t) = v(x, y, t)$, which can transform Eq (3.1) with only a single equation and a single unknown function $v(x, y, t)$ into a coupled Burgers' equation, analogous to Eq (1.1), and then apply the method proposed in this paper. The corresponding differential matrices of Burgers' equation (3.1) can be expressed as

The numerical results obtained through the barycentric rational interpolation method, utilizing various values of $R$, are presented in Figure 9. As $R$ increases, the solutions exhibit steeper gradients. For example, when $R = 40, 50$, the numerical solutions illustrated in Figure 9(c), (d) exhibit steeper gradients compared with those shown in Figure 9(a), (b). Additionally, the region where the error exceeds $10^{-3}$ is confined to a small area near the line defined by $x+y = t$. To further elucidate the long-term behavior of this scheme, we present the numerical solutions and the corresponding errors for various values of $t$ in Figure 10. It is apparent that the region characterized by a steep slope shifts as $t$ increases, reflecting the behavior of the exact solutions depicted in Figure 11.

Example 4. We consider the Eq (3.1) with the exact solution^[32]

over a domain $(x, y)\in [0, 1]\times [0, 1]$ and $t\geq 0$, with the initial and boundary conditions specified by the given exact solution.

In contrast to Example 3, the present example demonstrates that even when the value of $R$ is relatively large and the number of computing nodes is configured to $m = n = l = 8$, the computational accuracy of this method remains consistently stable, achieving an accuracy level of up to $10 \times 10^{-6}$. To facilitate a clearer understanding of the differences between the numerical and exact solutions, we present both the exact and numerical solutions for $R = 50,100,500, 1000$ at $t = 1$ in Figure 12. It is evident that, despite the relatively large value of $R$, the solutions remain smooth and devoid of shocks. The minimal discrepancies observed between the exact and numerical solutions further underscore the high accuracy of the method.

4.   Conclusions

This paper presents an efficient methodology for the numerical solution of a $2$D coupled Burgers' equation. Specifically, the barycentric interpolation collocation method is employed to approximate the unknown function, while a specially constructed differential matrix is utilized to derive the numerical solution of the $2$D coupled Burgers' equation. The results from numerical experiments demonstrate the effectiveness and accuracy of the proposed scheme in solving various instances of Burgers' equations. Furthermore, it is shown that the current algorithm can yield satisfactory numerical results for larger parameter values of $R$, specifically for $R \geq 1000$. It is important to acknowledge that the method proposed in this article has certain limitations and challenges when calculating the coupled equations in irregular regions. In the future, we will continue to explore the application of this method for solving the equations defined in these irregular regions and will conduct further theoretical analyses of coupled equations, including convergence analysis. Additionally, we plan to extend this method to tackle equations that have not been investigated in other fields, as well as coupled equations of other different categories.

Use of AI tools declaration

The authors declare they have not used Artificial Intelligence (AI) tools in the creation of this article.

Acknowledgments

This work was supported by the Natural Science Foundation of Shandong Province (Grant No. ZR2022MA003).

Conflict of interest

The authors declare there is no conflicts of interest.