Numerical simulation of Burgers’ equation using cubic B-splines

1 Introduction

Consider the nonlinear Burgers’ equation given by

where

with initial condition

and boundary conditions

where ν > 0 is the kinematic viscosity parameter and g, f₁ and f₂ are given functions of the arguments which are sufficiently smooth.

Burgers’ equation is one of the most fundamental nonlinear partial differential equations which has its impact in numerous areas of scientific research such as viscous flow and turbulence, gas dynamics, shock theory, traffic flow, stochastic process, cosmology, etc. Burgers’ equation was introduced in 1915 by a British-American mathematician Harry Bateman [6]. He was interested in vanishing the viscosity limit, i.e. ν →0. The applications of Burgers’ equation began when a Dutch physicist J. M. Burger [7] made use of it to elucidate some important characteristics of turbulence such as balance of energy and appearance of dissipation layers. He did immense research in this equation and to honour his contributions, this equation is well known as Burgers’ equation. Research in Burgers’ equation became more rigorous when Hopf [8] and Cole [9] in 1950 and 1951 respectively, independently developed a nonlinear transformation which could map the nonlinear Burgers’ equation to a linear Heat equation. Thus the analytical solution of Burgers’ equation could be obtained for a restricted set of initial conditions. Burgers’ equation, due to the simultaneous presence of convective term (u∂u∂x) and diffusive term (v∂2u∂x2), can be employed to study the opposite effects of convection and diffusion at a basic level, thus making the equation all the more celebrated. When ν approaches zero, Eq. (1.1a) becomes the inviscid Burgers’ equation which is a model for nonlinear wave propagation. When u approaches zero, Eq. (1.1a) becomes the Heat equation. Thus, it can be thought of as a hyperbolic problem with artificial diffusion for small kinematic viscosity ν and Heat equation for very small u.

Burgers’ equation is basically the Navier Stokes equation without the pressure term. For low kinematic viscosity there can be velocity discontinuities, i.e. shocks. Kreiss and Kreiss [10] considered the viscous Burgers’ equation with the initial and non-homogeneous Dirichlet boundary conditions to study the effect of shock on convergence of steady state solution. They proved that to speed up the convergence to the steady state, the shock should be located at the boundary. On the other hand, location of shock in the interior leads to very slow convergence. The speed of convergence was studied by analyzing the corresponding eigenvalue problem. Bec and Khanin [11] explained Burgers’ turbulence as the study of solutions to the one or multi-dimensional Burgers’ equation with random initial conditions or random forcing. Recently, it has become a point of attraction to many researchers due to the new emerging application of Burgers’ model to statistical physics, cosmology and fluid dynamics.

One of the earliest attempts to solve Burgers’ equation was by Benton and Platzman [12] who presented 35 distinct analytic solutions of the equation with different initial conditions. Rodin [13] studied some approximate and exact solution of boundary value problems for Burgers’ equation with the help of Hopf-Cole transformation. Varoglu and Finn [14] developed a finite element method based on weighted residual method where the method of characteristic was employed using space time elements. Cadwell et al. [15] attempted a piecewise approximation method (finite element method) using two elements with the aim of ‘chasing the peak’ by altering the size of the elements at each stage by using the information at the previous step. Cadwell and Smith [16] extended this method to the general case of n elements. They observed that for large Reynold’s number, the finite element results were much superior to the finite difference results. A number of other authors used moving node finite element method to solve Burgers’ equation. Gelinas et al. [17] presented a node moving finite element method which can be applied to large gradients or shocks with high resolution and accuracy. In their system, the nodes move systematically and continuously to those regions where they are required the most. Cadwell et al. [18] further developed a moving node finite element method by using an algorithm which was a generalization of that considered by Caldwell et al. [15]. Several other finite element techniques were employed for solving Burgers’ equation [19], [20], [2]. Asaithambi [21] presented a simple numerical method based on automatic differentiation or algorithmic differentiation that used second-order finite-differences for the spatial derivatives and marched the solution in time using a Taylor series expansion. Kadalbajoo et al. [22] developed an implicit scheme for solving Burgers’ equation. They used a standard backward Euler scheme with constant time step to discretize in temporal direction and a standard upwind finite difference scheme to discretize in spatial direction on piecewise uniform mesh. The quasi-linearization process was used to tackle the nonlinearity. A numerical method based on Crank-Nicolson scheme was put forward by Kadalbajoo and Awasthi [23], where they first reduced the Burgers’ equation to a linear Heat equation using Hopf-Cole transformation and then dicsretized using Crank-Nicolson scheme. The mesh size could be chosen without any restriction. It was shown that the scheme was second order accurate in both space and time and also unconditionally stable. Rashidi et al. [24] proposed the homotopy perturbation method (HPM) for solving generalized Burgers’ and Burger-Fisher equations. In most cases, the method gave good results even with a few terms. Rashidi et al. [25] also employed the homotopy analysis method (HAM) effectively to Burgers’ equation and regularized long wave (RLW) equations and made a comparative study of their results with the HPM to establish the effectiveness of the former method over the latter. Rashidi and Erfani [26] then extended the differential transform method (DTM) to develop approximate analytical solutions of the Burgers’ equation and found that the results were more accurate as compared to HAM. Adomian’s decomposition method was tried by Abbasbandy and Darvishi [27], [28]. Several mesh free methods [29], [30], [31], [32] have also been developed for solving Burgers’ equation. In 2015, Vijitha Mukundan and Ashish Awasthi [33] presented new efficient numerical techniques for solving Burgers’ equation. They used Hopf-Cole transformation to get one-dimensional diffusion equation which was semi-discretized by using method of lines (MOL). Resulting system of ODEs was solved by backward differentiation formulas (BDF) of order one, two and three.

