Finite Difference Methodfor(2+1)-Kuramoto-Sivashinsky Equation

BEZIA Abdelhamidand MABROUK Anouar Ben

1Algebra and Number Theory Laboratory,Faculty of Mathematics,University of Sciences and Technology Houari Boumediene,BP 32 EL-Alia 16111,Bab Ezzouar,Algiers,Algeria.

2D′epartement de Math′ematiques,Institut Sup′erieur de Math′amatiques Appliqu′ees et Informatique de Kairouan,Avenue Assad Ibn Al-Fourat,Kairouan 3100,Tunisia.

3Laboratory of Algebra,Number Theory and Nonlinear Analysis.Department of Mathematics,Faculty of Sciences,University of Monastir,5000 Monastir,Tunisia.

4Department of Mathematics,College of Sciences,Tabuk University,Saudi Arabia Kingdom.

Abstract.This paper investigates a solution technique for solving a two-dimensional Kuramoto-Sivashinsky equation discretized using a finite difference method.It consists of an order reduction method into a coupled system of second-order equations,and to formulate the fully discretized,implicit time-marched system as a Lyapunov-Sylvester matrix equation.Convergence and stability is examined using Lyapunov criterion and manipulating generalized Lyapunov-Sylvester operators.Some numerical implementations are provided at the end to validate the theoretical results.

Key Words:Kuramoto-Sivashinsky equation,Finite difference method,Lyapunov-Sylvester operators.

1 Introduction

Kuramoto-Sivashinsky(KS)equation is one of the well known models for for chaotic spatially extended systems[9].The KS equation arises in the description of stability of flame fronts,reaction-diffusion systems and many other physical settings[10,12].Similarly to[20],in the context of the present paper the two-dimensional generalized KS equation describes the evolution of a(2+1)-dimensional surface defined as a function on a two-dimensional plane and is growing in the direction perpendicular to that plane.Therefore,the present paper is devoted to the development of a computational method based on two-dimensional finite difference scheme to approximate the solution of the nonlinear KS equation

with initial conditions

and boundary conditions

on a rectangular domain ?=[L0,L1]×[L0,L1]in R2,t0≥0 is a real parameter fixed as the initial time.is the time derivative,? is the space gradient operator and Δ =+is the Laplace operator in R2,q,κ,λ are real parameters. ? and ψ are twice differentiable real valued functions on ?.

We propose to apply an order reduction of the derivation and thus to solve a coupled systemof equation involving second order differential operators.We set v=qu?κΔu and thus we have to solve the system

The Kuramoto-Sivashinsky equation(KS)is one of the most famous equations in mathphysicsformany decades.Ithas itsoriginintheworkofKuramotosince the70-thdecade of the 20-th century in his study of reaction-diffusion equation[21].The equation was then considered by Sivashinsky in modeling small thermal diffusion instabilities for laminar flames and modeling the reference flux of a film layer on an inclined plane[32,33].Since then the KS equation has experienced a growing developmentin theoretical mathematics,numerical as well as physical mechanics,nonlinear physics,hydrodynamics[28],incombustiontheory,chemistry,plasma physics,particle distributionsadvection,surface morphology,...etc.

Forexample,in[1],an anisotropic versionof the KS equationhas beenproposedleading to global resolutions of the equation on rectangular domains.Sufficient conditions were given for the existence of global solution.See for example[6–13].

In[14],author study the motion of passive tracers in a two-dimensional turbulent velocity field generated by the Kuramoto-Sivashinsky equation.In[15],using the numerical simulation approach we study the evolution of(2+1)-dimensional surface morphology in KS model.

The long-wavelength properties of the Kuramoto-Sivashinsky equation are studied in(2+1)dimensions using numerical and analytic techniques.See[11]and[18].using numerical and analytic techniques Nadja fikhah and Ahangari[23],based on the theory of Lie algebras,they have studied the problem of determining the largest possible set of symmetries(K-S)model in two spatial and one temporal dimensions.

From the dimensional point of view,KS equation has been widely studied in one dimension,[10,12,16,24,25,34].However,this equation since its appearance is related to the modeling of flame spread which is a two-dimensionalproblem.Inthis context,an important result was developed in[30]where the authors showed by adapting the method developed in[26]the existence of a set of bounded solutions on a rectangular domain.This major importance of the two-dimensional model was a main motivation behind this work.A model representing a nonlinear dynamical systemdefined in a two-dimensional space is considered,where the solution u(x,y,t)satisfies a fourth order partial differential equation of the form(1.1),where u is the height of the interface,q is a pre-factor proportional to the coefficient of surface tension expressed by q?2u.The quantity κ?4u represents the result of the diffusion surface due to the chemical potential gradient induced curvature.The pre-factor κ represents the surface diffusion term.The quantity λ|?u|2is due to the existence of overhangs and vacancies during the deposition process.Finally,the quantity q?2u+λ‖?u‖2is referred to as modeling the effect of deposited atoms[15].

