Open access peer-reviewed chapter

# Numerical Simulation of Wave (Shock Profile) Propagation of the Kuramoto-Sivashinsky Equation Using an Adaptive Mesh Method

By Denson Muzadziwa, Stephen T. Sikwila and Stanford Shateyi

Submitted: May 17th 2017Reviewed: October 23rd 2017Published: February 28th 2018

DOI: 10.5772/intechopen.71875

## Abstract

In this paper, the Kuramoto-Sivashinsky equation is solved using Hermite collocation method on an adaptive mesh. The method uses seventh order Hermite basis functions on a mesh that is adaptive in space. Numerical experiments are carried out to validate effectiveness of the method.

### Keywords

• Kuramoto-Sivashinsky equation
• collocation method
• moving mesh partial differential equation
• numerical solution

## 1. Introduction

The Kuramoto-Sivashinsky equation (KSe) is a non-linear fourth order partial differential equation (PDE) discovered separately by Kuramoto and Sivashinsky in the study of non-linear stability of travelling waves. Sivashinsky [1] came up with the equation while modelling small thermal diffusive instabilities in laminar flame fronts. Kuramoto [2, 3, 4, 5] derived the equation in the study of the Belousov-Zhabotinsky reaction as a model of diffusion induced chaos. The KSe is of interest to many researchers because of its ability to describe several physical contexts such as long waves on thin films or on the interface between two viscous fluids [6] and unstable drift waves in plasmas. The equation is also used as a model to describe spatially uniform oscillating chemical reaction in a homogeneous medium and fluctuations in fluid films on inclines [7]. In one dimension, consider the KSe of the form

ut+uux+2ux2+4ux4=0,t>0.E1

The second derivative term is an energy source and thus has a distributing effect. The non-linear term is a correction to the phase speed and responsible for transferring energy. The fourth derivative term is the dominating term and is responsible for stabilising the equation. Several methods have been used to solve the KSe numerically and these include Chebyshev spectral collocation method [8], Quintic B-spline collocation method [9], Lattice Boltzmann method [10], meshless method of lines [11], Fourier spectral method [12] and septic B-spline collocation method [13].

## 2. Grid generation

Generation of an adaptive mesh in the spatial domain is based on the r-refinement technique [14] which relocates a fixed number of nodal points to regions which need high spatial resolution in order to capture important characteristics in the solution. This has the benefit of improving computational effort in those regions of interest whilst using a fixed number of mesh points. The relocation of the fixed number of nodal points at any given time is achieved by solving Moving Mesh Partial Differential Equations (MMPDEs) [15, 16] derived from the Equidistribution Principle (EP). The EP [17] makes use of a measure of the solution error called a monitor function, denoted by M which is a positive definite and user defined function of the solution and/or its derivatives. Mesh points are then chosen by equally distributing the error in each subinterval. In this paper, MMPDE4 [15] is chosen to generate the adaptive mesh because of its ability to stabilise mesh trajectories and ability to give unique solutions for the mesh velocities with Dirichlet boundary conditions. MMPDE4 is given by

ξMẋξ=1τξMxξE2

where τis the relaxation parameter and it plays the role of driving the mesh towards equidistribution. Central finite difference approximation of MMPDE4 in space on the interval axbgives

Mi+1+Mi21N2ẋi+1ẋiMi+Mi121N2ẋiẋi1=Eiτ,E3

where

Ei=Mi+1+Mi21N2xi+1xiMi+Mi121N2xixi1,i=2,,NE4
x1=axN+1=b.E5

The modified monitor function given by

Mxt=1+α2ux2+α22ux2212E6

is used. It is composed of the standard arc-length monitor and the curvature monitor functions. Smoothing on the monitor function is done as described in [15]. Values of the smoothed monitor function M˜at the grid points are given by

M=k=ipi+pMk2γ1+γkik=ipi+pγ1γkiE7

where the parameter pis called the smoothing index which determines the extent of smoothing and is non-negative. γis non-negative and is called the smoothing index and determines the rigidity of the grid.

## 3. Discretization in time

The Crank-Nicolson scheme for the KSe is

un+1unδt+uuxn+1+uuxn2+uxxn+1+uxxn2+uxxxxn+1+uxxxxn2=0E8

where δtis the time step. Rubin and Graves [18] suggested the expression

uuxn+1=un+1uxn+unuxn+1uuxnE9

for the linearization of the non-linear term uuxn+1. Expression (9) is substituted into (1) and the terms are rearranged to give

un+1+δt2un+1uxn+unuxn+1+uxxn+1+uxxxxn+1=unδt2uxxn+uxxxxnE10