Inan and Bahadir [34] solved the nonlinear Burgers’ equation using an explicit exponential finite difference method. They employed Hopf-Cole transformation to convert the nonlinear Burgers’ equation to linear Heat equation. Inan and Bahadir [35] further illustrated an implicit and a fully implicit exponential finite difference method for solving Burgers’ equation. They used Newton’s method to solve the nonlinear system at each time step. Inan and Bahadir [36] improved on their exponential method which lead to Crank-Nicolson exponential finite difference solution of Burgers’ equation. In this method, they used Crank-Nicolson scheme for discretization and Newton Raphson method to solve the system at each time step. Recently in 2015, Inan and Bahadir [37] came up with an alternate implicit and Crank-Nicolson exponential finite difference method where they first reduced the Burgers’ equation to linear Heat equation by Hopf-Cole transformation and then solved the same using implicit and Crank-Nicolson exponential finite difference method.

A remarkable method to solve Burgers’ equation is by using splines. The term ‘splines’ is adopted from the name of a flexible strip of metal used by drafters in olden times to assist in drawing curves. Usage of spline functions has great advantages. They are constructed in such a way that they possess high degree of smoothness even at the juncture between two polynomial pieces, i.e. n^th degree spline function is continuous up to (n - 1)^th derivative. They are easy to store, manipulate and evaluate on a digital computer. The matrices arising naturally on using splines have convenient determinantal properties. Moreover, every continuous function on the interval [a, b]can be approximated arbitrarily well by polynomial splines, provided a sufficient number of knots are allowed. These desirable properties of splines have driven researchers to develop different types of spline functions of various degree. B-spline or Basis spline is a spline function that has minimal support with respect to given degree of smoothness and domain partition. Any spline function of a particular degree and smoothness can be written as a linear combination of B-splines of the same degree and smoothness. A remarkable property of B-spline that it has minimal support bestows it with a computational edge over finite element method, while its fourth order accuracy marks it superior to finite difference methods.

Splines was used in a mathematical context for the first time in 1946 by I. J. Schoenberg [38], in connection with smooth piecewise polynomial approximation. Since then, it has come a long way in solving Burgers’ equation. B-spline was also introduced by Schoenberg to approximate equidistant data by analytic functions by means of integrals. B-splines for arbitrarily spaced data was later given by Curry and Schoenberg [39]. Roach and Yanping [40] used linear exponential spline function as an interpolant to generate finite difference representation of Burgers’ equation. The time marching was performed using backward Euler method. In 1992, L.R.T. Gardner et al. [41] developed a collocation solution of Burgers’ equation using cubic B-spline finite element method. Kutluay et al. [42] presented a least-squares quadratic B-spline finite element method for Burgers’ equation by applying the classical weighted residual method over the finite elements. Dag et. al [43] solved time split and space split Burgers’ equation using B-spline collocation method. Aksan [44] brought forward another method in which Burgers’ equation was converted to a set of nonlinear ordinary differential equations by the method of discretization in time and then each of them was solved by applying the quadratic B-spline finite element method. A numerical solution based on collocation using septic splines was developed by Ramadan et al. [45]. Here, the time discretization was done using Crank-Nicolson scheme. The method proved to be unconditionally stable. Ramadan et al. [46] also proposed a non-polynomial spline solution of Burgers’ equation. The non-polynomial spline function in this work has a trigonometric part and a polynomial part of first degree. The C^∞ differentiability of the trigonometric part of non-polynomial spline compensates for the loss of smoothness inherited by polynomial splines. Moving boundary conditions were used to improve the accuracy. In another study of Burgers’ equation by Saka et al. [47], time and space splitting techniques were applied to the Burgers’ equation and the modified Burgers’ equation and then collocation procedure using quintic B-spline was employed to approximate the resulting systems. The method was proven to be better than quadratic B-spline F.E.M and quartic B-spline collocation method. Many other numerical methods for solving Burgers’ equation using splines were developed such as that by Jiang and Wang[48], etc. In 2012, R.C. Mittal and R.K. Jain [49] developed a method for solving Burgers’ equation based on collocation of modified cubic B-splines over finite elements. In 2014, Geetha Arora and Brajesh Kumar [50] used modified cubic B-splines in differential quadrature method to give a system of ordinary differential equation, which was solved by strong stability preserving time stepping Runge Kutta (SSP-RK43) scheme.

