# A Composite Algorithm for Numerical Solutions of Two-Dimensional Coupled Burgers’ Equations

Authors: Vikas Kumar; Sukhveer Singh; Mehmet Emir Koksal
Journal: Journal of Mathematics
 (2021)

---

## Abstract

In this study, a new composite algorithm with the help of the finite difference and the modified cubic trigonometric B-spline differential quadrature method is developed. The developed method was applied to two-dimensional coupled Burgers’ equation with initial and Dirichlet boundary conditions for computational modeling. The established algorithm is better than the traditional differential quadrature algorithm proposed in literature due to more smoothness of cubic trigonometric B-spline functions. In the development of the algorithm, the first step is semidiscretization in time with the forward finite difference method. Furthermore, the obtained system is fully discretized by the modified cubic trigonometric B-spline differential quadrature method. Finally, we obtain coupled Lyapunov systems of linear equations, which are analyzed by the MATLAB solver for the system. Moreover, comparative study of these solutions with the numerical and exact solutions which are appeared in the literature is also discussed. Finally, it is found that there is good suitability between exact solutions and numerical solutions obtained by the developed composite algorithm. The technique can be extended for various multidimensional Burgers’ equations after some modifications.


---

## Body

## 1. Introduction

In this paper, the authors considered the following dimensionless form of two-dimensional (2D) coupled Burgers’ equation:(1a)1Re∇2u=∂u∂t+u∂u∂x+v∂u∂y,(1b)1Re∇2v=∂v∂t+u∂v∂x+v∂v∂y,with initial conditions (ICs),(2a)ux,y,0=ψ1x,y,x,y∈α,β×γ,δ,(2b)vx,y,0=ψ2x,y,x,y∈α,β×γ,δ,and Dirichlet boundary conditions (BCs),(3)uα,y,t=h1y,t,uβ,y,t=h2y,t,ux,γ,t=h3x,t,ux,δ,t=h4x,t,where ∇2=∂2/∂x2+∂2/∂y2 is Laplace operator, ux,y,t and vx,y,t are velocity components to be determined. Also, hi,i=1,2,…,4 are known smooth functions and Re=ρu⟶L/μ is the Reynolds number with density ρ, viscosity μ, characteristic length L, and u⟶=u,vT.The nonlinear convection-diffusion model is simply represented by Burgers’ equation [1]. This famous equation describes the flow theory through a shockwave moving in viscous liquid [2], phenomena of turbulence [3], and various other kinds of phenomena in aerodynamics.Due to its extensive scope of applicability, various numerical schemes have been constructed to study its numerical solutions. Moreover, due to its application in various fields of science and technology, researchers and scientists are still interested in developing algorithms to find their numerical and exact solutions. A great number of works has been studied for finding approximate solutions of Burgers’ equation, for example, cubic spline method [4], finite element and difference methods [5–8], multilevel alternating direction implicit schemes [9], and various explicit and implicit methods [10, 11]. Furthermore, the decomposition method [12], spectral method [13], Chebyshev collocation method [14], and local discontinuous radial basis function collocation method [15, 16] are investigated in literature. Also, Haar wavelet quasilinearization approach [17] and differential quadrature methods (DQMs) [18–23] have been developed. In recent years, new meshless methods [24–26] for various types of Burgers’ equations have been developed.In Lagrange interpolation-based DQMs [18–22], Lagrange’s fundamentals are used to compute the weighting coefficients. In these cases, as the number of grid points increases the weights become unstable. Herein, to reduce this instability, the modified cubic trigonometric B-spline functions are used to the weighting coefficients of DQMs.In this article, a new numerical algorithm is developed based on the finite difference and the modified cubic trigonometric B-spline (CTBS) DQMs for approximate solutions of coupled two-dimensional Burgers’ equations’ weighting coefficients (WCs) of DQM are calculated by using the modified CTBS functions as test functions which are different from the conventional technique of Lagrange interpolation [27]. Some well-known test problems are worked out to inspect the correctness and competence of the planned approach. The techniques lead to correct results with insignificant L∞,RMS and L2 errors.


## 2. Differential Quadrature Method

Recently, DQMs have become popular for solving nonlinear partial differential equations (PDEs) arising in nonlinear phenomena. DQMs discretize the first and second derivatives over 1D domainΩ=α,β as follows:(4)uxxi,t=∑j=1Nαij1uxj,t,i=1,2,…,N,(5)uxxxi,t=∑j=1Nαij2uxj,t,j=1,2,…,N,where αij1 and αij2 are unknown coefficients weighting the first and second derivatives, respectively, and xi,i=1,2,…,N, are uniform grids as well as nonuniform grids that exist in the domain. Bellman et al. [28] introduced two approaches to calculate WCs. Furthermore, to modify Bellman’s approaches for finding WCs, many efforts have been carried out such as Lagrange interpolated cosine functions, spline functions, Legendre polynomials, Lagrange interpolation polynomials, and radial basis functions (see [19, 29–35] and the references therein) to determine these coefficients. In this study, we determine WCs with the use of CTBS functions after some modifications.