In the present work we propose to serve of algebraic operators to approximate the solutions ofthe Kuramoto-Sivashinsky(KS)equationin the two spatial and one temporal dimension without adapting classical developments based on separation of variables,radial solutions,tridiagonal operators,...etc.

This was onefold of the present paper.A second crucial idea is to transform the continuous K-S equation into a generalized Lyapunov-Sylvester equation of the form

where Aiand Biare appropriate matrices depending on the discretization procedure and the problem parameters.Xnrepresents the numerical solution at time n and Cnis usually depending on the past values Xk,k≤n?1 of the the solution.The equation(1.5)is known as generalized Lyapunov-Sylvester equation.Such equations have their origin in the work of Sylvester on classical matrices equations.In the particular case

The equation is known restrictively as Lyapunov one[31].Generally speaking,the equation

is very difficult to be inverted and remains an open problem in algebra.Nevertheless,some works have been developed and proved that under suitable conditions on the coefficient matrices,one may get a unique solution,but it’s exact computation remains hard.It necessitates to compute eigenvalues and precisely bounds/estimates of eigenvalues or direct inverses of big matrices which remains a complicated problem and usually inappropriate.

In[22],a native method to solve(1.7)is investigated based on Kronecker product and equivalent matrix-vector equation.The Sylvester’s equation is transformed into a linear equivalent one on the form Gx=c,with a matrix G obtained by tensor products issued from theand the.However,the general case remained already complicated.The authors have been thus restricted to special cases where the matrices Ajand Bjare scalar polynomials based on spacial and fixed matrices A and B.Denote by σ(A)the spectrum of A and σ(B)the one of B,the spectrum σ(G)may then be determined in terms these spectra.Indeed,with the assumptions theand the,the equation(1.7)can be written as

where αjkare complex numbers.Hence,the tensor matrix G will be written on the form G=φ(A,B),where φ is the 2-variables polynomial φ(x,y)=∑j,kαj,kxiyk.Thus,G is singular if and only if φ(λ,μ)=0 for some eigenvalues λ and μ of A and B respectively.

In the case of square matrices an other criterion of existence of the solution of(1.7)was pointed out by Roth[29].It was shown that the solution is unique if and only if the matrices

are similar.

However,innumerical studiesofPDE’s we may be confrontedwithmatrices G where the computation of spectral properties are not easy and necessitate enormous calculus and sometimes,induces slow algorithms and bed convergence rates.Thus,one motivation here and issued from[6]consists to resolve such problem and prove the invertibility of the algebraic operator yielded in the numerical scheme by appalling topological method Instead of using classical ones such as tri-diagonal transformations.We thus aim to prove that generalized Lyapunov-Sylvester operators can be good candidates for investigating numerical solutions of PDEs in multi-dimensional spaces.

The paper is organized as follows.In next section the discretization of(1.4)is developed,the solvability of the scheme is analyzed.In section 4,the consistency of the method is shown and next,the stability and convergence are proved based on Lyapunov method and Lax-Richtmyer theorem.Finally,a numerical implementation is provided leading to the computation of the numerical solution and error estimates.

2 The numerical scheme

Let J∈N?and h=be the space step.Denote for(j,m)∈I2={0,1,...,J}2,xj=L0+jh and ym=L0+mh.Next,let l=Δt be the time step and tn=t0+nl,n∈N be the time grid.We denote alsothe associated discrete space.Finally,for(j,m)∈I2and n≥0,denotes the net function u(xj,ym,tn)andis the numerical solution.

The following discrete approximations will be applied for the different differential operators involved in the problem.For time derivatives,we set

and for space derivatives,we shall use for the one order

and for the second order ones,we set

Finally,for n∈N?and α∈R,we denote

to designate the estimation ofwith an α-extrapolation/interpolation barycenter method.By applying these discrete approximations,we obtain

where α,β fixed real numbers.In what follows,we set

We thus get

Taking into account the boundary conditions,we obtain the full matrix expression

Next,for Q∈M(J+1)2(R),we denoteLQthe linear operator which associates to X∈M(J+1)2(R)its imageLQ(X)=QX+XQT.Thus,equation(2.1)will be written as