In this paper, we reduce the Burgers’ equation to a linear Heat equation using Hopf-Cole transformation. Linear Heat equation thus obtained is further solved using cubic B-splines in space. Time discretization is done using both Crank-Nicolson scheme as well as backward Euler scheme. We have considered various test examples and computed L₂ and L_∞ error norms to establish the accuracy of the proposed methods. Stability check using Von Neumann stability analysis has confirmed that the method is unconditionally stable. We have also compared our numerical results with those of Kutluay et al. [1], Ozis et al. [2], Dag et al. [3], Salkuyeh et al. [4] and Korkmaz et al. [5] and established that our results are more accurate as compared to the above mentioned works.

2 Hopf-Cole Transformation

Hopf [8] and Cole [9] in 1950 and 1951 respectively came up with a nonlinear transformation which could reduce the Burgers’ equation to a linear Heat equation for a restricted set of initial conditions. The transformation was thus named after them as Hopf-Cole transformation. The concept of Hopf-Cole transformation is of prime significance in the theory of conservation laws. Hopf and Cole proved that in the context with initial and boundary conditions of Eq. (1.1a), the nonlinear transformation

satisfies the viscous Burgers’ equation

where ϕ(x, t) is a solution to the linear Heat equation

Therefore, by using the nonlinear Hopf-Cole transformation, the nonlinear Burgers’ equation Eq. (1.1a) with initial condition Eq. (1.1b) and Dirichlet boundary conditions (f₁(t) = f₂(t) = 0, 0 ≤ t ≤ Τ ) is reduced to linear Heat equation

with initial condition

and boundary conditions

3 Exact solution of Burgers’ equation

The Fourier series solution to the linearized Heat equation is

with Fourier coefficients at t = 0 as

Using the Hopf-Cole transformation the exact solution to Burgers’ equation Eq. (1.1a) is

4 Numerical θ scheme based on cubic B-splines

Consider a uniform partition of the solution domain [a, b] with nodal points x_j as a = x₀ < x₁ < ... < x_N = b. The cubic B-spline at the nodal points is given by

where, B_-1, B₀, ..., B_N-1, B_N, B_N+1 forms a basis of all cubic continuous functions over [a, b]. Thus, any cubic continuous function on [a, b] can be written as a linear combination of the cubic B-splines as

where δ_m(t), m = -1,..., N + 1 are the time dependent parameters to be obtained from boundary conditions and B-spline formulation. It may be noted that each cubic B-spline is non-zero over four adjacent intervals.

The values of B_j(x)and its derivatives obtained from Eq. (4.1) are tabulated in Table 1.

Using approximation function Eq. (4.2) and Table 1, we have

The linear Heat equation Eq. (2.2) with initial condition Eq. (2.3) and boundary condition Eq. (2.4) is to be solved by making use of cubic B-spline.

To begin with, consider the weighted average of linear Heat equation as follows:

In this paper, we make use of Crank-Nicolson scheme, i.e. θ=12and backward Euler scheme, i.e.θ =1.

Discretizing time derivative of Eq. (4.6) by finite difference

Substituting Eq. (4.3) and Eq. (4.5) in Eq. (4.7), we have

The parameters δ_-1 and δ_Ν+1in Eq. (4.8) are eliminated using boundary condition as follows:

From Eq. (4.4) and boundary conditions Eq. (2.4), we have

The initial vector δ⁰, i.e. δ00,δ10,…,δN0 are obtained from Eq. (4.3) for n = 0. The solution of linear Heat equation is thus obtained by substituting values of δmnin Eq. (4.3). We thereby find the solution of Burgers’ equation Eq. (1.1a) in terms of the solution of linear Heat equation Eq. (2.2) using Hopf-Cole transformation Eq. (2.1)

The accuracy in temporal direction is improved by Richardson extrapolation as follows:

where h and k are the grid lengths in x and t directions respectively.

5 Stability Analysis

The stability analysis is investigated using a versatile technique known as Fourier stability method. Using this method, we assume

where ζ is known as amplification factor. For stability we require ζ < 1.

Using Eq. (5.1) it follows that

Substituting Eq. (5.1) to Eq. (5.4) in Eq. (4.8), we get

Dividing through out by ζⁿe^iαm