### 2.1. Cubic Trigonometric B-Spline Functions

In this section, we mesh the solution domainα≤x≤β into N subintervals xi,xi+1,i=0,1,…,N−1 with the help of knots xi such that α=x0<x1,…,<xN=β is a uniform partition with step length a=xi+1−xi=β−α/N,i=0,1,…,N−1.Now, the piecewise CTBS basis functionsτBjx over the uniform mesh are defined as follows [36, 37]:(6)Bjx=1ωr3xi,x∈xi−2,xi−1,rxirxisxi+2+sxi+3rxi+1+sxi+4r2xi+1,x∈xi−1,xi,sxi+4rxi+1sxi+3+sxi+4rxi+2+rxis2xi+3,x∈xi,xi+1,s3xi+4,x∈xi+1,xi+2,0,otherwise,where(7)rxi=sinx−xi2,sxi=sinxi−x2,ω=sina2sinasin3a2.The basis over the regionα≤x≤β is formed by the set(8)B−1x,B0x,…,BNx,BN+1x.Every CTBS covers four elements. Now, with the help of Table1, we have tabulated the values of Bjx and its derivatives as follows:Table 1
Coefficients of CTBSs and their derivatives at knotsxj.
xxj−2xj−1xjxj+1xj+2Bjx0α1α2α10Bj′x0α30α40Bj″x0α5α6α50(9)α1=sin2a/2sinasin3a/2,α2=21+2cosa,α3=−34sin3a/2,α4=34sin3a/2,α5=31+3cosa16sin2a/22cosa/2+cos3a/2,α6=3cos2a/22sin2a/21+2cosa.


### 2.2. Modified Cubic Trigonometric B-Spline Functions

In this work, we compute WCs of DQM with the help of modified CTBS function defined in (6) as follows:(10)B˜0x=B0x+2B−1x,j=0,B˜1x=B1x−B−1x,j=1,B˜jx=Bjx,j=2,3,…,N−2,B˜N−1x=BN−1x−BN+1x,j=N−1,B˜Nx=BNx+2BN+1x,j=N.It is worth mentioning that the modified functionsτB˜jx,j=0,1,…,N, are linearly independent. On the solution domain α,β, these functions create a family of basis functions.


### 2.3. Weighting Coefficients for Modified Cubic Trigonometric B-Spline Differential Quadrature Method

Now, substitute the modified functionsτB˜jx,forj=0,1,…,N, into equation (4). The matrix form of the equation is as follows:(11)Aw01=B0Aw11=B1⋮AwN−11=BN−1AwN1=BN,where A is N+1×N+1 coefficient matrix:(12)A=α2+2α1α10α2α1α1α2α1.........α1α20α1α2+2α1.The matrixwk1=αk01,αk11,…,αkN1T,fork=0,1,…,N, and Bk at xk,fork=0,1,…,N, are as follows:(13)B0=−2α42α400⋮00,B1=−α40α40⋮00,…,BN−1=000⋮−α40α4,BN=000⋮0−2α42α4.Furthermore, with the help of Thomas algorithm WCs,αij1 are achieved as solutions of tridiagonal systems of equation (11). Similarly, with the help of the above method, it is easy to calculate second-order WCs βij2.


### 2.4. Two-Dimensional Modified Cubic Trigonometric B-Spline Differential Quadrature Method