Now,applying similar techniques as previously we obtain

or equivalently,

Now,applying the boundary conditions,we obtain

As a result,we obtain finally the following discrete coupled system.

3 Solvability of the discrete method

In this section we will examine the solvability of the discrete scheme.the main idea is by transforming the system(2.4)-(2.5)to an equality of the form

with an appropriate function φ.We prove precisely that φ can be expressed with general Lyapunov-Sylvesterform.Next using general properties of such operators we prove that φ is an isomorphism.In most studies,even recent ones such as[2]the authors proved that φ or some modified versions may be contractive by inserting translation-dilation parameters leading to fixed point theory.In the present context it seams that such transformation is not possibleas‖φ‖may be greaterthan 1 and thus no contraction may occur.To overcome this problem we come back to differential calculus and topological properties.

Theorem 3.1.The system(2.4)-(2.5)is uniquely solvable whenever U0and U1are known.

The proof of this result is based on the following preliminary lemma.

Lemma 3.1.Let E be a finite dimensional vector space on R(or C),and(φn)nbe a sequence of endomorphisms converging uniformly to an invertible endomorphism φ.Then there exist n0such that the endomorphism φnis invertible for all n≥n0.

Proof.Consider the endomorphism φ on M(J+1)2(R)×M(J+1)2(R)defined by

where Γ(X)=qX?δκLA(X).To prove Theorem 3.1.,we show that φ is a one to one,and so kerφ is reduced to 0.Indeed,

is equivalent to

which is equivalent to

For β/=0,we get

Choosing l=o?hs+4?(which is always possible),the operatorK=I?σαΓLAtends uniformly to I whenever h tends to zero.Indeed,denote

We have

Thus,

Since we have ‖W‖≤4σ|α|[|q|+4δ|κ|],we obtain

For l=o?h4+s?,this implies that

Consequently the operatorKconverge uniformly to the identity whenever h tends towered 0 and??,with s＞0.Thus,using Lemma 3.1 φ is invertible for l,h small enough with.For β=0,we obtain the system

and thus,

For the same assumption on l and h as above the same operatorK(X)=X?σαΓLA(X)tends toward the identity as h tends to 0.

4 Consistency

Theconsistencyoftheproposedmethodwill be provedbyevaluating thelocal truncation error arising from the scheme introduced for the discretization of the system(1.4).

Applying Taylor Taylor’s expansion,and assuming that u and v to be sufficiently differentiable,we get

Similarly,

Hence,

Next,we get also

and

and finally,

Thus,

Similarly,

We have also

Finally,

Now,

We now examine the second equation in(1.4).Applying the same calculus as above,we get

It results from above that the principal part of the first equation in system(1.4)is

The principal part of the local error truncation due to the second equation is

As a result,we get the following lemma.

Lemma 4.1.The numerical method is consistent with an order 2 in space and time.

Proof.It is clear that the two operatorsandtend towards 0 as l and h tend to 0,which ensures the consistency of the method.Furthermore,the method is consistent with an order 2 in time and space.

5 Stability and convergence

The stability of the discrete scheme will be evaluated using Lyapunov criterion which states that a linear systemL(xn+1,xn,xn?1,...)=0 is stable in the sense of Lyapunov if for any bounded initial solution x0,the solution xnremains bounded for all n≥0.In this section we prove precisely the following result.

Lemma 5.1.The solution(Un,Vn)is bounded independently of n whenever the initial solution?U0,V0?is.

Proof.We proceed by recurrence on n.Assume ‖(U0,V0)‖≤η for some positive η.The system(2.4)-(2.5)can be written on the form

The last equation gives

Substituting in the first one,we obtain

Next,recall that

Thus,

and consequently,

Finally,(5.2)yields that

Setting ω=|q|+8δ|κ|,we obtain

We now evaluate ‖Vn+1‖.Applying Γ for the first equation in the system(5.1),we get

By replacing in the second equation of(3.2),we obtain

We get from(5.5)and(5.3),

Now coming back to(1.4)and applying boundary conditions,we get

where,

and

Hence,

Now,the Lyapunov criterion for stability states exactly that?ε＞0,there exists η＞0 such that whenever‖(U0,V0)‖≤η,we have ‖(Un,Vn)‖≤ε,?n≥0.For n=1,and any ε given such that‖(U1,V1)‖≤ε,we seek an η＞0 for which

By direct substitution in(5.4),for n=0,we obtain

From(5.8),we obtain

Observing that

we get