Putting e^iα = cos α+ i sin αand e^-iα = cos α - i sin α

Simplifying Eq. (5.7) we get;

By Fourier stability method, solution is stable for ζ < 1.

For simplifying the notation, we put

(2cosα)(1+6h2Δt(1−θ)να2)+(4−12h2Δt(1−θ)να2)=a and (2cosα)[1−6h2Δtθνα2]+(4+12h2Δtθνα2)=b so that ζ=ab

To prove that ζ=ab≤1, it is enough to prove that a ≤ b or a – b <0

When θ = 1;

a=(2cosα+4)andb=(2cosα)[1−6h2Δtνα2]+(4+12h2Δtνα2)

so that Eq. (5.9) reduces to

Here,

Hence a ≤ b and therefore ζ=ab≤1

When θ=12,

a=(2cosα)(1+3h2Δtνα2)+(4−6h2Δtνα2) and

b=(2cosα)[1−3h2Δtνα2]+(4+6h2Δtνα2) so that Eq. (5.9) becomes

In this case,

Hence a ≤ b and therefore ζ ≤ 1.

Hence, for both values of Θ, the scheme is unconditionally stable.

6 Numerical Experiments and Discussion

In this section, we consider some test examples to demonstrate the efficiency and accuracy of the proposed method. The numerical solutions obtained by proposed method at various times and for different values of coefficients of viscosity are plotted in the cases when θ = 0.5 and θ = 1 and compared with exact solution. The accuracy and efficiency of the method for each example considered is examined by calculating the L₂ and L_∞ error norms which are defined as follows:

where ui⋆is the numerical solution obtained by the proposed method and uiexact is the exact solution.

The 2-D plots and surface plots have been plotted to confirm the numerical results. The numerical results have also been compared with those of Kutluay et al. [1], Ozis et al. [2], Dag et al. [3], Salkuyeh et al. [4] and Korkmaz et al. [5] When θ=12, the system of equation Eq. (4.8) becomes

and when θ = 1, the system of equation Eq. (4.8) reduces to

Example 3

Consider the Burgers’ equation Eq. (1.1a) with the initial condition

u(x,0)=2νπsin⁡(πx)σ+cos(πx)),0<x<1,(6.19)

and homogeneous boundary conditions

u(0,t)=u(1,t)=0,0≤t≤T.(6.20)

where, σ > 0 is a parameter.

By Hopf-Cole transformation, the Burgers’ equation Eq. (1.1a) is transformed to the linear Heat equation Eq. (2.2) with the initial condition

ϕ(x,0)=σ+cos⁡(πx)σ+1,0≤x≤1,(6.21)

and Neumann type boundary conditions

ϕx(0,t)=ϕx(1,t)=0,0≤t≤T.(6.22)

The exact solution of the above equation is given by

u(x,t)=2πνe−π2νtsin⁡(πx)σ+e−π2νtcos⁡(πx);0≤x≤1(6.23)

To obtain the numerical solution, the required starting value of δ⁰, i.e. δ00,δ10,…,δN0is obtained from Eq. (4.3) by substituting n = 0, i.e.

δm−10+4δm0+δm+10=σ+cos⁡(πxm0)σ+1;m=−1,0,...N+1(6.24)

In the tables that follow, we have given the L₂and L_∞ error norms of the examples discussed for different values of coefficient of viscosity (v), time step (Δt) and space step (h) for θ = 0.5 and θ = 1 . We have also presented the comparison of our numerical results with those of Kutluay et al. [1], Ozis et al. [2], Dag et al. [3], Salkuyeh et al. [4] and Korkmaz et al. [5].

Tables 2 to 4 give the L₂ and L_∞ error norms of Example l for v = 1, 0.1 and 0.01 respectively. The L₂ and L_∞ error norms of Example 2 are given in Tables 5 to 7 and those of Example 3 are shown in tables 8 to 10. The L₂ and L_∞ error norms are very small establishing that numerical solution and analytical solution are very close. We have also presented the graphs for numerical solution in comparison with the exact solution for all the examples.

Table 2

L₂ and L_∞error norms of Example 1 for ν = 1, ∆t = 0.001 and h = 0.0125

t	θ = 0.5		θ=1
	L2	L∞	L2	L∞
0.1	3.26483×10−4	2.06493×10−3	3.30030×10−4	2.08277×10−3
0.2	1.39087×10−4	8.79499×10−4	7.25395×10−5	4.58657×10−4
0.3	5.83493×10−5	3.69059×10−4	8.44307×10−6	5.33990×10−5
0.4	2.41857×10−5	1.52965×10−4	3.77824×10−6	2.38963×10−5
0.5	9.92293×10−6	6.27582×10−5	3.99013×10−6	2.52359×10−5

Table 3

L₂ and L_∞error norms of Example 1 for ν = 0.1, ∆t = 0.0001 and h = 0.01