## 4. Septic Hermite collocation method

Consider the mesh on the domain abwhich is a solution of MMPDE4 given by

a=X1t<X2t<<XN+1t=bE11

The variable spatial length of each interval is given by Hiwhere Hi=Xi+1tXitfor i=1,,N. For some xϵXitXi+1t, define the local variable sas

s=xXitHitE12

such that 01for every subinterval of the mesh (11). Define the septic Hermite basis functions with the local variables sas

L0,0=20s3+10s2+4s+1s14L0,1=s10s2+4s+1s14L0,2=s224s+1s14L0,3=s36s14L1,0=20s370s2+84s35s4L1,2=s42s124s5L1,3=s46s13E13

For l=0,1,2,3the functions L0,lsand L1,lsyield the following conditions

dkdskL0,l0=δk,l,dkdskL0,l1=0,k,l=0,1,2,3
dkdskL0,l0=0,dkdskL1,l1=δk,l,k,l=0,1,2,3

where δk,ldenotes the Kronecker delta. The physical solution uxton the mesh (11) is approximated by the piecewise Hermite polynomial [19]

U(x,t)=Ui(t)L0,0(s)+Ux,iHi(t)L0,1(s)+Uxx,i(t)Hi2(t)L0,2(s)+Uxxx,i(t)Hi3(t)L0,3(s)+Ui+1(t)L1,0(s)+Ux,i+1Hi(t)L1,1(s)+Uxx,i+1(t)Hi2(t)L1,2(s)+Uxxx,i+1(t)Hi3(t)L1,3(s),E14

Where Uit,Ux,it,Uxx,itand Uxxx,itare the unknown variables. Derivatives of Uxtwith respect to the spatial variable xfor xXitXi+1tare obtained by direct differentiation of (14) to give

(l)U(x,t)x(l)=1Hi(t)(l)[Ui(t)d(l)L0,0ds(l)+Ux,i(t)Hi(t)d(l)L0,1ds(l)+Uxx,i(t)Hi2(t)d(l)L0,2ds(l)+Uxxx,i(t)Hi3(t)d(l)L0,3ds(l)+Ui+1(t)d(l)L1,0ds(l)+Ux,i+1(t)Hi(t)d(l)L1,1ds(l)+Uxx,i+1(t)Hi2(t)d(l)L1,2ds(l)+Uxxx,i+1(t)Hi3(t)d(l)L1,3ds(l)]E15

for l=1,2,3,4.In each subinterval XitXi+1tof the mesh (11), define four Gauss-Legendre points

0<ρ1<ρ2<ρ3<ρ4<1

which are given by

ρ1=12525+703070
ρ2=12525703070
ρ3=1ρ1
ρ4=1ρ2

One regards these points as the collocation points in each subinterval of the mesh (11). Scaling of the Gauss-Legendre points into subsequent intervals is done by defining the collocation points as

Xij=Xi+Hiρj,i=1,.,N,j=1,2,3,4.E16

and redefining the local variable sas

sji=XijXiHiE17

for i=1,,Nand j=1,2,3,4. Evaluation of the Hermite polynomial approximation (14), its first, second and fourth derivatives (15) is then done at the four internal collocation points in each subinterval XiXi+1and substitution of the expressions into (10) gives the difference equation

βj1iUin+1+βj2iUx,in+1+βj3iUxx,in+1+βj4iUxxx,in+1+βj5iUi+1n+1+βj6iUx,i+1n+1+βj7iUxx,i+1n+1+βj8iUxxx,i+1n+1=ψijnE18

where

ψijn=UintL0,0sj+Ux,inHitL0,1sj+Uxx,intHi2tL0,2sj+Uxxx,intHi3tL0,3sj+Ui+1ntL1,0sj+Ux,i+1nHitL1,1sj+Uxx,i+1nHitL1,1sj+Uxxx,i+1nHitL1,1sjδt2Hi2[UintL0,0sj+Ux,inHitL0,1sj+Uxx,intHi2tL0,2sj+Uxxx,intHi3tL0,3sj+Ui+1ntL1,0sj+Ux,i+1nHitL1,1sj+Uxx,i+1nHitL1,1sj+Uxxx,i+1nHitL1,1sj]δt2Hi4[UintL0,0ivsj+Ux,inHitL0,1ivsj+Uxx,intHi2tL0,2ivsj+Uxxx,intHi3tL0,3ivsj+Ui+1ntL1,0ivsj+Ux,i+1nHitL1,1ivsj+Uxx,i+1nHitL1,1ivsj+Uxxx,i+1nHitL1,1ivsj]E19