In order to apply this method to 2D nonlinear problems, first of all, decompose the domainΩ=x,y:α1≤x≤β1;α2≤y≤β2π/6 as Ω1=xi,yj,i=1,2,…,N;j=1,2,…,M by adopting step length Δx=xi−xi−1 and Δy=yj−yj−1 in x and y direction, respectively. This modified technique helps to estimate the 1st order partial derivatives of ux,y,t at a point as follows:(14)uxxi,yj,t=∑k=1Nαik1uxk,yj,t,i=1,2,…,N,(15)uyxi,yj,t=∑k=1Mβjk1uxi,yk,t,j=1,2,…,M,where αij1 is WCs for the 1st order derivatives w.r.t. x. Similarly, βjk1 are coefficients w.r.t. y.In order to compute the 2D WCs, we can define the functionsτB˜jy,j=0,1,…,N, as in equation (10). Furthermore, take the test functions as Tijx,y=τB˜ixτB˜jy. Now, with the help of the axioms of vector space and substituting the value of Tijx,y into equations (14) and (15), we have(16)B˜j′xi=∑k=1Nαik1B˜jxk,j,i=1,2,…,N,B˜j′yi=∑k=1Mβik1B˜jyk,j,i=1,2,…,M.Furthermore, applying the well-known algorithm “Thomas algorithm” and proceeding with the same methods as in the case of equation (11), the solutions of the systems give the value of αij1 and βjk1. In 2D case, the WCs in higher-order derivatives can be considered as follows:(17)αijr=rαij1αiir−1−αijr−1xi−xj,fori,j=1,2,…,N;i≠j;r=2,3,…,N−1,αiir=−∑j=1,j≠iNαijr,fori=j,βijr=rβij1βiir−1−βijr−1yi−yj,fori,j=1,2,…,N;i≠j;r=2,3,…,N−1,βiir=−∑j=1,j≠iMβijr,fori=j,where αijr and βijr are WCs for rth order partial derivatives w.r.t. x and y, respectively.


## 2.1. Cubic Trigonometric B-Spline Functions

In this section, we mesh the solution domainα≤x≤β into N subintervals xi,xi+1,i=0,1,…,N−1 with the help of knots xi such that α=x0<x1,…,<xN=β is a uniform partition with step length a=xi+1−xi=β−α/N,i=0,1,…,N−1.Now, the piecewise CTBS basis functionsτBjx over the uniform mesh are defined as follows [36, 37]:(6)Bjx=1ωr3xi,x∈xi−2,xi−1,rxirxisxi+2+sxi+3rxi+1+sxi+4r2xi+1,x∈xi−1,xi,sxi+4rxi+1sxi+3+sxi+4rxi+2+rxis2xi+3,x∈xi,xi+1,s3xi+4,x∈xi+1,xi+2,0,otherwise,where(7)rxi=sinx−xi2,sxi=sinxi−x2,ω=sina2sinasin3a2.The basis over the regionα≤x≤β is formed by the set(8)B−1x,B0x,…,BNx,BN+1x.Every CTBS covers four elements. Now, with the help of Table1, we have tabulated the values of Bjx and its derivatives as follows:Table 1
Coefficients of CTBSs and their derivatives at knotsxj.
xxj−2xj−1xjxj+1xj+2Bjx0α1α2α10Bj′x0α30α40Bj″x0α5α6α50(9)α1=sin2a/2sinasin3a/2,α2=21+2cosa,α3=−34sin3a/2,α4=34sin3a/2,α5=31+3cosa16sin2a/22cosa/2+cos3a/2,α6=3cos2a/22sin2a/21+2cosa.


## 2.2. Modified Cubic Trigonometric B-Spline Functions

In this work, we compute WCs of DQM with the help of modified CTBS function defined in (6) as follows:(10)B˜0x=B0x+2B−1x,j=0,B˜1x=B1x−B−1x,j=1,B˜jx=Bjx,j=2,3,…,N−2,B˜N−1x=BN−1x−BN+1x,j=N−1,B˜Nx=BNx+2BN+1x,j=N.It is worth mentioning that the modified functionsτB˜jx,j=0,1,…,N, are linearly independent. On the solution domain α,β, these functions create a family of basis functions.


## 2.3. Weighting Coefficients for Modified Cubic Trigonometric B-Spline Differential Quadrature Method

Now, substitute the modified functionsτB˜jx,forj=0,1,…,N, into equation (4). The matrix form of the equation is as follows:(11)Aw01=B0Aw11=B1⋮AwN−11=BN−1AwN1=BN,where A is N+1×N+1 coefficient matrix:(12)A=α2+2α1α10α2α1α1α2α1.........α1α20α1α2+2α1.The matrixwk1=αk01,αk11,…,αkN1T,fork=0,1,…,N, and Bk at xk,fork=0,1,…,N, are as follows:(13)B0=−2α42α400⋮00,B1=−α40α40⋮00,…,BN−1=000⋮−α40α4,BN=000⋮0−2α42α4.Furthermore, with the help of Thomas algorithm WCs,αij1 are achieved as solutions of tridiagonal systems of equation (11). Similarly, with the help of the above method, it is easy to calculate second-order WCs βij2.


## 2.4. Two-Dimensional Modified Cubic Trigonometric B-Spline Differential Quadrature Method

