Lyapunov-Sylvester operators for Kuramoto-Sivashinsky Equation
aa r X i v : . [ m a t h . NA ] N ov Generalized Lyapunov-Sylvester operators forKuramoto-Sivashinsky Equation
Abdelhamid BEZIA Anouar Ben MabroukOctober 16, 2018
Abstract
A numerical method is developed leading to algebraic systems based on generalizedLyapunov-Sylvester operators to approximate the solution of two-dimensional Kuramoto-Sivashinsky equation. It consists of an order reduction method and a finite difference dis-cretization which is proved to be uniquely solvable, stable and convergent by using Lyapunovcriterion and manipulating generalized Lyapunov-Sylvester operators. Some numerical im-plementations are provided at the end to validate the theoretical results.
Mathematics Subject Classification (2010) . 65M06, 65M12, 65M22, 15A30, 37B25.
Key words : Kuramoto-Sivashinsky equation, Finite difference method, LyapunovSylvesteroperators.
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 Kuramoto-Sivashinsky equation ∂u∂t = q ∆ u − κ ∆ u + λ |∇ u | , (( x, y ) , t ) ∈ Ω × ( t , + ∞ ) , (1)with initial conditions u ( x, y, t ) = ϕ ( x, y ) ; ( x, y ) ∈ Ω (2)and boundary conditions ∂u∂η ( x, y, t ) = 0 ; (( x, y ) , t ) ∈ ∂ Ω × ( t , + ∞ ) , (3)on a rectangular domain Ω = [ L , L ] × [ L , L ] in R , t ≥ ∂∂t is the time derivative, ∇ is the space gradient operator and ∆ = ∂ ∂x + ∂ ∂y is the Laplaceoperator in R , q , κ , λ are real parameters. ϕ and ψ are twice differentiable real valued functionson Ω.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 have1o solve the system ∂u∂t = ∆ v + λ |∇ u | , ( x, y, t ) ∈ Ω × ( t , + ∞ ) v = qu − κ ∆ u, ( x, y, t ) ∈ Ω × ( t , + ∞ )( u, v ) ( x, y, t ) = ( ϕ, ψ ) ( x, y ) , ( x, y ) ∈ Ω −→∇ ( u, v )( x, y, t ) = 0 , ( x, y, t ) ∈ ∂ Ω × ( t , + ∞ ) (4)The Kuramoto-Sivashinsky equation (KS) is one of the most famous equations in math-physicsfor many decades. It has its origin in the work of Kuramoto since the 70-th decade of the 20-thcentury in his study of reaction-diffusion equation. The equation was then considered by Sivashin-sky in modeling small thermal diffusion instabilities for laminar flames and modeling the referenceflux of a film layer on an inclined plane. Since then the KS equation has experienced a growing de-velopment in theoretical mathematics, numerical as well as physical mechanics, nonlinear physics,hydrodynamics, chemistry, plasma physics, particle distributions advection, surface morphology,...etc. See [1], [11], [13], [16], [18], [20], [23], [25], [28], [29]. In [20], the symmetry problem of themodel was studied based on the theory of Lie algebras. In [1], an anisotropic version of the KSequation has been proposed leading to global resolutions of the equation on rectangular domaines.Sufficient conditions were given for the existence of global solution.From the dimensional point of view, KS equation has been widely studied in one dimension,[8], [10], [14], [21], [22], [32]. However, this equation since its appearance is related to the mod-eling of flame spread which is a two-dimensional problem. In this context, an important resultwas developed in [30] where the authors showed by adapting the method developed in [24] theexistence of a set of bounded solutions on a rectangular domain. This major importance of thetwo-dimensional model was a main motivation behind this work. A model representing a nonlineardynamical system defined 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), where u is the height of theinterface, q is a pre-factor proportional to the coefficient of surface tension expressed by q ∇ u .The quantity κ ∇ u represents the result of the diffusion surface due to the chemical potentialgradient induced curvature. The pre-factor κ represents the surface diffusion term. The quantity λ |∇ u | is due to the existence of overhangs and vacancies during the deposition process. Finally,the quantity ν ∇ u + λ ∇ u | u is referred to as modeling the effect of deposited atoms. cf. [13].In the present work we propose to serve of algebraic operators to approximate the solutionsof the Kuramoto-Sivashinsky(K-S) equation in the two patial and one temporal dimesion with-out adapting classical developments based on separation of variables, radial solutions, tridiagonaloperators, ... etc.This was onefold of the present paper. A second crucial idea is to transform the continuousK-S equation into a generalized Lyapunov-Sylvester equation of the form X i A i X n B i = C n (5)where A i and B i are appropriate matrices depending on the discretization procedure and theproblem parameters. X n represents the numerical solution at time n and C n is usually dependingon the past values X k , k ≤ n − X A i XA Ti = C (6)2he equation is knwon restrictively as Lyapunov one [27]. Generally speaking, the equation X i A i XB i = C (7)is very difficult to be inverted and remains an open problem in algebra. Nevertheless, some workshave been developed and proved that under suitable conditions on the coefficient matrices, onemay get a unique solution, but it’s exact computation remains hard. It necessitates to computeeigenvalues and precisely bounds/estimates of eigenvalues or direct inverses of big matrices whichremains a complicated problem and usually inappropriate.In [19], a native method to solve (7) is investigated based on Kronecker product and equivalentmatrix-vector equation. The Sylvester’s equation is transformed into a linear equivalent one onthe form Gx = c, with a matrix G obtained by tensor products issued from the A ′ j s and the B ′ j s .However, the general case remained already complicated. The authors have been thus restrictedto special cases where the matrices A j and B j are scalar polynomials based on spacial and fixedmatrices 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 the A ′ j s and the B ′ j s, the equation (7) can be writen as X j,k α j,k A j XB k = C, where α jk are complexe numbers. Hence, the tensor matrix G will be written on the form G = φ ( A, B ) , where φ is the 2-variables polynomial φ ( x, y ) = P j,k α j,k x i y k . Thus, G is singular if andonly 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 (7) was pointedout by Roth [31]. It was shown that the solution is unique if and only if the matrices (cid:20) A C B (cid:21) , and (cid:20) A B (cid:21) are similar.However, in numerical stydies of PDE’s we may be confronted with matrices G where thecomputation of spectral properties are not easy and necessitate enormes 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 yieldedin the numerical scheme by applaying topological method Instead of using classical ones such astri-diagonal transformations. We thus aim to prove that generilazed Lyapunov-Sylvester operatorscan 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 (4) is developed, thesolvability of the scheme is analyzed. In section 4, the consistency of the method is shown andnext, the stability and convergence are proved based on Lyapunov method and Lax-Richtmyertheorem. Finally, a numerical implementation is provided leading to the computation of thenumerical solution and error estimates. 3 The numerical scheme
Let J ∈ N ∗ and h = L − L J be the space step. Denote for ( j, m ) ∈ I = { , , ..., J } , x j = L + jh and y m = L + mh . Next, let l = ∆ t be the time step and t n = t + nl, n ∈ N be the time grid.We denote also e Ω = { ( x j , y m , t n ) ; ( x j , y m , t n ) ∈ I × N } the associated discrete space. Finally,for ( j, m ) ∈ I and n ≥ u nj,m denotes the net function u ( x j , y m , t n ) and U nj,m is the numericalsolution.The following discrete approximations will be applied for the different differential operatorsinvolved in the problem. For time derivatives, we set ∂u∂t U n +1 j,m − U n − j,m l and for space derivatives, we shall use for the one order ∂u∂x U nj +1 ,m − U nj − ,m h , and ∂u∂y U nj,m +1 − U nj,m − h and for the second order ones, we set ∂ u∂x U n,αj +1 ,m − U n,αj,m + U n,αj − ,m h , and ∂ u∂y U n,αj,m +1 − U n,αj,m + U n,αj,m − h . Finally, for n ∈ N ∗ and α ∈ R , we denote U n,α = αU n − + (1 − α ) U n + αU n +1 to designate the estimation of U nj,m with an α -extrapolation/interpolation barycenter method. Byapplying these discrete approximations, we obtain U n +1 j,m − U n − j,m = 2 lh h V n,βj − ,m − V n,βj,m + V n,βj +1 ,m + V n,βj,m − − V n,βj,m + V n,βj,m +1 i + 2 lh λ (cid:20) (cid:0) U nj +1 ,m − U nj − ,m (cid:1) + 14 (cid:0) U nj,m +1 − U nj,m − (cid:1) (cid:21) . In what follows, we set F nj,m = 14 h(cid:0) U nj +1 ,m − U nj − ,m (cid:1) + (cid:0) U nj,m +1 − U nj,m − (cid:1) i We thus get U n +1 j,m − U n − j,m = σ [ βV n +1 j − ,m + (1 − β ) V nj − ,m + βV n − j − ,m − βV n +1 j,m − − β ) V nj,m − βV n − j,m + βV n +1 j +1 ,m + (1 − β ) V nj +1 ,m + βV n − j +1 ,m + βV n +1 j,m − + (1 − β ) V nj,m − + βV n − j,m − − βV n +1 j,m − − β ) V nj,m − βV n − j,m + βV n +1 j,m +1 + (1 − β ) V nj,m +1 + βV n − j,m +1 ] + σλF nj,m σ = lh . Otherwise, this can be written as U n +1 j,m − σβ (cid:2) V n +1 j − ,m − V n +1 j,m + V n +1 j +1 ,m + V n +1 j,m − − V n +1 j,m + V n +1 j,m +1 (cid:3) = U n − j,m + σ (1 − β ) (cid:2) V nj − ,m − V nj,m + V nj +1 ,m + V nj,m − − V nj,m + V nj,m +1 (cid:3) + σβ (cid:2) V n − j − ,m − V n − j,m + V n − j +1 ,m + V n − j,m − − V n − j,m + V n − j,m +1 (cid:3) + σλF nj,m . Taking into account the boundary conditions, we obtain the full matrix expression U n +1 − σβ (cid:0) AV n +1 + V n +1 A T (cid:1) = U n − + σ (1 − β ) (cid:0) AV n + V n A T (cid:1) + σβ (cid:0) AV n − + V n − A T (cid:1) + σλF n (8)where A = − · · · · · · − · · ·
00 . . . . . . . . . . . . ...... . . . . . . . . . . . . 00 · · · − · · · · · · − U n = (cid:0) U nj,m (cid:1) ≤ j,m ≤ J and V = (cid:0) V nj,m (cid:1) ≤ j,m ≤ J .Next, for Q ∈ M ( J +1) ( R ) , we denote L Q the linear operator which associates to X ∈ M ( J +1) ( R ) its image L Q ( X ) = QX + XQ T . Thus, equation ( ?? ) will be written as U n +1 − σβ L A (cid:0) V n +1 (cid:1) = U n − + σ (1 − β ) L A ( V n ) + σβ L A (cid:0) V n − (cid:1) + σλF n . (9)Now, applying similar techniques as previously we obtain V n,βj,m = qU n,αj,m − κh (cid:0) U n,αj − ,m − U n,αj,m + U n,αj +1 ,m + U n,αj,m − − U n,αj,m + U n,αj,m +1 (cid:1) . Let δ = h . This results in the following equation βV n +1 j,m + (1 − β ) V nj,m + βV n − j,m = q (cid:2) αU n +1 j,m + (1 − α ) U nj,m + αU n − j,m (cid:3) − δκ [ αU n +1 j − ,m + (1 − α ) U nj − ,m + αU n − j − ,m − αU n +1 j,m − − α ) U nj,m − αU n − j,m + αU n +1 j +1 ,m + (1 − α ) U nj +1 ,m + αU n − j +1 ,m + αU n +1 j,m − + (1 − α ) U nj,m − + αU n − j,m − − αU n +1 j,m − − α ) U nj,m − αU n − j,m + αU n +1 j,m +1 + (1 − α ) U nj,m +1 + αU n − j,m +1 ] , or equivalently, βV n +1 j,m − qαU n +1 j,m + δκα (cid:2) U n +1 j − ,m − U n +1 j,m + U n +1 j +1 ,m + U n +1 j,m − − U n +1 j,m + U n +1 j,m +1 (cid:3) = − (1 − β ) V nj,m − βV n − j,m + q (cid:2) (1 − α ) U nj,m + αU n − j,m (cid:3) − δκ [(1 − α ) (cid:2) U nj − ,m − U nj,m + U nj +1 ,m + U nj,m − − U nj,m + U nj,m +1 (cid:3) + (cid:2) U n − j − ,m − U n − j,m + U n − j +1 ,m + U n − j,m − − U n − j,m + U n − j,m +1 (cid:3) ] . βV n +1 − α (cid:2) qU n +1 − δκ L A (cid:0) U n +1 (cid:1)(cid:3) = (1 − α ) ( qU n − δκ L A ( U n )) + α (cid:0) qU n − − δκ L A (cid:0) U n − (cid:1)(cid:1) − (1 − β ) V n − βV n − . (10)As a result, we obtain finally the following discrete coupled system. U n +1 − σβ L A ( V n +1 ) = U n − + σ (1 − β ) L A ( V n ) + σβ L A ( V n − ) + σλF n βV n +1 − α [ qU n +1 − δκ L A ( U n +1 )] = (1 − α ) ( qU n − δκ L A ( U n ))+ α ( qU n − − δκ L A ( U n − )) − (1 − β ) V n − βV n − (11) In this section we will examine the solvability of the discrete scheme. the main idea is by trans-forming the system (11) to an equality of the forme ( U n +1 , V n +1 ) = φ ( U n , V n , U n − , V n − ) withan appropriate function φ . We prove precisely that φ can be expressed with general Lyapunov-Sylvester form. 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 versionsmay be contractive by inserting translation-dilation parameters leading to fixed point theory. Inthe present context it seams that such transformation is not possible as k φ k may be greater than1 and thus no contraction may occur. To overcome this problem we come back to differentialcalculus and topological properties. Theorem 1
The system (11) is uniquely solvable whenever U and U are known. The proof of this result is based on the following preliminary lemma.
Lemma 2
Let E be a finite dimensional vector space on R (or C ), and ( φ n ) n be a sequence ofendomorphisms converging uniformly to an invertible endomorphism φ . Then there exist n suchthat the endomorphism φ n is invertible for all n ≥ n . Proof.
Consider the endomorphism φ on M ( J +1) ( R ) × M ( J +1) ( R ) defined by φ ( X, Y ) = ( X − σβ L A ( Y ) , βY − α Γ ( X )) (12)where Γ ( X ) = qX − δκ L A ( X ) . To prove Theorem 2, we show that φ is a one to one, and so ker φ is reduced to 0 . Indeed, φ ( X, Y ) = 0 ⇔ ( X − σβ L A ( Y ) , βY − α [ qX − δκ L A ( X )]) = (0 , . Which is equivalent to (cid:26) X = σβ L A ( Y ) βY = α Γ ( X ) (13)for β = 0 we get Y = σα Γ L A ( Y ) (14)6hoosing l = o ( h s +4 ) (which is always possible), the operator K = I − σα Γ L A tends uniformly to I whenever h tends to zero. Indeed, denote W = σα (cid:0) qA − δκA (cid:1) = 2 lh α (cid:18) qA − h κA (cid:19) .We have K ( X ) = X − σα Γ L A ( X ) = Y − W X − XW T + 2 σαδκAXA T . thus, k ( K − I ) X k = (cid:13)(cid:13) W X + XW T + 2 ( σαδκ ) AXA T (cid:13)(cid:13) ≤ k W k k X k + 32 σαδ | κ | k X k Since we have k W k ≤ σ | α | [ | q | + 4 δ | κ | ] , we obtain, k ( K − I ) X k ≤ σα [ | q | + 8 δ | κ | ] k X k For l = o ( h s ) this implies that, k ( K ( X ) − I ( X )) k ≤ α (cid:2) | q | h s + 8 h s | κ | (cid:3) k X k Consequentely the operator K converge uniformly to the identity whenevr h tends towred 0 and l = o ( h s ) , withs s > . Thus, using Lemma 1 φ is invertible for l , h small enough with l = o ( h s ) . For β = 0 we obtain the system (cid:26) U n +1 = U n − + σ L A ( V n ) + σλF n V n = α Γ ( U n +1 ) + (1 − α ) Γ ( U n ) + α Γ ( U n − ) (15)and thus, U n +1 − σα Γ L A (cid:0) U n +1 (cid:1) = σ [(1 − α ) Γ L A ( U n ) + U n − + α Γ L A (cid:0) U n − (cid:1) + λF n ] (16)For the same assumption on l and h as above the same operator K ( X ) = X − σα Γ L A ( X ) tendstoward the identity as h tends to 0. The consistency of the proposed method will be proved by evaluating the local truncation errorarising from the scheme introduced for the discretization of the system (4)Applying Taylor Taylor’s expansion, and assuming that u and v to be sufficiently differentiable,we get U n +1 j,m = u + l ∂u∂t + l ∂ u∂t + l ∂ u∂t + l ∂ u∂t + ... (17)Similarly, U n − j,m = u − l ∂u∂t + l ∂ u∂t − l ∂ u∂t + l ∂ u∂t + ... (18)7ence, U n +1 j,m − U n − j,m l = ∂u∂t + l ∂ u∂t + ... (19)Next, we get also V n − j − ,m = (cid:20) v − l ∂v∂t + l ∂ v∂t − l ∂ v∂t + l ∂ v∂t (cid:21) − h (cid:20) ∂v∂x − l ∂ v∂t∂x + l ∂ v∂t ∂x − l ∂ v∂t ∂x + l ∂ v∂t ∂x (cid:21) + h (cid:20) ∂ v∂x − l ∂ v∂t∂x + l ∂ v∂t ∂x − l ∂ v∂t ∂x + l ∂ v∂t ∂x (cid:21) − h (cid:20) ∂ v∂x − l ∂ v∂t∂x + l ∂ v∂t ∂x − l ∂ v∂t ∂x + l ∂ v∂t ∂x (cid:21) + h (cid:20) ∂ v∂x − l ∂ v∂t∂x + l ∂ v∂t ∂x − l ∂ v∂t ∂x + l ∂ v∂t ∂x (cid:21) + ... (20)and V nj − ,m = v − h ∂v∂x + h ∂ v∂x − h ∂ v∂x + h ∂ v∂x + ... (21)and finally, V n +1 j − ,m = (cid:20) v + l ∂v∂t + l ∂ v∂t + l ∂ v∂t + l ∂ v∂t (cid:21) − h (cid:20) ∂v∂x + l ∂ v∂t∂x + l ∂ v∂t ∂x + l ∂ v∂t ∂x + l ∂ v∂t ∂x (cid:21) + h (cid:20) ∂ v∂x + l ∂ v∂t∂x + l ∂ v∂t ∂x + l ∂ v∂t ∂x + l ∂ v∂t ∂x (cid:21) − h (cid:20) ∂ v∂x + l ∂ v∂t∂x + l ∂ v∂t ∂x + l ∂ v∂t ∂x + l ∂ v∂t ∂x (cid:21) + h (cid:20) ∂ v∂x + l ∂ v∂t∂x + l ∂ v∂t ∂x + l ∂ v∂t ∂x + l ∂ v∂t ∂x (cid:21) + ... (22)Thus, V n,βj − ,m = v + βl ∂ v∂t + β l ∂ v∂t + h (cid:20) − ∂v∂x − βl ∂ v∂t ∂x − β l ∂ v∂t ∂x (cid:21) + h (cid:20) ∂ v∂x + βl ∂ v∂t ∂x + β l ∂ v∂t ∂x (cid:21) + h (cid:20) − ∂ v∂x − βl ∂ v∂t ∂x − β l ∂ v∂t ∂x (cid:21) + h (cid:20) ∂ v∂x + βl ∂ v∂t ∂x + β l ∂ v∂t ∂x (cid:21) + ... (23)8imilarly, V n,βj,m = v + βl ∂ v∂t + β l ∂ v∂t + ... (24)We have also V n,βj +1 ,m = v + βl ∂ v∂t + β l ∂ v∂t + h (cid:20) ∂v∂x + βl ∂ v∂t ∂x + β l ∂ v∂t ∂x (cid:21) + h (cid:20) ∂ v∂x + βl ∂ v∂t ∂x + β l ∂ v∂t ∂x (cid:21) + h (cid:20) ∂ v∂x + βl ∂ v∂t ∂x + β l ∂ v∂t ∂x (cid:21) + h (cid:20) ∂ v∂x + βl ∂ v∂t ∂x + β l ∂ v∂t ∂x (cid:21) + ... (25)Finally, V n,βj,m − − V n,βj,m + V n,βj,m +1 h = (cid:20) ∂ v∂y + βl ∂ v∂t ∂y + β l ∂ v∂t ∂y (cid:21) + h (cid:20) ∂ v∂y + βl ∂ v∂t ∂y + β l ∂ v∂t ∂y (cid:21) + ... (26)Now, U n +1 j,m − U n − j,m l − " V n,βj − ,m − V n,βj,m + V n,βj +1 ,m h + V n,βj,m − − V n,βj,m + V n,βj,m +1 h = ∂u∂t − ∆ v + l ∂ u∂t − βl ∂ ∂t (∆ v ) − h (cid:18) ∂ v∂x + ∂ v∂y (cid:19) + βl ∂ ∂t (cid:18) ∂ v∂x + ∂ v∂y (cid:19) − β l (cid:20) ∂ v∂t ∂x + ∂ v∂t ∂y (cid:21) − β l h (cid:20) ∂ v∂t ∂x + ∂ v∂t ∂y (cid:21) + ... (27)We now examine the second equation in (4). Applying the same calculus as above, we get v = ( qv − κ (∆ u )) − βl ∂ v∂t + qαl ∂ u∂t − καl ∂ (∆ u ) ∂t − κ h (cid:18) ∂ u∂x + ∂ u∂y (cid:19) + o (cid:0) l + h (cid:1) . It results from above that the principal part of the first equation in system 4 is L u,v ( t, x, y ) = βl ∂ v∂t − αql ∂ u∂t − ακl ∂ (∆ v ) ∂t − κ l (cid:18) ∂ v∂x + ∂ v∂y (cid:19) + o (cid:0) l + h (cid:1) . (28)9he principal part of the local error truncation due to the second equation is L u,v ( t, x, y ) = β l ∂ ∂t ( v − qu − κα ∆ u ) − κ h (cid:18) ∂ v∂x + ∂ v∂y (cid:19) + o (cid:0) l + h (cid:1) . (29)As a result, we get the following lemma. Lemma 3
The numerical method is consistent with an order in space and time. Proof.
It is clear that the two operators L u,v and L u,v tend towards 0 as l and h tend to 0, whichensures the consistency of the method. Furthermore, the method is consistent with an order 2 intime and space. The stability of the discrete scheme will be evaluated using Lyapunov criterion which states that alinear system L ( x n +1 , x n , x n − , ... ) = 0 is stable in the sense of Lyapunov if for any bounded initialsolution x , the solution x n remains bounded for all n ≥ . In this section we prove precisely thefollowing result.
Lemma 4
The solution ( U n , V n ) is bounded independently of n whenever the initial solution ( U , V ) is. Proof.
We proceed by recurrence on n. Assume k ( U , V ) k ≤ η for some positive η . The system(11) can be written on the form (cid:26) U n +1 − σβ L A ( V n +1 ) = U n − + σ (1 − β ) L A ( V n ) + σβ L A ( V n − ) + σλF n .βV n +1 = α Γ ( U n +1 ) + (1 − α ) Γ ( U n ) + α Γ ( U n − ) − (1 − β ) V n − βV n − . (30)The last equation gives, β L A (cid:0) V n +1 (cid:1) = α L A Γ (cid:0) U n +1 (cid:1) + (1 − α ) L A Γ ( U n )+ α L A Γ (cid:0) U n − (cid:1) − (1 − β ) L A ( V n ) − β L A (cid:0) V n − (cid:1) . Substituting in the first one, we obtain K (cid:0) U n +1 (cid:1) = σ (1 − α ) L A Γ ( U n ) + σα L A Γ (cid:0) U n − (cid:1) + U n − + σλF n . (31)Next, recall that (cid:12)(cid:12) F nj,m (cid:12)(cid:12) = 14 (cid:12)(cid:12)(cid:12)(cid:0) U nj +1 ,m − U nj − ,m (cid:1) + (cid:0) U nj,m +1 − U nj,m − (cid:1) (cid:12)(cid:12)(cid:12) . Thus, (cid:12)(cid:12) F nj,m (cid:12)(cid:12) ≤ h(cid:0) U nj +1 ,m (cid:1) + (cid:0) U nj − ,m (cid:1) + (cid:0) U nj,m +1 (cid:1) + (cid:0) U nj,m − (cid:1) i and consequently, k F n k ≤ k U n k . (32)10inally, (31) yields that (cid:13)(cid:13) K (cid:0) U n +1 (cid:1)(cid:13)(cid:13) ≤ σ | − α | kL A Γ ( U n ) k + σ | α | (cid:13)(cid:13) L A Γ (cid:0) U n − (cid:1)(cid:13)(cid:13) + (cid:13)(cid:13) U n − (cid:13)(cid:13) + σ | λ | k F n k . Setting ω = | q | + 8 δ | κ | , we obtain (cid:13)(cid:13) K (cid:0) U n +1 (cid:1)(cid:13)(cid:13) ≤ ωσ | − α | k U n k + [1 + 8 ωσ | α | ] (cid:13)(cid:13) U n − (cid:13)(cid:13) + 2 σ | λ | k U n k . (33)We now evaluate k V n +1 k . Applying Γ for the first equation in the system (30), we getΓ (cid:0) U n +1 (cid:1) = σβ Γ L A (cid:0) V n +1 (cid:1) + Γ (cid:0) U n − (cid:1) + σ (1 − β ) Γ ( L A ( V n ))+ σβ Γ (cid:0) L A (cid:0) V n − (cid:1)(cid:1) + σλ Γ ( F n ) . By replacing in the second equation of (13) we obtain β K (cid:0) V n +1 (cid:1) = (1 − β ) [ σα Γ ( L A ( V n )) − V n ]+ β (cid:2) σα Γ (cid:0) L A (cid:0) V n − (cid:1)(cid:1) − V n − (cid:3) +2 α Γ (cid:0) U n − (cid:1) + (1 − α ) Γ ( U n ) + σαλ Γ ( F n ) , (34)We get from (34) and (32), (cid:13)(cid:13) β K (cid:0) V n +1 (cid:1)(cid:13)(cid:13) ≤ | − β | [8 σ | α | ω + 1] k V n k (35)+ | β | (cid:2) [8 σ | α | ω + 1] (cid:13)(cid:13) V n − (cid:13)(cid:13)(cid:3) + 2 | α | ω (cid:13)(cid:13) U n − (cid:13)(cid:13) + (1 − α ) ω k U n k + 2 σαλω k U n k . Now comming back to (4) and appalying boundry conditions, we get U − = U + l e ϕ and V − = qU + e ψ (36)where, e ϕ = − q ∆ ϕ + κ ∆ ϕ − λ |∇ ϕ | and e ψ = − (cid:0) lq + κ (cid:1) ∆ ϕ + 2 lκq ∆ ϕ − lκ ∆ ϕ − λlq |∇ ϕ | + λlκ ∆ (cid:0) |∇ ϕ | (cid:1) . Hence, (cid:13)(cid:13) U − (cid:13)(cid:13) ≤ (cid:13)(cid:13) U (cid:13)(cid:13) + l k e ϕ k and (cid:13)(cid:13) V − (cid:13)(cid:13) ≤ | q | (cid:13)(cid:13) U (cid:13)(cid:13) + (cid:13)(cid:13)(cid:13) e ψ (cid:13)(cid:13)(cid:13) . (37)Now, the Lyapunov criterion for stability states exactly that ∀ ε > , ∃ η > (cid:13)(cid:13)(cid:0) U , V (cid:1)(cid:13)(cid:13) ≤ η ⇒ k ( U n , V n ) k ≤ ε, ∀ n ≥ . For n = 1 , and any ε given such that k ( U , V ) k ≤ ε, we seek an η > k ( U , V ) k < η .By direct substitution in (33), for n = 0 , we obtain (cid:13)(cid:13) K (cid:0) U (cid:1)(cid:13)(cid:13) ≤ ωσ | − α | (cid:13)(cid:13) U (cid:13)(cid:13) + [1 + 8 ωσ | α | ] (cid:13)(cid:13) U − (cid:13)(cid:13) + 2 σ | λ | (cid:13)(cid:13) U (cid:13)(cid:13) . , we obtain (cid:13)(cid:13) K (cid:0) U (cid:1)(cid:13)(cid:13) ≤ σ | λ | (cid:13)(cid:13) U (cid:13)(cid:13) + (1 + 8 ωσ [ | − α | + | α | ]) (cid:13)(cid:13) U (cid:13)(cid:13) + l (1 + 8 σ | α | ω ) k e ϕ k . Observing that, | − α | + | α | ≤ (1 + 3 | α | ) , we get (cid:13)(cid:13) K (cid:0) U (cid:1)(cid:13)(cid:13) ≤ σ | λ | (cid:13)(cid:13) U (cid:13)(cid:13) + 8 ωσ (1 + 3 | α | ) (cid:13)(cid:13) U (cid:13)(cid:13) + l (1 + 8 σ | α | ω ) k e ϕ k . Next choosing l = o ( h s ) small enough, we obtain (cid:13)(cid:13) U (cid:13)(cid:13) ≤ σ | λ | (cid:13)(cid:13) U (cid:13)(cid:13) + 16 ωσ (1 + 3 | α | ) (cid:13)(cid:13) U (cid:13)(cid:13) + 2 l (1 + 8 σ | α | ω ) k e ϕ k . Now, for ε >
0, we seek η > σ | λ | η + 16 ωσ (1 + 3 | α | ) η + 2 l (1 + 8 σ | α | ω ) k e ϕ k < ε (38)or otherwise, 4 σ | λ | η + 16 ωσ (1 + 3 | α | ) η + 2 l (1 + 8 σ | α | ω ) k e ϕ k − ε < . The discriminant of the last inequatlity is∆ ′ = 64 ( ωσ (1 + 3 | α | )) − σ | λ | (2 l (1 + 8 σ | α | ω ) k e ϕ k − ε ) . For the same assumption on l and h as above,∆ ′ ∼
64 ( ωσ (1 + 3 | α | )) + 4 σ | λ | ε > . Consenquentely, there are two zeros η < η of the inequality above. Furthermore, replacing η with 0 we get a negative quantity, thus η < < η . As a result, η is the good candidate. Now,choosing k ( U , V ) k ≤ η , we get immediately k U k < ε. Next, already with n = 0 , we get similarly to the previous case (cid:13)(cid:13) β K (cid:0) V (cid:1)(cid:13)(cid:13) ≤ | − β | [8 σ | α | ω + 1] (cid:13)(cid:13) V (cid:13)(cid:13) + | β | (cid:2) [8 σ | α | ω + 1] (cid:13)(cid:13) V − (cid:13)(cid:13)(cid:3) +2 | α | ω (cid:13)(cid:13) U − (cid:13)(cid:13) + | − α | ω (cid:13)(cid:13) U (cid:13)(cid:13) + 2 σαλω (cid:13)(cid:13) U (cid:13)(cid:13) Choosing l = o ( h s ) small enough as above, and µ = 8 σ | α | ω + 1, we obtain | β | (cid:13)(cid:13) V (cid:13)(cid:13) ≤ (cid:13)(cid:13) β K (cid:0) V (cid:1)(cid:13)(cid:13) ≤ µ | − β | (cid:13)(cid:13) V (cid:13)(cid:13) + µ | β | (cid:2)(cid:13)(cid:13) V − (cid:13)(cid:13)(cid:3) +2 | α | ω (cid:13)(cid:13) U − (cid:13)(cid:13) + | − α | ω (cid:13)(cid:13) U (cid:13)(cid:13) + 2 σαλω (cid:13)(cid:13) U (cid:13)(cid:13) Next, recall that (cid:13)(cid:13) U − (cid:13)(cid:13) ≤ (cid:13)(cid:13) U (cid:13)(cid:13) + k e ϕ k and (cid:13)(cid:13) V − (cid:13)(cid:13) ≤ | q | (cid:13)(cid:13) U (cid:13)(cid:13) + l (cid:13)(cid:13)(cid:13) e ψ (cid:13)(cid:13)(cid:13) , we get | β | (cid:13)(cid:13) V (cid:13)(cid:13) ≤ µ | − β | (cid:13)(cid:13) V (cid:13)(cid:13) + 2 µ | β | (cid:16) | q | (cid:13)(cid:13) U (cid:13)(cid:13) + (cid:13)(cid:13)(cid:13) e ψ (cid:13)(cid:13)(cid:13)(cid:17) +4 | α | ω (cid:0)(cid:13)(cid:13) U (cid:13)(cid:13) + l k e ϕ k (cid:1) + 2 | − α | ω (cid:13)(cid:13) U (cid:13)(cid:13) + 4 σαλω (cid:13)(cid:13) U (cid:13)(cid:13) . | β | (cid:13)(cid:13) V (cid:13)(cid:13) ≤ σαλω (cid:13)(cid:13) U (cid:13)(cid:13) + 2 µ | − β | (cid:13)(cid:13) V (cid:13)(cid:13) +2 ( µ | q | | β | + ω (2 | α | + | − α | )) (cid:13)(cid:13) U (cid:13)(cid:13) +2 µ | β | (cid:13)(cid:13)(cid:13) e ψ (cid:13)(cid:13)(cid:13) + 4 | α | ωl k e ϕ k . Now, proceeding as for U , we prove that for all ε > , there is an η ′ > k V k < ε whenever k ( U , V ) k ≤ η ′ . Finally η = min ( η , η ′ ) answers the question.Assume now that (cid:0) U k , V k (cid:1) is bounded for k = 1 , , ..., n by ε whenever ( U , V ) is bounded by η and let ε > . We shall prove that it is possible to choose η satisfying k ( U n +1 , V n +1 ) k ≤ ε. Lemma 5 (Banach) Let
E, F be two Banach spaces and φ : E → F be linear. If φ is continueand bijective then φ is a homeomorphism.Consider as above the endomorphism φ on M ( J +1) ( R ) × M ( J +1) ( R ) defined by φ ( X, Y ) = ( X − σβ L A ( Y ) , βY − α [ qX − δκ L A ( X )]) . Consider also f ( X, Y, Z, W ) = X + αβ L A ( Y ) + σ (1 − β ) + σβ L A ( Z ) + σλW and f ( X, Y, Z, T ) = α Γ ( X ) − βY + (1 − α ) Γ ( Z ) − (1 − β ) T. We obtain thus φ (cid:0) U n +1 , V n +1 (cid:1) = (cid:0) U n +1 − σβ L A (cid:0) V n +1 (cid:1) , βV n +1 − α Γ (cid:0) U n +1 (cid:1)(cid:1) = (cid:0) f (cid:0) U n − , V n − , V n , F n (cid:1) , f (cid:0) U n − , V n − , U n , V n (cid:1)(cid:1) . We have already proved that φ is one to one for l and h small enough. Since φ is linear andcontinuous, then φ has a continuous inverse function. So, φ is a homeomorphism on M ( J +1) ( R ) × M ( J +1) ( R ) . Furthermore (cid:0) U n +1 , V n +1 (cid:1) = φ − (cid:0) f (cid:0) U n − , V n − , V n , F n (cid:1) , f (cid:0) U n − , V n − , U n , V n (cid:1)(cid:1) . As φ − is continuous, there exists C > such that (cid:13)(cid:13)(cid:0) U n +1 , V n +1 (cid:1)(cid:13)(cid:13) = (cid:13)(cid:13) φ − (cid:0) f (cid:0) U n − , V n − , V n , F n (cid:1) , f (cid:0) U n − , V n − , U n , V n (cid:1)(cid:1)(cid:13)(cid:13) ≤ C (cid:13)(cid:13)(cid:0) f (cid:0) U n − , V n − , V n , F n (cid:1) , f (cid:0) U n − , V n − , U n , V n (cid:1)(cid:1)(cid:13)(cid:13) . By choosing (cid:13)(cid:13)(cid:0) U k , V k (cid:1)(cid:13)(cid:13) ≤ η , for all k = 0 , , ...n , we get (cid:13)(cid:13) f (cid:0) U n − , V n − , V n , F n (cid:1)(cid:13)(cid:13) ≤ σ | λ | η + (1 + 8 σ (1 + 3 | β | )) η. Now, it suffices to prove that there exists η > σ (1 + 3 | β | )) η + 2 σ | λ | η ≤ ε ⇔ σ | λ | η + (1 + 8 σ (1 + 3 | β | )) η − ε ≤ . σ (1 + 3 | β | )) + 8 σ | λ | ε > . Hence, there are two zeros of the last equality η = − (1 + 8 σ (1 + 3 | β | )) − √ ∆4 σ | λ | and η = − (1 + 8 σ (1 + 3 | β | )) + √ ∆4 σ | λ | . It is straightforward that 0 ∈ ] η , η [ . Hence η >
0. Now for f we get (cid:13)(cid:13) f (cid:0) U n − , V n − , U n , V n (cid:1)(cid:13)(cid:13) ≤ [ ω (1 + 3 | α | ) + (1 + 3 | β | )] η. For η ≤ η = εω (1 + 3 | α | ) + (1 + 3 | β | ) , we obtain k V n +1 k ≤ ε. Finally, choosing η the minimum between η and η the criterion is proved. Lemma 6
As the numerical scheme is consistent and stable, it is then convergent.
We propose in this section to present some numerical examples to validate the theoretical resultsdeveloped previously. The error between the exact solutions and the numerical ones via an L discrete norm will be estimated. The matrix norm used will be k X k = N X i,j =1 | X ij | ! / for a matrix X = ( X ij ) ∈ M N +2 C . Denote u n the net function u ( x, y, t n ) and U n the numericalsolution. We propose to compute the discrete error Er = max n k U n − u n k (39)on the grid ( x i , y j ), 0 ≤ i, j ≤ J + 1 and to validate the convergence rate of the proposed schemeswe propose to compute the proportion C = Erl + h . We consider the inhomogeneous problem ∂u∂t = ∆ v + λ |∇ u | + g ( x, y, t ) , ( x, y, t ) ∈ Ω × ( t , + ∞ ) v = qu − κ ∆ u, ( x, y, t ) ∈ Ω × ( t , + ∞ )( u, v ) ( x, y, t ) = ( ϕ, ψ ) ( x, y ) , ( x, y ) ∈ Ω −→∇ ( u, v )( x, y, t ) = 0 , ( x, y, t ) ∈ ∂ Ω × ( t , + ∞ ) (40)on the rectangular domain Ω = [ − , × [ − , g ( x, y, t ) = Ke ( t )[ e ( t )[ C ( x ) S ( x ) C ( y ) + C ( x ) C ( y ) S ( y )] (41) − (cid:2) C ( x ) S ( x ) C ( y ) S ( y ) (cid:3) ] (42)14 N Er
210 640 1 , . −
12 1320 2 , . −
14 2450 6 , . −
16 4183 1 , . −
18 6707 6 , . −
20 10233 2 , . −
22 14997 1 , . −
24 21258 4 , . −
25 25039 3 , . −
30 52015 6 , . − Table 1: Error estimates for l = o ( h . )and the exact solution ( u, v )( x, y, t ) = ( e ( t ) C ( x ) C ( y ) , ... ) , with C ( x ) = cos( πx S ( x ) = sin( πx e ( t ) = e − π t/ and K = 9 π . In the following tables, numerical results are provided. We computed for different space and timesteps the discrete L -error estimates defined by (39). The time interval is [0 ,
1] for a choice t = 0and 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 thatsome variations done on these latter parameters have induced an important variation in the errorestimates which explains their effect as they calibrates the position of the approximated solutionaround the exact one. The parameters q, λ and κ have the role of viscosity-type coefficients andfixed to the values κ = 1 , λ = − π and q = π . The following tables outputs the error estimatesrelatively to the discrete L − norm defined above for different values of the space and time steps.First, we provide estimates when the optimal condition l = o ( h s ) , s > s is fixedto 0 .
01 for the first table. Next, to control more the effect of such assumption which is due tothe presence of a second order Laplacian in the original problem, we tested the convergence ofthe scheme at some orders less than the optimal power fixed to 4 . The second table provides theestimates with a slightly sub-critical power 4 − s, s > s is fixed to 0 . 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 equa-tion is a 4-th order partial differential equation. Thus, in a first step it was recasted into a systemof second order partial differential equations by applying a reduction order. Next, the continuous15
N Er
210 616 1 , . −
12 1273 2 , . −
14 2355 6 , . −
16 4012 2 , . −
18 6419 6 , . −
20 7973 4 , . −
22 14295 1 , . −
24 20228 5 , . −
25 23806 3 , . −
30 49273 7 , . − Table 2: Error estimates for l = o ( h . ) J N Er
210 128 3 , . −
12 220 8 , . −
14 350 2 , . −
16 523 1 , . −
18 746 9 , . −
20 1024 2 , . −
22 1364 1 , . −
24 1772 6 , . −
25 2004 5 , . −
30 3468 1 , . −