and

βj1i=L0,0sj+δt2Ux,inL0,0sj+δt2HitUinL0,0sj+δt2Hi2tL0,0sj+δt2Hi2tL0,0ivsjβj2i=HitL0,1sj+δt2Ux,inHitL0,1sj+δt2HitUinHitL0,1sj+δt2Hi2tHitL0,1sj+δt2Hi2tHitL0,1ivsjβj3i=Hi2tL0,2sj+δt2Ux,inHi2tL0,2sj+δt2HitUinHi2tL0,2sj+δt2Hi2tHi2tL0,2sj+δt2Hi4tHi2tL0,2ivsjβj4i=Hi3tL0,3sj+δt2Ux,inHi3tL0,3sj+δt2HitUinHi3tL0,3sj+δt2Hi2tHi3tL0,3sj+δt2Hi4tHi3tL0,3ivsjβj5i=L1,0sj+δt2Ux,inL1,0sj+δt2HitUinL1,0sj+δt2Hi2tL1,0sj+δt2Hi4tL1,0ivsjβj6i=HitL1,1sj+δt2Ux,inHitL1,1sj+δt2HitUinHitL1,1sj+δt2Hi2tHitL1,1sj+δt2Hi4tHitL1,1ivsjβj7i=Hi2tL1,2sj+δt2Ux,inHi2tL1,2sj+δt2HitUinHi2tL1,2sj+δt2Hi2tHi2tL1,2sj+δt2Hi4tHi2tL1,2ivsjβj8i=Hi3tL1,3sj+δt2Ux,inHi3tL1,3sj+δt2HitUinHi3tL1,3sj+δt2Hi2tHi3tL1,3sj+δt2Hi4tHi3tL1,3ivsjE20

From the boundary conditions (28) and (29), one gets

Ux1=σUxx1=βUxN+1=ωUxxN+1=ζE21

which results in a consistent system of 4N+4equations in 4N+4unknowns.

## 5. Solution approach for the PDE system

The PDE system is solved using the rezoning approach which works best with the decoupled solution procedure [20]. The rezoning approach allow varying criteria of convergence for the mesh and physical equation since in practice the mesh does not require the same level of accuracy to compute as compared to the physical solution. The algorithm for the rezoning approach is as follows:

1. Solve the given physical PDE on the current mesh.

2. Use the PDE solution obtained to calculate the monitor function.

3. Find the new mesh by solving a MMPDE.

4. Adjust the current PDE solution to suite the new mesh by interpolation.

5. Solve the physical PDE on the new mesh for the solution in the next time.

## 6. Solution adjustment by interpolation

Discretization of the time domain tatbis done using the following finite sequence

ta=t0<<tn<<tk=tbE22

At each time t=tn=n×dt, consider a non-uniform spatial mesh Xini=1N+1given by

a=X1n<<XN+1n=bE23

where Xin=Xitnwith Hin=Xi+1nXinbeing a non-uniform spatial step for i=1,,N. At the same time step t=tnone also considers the approximations to the exact solution uxtand its derivatives given by Uini=1N+1and Uilni=1N+1respectively where Uilnrepresents the lthderivative approximation with respect to the variable xat the time t=tnFor l=1,2,3. A new mesh X˜ini=1N+1is generated by (2) at each current time step tn. The goal is to determine the new approximations U˜ini=1N+1and U˜ilni=1N+1which are related to the new mesh X˜ini=1N+1in a similar manner the approximations Uini=1N+1and Uilni=1N+1are related to the old mesh Xini=1N+1in each subinterval XiXi+1. This process of updating the solution and its derivatives from the old mesh to the new mesh is achieved by interpolation. One considers the septic Hermite interpolating polynomial, a piecewise polynomial which allows the function values and its three consecutive derivatives to be satisfied in each subinterval XiXi+1. The Hermite polynomial (14) is written in compact form as

l=03hllUilL0,ls+l=03hllUi+1lL1,lsE24

where the 4N+1unknowns are given by

Uii=luxlXittUi+1lluxlXi+1tt,l=0,1,2,3.

Given the partition (23) and approximations Uilnfor l=0,1,2,3, suppose interpolation of Ulxis required at x=Xi˜nwhere Xi˜nXinXi+1nfor i=1,,N. Firstly, the local coordinate sof Xi˜nis defined as

s=Xi˜nXinHinE25

U˜lX˜inis then defined as

U˜lX˜in=i=03HilpUildlL0,lsdsl+i=03HilpUi+1ldlL0,lsdslE26