t	θ = 0.5		θ=1
	L2	L∞	L2	L∞
1	2.05415×10−4	1.94686×10−3	2.95631×10−4	2.39409×10−3
2	6.52273×10−5	6.41337×10−4	1.22403×10−5	1.11272×10−4
3	2.41200×10−5	2.40547×10−4	1.50685×10−5	1.50261×10−4
4	9.23729×10−6	9.23372×10−5	1.25801×10−5	1.25745×10−4
5	3.55340×10−6	3.5532×10−5	7.29220×10−6	7.29177×10−5

Table 4

L₂ and L_∞ error norms of Example 1 for ν = 0.01, ∆t = 0.005 and h = 0.0125

t	θ = 0.5		θ = 1
	L2	L∞	L2	L∞
5	1.17747×10-4	5.98388×10-4	4.00832×10-4	1.72249×10-3
10	2.79743×10-5	2.05344×10-4	5.50154×10-5	2.98833×10-4
15	1.74834×10-5	1.20216×10-4	1.60754×10-5	1.00908×10-4
20	1.05552×10-5	7.29685×10-5	6.28961×10-6	4.26933×10-5
25	6.57222×10-6	4.59493×10-5	2.68364×10-6	1.87708×10-5

Table 5

L₂ and L_∞ error norms of Example 2 for ν = 1, ∆t = 0.0001, h = 0.01

t	θ = 0.5		θ = 1
	L2	L∞	L2	L∞
0.1	2.83954×10-4	2.00787×10-3	2.88138×10-4	2.02929×10-3
0.2	1.17263×10-4	8.28997×10-4	4.85868×10-5	3.43450×10-4
0.3	4.79970×10-5	3.39416×10-4	3.50799×10-6	2.47984×10-5
0.4	1.95002×10-5	1.37889×10-4	9.36010×10-6	6.61869×10-5
0.5	7.86838×10-6	5.56379×10-5	6.49135×10-6	4.59008×10-5

Table 6

L₂ and L_∞ error norms of Example 2 for ν = 0.1, ∆t = 0.001, h = 0.01

t	θ = 0.5		θ = 1
	L2	L∞	L2	L∞
1	2.60643×10−4	1.80635×10−3	3.51037×10-4	2.10778×10−3
2	9.44114×10-5	6.58582×10−4	3.97978×10-5	2.74071×10-4
3	3.72652×10-5	2.62888×10-4	2.78606×10−6	1.94185×10-5
4	1.50501×10-5	1.06392×10−4	7.22900×10-6	5.10839×10-5
5	6.06745×10-6	4.29058×10-5	5.00608×10-6	3.53980×10-5

Table 7

L₂ and L_∞ error norms of Example 2 for ν = 0.01, ∆t = 0.005 and h = 0.005

t	θ = 0.5		θ = 1
	L2	L∞	L2	L∞
5	5.05467×10-5	5.46951×10-4	5.76664×10-4	3.39188×10-3
10	2.35883×10-5	2.16993×10-4	5.54081×10-5	3.87510×10-4
15	1.14067×10-5	1.06927×10-4	1.05448×10-5	8.37374×10-5
20	6.22943×10-6	6.01938×10-5	2.03732×10-6	1.76742×10-5
25	3.66978×10-6	3.61447×10-5	2.49273×10-7	2.50222×10-6

Table 8

L₂ and L_∞ error norms of Example 3 for ν = l, ∆t = 0.0001 and h = 0.0125

t	θ = 0.5		θ = 1
	L2	L∞	L2	L∞
0.1	1.08102×10-3	6.72107×10-3	1.08885×10-3	6.76936×10-3
0.2	4.41286×10-4	2.78398×10-3	2.29689×10-4	1.45036×10-3
0.3	1.84079×10-4	1.16382×10-3	2.66087×10-5	1.68287×10-4
0.4	7.62220×10-5	4.82117×10-4	1.19130×10-5	7.53275×10-5
0.5	3.12724×10-5	1.97787×10-4	1.25750×10-5	7.95332×10-5

Table 9

L₂ and L_∞ error norms of Example 3 for ν = 0.1, ∆t = 0.001 and h = 0.01

t	θ = 0.5		θ = 1
	L2	L∞	L2	L∞
1	9.14537×10-5	6.34567×10-4	9.22393×10-5	6.39987×10-4
2	3.60692×10-5	2.54355×10-4	1.49080×10-5	1.05242×10-4
3	1.46735×10-5	1.03730×10-4	1.07794×10-6	7.59460×10-6
4	5.95605×10-6	4.21173×10-5	2.85934×10-6	2.02172×10-5
5	2.40313×10-6	1.69930×10-5	1.98256×10-6	1.40191×10-5

Table 10