In order to apply this method to 2D nonlinear problems, first of all, decompose the domainΩ=x,y:α1≤x≤β1;α2≤y≤β2π/6 as Ω1=xi,yj,i=1,2,…,N;j=1,2,…,M by adopting step length Δx=xi−xi−1 and Δy=yj−yj−1 in x and y direction, respectively. This modified technique helps to estimate the 1st order partial derivatives of ux,y,t at a point as follows:(14)uxxi,yj,t=∑k=1Nαik1uxk,yj,t,i=1,2,…,N,(15)uyxi,yj,t=∑k=1Mβjk1uxi,yk,t,j=1,2,…,M,where αij1 is WCs for the 1st order derivatives w.r.t. x. Similarly, βjk1 are coefficients w.r.t. y.In order to compute the 2D WCs, we can define the functionsτB˜jy,j=0,1,…,N, as in equation (10). Furthermore, take the test functions as Tijx,y=τB˜ixτB˜jy. Now, with the help of the axioms of vector space and substituting the value of Tijx,y into equations (14) and (15), we have(16)B˜j′xi=∑k=1Nαik1B˜jxk,j,i=1,2,…,N,B˜j′yi=∑k=1Mβik1B˜jyk,j,i=1,2,…,M.Furthermore, applying the well-known algorithm “Thomas algorithm” and proceeding with the same methods as in the case of equation (11), the solutions of the systems give the value of αij1 and βjk1. In 2D case, the WCs in higher-order derivatives can be considered as follows:(17)αijr=rαij1αiir−1−αijr−1xi−xj,fori,j=1,2,…,N;i≠j;r=2,3,…,N−1,αiir=−∑j=1,j≠iNαijr,fori=j,βijr=rβij1βiir−1−βijr−1yi−yj,fori,j=1,2,…,N;i≠j;r=2,3,…,N−1,βiir=−∑j=1,j≠iMβijr,fori=j,where αijr and βijr are WCs for rth order partial derivatives w.r.t. x and y, respectively.


## 3. Numerical Algorithm for Two-Dimensional Coupled Burgers’ Equation

In this section, the numerical algorithm is developed in the following sections.


### 3.1. Semidiscretization in Time

Applying forward difference on time derivatives and weighted average on spatial derivatives, we have(18a)ux,y,tn+1−ux,y,tnΔt=1Reθ∇2ux,y,tn+1+1+θ∇2ux,y,tn−uuxn+1−uuyn+1,n=0,1,…,K,(18b)vx,y,tn+1−vx,y,tnΔt=1Reθ∇2vx,y,tn+1+1+θ∇2vx,y,tn−vvxn+1−vvyn+1,n=0,1,…,K,where ux,y,tn+1=ux,y,t+nΔt,vx,y,tn+1=vx,y,t+nΔt,Δt step length in time direction, and 0≤θ≤1. The nonlinear term is linearized in the following manner:(19a)uuxn+1=unuxn,vuyn+1=vnuyn,(19b)uvxn+1=unvxn,vvyn+1=vnvyn,with ICs(20a)ux,y,t0=ψ1x,y,x,y∈α,β×γ,δ,(20b)vx,y,t0=ψ2x,y,x,y∈α,β×γ,δ,and prescribed BCs (3).After simplification, equations (18a) and (18b) can be written as follows:(21a)un+1−ΔtReθ∇2un+1=1−Δtuxnun+Δt1−θΔ2un−Δtvnuyn,(21b)vn+1−ΔtReθ∇2vn+1=1−Δtvynvn+Δt1−θΔ2vn−Δtunvxn,which is a system of second-order differential equations, where un+1x,y=ux,y,tn+1 and equations (21a) and (21b) are a system of second-order differential equations.


### 3.2. Fully Discretization in Space

In this section, spatial derivatives that occur in equations (21a) and (21b) are discretized by modified CTBS DQM over the given domain. After spatial discretization, equations (21a) and (21b) convert into a system of linear equations for each n in the following form:(22a)uijn+1−ΔtθRe∑k=1Nαik2ukjn+1+∑k=1Mβjk2uikn+1=1−Δt∑k=1Nβik1ukjnuijn+Δt1−θ∑k=1Nαik2ukjn+∑k=1Mβjk2uikn−Δtuijn∑k=1Mβjk1uikn,(22b)vijn+1−ΔtθRe∑k=1Nαik2vkjn+1+∑k=1Mβjk2vikn+1=1−Δt∑k=1Nαik1vkjnvijn+Δt1−θ∑k=1Nαik2vkjn+∑k=1Mβjk2vikn−Δtvijn∑k=1Mαjk1vikn,where uijn=unxi,yi and αik2andβik2 are WCs of 2nd order partial derivatives w.r.t. x and y.


### 3.3. Implementation of Dirichlet Boundary Conditions