for l=0,1,2,3to give the interpolated values of U˜and the first three consecutive derivatives on the new subinterval X˜inX˜i+1n. In order to compute the approximations of Uat the next time step t=tn+1denoted by Uini=1N+1, the values of the new mesh X˜ini=1N+1and the updated approximations U˜ini=1N+1are used in a septic Hermite collocation numerical scheme. The new approximations Uin+1i=1N+1and the new mesh X˜in+1i=1N+1become the starting conditions for repeating the whole adaptive process.

## 7. Numerical results

Consider the KSe

ut+uux+2ux2+4ux4=0,t>0E27

in the domain 3030,t>0with boundary conditions

u30t=σ,ux30t=βE28
u30t=ω,ux30t=ζE29

Where σ,β,ωand ζare obtained from the exact solution

uxt=c+151911199tanh3kxctx0+11tanhkxctx0E30

With c=0.1,x0=12and k=121119.

Figures 1 and 2 show the behaviour of the numerical solution and the absolute error, respectively of the KSe equation on a stationary mesh using Hermite collocation method at t=4with N=100and δt=0.001. In Figure 1, one observes that the numerical solution tracks the exact solution with the absolute error variation as shown in Figure 2.

Figure 3 shows the solution obtained by the collocation method on a stationary mesh for time t=0,1,2,3,4. The movement of the solution is from left to right as time increases and the solution tracks the exact solution with no oscillations. One also observes that the concentration of mesh points is higher in the flatter regions of the solution profile in comparison to the concentration in the steeper region.

Figures 4 and 5 show the numerical solution profile and the behaviour of the maximum absolute error, respectively at t=4with N=100,δt=0.001and α=8on an adaptive mesh. In Figure 4, one observes that the numerical solution is able to track the exact solution and the distribution of mesh points is almost equal along the solution profile which enables resolution of the solution with minimum errors.

Figure 6 shows the numerical solution profiles produced by the adaptive collocation method for time t=0,1,2,3,4. One observes that the solution moves from left to right as time progresses. The mesh points at different times keep on tracking the solution profile and maintain an almost equal distribution along the profile up to final time T=4. Figure 7 shows the paths taken by the mesh points in tracking the solution profile. In Table 1, the infinity norm error for an adaptive collocation method is calculated and results are compared with the method in [13]. Results show improvements in the maximum point wise errors when an adaptive Hermite collocation method is used.

TimeHermite collocationMethod in [19]
0.59.0×1041.03619×103
11.4×1031.63762×103
1.51.9×1032.07273×103
21.7×1032.48375×103
2.52.0×1032.79434×103
32.1×1033.00439×103
3.52.1×1033.16038×103
42.1×1033.43704×103

### Table 1.

Comparison of maximum pointwise errors in the numerical solution of the KSe on an adaptive mesh at different times with δt=0.001and N=100.

## 8. Conclusions

The KSe is solved using an adaptive mesh method with discretization in the spatial domain done using seventh order Hermite basis functions. Numerical results show that Hermite collocation method on a non-uniform adaptive mesh is able to improve the accuracy of the numerical solution of the KSe. The method is able to keep track of the region of rapid solution variation in the KSe, which is one of the desired properties of an adaptive mesh method.

chapter PDF
Citations in RIS format
Citations in bibtex format

## More

© 2018 The Author(s). Licensee IntechOpen. This chapter is distributed under the terms of the Creative Commons Attribution 3.0 License, which permits unrestricted use, distribution, and reproduction in any medium, provided the original work is properly cited.

## How to cite and reference

### Cite this chapter Copy to clipboard

Denson Muzadziwa, Stephen T. Sikwila and Stanford Shateyi (February 28th 2018). Numerical Simulation of Wave (Shock Profile) Propagation of the Kuramoto-Sivashinsky Equation Using an Adaptive Mesh Method, Finite Element Method - Simulation, Numerical Analysis and Solution Techniques, Răzvan Păcurar, IntechOpen, DOI: 10.5772/intechopen.71875. Available from:

### Related Content

Next chapter

#### Numerical Analysis on the Simulated Heavy Rainfall Event of Tropical Cyclone Fung-Wong

By Lei-Ming Ma and Xu-Wei Bao

First chapter

#### Finite Element and Finite Difference Methods for Elliptic and Parabolic Differential Equations

By Aklilu T. G. Giorges

We are IntechOpen, the world's leading publisher of Open Access books. Built by scientists, for scientists. Our readership spans scientists, professors, researchers, librarians, and students, as well as business professionals. We share our knowledge and peer-reveiwed research papers with libraries, scientific and engineering societies, and also work with corporate R&D departments and government entities.

View all Books