L₂ and L_∞ error norms of Example 3 for ν = 0.01, ∆t = 0.005 and h = 0.0125

t	θ = 0.5		θ = 1
	L2	L∞	L2	L∞
5	1.22382×10-5	7.46116×10-5	1.51129×10-5	9.16912×10-5
10	7.71923×10-6	4.81380×10-5	7.75865×10-6	4.83819×10-5
15	5.04586×10-6	3.17368×10-5	4.17679×10-6	2.62915×10-5
20	3.32657×10-6	2.09931×10-5	2.26833×10-6	1.43235×10-5
25	2.19179×10-6	1.38506×10-5	1.22415×10-6	7.73911×10-6

In Tables 11 to 13, we have given the comparison of Example 1 at θ = 0.5 and θ = 1 with those of Kutluay et al. [1], Ozis et al [2], Dag et al [3], Salkuyeh et al [4] and Korkmaz et al [5] for Example 1. Tables 14 to 16 give the comparison of numerical results of Example 2 by our method for θ = 0.5 and θ = 1 with those of Kutluay et al [1], Ozis et al [2] and Salkuyeh et al [4]. It can be seen that results by the proposed method are more accurate as compared to the results presented in the above mentioned papers.

Table 11

Comparison of the numerical solutions of Example 1 for θ = 0.5 and θ = 1 when ν = 1 in comparison with the numerical solutions obtained by Kutluay et al. [1], Korkmaz [5] and Dag [3] for h = 0.0125 and ∆t = 0.0001 and exact solutions

x	t	[1]	[5]	[3]	θ = 1	θ = 0.5	Exact solution
0.25	0.4	0.01359	0.01363	0.01355	0.013556	0.013566	0.01357
	0.6	0.00189	0.00190	0.00188	0.001886	0.001881	0.00189
	0.8	0.00026	0.00026	0.00026	0.000262	0.000262	0.00026
	1.0	0.00004	0.00003	0.00004	0.000036	0.000036	0.00004
0.5	0.4	0.01927	0.01932	0.01920	0.019226	0.019227	0.01924
	0.6	0.00268	0.00269	0.00266	0.002671	0.002671	0.00267
	0.8	0.00037	0.00037	0.00037	0.000372	0.000371	0.00037
	1.0	0.00005	0.00005	0.00005	0.000051	0.000051	0.00005
0.75	0.4	0.01365	0.01363	0.01361	0.0136250	0.0136250	0.01363
	0.6	0.00189	0.00190	0.00188	0.0018993	0.0018892	0.00189
	0.8	0.00026	0.00026	0.00026	0.0002624	0.0002623	0.00026
	1.0	0.00004	0.00004	0.00004	0.0000364	0.0000365	0.00004

Table 12

Comparison of the numerical solutions of Example 1 for θ = 0.5 and θ = 1 when ν = 1, h = 0.0125, and ∆t = 0.0001 in comparison with the numerical solutions obtained by Ozis et al. [2] and exact solutions

x	t	[2]	θ = 1	θ = 0.5	Exact solution
0.25	0.1	0.25434	0.25350	0.25350	0.25364
	0.15	0.15641	0.15651	0.15652	0.15660
	0.2	0.09621	0.09638	0.09639	0.09644
	0.25	0.05904	0.05917	0.05919	0.05921
0.5	0.1	0.37087	0.371372	0.371373	0.37158
	0.15	0.22632	0.226693	0.226703	.226820
	0.2	0.13812	0.138386	0.138402	0.13847
	0.25	0.08430	0.084489	0.084496	0.08453
0.75	0.1	0.27073	0.272436	0.272437	0.27258
	0.15	0.16382	0.164279	0.164281	0.16437
	0.2	0.09915	0.099378	0.099384	0.09943
	0.25	0.06016	0.060307	0.060317	0.06034

Table 13

Comparison of the numerical solutions of Example 1 for θ = 0.5 and θ = lwhenν = 0.1 in comparison with the numerical solutions obtained by Kutluay [1], Ozis [2], Korkmaz [5] and Salkuyeh [4] and exact solutions

x	t	[1]	[2]	[5]	[4]	θ = 1	θ = 0.5	Exact solution
0.25	0.4	0.31215	0.31420	0.30910	0.30704	0.30867	0.30867	0.30889
	0.6	0.24360	0.24368	0.24093	0.23898	0.24059	0.24059	0.24074
	0.8	0.19815	0.19753	0.19568	-	0.19557	0.19557	0.19568
	1.0	0.16473	0.16387	0.16274	-	0.16248	0.16248	0.16256
0.5	0.4	0.57293	0.57630	0.56973	0.56867	0.56944	0.56944	0.56963
	0.6	0.45088	0.45165	0.44736	0.44578	0.44703	0.44703	0.44721
	0.8	0.36286	0.36241	0.35943	-	0.35910	0.35909	0.35924
	1.0	0.29532	0.29433	0.29213	-	0.29179	0.29179	0.29192