The Dirichlet BCs given in equation (3) as α,y,t=h1y,t,uβ,y,t=h2y,t,ux,γ,t=h3x,t,andux,δ,t=h4x,t can be implemented directly as follows:(23a)u1,j=h1yj,t,uN,j=h2yj,t,yj∈γ,δ,j=1,2,…,M,(23b)ui,1=h3xi,t,ui,M=h4xi,t,xi∈α,β,i=1,2,…,N.As a result of applying the BCs on systems (23a) and (23b), the system can be written as follows:(24a)uijn+1−ΔtθRe∑k=1Nαjk2ukjn+1+∑k=1Mβjk2uikn+1=Sij,(24b)vijn+1−ΔtθRe∑k=1Nαik2ukjn+1+∑k=1Mβjk2vikn+1=Tij,Sij=1−Δt∑k=1Nαjk1ukjnuijn+Δt1−θ∑k=1Nαik2ukjn+∑k=1Mβjk2uikn−Δtuijn∑k=1Mβjk1uikn,Tij=1−Δt∑k=1Mβik1vkjnvijn+Δt1−θ∑k=1Nαik2vkjn+∑k=1Mβjk2vikn−Δtuijn∑k=1Nαik1vkjn.The system of equations (18a) and (18b) is a Lyapunov system of equations of the form(25a)A1U+UB1+C1=0,(25b)A2U+UB2+C2=0,where(26)A1=11⋱11N−2×N−2−ΔtθReα222α232⋯α2N−22α2N−12α322α332⋯α3N−22α3N−12⋮⋮⋯⋮⋮αN−222αN−232⋯αN−2N−22αN−2N−12αN−122αN−232⋯αN−1N−22αN−1N−12N−2×N−2,B1=−ΔtθRβ222β232⋯β2M−22β2M−12β322β332⋯β3M−22β3M−12⋮⋮⋯⋮⋮βM−222βM−232⋯βM−2M−22βM−2M−12βM−122βM−232⋯βM−1M−22βM−1M−12M−2×M−2,C1=S22S23⋯S2M−2S2M−1S32S33⋯S3M−2S3M−1⋮⋮⋯⋮⋮SN−22SN−23⋯SN−2M−2SN−2M−1SN−12SN−13⋯SN−1M−2SN−1M−1N−2×M−2.Also,(27)A2=11⋱11N−2×N−2−ΔtθReα222α232⋯α2N−22α2N−12α322α332⋯α3N−22α3N−12⋮⋮⋯⋮⋮αN−222αN−232⋯αN−2N−22αN−2N−12αN−122αN−232⋯αN−1N−22αN−1N−12N−2×N−2,V=v22n+1v23n+1⋯v2M−2n+1v2M−1n+1v32n+1v33n+1⋯v3M−2n+1v3M−1n+1⋮⋮⋯⋮⋮vN−22n+1vN−23n+1⋯bN−2M−2n+1bN−2M−1n+1vN−12n+1bN−23n+1⋯bN−1M−2n+1vN−1M−1n+1N−2×M−2,B2=−ΔtθRβ222β232⋯β2M−22β2M−12β322β332⋯β3M−22β3M−12⋮⋮⋯⋮⋮βM−222βM−232⋯βM−2M−22βM−2M−12βM−122βM−232⋯βM−1M−22βM−1M−12M−2×M−2,C2=−T22T23⋯T2M−2T2M−1T32T33⋯T3M−2T3M−1⋮⋮⋯⋮⋮TN−22TN−23⋯TN−2M−2TN−2M−1TN−12TN−13⋯TN−1M−2TN−1M−1N−2×M−2.Equations (19a) and (19b) are the coupled Lyapunov system, first solved for n=0 and then solved simultaneously for n=1,2,…,K by developing code in MATLAB 7.


## 3.1. Semidiscretization in Time

Applying forward difference on time derivatives and weighted average on spatial derivatives, we have(18a)ux,y,tn+1−ux,y,tnΔt=1Reθ∇2ux,y,tn+1+1+θ∇2ux,y,tn−uuxn+1−uuyn+1,n=0,1,…,K,(18b)vx,y,tn+1−vx,y,tnΔt=1Reθ∇2vx,y,tn+1+1+θ∇2vx,y,tn−vvxn+1−vvyn+1,n=0,1,…,K,where ux,y,tn+1=ux,y,t+nΔt,vx,y,tn+1=vx,y,t+nΔt,Δt step length in time direction, and 0≤θ≤1. The nonlinear term is linearized in the following manner:(19a)uuxn+1=unuxn,vuyn+1=vnuyn,(19b)uvxn+1=unvxn,vvyn+1=vnvyn,with ICs(20a)ux,y,t0=ψ1x,y,x,y∈α,β×γ,δ,(20b)vx,y,t0=ψ2x,y,x,y∈α,β×γ,δ,and prescribed BCs (3).After simplification, equations (18a) and (18b) can be written as follows:(21a)un+1−ΔtReθ∇2un+1=1−Δtuxnun+Δt1−θΔ2un−Δtvnuyn,(21b)vn+1−ΔtReθ∇2vn+1=1−Δtvynvn+Δt1−θΔ2vn−Δtunvxn,which is a system of second-order differential equations, where un+1x,y=ux,y,tn+1 and equations (21a) and (21b) are a system of second-order differential equations.


