Numerical Solutions and Stability Analysis of Modified Burgers Equation Via Modified Cubic B-Spline

Abstract

Keywords: Differential quadrature method; modified Burgers’ equation; B-spline; Runge-Kutta method.

1. Introduction

The nonlinear Burgers’ equation (BE) first was first launched by Bateman [1] and then handled by Burgers’ [2], has the following form 𝑢𝑡 + 𝑢𝑢𝑥 − 𝜗𝑢𝑥𝑥 = 0,

in which 𝜗 is a positive parameter which may be indicate the viscosity and the subscripts 𝑡 and 𝑥 show time and space derivatives, respectively. The viscosity parameter serves as a regulator the balance between the convection and diffusion terms. When the viscosity parameter is taken zero (𝜗 = 0) the Burgers’ equation is turns into inviscid Burgers’ equation, which is the model for evaluation of shock waves that has a great importance in physical applications [3]. Seydaoğlu [4] is presented numerical solution of the Burgers’ equation via splitting methods for small value of viscosity parameter. BE is very important in fluid dynamics especially for turbulence problems, gas dynamics, heat conduction, continuous stochastic processes and theory of shock waves [5]. The modified Burgers’ equation (MBE) discussed in the present study depends on the Burgers’ equation (BE) of the following form 𝑢𝑡 + 𝑢2 𝑢𝑥 − 𝜗𝑢𝑥𝑥 = 0. *

Corresponding Author: e-mail: alibashan@gmail.com, tel: (372) 257 40 10 129

The MBE has the strong non-linear behaviours and also has widely been utilized in physical problems, for example for example turbulence transport, non-linear waves in a medium having low-frequency pumping or absorption, wave process in thermoelastic media, dispersion and transport of pollutants inside rivers and also sediment transportation, ion reflection available at quasi-perpendicular shocks [6]. Recently, some numerical studies of the equation have been presented: Bratsos [6] has used a finite difference scheme for calculating the numerical solution of the MBE. Ramadan and Danaf [8] and Ramadan et al. [7] have used the quintic and the septic B-spline collocation methods, respectively. A special lattice Boltzmann model is developed by Duan et al. [9]. Saka et al. [10] have developed a Galerkin FEM solution of the MBE. A solution bas ed on sextic B-spline collocation method is proposed by Irk [11]. Roshan et al. [12] has applied a Petrov-Galerkin method. A discontinuous Galerkin method is presented by Zhang Rong-Pei et al. [13]. Başhan et al. [14] have utilized quintic B-spline DQM to solve it numerically and Karakoç et al. [15] have solved the problem via two different methods namely, subdomain FEM and DQM and used quartic B-spline base functions at the both of the methods. Aswin and Awasthi [16] are investigated numerical solution of the MBE via iterative differential quadrature algorithms. Lakshmi and Awasthi [17] researched solution of the MBE numerically via quintic spline collocation method and Crank-Nicolson scheme. Bellman et al. [18] has firstly introduced DQM for solving the differential equations easily, in 1972 where partial derivative of a function in terms of a space direction is mentioned as a linear weighted sum of all of the values of the function for every grid point on the space direction [19]. DQM has widely been popular in recent years thanks to the easy usage of the method for many applications. Several authors have used different types of DQMs by using different test problems [20–33]. In the past few decades, DQM has become a very efficient and effective method for obtaining approximate solutions for different types of partial differential equations due to its simplicity for application. The DQM has several advantages over the other conventional techniques, basically, it prevents linearization and perturbation for finding better solutions of present nonlinear equations. Since MCB-DQM do not need transformation for solving the equation, the method has been preferred. In the present work, we have applied a MCB-DQMs for obtaining numerical solutions of the MBE. To illustrate the outstanding advantages of the MCB-DQMs and make a comparison of numerical solutions we have taken different values of 𝜗. A stability analysis has also been given.

2. Modified Cubic B-Spline Based Differential Quadrature Method

First, we will handle the Eq.(2) having the following initial condition 𝑈(𝑥, 0) = 𝑓(𝑥), 𝑎 ≤ 𝑥 ≤ 𝑏, and boundary conditions 𝑈(𝑎, 𝑡) = ℎ1 (𝑡), 𝑈(𝑏, 𝑡) = ℎ2 (𝑡),

in which ℎ1 (𝑡) and ℎ2 (𝑡) are constants. Let us take the grid distribution 𝑎 = 𝑥1 < 𝑥2 < ⋯ < 𝑥𝑁 = 𝑏 of the closed interval [𝑎, 𝑏]. In case a predefined function 𝑈(𝑥) is smooth enough on the solution region of the problem, its derivatives with regard to 𝑥 at a nodal point 𝑥𝑖 may also be approximated by a linear combination the values of the function over the solution region, that is, (𝑟)