0.75	0.4	0.63038	0.62598	0.62573	0.62284	0.62557	0.62557	0.62544
	0.6	0.49268	0.49038	0.48760	0.48380	0.48719	0.48719	0.48721
	0.8	0.37912	0.37713	0.37434	-	0.37384	0.37384	0.37392
	1.0	0.29204	0.29014	0.28788	-	0.28738	0.28738	0.28747

Table 14

Comparison of the numerical solutions of Example 2 for θ = 0.5 and θ = 1 when ν = 1, h = 0.0125 and ∆t = 0.0001, in comparison with the numerical solutions obtained by Kutluay et al. [1], Salkuyeh [4] and exact solutions

x	t	[1]	[4]	θ = 1	θ = 0.5	Exact solution
0.25	0.4	0.01403	0.01344	0.013991	0.013997	0.01400
	0.6	0.00195	0.00183	0.001946	0.001948	0.00195
	0.8	0.00027	-	0.000270	0.000271	0.00027
	1.0	0.00004	-	0.000037	0.000038	0.00004
0.5	0.4	0.01988	0.01912	0.019830	0.019840	0.01985
	0.6	0.00276	0.00260	0.002755	0.002756	0.00276
	0.8	0.00038	-	0.000382	0.000383	0.00038
	1.0	0.00005	-	0.0000531	0.000532	0.00005
0.75	0.4	0.01409	0.01350	0.014053	0.014060	0.01407
	0.6	0.00195	0.00183	0.001948	0.001950	0.00195
	0.8	0.00027	-	0.000270	0.000271	0.00027
	1.0	0.00004	-	0.000038	0.000038	0.00004

Table 15

Comparison of the numerical solutions of Example 2 for θ = 0.5 and θ = 1 when ν = 1, h = 0.0125 and∆t = 0.0001, in comparison with the numerical solutions obtained by Ozis [2],Xu and exact solutions

x	t	[2]	θ = 1	θ = 0.5	Exact solution
0.25	0.1	0.26245	0.261334	0.261337	0.26148
	0.15	0.16157	0.16138	0.16140	0.16148
	0.2	0.09948	0.099403	0.099419	0.09947
	0.25	0.06111	0.061043	0.061057	0.06108
0.5	0.1	0.38266	0.383211	0.383212	0.38342
	0.15	0.23351	0.233912	0.233931	0.23406
	0.2	0.14251	0.142791	0.142814	0.14289
	0.25	0.08698	0.087168	0.087189	0.08723
0.75	0.1	0.27974	0.281420	0.281416	0.28157
	0.15	0.16917	0.169634	0.169647	0.16974
	0.2	0.10236	0.102586	0.102602	0.10266
	0.25	0.06211	0.062244	0.062259	0.06229

Table 16

Comparison of the numerical solutions of Example 2 for θ = 0.5 and θ = 1 when ν = 0.1 in comparison with the numerical solutions obtained by Kutluay [1], Ozis [2] and Salkuyeh [4] when h = 0.0125 and ∆t = 0.0001 and exact solutions

x	t	[1]	[2]	[4]	θ = 1	θ = 0.5	Exact solution
0.25	0.4	0.32091	0.32681	0.31247	0.317281	0.317285	0.31752
	0.6	0.24910	0.25120	0.24148	0.245980	0.245983	0.24614
	0.8	0.20211	0.20275	0.19524	0.199443	0.199446	0.19956
	1.0	0.16782	0.16786	0.16153	0.165514	0.165516	0.16560
0.5	0.4	0.58788	0.59661	0.58176	0.584345	0.584338	0.58454
	0.6	0.46174	0.46582	0.45414	0.457800	0.457977	0.45798
	0.8	0.37111	0.37295	0.36283	0.367250	0.367249	0.36740
	1.0	0.30183	0.30256	0.29336	0.298218	0.298219	0.29834
0.75	0.4	0.65054	0.64685	0.63858	0.645764	0.645737	0.64562
	0.6	0.50825	0.50855	0.49362	0.502660	0.502644	0.50268
	0.8	0.39068	0.39120	0.37570	0.385259	0.385250	0.38534
	1.0	0.30057	0.30068	0.28663	0.295764	0.295761	0.29586

Figure 1 gives the graph of numerical solution of Example 1 for different values of v, h and Δt. In Figure 2, we have given the absolute error of our numerical method at θ = 0.5 and θ = 1. Solution profile of Example 2 is given in Figure 3 and absolute error of Example 2 for θ = 0.5 and θ = 1 is presented in Figure 4. Figures 5 and 6 respectively represent the numerical solution and the absolute error at θ = 0.5 and θ = 1 of Example 3. Figures confirm that numerical solutions are in close agreement with the analytical solution. The L₂and L_∞ error norms as well as the plots of absolute errors convey that θ = 0.5 (i.e. Crank-Nicolson scheme) has better accuracy than θ = 1 (i.e. backward Euler scheme)