## 3.2. Fully Discretization in Space

In this section, spatial derivatives that occur in equations (21a) and (21b) are discretized by modified CTBS DQM over the given domain. After spatial discretization, equations (21a) and (21b) convert into a system of linear equations for each n in the following form:(22a)uijn+1−ΔtθRe∑k=1Nαik2ukjn+1+∑k=1Mβjk2uikn+1=1−Δt∑k=1Nβik1ukjnuijn+Δt1−θ∑k=1Nαik2ukjn+∑k=1Mβjk2uikn−Δtuijn∑k=1Mβjk1uikn,(22b)vijn+1−ΔtθRe∑k=1Nαik2vkjn+1+∑k=1Mβjk2vikn+1=1−Δt∑k=1Nαik1vkjnvijn+Δt1−θ∑k=1Nαik2vkjn+∑k=1Mβjk2vikn−Δtvijn∑k=1Mαjk1vikn,where uijn=unxi,yi and αik2andβik2 are WCs of 2nd order partial derivatives w.r.t. x and y.


## 3.3. Implementation of Dirichlet Boundary Conditions

The Dirichlet BCs given in equation (3) as α,y,t=h1y,t,uβ,y,t=h2y,t,ux,γ,t=h3x,t,andux,δ,t=h4x,t can be implemented directly as follows:(23a)u1,j=h1yj,t,uN,j=h2yj,t,yj∈γ,δ,j=1,2,…,M,(23b)ui,1=h3xi,t,ui,M=h4xi,t,xi∈α,β,i=1,2,…,N.As a result of applying the BCs on systems (23a) and (23b), the system can be written as follows:(24a)uijn+1−ΔtθRe∑k=1Nαjk2ukjn+1+∑k=1Mβjk2uikn+1=Sij,(24b)vijn+1−ΔtθRe∑k=1Nαik2ukjn+1+∑k=1Mβjk2vikn+1=Tij,Sij=1−Δt∑k=1Nαjk1ukjnuijn+Δt1−θ∑k=1Nαik2ukjn+∑k=1Mβjk2uikn−Δtuijn∑k=1Mβjk1uikn,Tij=1−Δt∑k=1Mβik1vkjnvijn+Δt1−θ∑k=1Nαik2vkjn+∑k=1Mβjk2vikn−Δtuijn∑k=1Nαik1vkjn.The system of equations (18a) and (18b) is a Lyapunov system of equations of the form(25a)A1U+UB1+C1=0,(25b)A2U+UB2+C2=0,where(26)A1=11⋱11N−2×N−2−ΔtθReα222α232⋯α2N−22α2N−12α322α332⋯α3N−22α3N−12⋮⋮⋯⋮⋮αN−222αN−232⋯αN−2N−22αN−2N−12αN−122αN−232⋯αN−1N−22αN−1N−12N−2×N−2,B1=−ΔtθRβ222β232⋯β2M−22β2M−12β322β332⋯β3M−22β3M−12⋮⋮⋯⋮⋮βM−222βM−232⋯βM−2M−22βM−2M−12βM−122βM−232⋯βM−1M−22βM−1M−12M−2×M−2,C1=S22S23⋯S2M−2S2M−1S32S33⋯S3M−2S3M−1⋮⋮⋯⋮⋮SN−22SN−23⋯SN−2M−2SN−2M−1SN−12SN−13⋯SN−1M−2SN−1M−1N−2×M−2.Also,(27)A2=11⋱11N−2×N−2−ΔtθReα222α232⋯α2N−22α2N−12α322α332⋯α3N−22α3N−12⋮⋮⋯⋮⋮αN−222αN−232⋯αN−2N−22αN−2N−12αN−122αN−232⋯αN−1N−22αN−1N−12N−2×N−2,V=v22n+1v23n+1⋯v2M−2n+1v2M−1n+1v32n+1v33n+1⋯v3M−2n+1v3M−1n+1⋮⋮⋯⋮⋮vN−22n+1vN−23n+1⋯bN−2M−2n+1bN−2M−1n+1vN−12n+1bN−23n+1⋯bN−1M−2n+1vN−1M−1n+1N−2×M−2,B2=−ΔtθRβ222β232⋯β2M−22β2M−12β322β332⋯β3M−22β3M−12⋮⋮⋯⋮⋮βM−222βM−232⋯βM−2M−22βM−2M−12βM−122βM−232⋯βM−1M−22βM−1M−12M−2×M−2,C2=−T22T23⋯T2M−2T2M−1T32T33⋯T3M−2T3M−1⋮⋮⋯⋮⋮TN−22TN−23⋯TN−2M−2TN−2M−1TN−12TN−13⋯TN−1M−2TN−1M−1N−2×M−2.Equations (19a) and (19b) are the coupled Lyapunov system, first solved for n=0 and then solved simultaneously for n=1,2,…,K by developing code in MATLAB 7.


