# On a seventh order convergent weakly L-stable Newton Cotes formula with application on Burger's equation

In this paper we derive 7^th order convergent integration formula in time which is weakly L-stable. To derive the method we use, Newton Cotes formula, fifth-order Hermite interpolation polynomial approximation (osculatory interpolation) and sixth-order explicit backward Taylor's polynomial approximation. The vector form of this formula is used to solve Burger's equation which is one dimensional form of Navier-Stokes equation. We observe that the method gives high accuracy results in the case of inconsistencies as well as for small values of viscosity, e.g., 10^-3. Computations are performed by using Mathematica 11.3. Stability and convergence of the schemes are also proved. To check the efficiency of the method we considered 6 test examples and several tables and figures are generated which verify all results of the paper.

• 4 publications
• 1 publication
• 1 publication
08/05/2021

### The Key Equation for One-Point Codes

For Reed-Solomon codes, the key equation relates the syndrome polynomial...
10/25/2020

### Backward difference formula: The energy technique for subdiffusion equation

Based on the equivalence of A-stability and G-stability, the energy tech...
11/13/2019

### Haar wavelets collocation on a class of Emden-Fowler equation via Newton's quasilinearization and Newton-Raphson techniques

In this paper we have considered generalized Emden-Fowler equation, y”(...
10/13/2017

### Compact Formulae in Sparse Elimination

It has by now become a standard approach to use the theory of sparse (or...
07/04/2019

### Second-order semi-implicit projection methods for micromagnetics simulations

Micromagnetics simulations require accurate approximation of the magneti...
04/06/2022

### The Challenge of Sixfold Integrals: The Closed-Form Evaluation of Newton Potentials between Two Cubes

The challenge of explicitly evaluating, in elementary closed form, the w...
04/20/2021

### Fully Implicit Spectral Boundary Integral Computation of Red Blood Cell Flow

An approach is presented for implicit time integration in computations o...

## 1 Introduction

Burgers’ equation with as coefficient of viscosity can be defined as

 ∂w∂t+w∂w∂x−νd2∂2w∂x2=0,(x,t)∈ΣT, (1.1)

where

 ΣT=(α0,α1)×(0,T], T>0,

with Dirichlet boundary conditions (BCs),

 w(αi,t)=0,  i=0,1  and  t∈(0,T], (1.2)

and initial conditions (ICs),

 w(x,0)=f(x),x∈(α0,α1). (1.3)

Linearized form of Burgers’ equation (by using Hopf-Cole transformation) is given as

 ∂ψ∂t=νd2∂2ψ∂x2 (1.4)

with the Neumann boundary conditions

 ψx(αi,t)=0,  i=0,1,

and the initial conditions

 ψ(x,0)=g(x).

The study of Burger’s equation is popular among the scientific community as it is very simple form of Naviers’ Stokes equation and due to it’s appearence in various field of applied mathematics and physics such as in the context of gas dynamics, in the theory of shock waves, traffic flow, mathematical modeling of turbulent fluid and in continuous stochastic processes. In , it was first introduced by Bateman [5]. Later in , it was introduced by Burger [6, 7] as a class of equation which delineate the mathematical model of turbulence. Recently in , Ryu etc. [36] propose some nowcasting rainfall models based on Burger’s equation. This equation has been solved analytically for some initial condition and solution is represented in the form of Fourier series expansion which converges slowly for small values of viscosity. Exact solution does not work very well for small values of viscosity and hence it always attracts researchers to test newly devloped numerical method on this nonlinear parabolic PDEs.

Recently, with development in computer speed several numerical schemes based on finite difference method, finite element method, spectral method, differential quadrature method, decomposition method, moving least squares particle method, Haar wavelet quasilinearization approach etc., have been developed to solve the Burger’s equation [8, 37, 2, 35, 31, 19, 47, 26, 1, 15, 3, 23, 24, 25, 33, 32, 17, 21, 22, 4, 18, 39, 40, 38, 16, 30, 45, 46].

Crank-Nicolson (CN) method [14, 41, 42, 28] is a second order method which is based on Trapezoidal formula which is A-stable but not L-stable. In the presence of inconsistencies [34] CN produces unwanted oscillations. Chawla etc. [11] produces generalised Trapezoidal formula (GTF()), where which is L-stable and gives a quite stable result. Chawla etc. [10] proposed a modified Simpson’s rule (ASIMP) which is A-stable and used to give fourth-order time integration formula but it produces unwanted oscillations like CN due to lack of L-stability. To remove this oscillation, Chawla etc. [12] produced L-stable version of Simpson’s rule and implemented it to derive a third-order time integration formula for the diffusion equation which gives stable and accurate result. Lajja [44] proposed L-stable derivative free error corrected Trapezoidal rule (LSDFECT). Verma etc. [43] developed a fifth order time integration formula for the diffusion equation which is weakly L-stable.

Here we derive order time integration formula which is weakly -stable and generalize above mentioned existing results. The issue of slow convergence of series solution for small forces analytical solution of Burgers’ equation to deviate from the exact solution. So, it is not easy to compute the solution for small values of . The newly developed method computes the solution even for small values of . To compute the numerical solution we use Mathematica 11.3 and find out that numerical solutions are in good agreement for small values of . The result are in good agreement with exact solution when inconsistencies are present in the initial and boundary condition.

The paper is organized as follows. In section 2, we give close form solution which we use to compute exact solution. In section 3, we derive higher order integration method in time for . In section 4, we use this technique combined with finite difference to solve Burgers’ equation and demonstrate the stability. In section 5, we illustrate the numerical results with tables and 2D-3D graphs.

## 2 Close Form Solution

Hopf [20] and Cole [13] gave idea that the equation (1.1) can be reduced to the following linear heat equation

 ∂ψ∂t=νd2∂2ψ∂x2, (2.1)

with the Neumann boundary condition (BC)

 ψx(αi,t)=0,  αi=i,  i=0,1, (2.2)

and the initial condition (IC)

 ψ(x,0)=g(x), (2.3)

by non-linear transformation

 ϕ=−νd(logψ),ϕ=ϕ(x,t), (2.4)

and

 w=ϕx. (2.5)

The analytical solution of the linearized heat equation (2.1) is given by

 ψ(x,t)=β0+∞∑l=1βlexp(−νdl2π2t2)cos(lπx), (2.6)

where and are Fourier coefficient and is given by

 β0 = ∫10exp(−1νd∫x0w0(ξ)dξ)dx, βl = 2∫10exp(−1νd∫x0w0(ξ)dξ)cos(lπx)dx,

where .

The analytical solution by Hopf-Cole transformation is

 w(x,t)=πνd∑∞l=1βlexp(−νdl2π2t2)lsin(lπx)β0+∑∞l=1βlexp(−νdl2π2t2)cos(lπx). (2.7)

## 3 Illustration of the proposed method

We consider the initial value problem

 u′(t)=f(t,u),u(t0)=η0. (3.1)

The Newton Cotes time integration formula is given by

 un+1=un+h840(41fn+216fn+1/6+27fn+2/6+272fn+3/6+27fn+4/6+216fn+5/6+41fn+1). (3.2)

Now, we use the fifth order Hermite approximation for , , , , which are given by

 un+1/6=115552(1500yn+552un+1+2250hu′n−210hu′n+1+125h2u′′n+25h2u′′n+1), (3.3) un+2/6=1243(192un+51un+1+48hu′n−18hu′n+1+4h2u′′n+2h2u′′n+1), (3.4) (3.5) un+4/6=1243(51un+192un+1+18hu′n−48hu′n+1+2h2u′′n+4h2u′′n+1), (3.6) un+5/6=115552(552un+15000un+1+210hu′n−2250hu′n+1+25h2u′′n+125h2u′′n+1), (3.7)

and sixth order Taylor’s approximation

 un=un+1−hu′n+1+h22u′′n+1−h36u′′′n+1+h424uivn+1−h5120uvn+1, (3.8)

to get

 ¯¯¯¯¯¯¯¯¯¯¯¯¯¯un+1/6 = 146656[44875un+1781un+1+6750hu′n+375h2u′′n−755hu′n+1+2752h2u′′n+1 (3.9) +125(−h36u′′′n+1−h46uivn+1−h56uvn+1)], ¯¯¯¯¯¯¯¯¯¯¯¯¯¯un+2/6 = 1729[568un+161un+1+144hu′n+12h2u′′n−62hu′n+1+10h2u′′n+1+ (3.10) 8(−h36u′′′n+1−h46uivn+1−h56uvn+1)], ¯¯¯¯¯¯¯¯¯¯¯¯¯¯un+3/6 = 164[31un+33un+1+10hu′n+h2u′′n−11hu′n+1+32h2u′′n+1+(−h36u′′′n+1− (3.11) h46uivn+1−h56uvn+1)],
 ¯¯¯¯¯¯¯¯¯¯¯¯¯¯un+4/6 = 1729[145un+584un+1+54hu′n+6h2u′′n−152hu′n+1+20h2u′′n+1+ (3.12) 8(−h36u′′′n+1−h46uivn+1−h56uvn+1)], ¯¯¯¯¯¯¯¯¯¯¯¯¯¯un+5/6 = 146656[1531un+45125un+1+630hu′n+75h2u′′n−6875hu′n+1+8752h2u′′n+1 (3.13) +125(−h36u′′′n+1−h46uivn+1−h56uvn+1)].

Now, we define

 ¯¯¯¯¯¯¯¯¯¯¯¯¯¯fn+1/6=f(xn+1/6,¯¯¯¯¯¯¯¯¯¯¯¯¯¯un+1/6), (3.14) ¯¯¯¯¯¯¯¯¯¯¯¯¯¯fn+2/6=f(xn+2/6,¯¯¯¯¯¯¯¯¯¯¯¯¯¯un+2/6), (3.15) ¯¯¯¯¯¯¯¯¯¯¯¯¯¯fn+3/6=f(xn+3/6,¯¯¯¯¯¯¯¯¯¯¯¯¯¯un+3/6), (3.16) ¯¯¯¯¯¯¯¯¯¯¯¯¯¯fn+4/6=f(xn+4/6,¯¯¯¯¯¯¯¯¯¯¯¯¯¯un+4/6), (3.17) ¯¯¯¯¯¯¯¯¯¯¯¯¯¯fn+5/6=f(xn+5/6,¯¯¯¯¯¯¯¯¯¯¯¯¯¯un+5/6). (3.18)

Now, the time integral formula (3.2) for the interval becomes

 un+1=un+h840(41fn+216¯¯¯¯¯¯¯¯¯¯¯¯¯¯fn+1/6+27¯¯¯¯¯¯¯¯¯¯¯¯¯¯fn+2/6+272¯¯¯¯¯¯¯¯¯¯¯¯¯¯fn+3/6+27¯¯¯¯¯¯¯¯¯¯¯¯¯¯fn+4/6+216¯¯¯¯¯¯¯¯¯¯¯¯¯¯fn+5/6+41fn+1). (3.19)

### 3.1 Local trunction error

Using Taylor’s series expansion, we have

 un+1/6 = 115552[1500un+552un+1+2250hu′n−210hu′n+1+125h2u′′n (3.20) +25h2u′′n+1]−25h66718464uvin−475h7282175488uviin+O(h8), un+2/6 = 1243[192un+51un+1+48hu′n−18hu′n+1+4h2u′′n+2h2u′′n+1] (3.21) −h665610uvin−h7137781uviin+O(h8), un+3/6 = 164[32un+32un+1+10hu′n−10hu′n+1+h2u′′n+h2u′′n+1] (3.22) −h646080uvin−h792160uviin+O(h8), un+4/6 = 1243[51un+192yn+1+18hu′n−48hu′n+1+2h2u′′n+4h2u′′n+1)] (3.23) −h665610uvin−11h71377810uviin+O(h8), un+5/6 = 115552[552un+15000un+1+210hu′n−2250hu′n+1+25h2u′′n (3.24) +125h2u′′n+1]−25h66718464uvin−575h7282175488uviin+O(h8), un = un+1−hu′n+1+h22u′′n+1−h36u′′′n+1+h424uivn+1−h5120uvn+1+h6720uvin (3.25) +h7840uviin+O(h8),

then it follows that

 un+1/6=¯¯¯¯¯¯¯¯¯¯¯¯¯¯un+1/6+425h7282175488uviin+O(h8), (3.26) un+2/6=¯¯¯¯¯¯¯¯¯¯¯¯¯¯un+2/6+4h7688905uviin+O(h8), (3.27)
 un+3/6=¯¯¯¯¯¯¯¯¯¯¯¯¯¯un+3/6+h7129024uviin+O(h8), (3.28) un+4/6=¯¯¯¯¯¯¯¯¯¯¯¯¯¯un+4/6+h7196830uviin+O(h8), (3.29) un+5/6=¯¯¯¯¯¯¯¯¯¯¯¯¯¯un+5/6+325h7282175488uviin+O(h8). (3.30)

Also, we have

 un+1 = un+h840[41fn+216fn+1/6+27fn+2/6+272fn+3/6+27fn+4/6 (3.31) +216fn+5/6+41fn+1]−h9u(9)1567641600.

From all of the above, we deduce that

 un+1 = un+h840(41fn+216¯¯¯¯¯¯¯¯¯¯¯¯¯¯fn+1/6+27¯¯¯¯¯¯¯¯¯¯¯¯¯¯fn+2/6+272¯¯¯¯¯¯¯¯¯¯¯¯¯¯fn+3/6+27¯¯¯¯¯¯¯¯¯¯¯¯¯¯fn+4/6 (3.32) +216¯¯¯¯¯¯¯¯¯¯¯¯¯¯fn+5/6+41fn+1)+tn(h),

where

 tn(h)=O(h8).

Thus the scheme (3.19) is seventh order convergent.

### 3.2 Stability of the formula (3.19)

Consider the test problem

 u′(t)=−λu(t),λ>0 (3.33)

and assume , then we have

 un+1=Ψ(s)un, (3.34)

where

 Ψ(s)=540(840−414s+82s2−7s3)453600+230040s+48600s2+5480s3+540s4+135s5+27s6.

From figure 1 it can be seen that and hence our scheme is not -stable. Since as and hence scheme is weakly -stable.

### 3.3 Stability Region for the formula (3.19)

We use the boundary locus method [29, p.64, chapter 7] and determine the boundary of the region. It can be easily seen that outside of the region it is unconditionally stable.

## 4 Solution of the Burgers’ equation

### 4.1 The final scheme

We discretize the solution space with uniform mesh expressed as For that, we partition the interval in to equal sub intervals with the spatial grid , where is a positive integer and is the spatial step.

Also partition the interval in to equal subintervals with the temporal grid where and is a positive integer.

Now define and consider linearized Burger’s equation (2.1) and compute the solution for a given and for on . Then we use (2.4)-(2.5) to deduce the following formula for computing the which is the solution of the nonlinear Burgers’ equation (1.1),

 w(xi,t)=(−νd2h)ψ(xi+h,t)−ψ(xi−h,t)ψ(xi,t).

Here we approximate second order spatial derivative by fourth order finite difference ratio which is given by

 ∂2ψ(x,t)∂2x≈16(ψ(x+h,t)+ψ(x−h,t))−30ψ(x,t)−(ψ(x+2h,t)+ψ(x−2h,t))12h2,

and convert the linearized Burgers’ equation into an initial value problem in vector form.

Now, we apply finite difference discretization on (2.1) with the Neumann boundary conditions

 ψx(αi,t)=0,  i=0,1,

we get the following equation

 ∂Ψ(t)∂t=−νd24h2DΨ(t), (4.1)

where and is the pentadiagonal matrix given by

 D=⎛⎜ ⎜ ⎜ ⎜ ⎜ ⎜ ⎜ ⎜ ⎜ ⎜ ⎜ ⎜⎝30−32200⋯00−1631−1610⋯001−1630−1610⋯0⋮⋮⋱⋱⋱⋱⋮0001−1630−16100001−1631−16000002−3230⎞⎟ ⎟ ⎟ ⎟ ⎟ ⎟ ⎟ ⎟ ⎟ ⎟ ⎟ ⎟⎠ (4.2)

and
Let , then applying the time integration formula on the initial value problem (4.1), we get

 Ψj+1 = Ψj−ρ840D(41Ψj+216¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯Ψj+1/6+27¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯Ψj+2/6+272¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯Ψj+3/6 (4.3) +27¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯Ψj+4/6+216¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯Ψj+5/6+41Ψj+1),

where

 ¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯Ψj+1/6 = 146656[(44875I−6750ρD+375ρ2D2)Ψj+(1781I+755ρD (4.4) ¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯Ψj+2/6 = 1729[(568I−144ρD+12ρ2D2)Ψj+(161I+62ρD+10ρ2D2 (4.5) +8(ρ3D36+ρ4D424+ρ5D5120))Ψj+1], ¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯Ψj+3/6 = 164[(31I−10ρD+ρ2D2)Ψj+(33I+11ρD+32ρ2D2+(ρ3D36 (4.6) +ρ4D424+ρ5D5120))Ψj+1], ¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯Ψj+4/6 = 1729[(145I−54ρD+6ρ2D2)Ψj+(584I+152ρD+20ρ2D2 (4.7) +8(ρ3D36+ρ4D424+ρ5D5120))Ψj+1], ¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯Ψj+5/6 = (4.8)

Now, use in Equation (4.3), we get

 (453600I+2300ρD+48600ρ2D2+5480ρ3D3+540ρ4D4+135ρ5D5+27ρ6D6) ¯¯¯¯¯¯¯¯¯¯¯Ψj+1 =540(840I−414ρD+84ρ2D2−7ρ3D3) ¯¯¯¯¯¯Ψj. (4.9)

This method is BY using (4.3) we can compute and hence is computed at different ’s for a given time level . Physical properties of the solutions are discussed later in form of figures and tables.

### 4.2 Unconditional stability of the scheme for the heat equation

Equation (4.3) can be written as

 Ψj+1=PΨj,

where

 P = 540(840I−414ρD+84ρ2D2−7ρ3D3)(453600I+230040ρD+48600ρ2D2+5480ρ3D3+540ρ4D4+135ρ5D5+27ρ6D6) (4.10) = L−11L2, (say),

where

 L1=(453600I+230040ρD+48600ρ2D2+5480ρ3D3+540ρ4D4+135ρ5D5+27ρ6D6) L2=540(840I−414ρD+84ρ2D2−7ρ3D3).
###### Lemma 4.1.

The matrix is similar to a symmetric matrix.

###### Proof.

Let us introduce a diagonal matrix

 Q=⎛⎜ ⎜ ⎜ ⎜ ⎜ ⎜ ⎜ ⎜⎝√21⋱1√2⎞⎟ ⎟ ⎟ ⎟ ⎟ ⎟ ⎟ ⎟⎠

such that

 ~D=Q−1DQ,

i.e., is similar to a symmetric matrix .
Now, we will show that is similar to symmetric matrix. Let

 ~P=Q−1PQ=Q−1L−11L2Q=[Q−1L−11Q][Q−1L2Q]=[Q−1L1Q]−1[Q−1L2Q]=~L−11~L2.
 ~L1 = (453600I+230040ρ~D+48600ρ2~D2+5480ρ3~D3+540ρ4~D4+135ρ5~D5+27ρ6~D6) ~L2 = 540(840I−414ρ~D+84ρ2~D2−7ρ3~D3)

but matrices and are symmetric and commute and therefore is similar to a symmetric matrix and therefore all the eigen values of the matrix are real. ∎

###### Lemma 4.2.

Eigen values of the matrix are

###### Proof.

Let be the eigen vectors of the matrix corresponding to the eigen value Then, we have

 (30−λl)v1−32v2+2v3=0,
 −16v1+(31−λl)v2−16v3+v4=0,
 vj−2−16vj−1+(30−λl)vj−16vj+1+vj+2=0,  j=2,3,4,...,N−1,
 vN−2−16vN−1+(31−λl)vN−16vN+1=0,
 2vN−1−32vN+(30−λl)vN+1=0.

We set then we get fourth order difference equation

 vj−2−16vj−1+(30−λl)vj−16vj+1+vj+2=0,  j=1,2,...,N+1, (4.11)

with the boundary conditions . The characteristc equation of the equation (4.11) is . Assume are the characteristc roots then we have

 m1+m2+m3+m4=16,
 m