in which 𝑟 is used to denote the order of derivative, 𝑤𝑖𝑗 are used to represent the weighting coefficients of the 𝑟 - th order derivative approximation, and 𝑁 is used to denote the number of

nodal points over the solution region. In this equation, the index 𝑗 is used to denote the fact that (𝑟) 𝑤𝑖𝑗 is the related weighting coefficients for the values of function 𝑈(𝑥𝑗 ). In the present study, we are going to need the first and the second order derivatives of he function 𝑈(𝑥). So, we will find the value of the equation (5) for the values of 𝑟 = 1, 2. Let Κ 𝑖 (𝑥) denote the cubic B-splines having meshes at the points 𝑥𝑖 in which the uniformly scattered 𝑁 nodal points are taken as 𝑎 = 𝑥1 < 𝑥2 < ⋯ < 𝑥𝑁 = 𝑏 on the real axis. Next, the cubic B-splines {Κ 0 , Κ1 , … , Κ 𝑁+1 } constitute a basis for functions described over [𝑎, 𝑏]. Those cubic B-spline basis functions Κ 𝑖 (𝑥) are given by the following equalities: (𝑥 − 𝑥𝑖−2 )3 , [𝑥𝑖−2 , 𝑥𝑖−1 ] (𝑥 − 𝑥𝑖−2 )3 − 4(𝑥 − 𝑥𝑖−1 )3 , [𝑥𝑖−1 , 𝑥𝑖 ] 1 Κ 𝑖 (𝑥) = 3 (𝑥𝑖+2 − 𝑥)3 − 4(𝑥𝑖+1 − 𝑥)3 , [𝑥𝑖 , 𝑥𝑖+1 ] ℎ (𝑥𝑖+2 − 𝑥)3 [𝑥 , 𝑖+1 , 𝑥𝑖+2 ] {0 , 𝑜𝑡ℎ𝑒𝑟𝑤𝑖𝑠𝑒 where ℎ = 𝑥𝑖 − 𝑥𝑖−1 for all 𝑖 [34]. Table 1 shows the values of cubic B-spline functions and their derivatives given at the nodal points. Table 1. The values of cubic B-spline functions and their derivatives. X Κ𝑖 Κ ′𝑖 Κ ′′ 𝑖

A matrix system of equations which is dominant in diagonal is obtained by using the modified cubic B-splines. This structure of diagonally dominant matrix is also great value to the stability analysis. The modification process of cubic B-spline functions can also be performed in a different way. We used the modification which is presented by Mittal and Jain [36] as given below Ψ1 (𝑥) = Κ1 (𝑥) + 2Κ 0 (𝑥) Ψ2 (𝑥) = Κ 2 (𝑥) − Κ 0 (𝑥) Ψ𝑠 (𝑥) = Κ 𝑠 (𝑥), for 𝑠 = 3,4, … , 𝑁 − 2 Ψ𝑁−1 (𝑥) = Κ 𝑁−1 (𝑥) − Κ 𝑁+1 (𝑥) Ψ𝑁 (𝑥) = Κ 𝑁 (𝑥) + 2Κ 𝑁+1 (𝑥)

where Ψ𝑡 , (𝑡 = 1, 2, … , 𝑁) constitutes a basis functions on the [a, b] domain.

3. Numerical Discretization Process

The modified Burgers’ equation of the form 𝑈𝑡 + 𝑈 2 𝑈𝑥 − 𝜗𝑈𝑥𝑥 = 0,

having the initial condition (3) and the boundary conditions (4) is rewritten as, 𝑈𝑡 = −𝑈 2 𝑈𝑥 + 𝜗𝑈𝑥𝑥 = 0.

After that the differential quadrature approach have been used in Eq. (18) 𝑑𝑈(𝑥𝑖 ) 𝑑𝑡

and ODE (19) is obtained. Next, ODE in Eq. (19) is integrated in terms of time via a suitable method. In this manuscript, we have chosen strong stability preserving low storage RungeKutta4(3)5[3S∗] method [37] thanks to its outstanding characteristics such as accuracy, stability and the allocated memory storage.

4. Numerical Results And Stability Analysis

In the present section, we will find out the numerical solutions of the MBE via MCB-DQM. The efficiency of the present method has been controlled via the error norms 𝐿2 and 𝐿∞ , respectively: 𝑁

𝐿∞ = ‖𝑈 𝑛 − 𝑢𝑛 ‖∞ ≅ max|𝑈𝑗 − 𝑢𝑗 |, 𝑗 = 1,2, … 𝑁 − 1. 𝑗

Calculation of the eigenvalues for the coefficient matrices is a requirement for the stability analysis of a numerical technique used for a differential equation which is nonlinear. With the application of DQM to the MBE, the equation becomes an ODE. The stability of a time-dependent problem: 𝜕𝑈 𝜕𝑡

having the appropriate initial and boundary conditions, in which 𝑙 is a spatial differential operator. When discretized via DQM, Eq. (20) turns into a group of ODEs in time as follows 𝑑{𝑈} 𝑑𝑡

in which {u} is a vector consisting of the unknown values of the function nodal points except for the both boundary points, {b} is a vector containing the non-homogenous part and the boundary conditions and A is the coefficient matrix. The stability of ODE given by Eq.(21) also determines that of numerical scheme for numerical integration of Eq.(21). When ODE given in Eq.(21) is not stable, it is possible that the numerical methods do not generate convergent solutions. It is obvious that the stability of Eq. (21) depends on the eigenvalues of matrix 𝐴, for its analytical solution is directly given by the eigenvalues of this matrix. When all 𝑅𝑒(𝜆𝑖 ) ≤ 0 for all 𝑖 is enough to illustrate the stability of the analytical solution of {𝑢} as 𝑡 → ∞ in which 𝑅𝑒(𝜆𝑖 ) denotes the real part of the eigenvalues 𝜆𝑖 of the matrix 𝐴. The matrix 𝐴 at Eq.(21) is evaluated as (1) (2) 𝐴𝑖𝑗 = −𝛼𝑖2 𝑤𝑖,𝑗 +𝜗𝑤𝑖,𝑗 where 𝛼𝑖 = 𝑈(𝑥𝑖 , 𝑡) [19]. It is necessary that the eigenvalues of matrix 𝐴 be inside the stability domain illustrated in Figure 1 [38].

Figure 1. SSP-Runge-Kutta method stability regions. The analytical solution of MBE is given in [35] as: 𝑈(𝑥, 𝑡) =

where 𝑐0 is a constant and 0 < 𝑐0 < 1 . For our numerical calculations, we take 𝑐0 = 0.5. We use the initial condition for the Eq.(22) evaluating at 𝑡 = 1 and the boundary conditions are taken as 𝑈(0, 𝑡) = 0 and 𝑈(1, 𝑡) = 0 . For the numerical simulation, we have chosen the various viscosity parameters ranging from 𝜗 = 0.01, to 𝜗 = 0.0005 and time step ∆𝑡 = 0.01, ∆𝑡 = 0.001 over the interval 0 ≤ 𝑥 ≤ 1. The behaviors of the numerical solutions for different values of viscosity are illustrated in Figure 2. As it is observed from the Figure 2, when we select the solution domain 0 ≤ 𝑥 ≤ 1 by the decreasing values of viscosity parameters 𝜗, the amplitude of waves have decreased, the waves became steep and both of the solutions with Method 1 and Method 2 have became indistinguishable. Also, we have seen from the Figure 2 that when the time increases the curve of the the numerical solution decay. With smaller viscosity value, the decay of numerical solution has become faster. In Figure 3 absolute errors are drawn at time 𝑡 = 10 for various viscosity parameters from 𝜗 = 0.01, to 𝜗 = 0.0005 by both methods, respectively. As it seen from Figure 3 error norms have similar values for Method 1 and Method 2. Maximum error occurs at the right hand boundary when greater value of viscosity 𝜗 = 0.01 is taken and with smaller value of viscosities 𝜗 = 0.005, 𝜗 = 0.001 and 𝜗 = 0.0005 maximum error has been monitored near the point at which the shock wave has its highest amplitude.

Figure 2. Behaviour of numerical solutions at various time for 𝜗 = 0.01, 𝜗 = 0.005, 𝜗 = 0.001 and 𝜗 = 0.0005, respectively. 136

The computed values of the error norms 𝐿2 and 𝐿∞ are presented at some selected times up to 𝑡 = 10. The obtained results are tabulated in Tables [2−5]. As it is seen from the Table 2, the error norms 𝐿2 and 𝐿∞ are sufficiently small and satisfactorily acceptable. Also, it is clearly observed from these tables that from comparison of present solutions with earlier works that MCB-DQMs using less number of grid points solutions are better than quintic B-spline DQM. [20] and septic collocation FEM [7]. When we have extended the solution domain from 0 ≤ 𝑥 ≤ 1 to 0 ≤ 𝑥 ≤ 1.3 present numerical results that compared and tabulated with earlier work [14] at Table 3 are acceptable. In Tables 4 and 5 clearly seen that if the value of viscosity 𝜗 decrease, the value of the error norms will decrease, too. We have obtained the better results if the value of viscosity parameter is get smaller. In Table 5 numerical results for 𝜗 = 0.005 and 𝜗 = 0.0005 are given. Also, at Table 5 the error norms of 𝐿2 and 𝐿∞ obtained by the present methods are enough small for acceptance. Table 2. 𝐿2 and 𝐿∞ error norms for 𝜗 = 0.01, ∆𝑡 = 0.01 at 0 ≤ 𝑥 ≤ 1. 𝑡 2 3 4 5 6 7 8 9 10

Method 1 N = 11 𝐿2 × 103 𝐿∞ × 103 0.520872 1.265996 0.547643 1.079945 0.535512 0.998953 0.519949 0.858252 0.521119 0.799671 0.548745 0.734032 0.598388 0.962144 0.659246 1.181440 0.722110 1.368425

Method 1 N = 11 𝐿2 × 103 𝐿∞ × 103 0.575145 1.399995 0.585931 1.087017 0.563132 1.027359 0.534899 0.905242 0.507374 0.782472 0.482391 0.731187 0.461106 0.675984 0.444932 0.620800 0.435297 0.569807

Method 1 N = 11 𝐿2 × 103 𝐿∞ × 103 0.127939 0.402897 0.122057 0.432444 0.112170 0.366878 0.103869 0.302373 0.096985 0.266505 0.091180 0.247102 0.086203 0.225710 0.081876 0.205377 0.078068 0.187063

Method 2 N = 11 𝐿2 × 103 𝐿∞ × 103 0.744386 1.497042 0.622754 1.210896 0.543888 0.988942 0.505238 0.888133 0.502769 0.772025 0.532253 0.762053 0.582727 0.995931 0.642355 1.203811 0.702863 1.380314

Quin. DQM.[20]N = 51 𝐿2 × 103 𝐿∞ × 103 0.688316 1.406116 0.611194 1.228470 0.551891 1.047041 0.524368 0.911470 0.536004 0.814737 0.583793 1.014095 0.652737 1.301495 0.727927 1.545607 0.800131 1.742584

Sep. Coll. [7] N = 51 𝐿2 × 103 𝐿∞ × 103 0.790430 1.703092 0.655193 1.183270 0.557679 0.996452 0.510562 0.856134 0.516723 0.761053 0.567744 1.065455 0.642754 1.358111 0.723643 1.604831 0.800256 1.802394

Table 3. 𝐿2 and 𝐿∞ error norms for 𝜗 = 0.01, ∆𝑡 = 0.01 at 0 ≤ 𝑥 ≤ 1.3. Method 2 N = 11 𝐿2 × 103 𝐿∞ × 103 0.855814 1.814488 0.667538 1.163823 0.545585 0.924348 0.475804 0.841016 0.432537 0.727781 0.403642 0.632487 0.384467 0.593127 0.373711 0.549625 0.371088 0.508189

Quin. DQM.[14]N = 51 𝐿2 × 103 𝐿∞ × 103 0.647514 1.418692 0.603831 1.248161 0.559737 1.074432 0.524882 0.939394 0.496279 0.834173 0.472949 0.751101 0.455623 0.685356 0.445747 0.631350 0.444390 0.587301

Quar. DQM. [14]N = 51 𝐿2 × 103 𝐿∞ × 103 0.888636 1.842686 0.783980 1.273734 0.684999 1.040082 0.594035 0.897385 0.507332 0.772407 0.427485 0.662413 0.357823 0.567976 0.301899 0.489850 0.263822 0.424509

Table 4. 𝐿2 and 𝐿∞ error norms for 𝜗 = 0.001, ∆𝑡 = 0.01 at 0 ≤ 𝑥 ≤ 1. Method 2 N = 11 𝐿2 × 103 𝐿∞ × 103 0.186973 0.651824 0.138003 0.457424 0.105227 0.303682 0.086320 0.241783 0.074980 0.207681 0.067703 0.193589 0.062680 0.176297 0.058969 0.159393 0.056064 0.144157

Quin. DQM.[20]N = 166 𝐿2 × 103 𝐿∞ × 103 0.127227 0.457137 0.110849 0.389233 0.098569 0.333200 0.090234 0.288512 0.084073 0.254679 0.079187 0.228346 0.075126 0.207124 0.071646 0.190023 0.068599 0.175928

Sep. Coll. [7] N = 201 𝐿2 × 103 𝐿∞ × 103 0.183549 0.818521 0.144142 0.523483 0.114411 0.356354 0.094787 0.254979 0.081418 0.213485 0.071898 0.188005 0.064837 0.168260 0.059412 0.152408 0.055115 0.139431

Figure 3. Absolute errors at region 0 ≤ 𝑥 ≤ 1, 𝑡 = 10, 𝑁 = 51 for various value of 𝜗 = 0.01, 𝜗 = 0.005, 𝜗 = 0.001 and 𝜗 = 0.0005, respectively.

Table 5. 𝐿2 and 𝐿∞ error norms at various time for ∆𝑡 = 0.001 at 0 ≤ 𝑥 ≤ 1. 𝑡 2 3 4 5 6 7 8 9 10

𝜗 = 0.005, N = 11 Method 1 Method 2 𝐿2 × 103 𝐿∞ × 103 𝐿2 × 103 𝐿∞ × 103 0.375651 1.019073 0.551859 1.366950 0.368724 0.886011 0.418319 0.890476 0.347760 0.694557 0.331494 0.644805 0.326879 0.655047 0.281799 0.558450 0.307983 0.587321 0.251519 0.506531 0.291229 0.520186 0.231408 0.447099 0.276497 0.469593 0.217136 0.392948 0.263755 0.441707 0.206833 0.364836 0.253180 0.412211 0.199821 0.342825

𝜗 = 0.0005, N = 11 Method 1 Method 2 𝐿2 × 103 𝐿∞ × 103 𝐿2 × 103 𝐿∞ × 103 0.056145 0.252287 0.083951 0.373187 0.059030 0.257393 0.067899 0.258396 0.057538 0.214512 0.057624 0.226150 0.055043 0.199802 0.051556 0.187015 0.052420 0.177984 0.047540 0.171625 0.049923 0.159320 0.044581 0.153390 0.047622 0.148779 0.042227 0.136386 0.045525 0.137607 0.040258 0.126996 0.043619 0.126915 0.038557 0.118341

An analysis about the stability of the matrix is also carried out for the MCB-DQM. For software, MATLAB R2013b program is used in order to calculate the eigenvalues of the coefficient matrix of single soliton problem. Eigenvalues of the suggested method for varying number of nodal points are demonstrated in Figure 4. The eigenvalues for N = 11, N = 21 and N = 31 have also imaginary parts for both methods. When grid points are chosen N = 41 and N = 51 eigenvalues that are obtained by Method 1 have only real parts but eigenvalues that are obtained by Method 2 have imaginary parts again. It is seen that all of the eigenvalues lie in the stability criteria region [38]. Also, maximum absolute value of the eigenvalues for different number of nodal points are tabulated in Table 6. Method 1 produces pure real eigenvalues when the grid numbers greater than 41 but method 2 produces pure real eigenvalues from 101. As it is seen from Table 6 when the number of the nodal points are taken higher since absolute value of eigenvalue grows, time step should be decreased to obtain the stable solution. In this study, an ordinary PC Intel(R) Core(TM) i7 − 4770 CPU @ 3.40 GHz Windows 10 Home is used for calculations. CPU time of the both methods are computed and reported at Table

7. Those results show that the present methods are practicable and speedy.

Table 6. Eigenvalues for different number of nodal points. Grid Number 𝑀𝑎𝑥|𝑅𝑒(𝜆)| 𝑀𝑎𝑥|𝐼𝑚(𝜆)|

Method 1: 𝜗 = 0.001 21 31 41 4.68 10.68 19.68 1.40 0.52 0.00 Method 2: 𝜗 = 0.001 21 31 41 2.54 5.84 10.45 2.15 3.23 3.90

Table 7. CPU time(sec.) for various value of 𝜗 at solution domain 0 ≤ 𝑥 ≤ 1. 𝜗 0.01 0.001 0.005 0.0005

Figure 4. Eigenvalues for various number of grid points: 𝜗 = 0.001.

5. Conclusion

Throughout this manuscript, we have tried to apply DQMs depending on the modified cubic B-spline basis functions for approximate solutions of the MBE. The accuracy and efficiency of the new method have been demonstrated by calculating the error norms 𝐿2 and 𝐿∞ . The outstanding feature of the methods is its ability to obtain acceptable good approximations despite of using even less number of nodal points. This can be easily seen from the tables present in this article. As it can be observed by the comparison between the calculated values of the error norms 𝐿2 and 𝐿∞ by this method and previous works, MCB-DQMs results are acceptable good. The

stability analysis of the numerical approximation by the eigenvalues has also been made. The obtained results show that MCB-DQMs may be used to obtain reasonable accurate numerical solutions of the MBE. So, MCB-DQMs is a reliable one to obtain the approximate solutions some physically important nonlinear partial differential equations.

Share and Cite