## 4. Numerical Experiments and Discussion

Under this heading, to check the correctness and competence of the algorithm modified CTBS DQM, two test problems have been considered, which are available in the literature. All the computation work is conducted by using MATLAB 7.0. The following formulas are used for computing maximum absolute errorL∞, root mean square (RMS) error, and L2 error, respectively:(28)L∞=max1≤i≤N1≤j≤Meij,RMS=1NM∑i=1N∑j=1Meij21/2,L2=∑i=1N∑j=1Meij21/2,where uijanduij¯ are approximate and exact solutions, respectively, and eij=uij−uij¯.Problem 1.
As the first problem, consider 2D Burgers’ equations (1a) and (1b). The exact solutions over the domain D=x,y:0≤x≤1,0≤y≤1 is generated by the Hopf-Cole transformation [12, 15, 17] and obtained as(29)ux,y,t=34−141+eRe4y−4x−t/32,(30)vx,y,t=34+141+eRe4y−4x−t/32.
ICs and BCs are taken from exact solutions (29) and (30). The numerical results are shown with the help of Tables 2 and 3 and Figures 1–4 in form of errors, three-dimension, and contour plots. Convection prevails the flow which causes the errors become larger and larger as we increase the value of Re. L∞ is smaller than [15] for T=2.0,Re=100 with less grid points N=M=40. The figures show that exact solutions and numerical solutions are well consistent in three-dimensional and contour form. Table 4 shows that, as we increase the values of M and N, the absolute errors decrease which shows the convergence of the method.Table 2
L∞,RMS, and L2 errors of Problem 1 at different times and Reynolds number for ux,y,t with N=M=40.
TRe = 50Re = 100L∞RMSL2L∞RMSL20.17.064E − 041.020E − 044.081E − 032.287E − 044.676E − 051.871E − 030.56.864E − 041.102E − 044.409E − 032.393E − 046.091E − 052.437E − 031.06.386E − 041.068E − 044.273E − 032.424E − 046.153E − 052.461E − 031.56.150E − 041.031E − 044.124E − 032.403E − 045.857E − 052.343E − 032.06.028E − 049.924E − 053.969E − 032.348E − 045.521E − 052.208E − 03Table 3
L∞,RMS, and L2 errors of Problem 1 at different times and Reynolds number for ux,y,t with N=M=40.
TRe = 200Re = 400L∞RMSL2L∞RMSL20.19.269E − 041.265E − 055.059E − 033.269E − 032.390E − 043.563E − 030.59.948E − 048.723E − 053.489E − 033.742E − 033.349E − 043.768E − 031.09.167E − 041.345E − 055.282E − 033.910E − 033.214E − 043.867E − 031.59.020E − 041.367E − 055.458E − 034.481E − 033.514E − 043.989E − 032.09.037E − 041.332E − 055.321E − 034.891E − 033.890E − 044.891E − 03Figure 1
NSs and ESs ofux,y,t in 3D form for T=1.0 of Problem 1.
(a)(b)Figure 2
NSs and ESs ofux,y,t in 3D form for T=2.0 of Problem 1.
(a)(b)Figure 3
NSs and ESs ofux,y,t in 3D form for T=3.0 of Problem 1.
(a)(b)Figure 4
NSs and ESs ofux,y,t in 3D form for T=5.0 of Problem 1.
(a)(b)Table 4
Maximum abosulte errorL∞ of the problems at different nodes for ux,y,t.
TProblem2Problem1Re = 100Re = 100N = M = 10N = M = 20N = M = 40N = M = 10N = M = 20N = M = 400.12.923E − 045.569E − 052.950E050.747E − 030.326E − 032.287E − 040.54.958E − 047.518E − 052.421E − 050.381E − 030.103E − 032.393E − 041.03.932E − 048.167E − 051.387E − 050.444E − 039.972E − 042.424E − 04Problem 2.
Consider 2D Burgers’ equations (1a) and (1b) over the computational domain D=x,y:0≤x≤1,0≤y≤1 with the ICs [12],(31)ux,y,0=−4πcos2πxsinπyRe2+sin2πxsinπy,x,y∈D,vx,y,0=−2πsin2πxcosπyRe2+sin2πxsinπy,x,y∈D,and BCs,(32)u0,y,t=−2πexp−5π2tResinπyRe,u1,y,t=−2πexp−5π2tResinπyRe,v0,y,t=0,v1,y,t=0,0≤y≤1,t≥0,ux,0,t=0,ux,1,t=0,vx,0,t=−πexp−5π2tResin2πxRe,vx,0.5,t=πexp−5π2tResin2πxRe,0≤x≤1,t≥0.
The exact solutions of the problem are given by(33)ux,y,t=−4πexp−5π2t/Recos2πxsinπyRe2+exp−5π2t/Resin2πxsinπy,vx,y,t=−2πexp−5π2t/Resin2πxcosπyRe2+exp−5π2t/Resin2πxsinπy.
Tables5 and 6 show L∞, RMS, and L2 errors for different values of Re and time, while Figures 5–8 show a comparison of numerical and exact solution in three-dimensional form. Convection prevails the flow which causes the errors become larger and larger as we increase the value of Re. The figures show that exact and numerical solutions are well consistent in three-dimensional and contour form. Table 4 shows that, as we increase the values of M and N, the absolute errors decrease which shows the convergence of the method.Table 5
L∞,RMS, and L2 errors of Problem 2 at different times and Reynolds number for ux,y,t with N=M=40.
TRe = 50Re = 100L∞RMSL2L∞RMSL20.11.123E − 042.437E − 054.874E − 042.950E058.381E − 061.676E − 050.57.569E − 053.865E − 057.730E − 042.421E − 051.759E − 053.519E − 041.04.621E − 053.983E − 057.966E − 041.892E − 052.120E − 054.240E − 041.52.818E − 053.509E − 057.018E − 041.478E − 052.227E − 054.454E − 042.01.271E − 052.949E − 055.898E − 041.154E − 052.190E − 054.380E − 04Table 6
L∞,RMS, and L2 errors of Problem 2 at different times and Reynolds number for ux,y,t with N=M=40.
TRe = 200Re = 400L∞RMSL2L∞RMSL20.17.557E − 062.561E − 065.122E − 051.915E − 063.575E − 071.422E − 050.56.847E − 066.965E − 061.393E − 044.103E − 068.296E − 073.318E − 051.06.637E − 069.197E − 061.839E − 044.844E − 061.089E − 064.359E − 051.56.944E − 061.040E − 052.079E − 045.187E − 061.262E − 064.985E − 052.06.938E − 061.112E − 052.224E − 045.378E − 061.353E − 065.413E − 05Figure 5
NSs and ESs ofux,y,t in 3D form for T=2.0 of Problem 2.
(a)(b)Figure 6
NSs and ESs ofux,y,t in contour form for T=2.0 of Problem 2.
(a)(b)Figure 7
NSs and ESs ofvx,y,t in 3D form for T=2.0 of Problem 2.
(a)(b)Figure 8
NSs and ESs ofvx,y,t in contour form for T=2.0 of Problem 2.
(a)(b)


## 5. Conclusion

In this study, a modified CTBS DQM and a new algorithm to reveal the computational modeling of 2D coupled Burgers’ equations are developed. The proposed algorithm is tested on two benchmark problems appearing in the literature. The main results of this study are summarized as follows:(i)
A different technique using modified CTBS functions is presented to determine the WCs of 2D DQM than Lagrange interpolation traditional technique [22].(ii)
CTBS DQ algorithm proposed in [33] has extended for 2D problems in different forms, and it has concluded the algorithm worked nicely for the same problems.(iii)
The developed algorithm is better than the DQ algorithms proposed in [31, 32, 34] due to more smoothness of CTBS functions.(iv)
The presented method leads to quite similar results to those treated in [12, 15, 17, 18] and good accuracy in the case of a small number of grid points.(v)
After some modifications, the presented method can be extended to solve 2D or higher-dimensional equations. In this way, it can be used to analyze many other biological, mechanical or physical events, such as reaction, linear diffusion, dispersion, and nonlinear convection.Comparison of numerical solutions (NSs) and exact solutions (ESs) are given on left and right sides, respectively, in Figures [1–8] for R=100andΔt=0.001.



---

*Source: 7240300-2021-07-01.xml*