Fig. 1

Numerical solution of Example lat θ = 0.5 and θ = 1 at various times for space step, h=0.025 and for different values of dt and (v), (a) ν = 1, Δt = 0.0001, (b) ν = 0.1, Δt = 0.001, (c) ν = 0.02, Δt = 0.005, (d) ν = 0.01, Δt = 0.005

Fig. 2

Absolute error of Example 1 at θ = 0.5 and θ = 1 for ν = 0.1, h = 0.0125 and Δt = 0.0001

Fig. 3

Numerical solution of Example 2 at θ = 0.5 and θ = 1 at various times for space step, h=0.025 and for different values of Δt and(v)(a)v = l, Δt = 0.0001, (b)v = 0.1, Δt = 0.001, (c)v = 0.01, Δt = 0.005, (d)v = 0.005, Δt = 0.01

Fig. 4

Absolute error of Example 2 at θ = 0.5 and θ = 1 for ν = 0.1, h = 0.0125 and Δt = 0.0001

Fig. 5

Numerical solution of Example 3 at θ = 0.5 and θ = 1 at various times for space step, h=0.025 and for different values of Δt and (v)(a)v = I, Δt = 0.0001, (b)v = 0.1, Δt = 0.001, (c)v = 0.01, Δt = 0.005, (d)v = 0.005, Δt = 0.01

Fig. 6

Absolute error of Example 3 at θ = 0.5 and θ = 1 for v = 0.005, h = 0.0125 and Δt = 0.0001

The physical behaviour of the solution of Example 1 is given in the surface plots 7 and 8. Figures 9 and 10 give the physical behaviour of the solution to Example 2. Surface plot of Example 3 is given in Figures 11 to 12. It can be seen that as time increases, the solution profile tries to attain steady state.

Fig. 7

Surface plots of Example 1 at θ = 0.5 for space step h = 0.025 and for different values of ν and Δt, (a)v = 1, Δt = 0.002, (b)v = 0.1, Δt = 0.02, (c)v = 0.01, Δt = 0.02, (d)v = 0.005, Δt = 0.01

Fig. 8

Surface plots of Example 1 at θ = 1 for space step h = 0.025 and for different values of ν and Δt, (a)v = l, Δt = 0.0001, (b)v = 0.1, Δt = 0.001, (c) ν = 0.01, Δt = 0.005, (d)v = 0.005, Δt = 0.01

Fig. 9

Surface plots of Example 2 at θ = 0.5 for space step h = 0.025 and for different values of ν and∆t, (a)ν = 1,∆t = 0.001, (b)ν = 0.1, ∆t = 0.01, (c)ν = 0.01, ∆t = 0.01, (d)ν = 0.005, ∆t = 0.01

Fig. 10

Surface plots of Example 2 at θ = 1 for space step h = 0.025 and for different values of ν and∆t, (a)ν = 1,∆t = 0.0001,(b)ν = 0.1, ∆t = 0.001, (c)ν = 0.01, ∆t = 0.005, (d)ν = 0.005, ∆t = 0.01

Fig. 11

Surface plots of Example 3 at θ = 0.5 for space step h = 0.025 and for different values of ν and Δt, (a)v = 1, Δt = 0.002, (b)v = 0.1, Δt = 0.02, (c)v = 0.01, Δt = 0.01, (d)v = 0.005, Δt = 0.01

Fig. 12

Surface plots of Example 3 at θ = 1 for space step h = 0.025 and for different values of ν and Δt, (a) ν = 1, Δt = 0.0001, (b) ν = 0.1, Δt = 0.001, (c)v = 0.01, Δt = 0.005, (d) ν = 0.005, Δt =0.01

7 Conclusion

Burgers’ equation is linearized by using Hopf-Cole transformation and the resultant linear Heat equation is solved by using cubic B-splines. For solving the linear Heat equation, time discretization is done using Crank-Nicolson scheme in first method and backward Euler scheme in second method. L₂ and L_∞error norms are computed and tabulated. The computed results are in good agreement with exact solution. The accuracy of the method is further improved by using Richardson extrapolation. It can be observed from the numerical results that Crank-Nicolson scheme is more accurate compared to backward Euler scheme. Stability analysis, which is performed using Von-Neumann stability method, proved that the proposed schemes are unconditionally stable. Comparison with the works of Kutluay et al. [1], Ozis et al. [2], Dag et al. [3], Salkuyeh et al. [4] and Korkmaz et al. [5] has validated the superiority of the proposed method over these methods. Numerical results establish that the proposed method is consistent and unconditionally stable.