Next choosing l=o?h4+s?small enough,we obtain

Now,for ε＞0,we seek η ＞0 such that

or otherwise,

The discriminant of the last inequality is

For the same assumption on l and h as above,

Consequently,there are two zeros η1＜η2of the inequality above.Furthermore,replacing η with 0 we get a negative quantity,thus η1＜0＜η2.As a result,η2is the good candidate.Now,choosing ‖(U0,V0)‖≤η2,we get immediately ‖U1‖＜ε.Next,already with n=0,we get similarly to the previous case

Choosing l=o?h4+s?small enough as above,and μ=8σ|α|ω+1,we obtain

Next,recall that

we get

Henceforth,

Now,proceeding as for U1,we prove that for all ε＞0,there is an＞0 satisfying‖V1‖＜ε wheneverFinallyanswers the questi?on.?Assume now thatUk,Vkis bounded for k=1,···,n by ε1wheneverU0,V0is bounded by η and let ε＞ 0.We shall prove that it is possible to choose η satisfying

Lemma 5.2.(Banach)Let E,F be two Banach spaces and φ:E→F be linear.If φ is continue and bijective then φ is a homeomorphism.

Consider as above the endomorphism φ on M(J+1)2(R)×M(J+1)2(R)defined by

Consider also

We obtain thus

We have already proved that φ is one to one for l and h small enough.Since φ is linear and continuous,then φ has a continuous inverse function.So,φ is a homeomorphism on M(J+1)2(R)×M(J+1)2(R).Furthermore

As φ?1is continuous,there exists C＞0 such that

By choosing ‖(Uk,Vk)‖≤η,for all k=0,1,...n,we get

Now,it suffices to prove that there exists η＞0 for which

The discernments is

Hence,there are two zeros of the last equality

It is straightforward that 0∈]η1,η2[.Hence η2＞0.Now for f2we get

For

we obtain ‖Vn+1‖ ≤ ε.Finally,choosing η the minimum between η2and η3the criterion is proved.

As the numerical scheme is consistent and stable,it is then convergent.

6 Numerical Implementations

We propose in this section to present some numerical examples to validate the theoretical results developed previously.The error between the exact solutions and the numerical ones via an L2discrete norm will be estimated.The matrix norm used will be

for a matrix X=(Xij)∈MN+2C.Denote unthe net function u(x,y,tn)and Unthe numerical solution.We propose to compute the discrete error

on the grid(xi,yj),0≤i,j≤ J+1 and to validate the convergence rate of the proposed schemes we propose to compute the proportion

We consider the inhomogeneous problem

on the rectangular domain ? =[?1,1]×[?1,1],where

and the exact solution

with

In Table 1,numerical results are provided.We computed for different space and time stepsthe discrete L2-error estimates defined by(6.1).The time interval is[0,1]for a choice t0=0 and T=1.The following results are obtained for different values of the parameters J(and thus h),l(and thus N).The parameters α and β are fixed to α =and β =.We just notice that some variations done on these latter parameters have induced an important variation in the error estimates which explains their effect as they calibrates the position of the approximated solution around the exact one.The parameters q,λ and κ have the role of viscosity-type coefficients and fixed to the values κ =1,λ=?and q=.The following table outputs the error estimates relatively to the discrete L2-norm defined above for different values of the space and time steps.First,we provide estimates when the optimal condition l=o(h4+s),s＞0 is ful filled.s is fixed to 0.01 for the first column of the table.Next,to control more the effect of such assumption which is due to the presence of a second order Laplacian in the original problem,we tested the convergence of the scheme at some orders less than the optimal power fixed to 4.The second column of the table provides the estimates with a slightly sub-critical power 4?s,s＞0 small enough.Here also s is fixed to 0.01,and finally in the third column of the table,we tested the discrete scheme for a strong sub-critical power.

7 Conclusion

This paper investigated the solution of the well-known Kuramoto-Sivashinsky equation in two-dimensional case by applying a two-dimensional finite difference discretization.The original equation is a 4-th order partial differential equation.Thus,in a first step it was recasted into a system of second order partial differential equations by applying a reduction order.Next,the continuous system of simultaneous coupled PDEs has been transformed into an algebraic discrete system involving a generalized Lyapunov-Syslvester type operators.Solvability,consistency,stability and convergence are thenestablished by applying well-known methods such as Lax-Richtmyer equivalence theorem and Lyapunov Stability and by examining the topological properties of the obtained Lyapunov-Sylvester type operators.The method was finally improved by developing a numerical example.

Table 1:Error estimates.