\documentclass[reqno]{amsart} \usepackage{hyperref} \AtBeginDocument{{\noindent\small \emph{Electronic Journal of Differential Equations}, Vol. 2016 (2016), No. 268, pp. 1--19.\newline ISSN: 1072-6691. URL: http://ejde.math.txstate.edu or http://ejde.math.unt.edu} \thanks{\copyright 2016 Texas State University.} \vspace{8mm}} \begin{document} \title[\hfilneg EJDE-2016/268\hfil Boussinesq Equation] {Lyapunov-Sylvesters operators for (2+1)-Boussinesq equation} \author[A. Bezia, A. Ben Mabrouk, K. Betina \hfil EJDE-2016/268\hfilneg] {Abdelhamid Bezia, Anouar Ben Mabrouk, Kamel Betina} \address{Abdelhamid Bezia \newline Algebra and Number Theory Laboratory, Faculty of Mathematics, University of Sciences and Technology Houari Boumediene, BP 32 EL-Alia 16111, Bab Ezzouar, Algiers, Algeria} \email{abdelhamid.bezia@gmail.com} \address{Anouar Ben Mabrouk \newline D\'epartement de Math\'ematiques, Institut Sup\'erieur de Math\'amatiques Appliqu\'ees et Informatique de Kairouan, Avenue Assad Ibn Al-Fourat, Kairouan 3100, Tunisia} \email{anouar.benmabrouk@fsm.rnu.tn} \address{Kamel Betina \newline Algebra and Number Theory Laboratory, Faculty of Mathematics, University of Sciences and Technology Houari Boumediene, BP 32 EL-Alia 16111, Bab Ezzouar, Algiers, Algeria} \email{kamelbetina@gmail.com} \thanks{Submitted February 8, 2016. Published October 7, 2016} \subjclass[2010]{65M06, 65M12, 65M22, 15A30, 37B25} \keywords{Boussinesq equation; finite difference method; Numerical solution; \hfill\break\indent Lyapunov-Sylvester operators} \begin{abstract} This article studies a technique for solving a two-dimensional Boussinesq 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. \end{abstract} \maketitle \numberwithin{equation}{section} \newtheorem{theorem}{Theorem}[section] \newtheorem{lemma}[theorem]{Lemma} \newtheorem{remark}[theorem]{Remark} \allowdisplaybreaks \section{Introduction} In this work we use Lyapunov-Sylvester algebraic operators to approximate the solutions of some PDEs such as Boussinesq one in higher dimensions without adapting classical developments based on separation of variables, radial solutions, etc. The crucial idea is to prove that simple methods of discretization of PDEs such as finite difference, finite volumes, can be transformed into well adapted algebraic systems such as Lyapunov-Sylvester ones leading to best algorithms when regarded for convergence rates, time execution and error estimates. In this article, fortunately, we are confronted with more complicated but fascinating forms to prove the invertibility of the algebraic operator appearing in the numerical scheme. Instead of using classical methods such as tri-diagonal transformations we applied a topological method to prove the invertibility. This is good as it did not necessitate to compute eigenvalues and precisely bounds/estimates of eigenvalues or direct inverses which remains a complicated problem in general linear algebra and especially for generalized Lyapunov-Sylvester operators. Recall that even though, bounds/estimates of eigenvalues can already be efficient in studying stability. Recall also that block tri-diagonal systems for classical methods can be already used here also and can be solved for example using iterative techniques, or highly structured bandwidth solvers, Kronecker-product techniques, etc. These methods have been subjects of more general discretization. See \cite{ElMikkawy1,ElMikkawy2,ElMikkawy3,ElMikkawy-Karawia,El-Mikkawy-Atlan,Jia-Li} for a review on tri-diagonal and block tri-diagonal systems, their advantages as well as their disadvantages. This article has many folds. One principal aim is to apply non tri-diagonal type algebraic methods to investigates numerical solutions for PDEs in multi-dimensional spaces. We aim to prove that Lyapunov-Sylvester operators can be good candidates for such aim and that they may give best solvers compared to tri-diagonal and/or block tri-diagonal ones. Recall that the later methods are unadvised because of many reasons. First they are costing methods from both the machine memory and time. In higher dimensions, they secondly need to transform the original problem into an external space of projection and thus solve an associated problem in the new space and next to left to the original one. These facts may induce as previously time and accuracy losing. This article is devoted to the development of a numerical method based on two-dimensional finite difference scheme to approximate the solution of the nonlinear Boussinesq equation in $\mathbb{R}^2$ written in the form \begin{equation}\label{eqn1-1} u_{tt}=\Delta u+qu_{xxxx}+(u^2)_{xx},\quad ((x,y),t)\in \Omega \times (t_{0},+\infty ) \end{equation} with initial conditions \begin{equation}\label{eqn1-3} u(x,y,t_{0})=u_{0}(x,y)\quad \text{and}\quad \frac{\partial u}{\partial t}(x,y,t_{0})=\varphi (x,y),\quad (x,y)\in \Omega \end{equation} and boundary conditions \begin{equation}\label{eqn1-4} \frac{\partial u}{\partial \eta }(x,y,t) =0,\quad ((x,y),t)\in\partial \Omega \times (t_{0},+\infty ). \end{equation} To reduce the derivation order, we set \begin{equation} v=qu_{xx}+u^2. \label{eqn1-2} \end{equation} We have to solve the system \begin{equation} \begin{gathered} u_{tt}=\Delta u+v_{xx},\quad (x,y,t)\in \Omega \times (t_{0},+\infty ) \\ v=qu_{xx}+u^2,\quad (x,y,t)\in \Omega \times (t_{0},+\infty ) \\ (u,v) (x,y,t_{0})=(u_{0},v_{0}) (x,y),\quad (x,y)\in\overline{\Omega } \\ \frac{\partial u}{\partial t}(x,y,t_{0})=\varphi (x,y),\quad (x,y)\in\overline{\Omega } \\ \frac{\partial }{\partial \eta }(u,v)(x,y,t)=0,\quad (x,y,t)\in \partial\Omega \times (t_{0},+\infty ) \end{gathered}\label{systemequation1} \end{equation} on a rectangular domain $\Omega =]L_{0},L_1[\times ]L_{0},L_1[$ in $\mathbb{R}^2$. $t_{0}\geq 0$ is a real parameter fixed as the initial time, $u_{tt}$ is the second order partial derivative in time, $\Delta =\frac{\partial ^2}{\partial x^2}+\frac{\partial ^2}{\partial y^2}$ is the Laplace operator in $\mathbb{R}^2$, $q$ is a real constant, $u_{xx}$ and $u_{xxxx}$ are respectively the second order and the fourth order partial derivative according to $x$. $\frac{\partial }{\partial \eta }$ is the outward normal derivative operator along the boundary $\partial \Omega $. Finally, $u$, $u_{0}$ and $\varphi $ are real valued functions with $u_{0}$ and $\varphi $ are $\mathcal{C}^2$ on $\overline{\Omega }$. $u$ (and consequently $v$) is the unknown candidates supposed to be $\mathcal{C}^{4}$ on $\overline{\Omega}$. Several papers have been devoted to the study of existence and uniqueness of solutions of problem \eqref{eqn1-1} and sometimes exact solutions are developed such as solitary, stationary, time-independent, one-dimensional ones. For example, in the case of a one-direction viscous fluid we may seek solutions of the form $u(x,y,t)=\alpha\psi(x)$. In this case, the problem is transformed into a one variable ordinary differential equation $$ q\psi''(x)+\psi(x)+\alpha\psi^2(x)=ax+b, $$ for some constants $a$ and $b$ depending on the initial-boundary conditions. Therefore, the existence and uniqueness problems are overcame using the well-known theory of ODEs. For more details on these facts, we may refer to \cite{BaoDan1,Bratsos1,Diaz1,Diaz2,Jafarizadeh1,Kano1,Kaya1,Liu1,Parlange1, Song1,Varlamov1,Wazwaz1,Yi1,Yi2}. The Boussinesq equation has a wide reputation in both theoretic and applied fields. It governs the flow of ground water, heat conduction, natural convection in thermodynamics for both volume and fluids in porous media, traveling-waves solutions, self-similar solutions, scattering method, mono and multi dimensional versions, reduction of multi dimensional equations with respect to algebras, etc. In \cite{Dehghan2006}, several finite difference schemes such as three fully implicit finite difference schemes, two fully explicit ones, an alternating direction implicit procedure and the Barakat and Clark type explicit formula are discussed and applied to solve a two-dimensional case. In \cite{Dehghan2008}, the solution of a generalized Boussinesq equation has been developed by means of the homeotypic perturbation method. It consisted in a technique method that avoids the discretization, linearization, or small perturbations of the equation and thus reduces the numerical computations. Next, \cite{Dehghan-Salehi}, a boundary-only meshfree method has been applied to approximate the numerical solution of the classical Boussinesq equation in one dimension. In \cite{Shokri-Dehghan}, a collocation and approximation of the numerical solution of the improved Boussinesq equation is obtained based on radial bases. A predictor-corrector scheme is provided and the Not-a-Knot method is used to improve the accuracy in the boundary. For this reason, many studies have been developed discussing the solvability of such equations. In \cite{Diaz1}, a Boussinesq system of hydrodynamics equations arising in a coupling between Navier–Stokes equations and thermodynamic ones in the the presence of density gradients and where thermodynamical coefficients such as viscosity, specific heat and thermal conductivity are not assumed to be constants and thus leading to a coupled system of quasi-parabolic equations. The authors studied the existence and uniqueness of weak solutions. In this model there are two paradigmatic situations as stated by the authors and related to the fast and the slow heat diffusion. In theoretical mathematical study of such systems, this may correspond to the singular or degenerate character of the heat equation which occur according to the relative behavior of the specific heat of the fluid and its thermal conductivity. By assuming local H\"older approximations the behaviour of the solution is studied near the origin. In \cite{Diaz2}, local strong solutions for a parabolic system based on Boussinesq equation are studied for buoyancy-driven flow with viscous heating. A modification of the classical Navier-Stokes-Boussinesq system motivated by unresolved issues regarding the global solvability of the classical system in situations where viscous heating cannot be neglected is developed. The authors applied a simple model to obtain a coupled system of two parabolic equations where a source term involving the square of the gradient of one of the unknowns appears. Local existence and uniqueness in time of strong solutions for the model problem are established. See for instance \cite{BaoDan1,Benmabrouk1,Bratsos1,Clarkson1,Diaz1,Diaz2,Jafarizadeh1, Kano1,Kaya1,Lai1,Liu1,Parlange1,Song1,Varlamov1,Wazwaz1,Yi1,Yi2} and the references therein for backgrounds on theses facts. For given matrices $A\in \mathbb{R}^{m\times m},B\in \mathbb{R}^{n\times n}$, and $C\in \mathbb{R}^{m\times n}$, the Sylvester equation is given by the form $AX+XB=C$. A brute force attack to obtain the the solution $X$ is to rewrite the Sylvester equation in standard\ $mn\times mn$ linear system $G\widetilde{x}=\widetilde{c}$ using the Kronecker Product,\cite{Jameson}. The Sylvester equation can be solved by Gaussian elimination with $O(m^{3}n^{3}) $ flops. This approach dramatically increases the complexity of the computation, and also cannot preserve the intrinsic properties of the problem in practice \cite{Simoncini}. We denote whatever the special structure of the large linear system $G\widetilde{x}=\widetilde{c}$ can be using only rational operation. In numerical analysis, for solving the Sylvester equation one using the Bratels-Stewart and the Golub-Nash-Van Loan algorithm use $O(m^{3}+n^{3})$ floating point operations, if one assume that an $M\times M$ matrix can be reduced to Schur form with $O(M^{3})$ operations. More precise details are given in \cite{Bartels} and \cite{Goloub}. In \cite{Kirrini} the author describe an algorithm that computes the solution $X$ over an arbitrary field $\mathbb{F}.$ The complexity of the algorithm for $A\in \mathbb{F}^{m\times m},\ \ B\in \mathbb{F}^{n\times n}$ and $m,n\leq N$ is $O(N^{\beta }.logN)$ arithmetic operations in $\mathbb{F}$, where $\beta>2$ is such that $M\times M$ matrices can be multiplied with $O(M^{\beta }) $ arithmetic operations. This algorithm is competitive in terms of arithmetic operation with and even faster than the classical algorithms, and by useful for generalizations for other field than $\mathbb{R}$ or $\mathbb{C}$. The method developed in this paper consists in replacing time and space partial derivatives by finite-difference approximations in order to transform the continuous problem into linear Lyapunov-Sylvester systems. An order reduction method is adapted leading to a system of coupled PDEs which is transformed by the next to a discrete algebraic one. The motivation behind the idea of applying Lyapunov operators was already evoked in our work \cite{Benmabrouk1}. We recall in brief that such a method leads to fast convergent and more accurate discrete algebraic systems without going back to the use of tri-diagonal and/or fringe-tridiagonal matrices already used when dealing with multidimensional problems especially in discrete PDEs. To recapitulate, the method developed here is favorable for many reasons \begin{itemize} \item The first motivation is the fact that it somehow does not change the geometric presentation of the problem as we propose to solve in the same two-dimensional space. We did not project the problem on tri-diagonal representations using the Kronecker product. Relatively to computer architecture, the process of projecting on different spaces and next lifting to the original one may induce degradation of error estimates and slow algorithms. \item The method developed is not just a resolution of a PDE. But, we recall that the resolution itself is not a negligible aim. Further, it proves the efficiency of algebraic operators other than classical tri-diagonal ones. \item We proved here that even when the two systems are equivalent in the sense that they present the same PDE, but with different forms and dimensions, such forms play a major role in the resolution. \item The fact obtaining fast algorithms is very important in computer sciences and makes itself a major aim in computer studies. Recall that the famous method known in mathematical studies of accelerating algorithms in the EM one (expectation-maximisation) which is based on more complicated theories. Here, we proved that we may obtain more rapid algorithms by using just a suitable representation and suitable discerete transformation of the PDE. We got faster algorithms without adding more parameters. \end{itemize} In the organization of this article, the next section is concerned with the introduction of the finite difference scheme. Section 3 is devoted to the discretization of the continuous reduced system obtained from \eqref{eqn1-1}-\eqref{eqn1-4} by the order reduction method. Section 4 deals with the solvability of the discrete Lyapunov equation obtained from the discretization method. In section 5, the consistency of the method is shown and next, the stability and convergence of are proved based on Lyapunov criterion. Finally, some numerical implementation is provided in section 6 leading to the computation of the numerical solution and error estimates. \section{Discrete two-dimensional Boussinesq equation} Consider the domain $\Omega =]L_{0},L_1[\times ]L_{0},L_1[\subset\mathbb{R}^2$ and an integer $J\in \mathbb{N}^{\ast }$. Denote $h=\frac{L_1-L_{0}}{J}$ for the space step, $x_{j}=L_{0}+jh$ and $y_{m}=L_{0}+mh$ for all $(j,m)\in {I}^2=\{0,1,\dots ,J\}^2$. Let $l=\Delta t$ be the time step and $t_n=t_{0}+nl$, $n\in \mathbb{N}$ for the discrete time grid. For $(j,m)\in {I}$ and $n\geq 0$, $u_{j,m}^n$ will be the net function $u(x_{j},y_{m},t_n)$ and $U_{j,m}^n$ the numerical solution. The following discrete approximations will be applied for the different differential operators involved in the problem. For time derivatives, we set as discrete initial condition \[ U^0=U^{-1}+l\varphi \] and for $n\geq1$, \[ u_t\rightsquigarrow \frac{U_{j,m}^{n+1}-U_{j,m}^{n-1}}{2l}\quad \text{and}\quad u_{tt}\rightsquigarrow \frac{U_{j,m}^{n+1}-2U_{j,m}^n+U_{j,m}^{n-1}}{l^2} \] and for space derivatives, we shall use \[ u_{x}\rightsquigarrow(U_x)_{j,m}=\frac{U_{j+1,m}^n-U_{j-1,m}^n}{2h} \quad\text{and}\quad u_{y}\rightsquigarrow(U_y)_{j,m}=\frac{U_{j,m+1}^n-U_{j,m-1}^n}{2h} \] for first order derivatives, and \begin{gather*} u_{xx}\rightsquigarrow(U_{xx})_{j,m} =\frac{U_{j+1,m}^{n,\alpha}-2U_{j,m}^{n,\alpha }+U_{j-1,m}^{n,\alpha }}{h^2},\\ u_{yy}\rightsquigarrow(U_{yy})_{j,m}=\frac{U_{j,m+1}^{n,\alpha}-2U_{j,m}^{n,\alpha } +U_{j,m-1}^{n,\alpha }}{h^2} \end{gather*} for second order ones, where for $n\in \mathbb{N}^{\ast }$ and $\alpha \in\mathbb{R}$, \[ u^{n,\alpha }=\alpha U^{n+1}+(1-2\alpha )U^n+\alpha U^{n-1}. \] Finally, we denote $\sigma =\frac{l^2}{h^2}$ and $\delta=\frac{q}{h^2}$. For $(j,m)\in \mathring{I}^2$ an interior point of the grid $I^2$, ($\mathring{I}=\{1,2,\dots ,J-1\}$), and $n\geq 1$, the following discrete equation is deduced from the first equation in the system \eqref{systemequation1}. \begin{equation} \begin{aligned} & U_{j,m}^{n+1}-2U_{j,m}^n+U_{j,m}^{n-1}\hfill \\ & = \sigma \alpha(U_{j-1,m}^{n+1}-4U_{j,m}^{n+1}+U_{j+1,m}^{n+1} +U_{j,m-1}^{n+1}+U_{j,m+1}^{n+1}) \\ &\quad +\sigma (1-2\alpha )(U_{j-1,m}^n-4U_{j,m}^n+U_{j+1,m}^n +U_{j,m-1}^n+U_{j,m+1}^n) \\ &\quad +\sigma \alpha (U_{j-1,m}^{n-1}-4U_{j,m}^{n-1}+U_{j+1,m}^{n-1} +U_{j,m-1}^{n-1}+U_{j,m+1}^{n-1}) \\ &\quad +\sigma \alpha (V_{j-1,m}^{n+1}-2V_{j,m}^{n+1} +V_{j+1,m}^{n+1}) \\ &\quad +\sigma (1-2\alpha )(V_{j-1,m}^n-2V_{j,m}^n+V_{j+1,m}^n)\\ &\quad +\sigma \alpha (V_{j-1,m}^{n-1}-2V_{j,m}^{n-1}+V_{j+1,m}^{n-1}) . \end{aligned} \label{eqn1-1discrete} \end{equation} Similarly, the following discrete equation is obtained from equation \eqref{eqn1-2}. \begin{equation} \begin{aligned} V_{j,m}^{n+1}+V_{j,m}^{n-1} & = 2\delta \alpha (U_{j-1,m}^{n+1}-2U_{j,m}^{n+1}+U_{j+1,m}^{n+1}) \\ &\quad +2\delta (1-2\alpha )(U_{j-1,m}^n-2U_{j,m}^n+U_{j+1,m}^n)\hfill \\ &\quad +2\delta \alpha (U_{j-1,m}^{n-1}-2U_{j,m}^{n-1}+U_{j+1,m}^{n-1}) +2\widehat{F(U_{j,m}^n)} \end{aligned} \label{eqn2-1discrete} \end{equation} where \[ F(u)=u^2,\;F^n=F(u^n)\quad \text{and}\quad \widehat{F^n}=\frac{F^{n-1}+F^n}{2}. \] The discrete boundary conditions are written for $n\geq 0$ as \begin{gather} \label{CBD1} U_{1,m}^n=U_{-1,m}^n\quad \text{and}\quad U_{J-1,m}^n=U_{J+1,m}^n, \\ \label{CBD2} U_{j,1}^n=U_{j,-1}^n\quad \text{and}\quad U_{j,J-1}^n=U_{j,J+1}^n. \end{gather} The parameter $q$ is related to the equation and has the role of a viscosity-type coefficient and thus it is related to the physical domain of the model. The barycenter parameter $\alpha $ is used to calibrates the position of the approximated solution around the exact one. Of course, these parameters affect surely the numerical solution as well as the error estimates. This fact will be recalled later in the numerical implementations part. In a future work in progress now, we are developing results on numerical solutions of 2D Schr\"{o}dinger equation on the error estimates as a function on the barycenter calibrations by using variable coefficients $\alpha _n$ instead of constant $\alpha $. The use of these calibrations permits the use of implicit/explicit schemes by using suitable values. For example for $\alpha =\frac{1}{2}$, the barycenter estimation \[ V^{n,\alpha }=\alpha V^{n+1}+(1-2\alpha )V^n+\alpha V^{n-1} =\frac{V^{n+1}+V^{n-1}}{2} \] which is an implicit estimation that guarantees an error of order $2$ in time. As motioned in the introduction, the main idea consists in applying Lyapunov-Sylvester operators to approximate the solution of the continuous problem \eqref{eqn1-1}-\eqref{eqn1-4} or its discrete equivalent system \eqref{eqn1-1discrete}-\eqref{CBD2}. Denote \begin{gather*} a_1=\frac{1}{2}+2\alpha \sigma ,\quad a_2=-\alpha \sigma , \\ b_1=1-2(1-2\alpha )\sigma ,\quad \,b_2=(1-2\alpha )\sigma , \\ c_1=(1-2\alpha )\delta \quad \text{and}\quad c_2=\alpha \delta . \end{gather*} Equation \eqref{eqn1-1discrete} becomes \begin{equation} \begin{aligned} &a_2U_{j-1,m}^{n+1}+a_1U_{j,m}^{n+1}+a_2U_{j+1,m}^{n+1} +a_2U_{j,m-1}^{n+1}+a_1U_{j,m}^{n+1}+a_2U_{j,m+1}^{n+1}\\ &\quad +a_2(V_{j-1,m}^{n+1}-2V_{j,m}^{n+1}+V_{j+1,m}^{n+1})\\ & = b_2U_{j-1,m}^n+b_1U_{j,m}^n+b_2U_{j+1,m}^n+b_2U_{j,m-1}^n +b_1U_{j,m}^n+b_2U_{j,m+1}^n \\ &\quad -a_2U_{j-1,m}^{n-1}-a_1U_{j,m}^{n-1}-a_2U_{j+1,m}^{n-1} -a_2U_{j,m-1}^{n-1}-a_1U_{j,m}^{n-1}-a_2U_{j,m+1}^{n-1}\\ &\quad +b_2(V_{j-1,m}^n-2V_{j,m}^n+V_{j+1,m}^n) -a_2(V_{j-1,m}^{n-1}-2V_{j,m}^{n-1}+V_{j+1,m}^{n-1}). \end{aligned} \label{eqn1-1discrete-forme-a} \end{equation} Equation \eqref{eqn2-1discrete} becomes \begin{equation} \begin{aligned} &V_{j,m}^{n+1}-2c_2 (U_{j-1,m}^{n+1}-2U_{j,m}^{n+1}+U_{j+1,m}^{n+1}) \\ & = 2c_1(U_{j-1,m}^n-2U_{j,m}^n+U_{j+1,m}^n) \\ &\quad +2c_2(U_{j-1,m}^{n-1}-2U_{j,m}^{n-1}+U_{j+1,m}^{n-1}) -V_{j,m}^{n-1}+2\widehat{F(U_{j,m}^n)}. \end{aligned} \label{eqn2-1discrete-forme-a} \end{equation} Denote \begin{gather*} A=\begin{pmatrix} a_1 & 2a_2 & 0 & \dots & \dots & 0\hfill \cr a_2 & a_1 & a_2 & \ddots & \ddots & \vdots \hfill \cr0 & \ddots & \ddots & \ddots & \ddots & \vdots \hfill \cr\vdots & \ddots & \ddots & \ddots & \ddots & 0\hfill \cr\vdots & \ddots & \ddots & a_2 & a_1 & a_2\hfill \cr0 & \dots & \dots & 0 & 2a_2 & a_1\hfill \end{pmatrix} ,\quad B= \begin{pmatrix} b_1 & 2b_2 & 0 & \dots & \dots & 0\hfill \cr b_2 & b_1 & b_2 & \ddots & \ddots & \vdots \hfill \cr0 & \ddots & \ddots & \ddots & \ddots & \vdots \hfill \cr\vdots & \ddots & \ddots & \ddots & \ddots & 0\hfill \cr\vdots & \ddots & \ddots & b_2 & b_1 & b_2\hfill \cr0 & \dots & \dots & 0 & 2b_2 & b_1 \end{pmatrix}, \\ R=\begin{pmatrix} -2 & 2 & 0 & \dots & \dots & 0\\ 1 & -2 & 1 & \ddots & \ddots & \vdots \\ 0 & \ddots & \ddots & \ddots & \ddots & \vdots\\ \vdots & \ddots & \ddots & \ddots & \ddots & 0\\ \vdots & \ddots & \ddots & 1 & -2 & 1\\ 0 & \dots & \dots & 0 & 2 & -2\ \end{pmatrix} \end{gather*} The system \eqref{CBD1}-\eqref{eqn2-1discrete-forme-a} can be written on the matrix form \begin{equation} \begin{gathered} \mathcal{L}_A(U^{n+1})+a_2RV^{n+1} =\mathcal{L}_B(U^n)-\mathcal{L}_A(U^{n-1})+R(b_2V^n-a_2V^{n-1}),\\ V^{n+1}-2c_2RU^{n+1} =2R(c_1U^n+c_2U^{n-1})-V^{n-1}+2\widehat{F^n} \end{gathered} \label{matrixform1} \end{equation} for all $n\geq 1$ where \[ U^n=(U_{j,m}^n)_{0\leq j,m\leq J},\quad V^n=(V_{j,m}^n)_{0\leq j,m\leq J}, \quad F^n=(F(U_{j,m}^n))_{0\leq j,m\leq J} \] and for a matrix $Q\in \mathcal{M}_{(J+1)^2}(\mathbb{R})$, $\mathcal{L}_{Q} $ is the Lyapunov operator defined by \[ \mathcal{L}_{Q}(X)=QX+XQ^{T},\quad \forall X\in \mathcal{M}_{(J+1)^2}(\mathbb{R}). \] Remark that $V$ is obtained from the auxiliary function $v$ that is applied to reduce the order of the original PDEs in $u$. This reduction yielded the Lyapunov-Sylvester system \eqref{matrixform1} above. A natural question that can be raised here turns around the ordering of $U$ and $V$. So, we stress the fact that no essential idea is fixed at advance but, this is strongly related to the system obtained. For example, in \eqref{matrixform1} above, it is easy to substitute the second equation into the first to omit the unknown matrix $V^{n+1}$ from the first equation. But in the contrary, it is not easier to do the same for $U^{n+1}$, due to the difficulty to substitute it from $\mathcal{L}_A(U^{n+1})$. It is also not guaranteed that the part $a_2RV^{n+1}$ in the first equation is invertible to substitute $V^{n+1}$. So, it is essentially the final system that shows the ordering in $U$ and $V$. \section{Solvability of the discrete problem} In \cite{Benmabrouk1}, the authors have transformed the Lyapunov operator obtained from the discretization method into a standard linear operator acting on one column vector by juxtaposing the columns of the matrix $X$ horizontally which leads to an equivalent linear operator characterized by a fringe-tridiagonal matrix. We used standard computation to prove the invertibility of such an operator. Here. we do not apply the same computations as in \cite{Benmabrouk1}, but we develop different arguments. The first main result is stated as follows. \begin{theorem}\label{theorem1} System \eqref{matrixform1} is uniquely solvable whenever $U^0$ and $U^1$ are known. \end{theorem} \begin{proof} It reposes on the inverse of Lyapunov operators. Consider the endomorphism $\Phi $ defined on $\mathcal{M}_{(J+1)^2}(\mathbb{R})\times \mathcal{M}_{(J+1)^2}(\mathbb{R})$ by $\Phi(X,Y)=(AX+XA^{T}+a_2RY,\frac{1}{2}Y-c_2RX)$. To prove Theorem \ref{theorem1}, it suffices to show that $ker\Phi $ is reduced to $0$. Indeed, \[ \Phi (X,Y)=0\Longleftrightarrow (AX+XA^{T}+a_2RY,\frac{1}{2} Y-c_2RX)=(0,0) \] or equivalently, \[ Y=2c_2RX\quad \text{and}\quad (A+2a_2c_2R^2)X+XA^{T}=0. \] So, the problem is transformed to the resolution of a Lyapunov type equation of the form \begin{equation} \mathcal{L}_{W,A}(X)=WX+XA^{T}=0 \label{Lyapunovequation} \end{equation} where $W$ is the matrix given by $W=A+2a_2c_2R^2$. Denoting \[ \omega =2a_2c_2,\quad \omega _1=a_1+6\omega ,\quad \overline{\omega }_1=\omega _1+\omega, \quad \omega _2=a_2-4\omega\,. \] The matrix $W$ is explicitly given by \[ W= \begin{pmatrix} \omega _1 & 2\omega _2 & 2\omega & 0 & \dots & \dots & \dots & 0 \\ \omega _2 & \overline{\omega }_1 & \omega _2 & \omega & \ddots & \ddots & \ddots & \vdots \\ \omega & \omega _2 & \omega _1 & \omega _2 & \omega & \ddots & \ddots & \vdots \\ 0 & \ddots & \ddots & \ddots & \ddots & \ddots & \ddots & \vdots \\ \vdots &\ddots & \ddots & \ddots & \ddots & \ddots & \ddots & 0 \\ \vdots &\ddots & \ddots & \omega & \omega _2 & \omega _1 & \omega _2 & \omega \\ \vdots & \ddots & \ddots & \ddots & \omega & \omega _2 & \overline{\omega }_1 & \omega _2 \\ 0 & \dots & \dots & \dots & 0 & 2\omega & 2\omega _2 & \omega _1 \end{pmatrix} \] \end{proof} Next, we use the following preliminary result of differential calculus (See \cite{HenriCartan} for example). \begin{lemma} \label{LemmeInversion} Let $E$ be a finite dimensional ($\mathbb{R}$ or $\mathbb{C}$) vector space and $(\Phi_n)_n$ be a sequence of endomorphisms converging uniformly to an invertible endomorphism $\Phi$. Then, there exists $n_{0}$ such that, for any $n\geq\,n_{0}$, the endomorphism $\Phi_n$ is invertible. \end{lemma} The proof is simple and can be found in any differential calculus references such as \cite{HenriCartan}. We recall it here for the convenience and clearness of the paper. Recall that the set $Isom(E)$ (the set of isomorphisms on $E$) is already open in $L(E)$ (the set of endomorphisms of $E$). Hence, as $\Phi\in\operatorname{Isom}(E)$ there exists a ball $B(\Phi,r)\subset \operatorname{Isom}(E)$. The elements $\Phi_n$ are in this ball for large values of $n$. So these are invertible. Assume now that $l=o(h^{2+s})$, with $s>0$ which is always possible. Then, the coefficients appearing in $A$ and $W$ will satisfy as $h\to 0 $ the following. \[ A_{i,i}=\frac{1}{2}+\varepsilon \,h^{2+2s}\to \frac{1}{2}. \] For $1\leq i\leq J-1$, \[ A_{i,i-1}=A_{i,i+1}=\frac{A_{0,1}}{2} =\frac{A_{J,J-1}}{2}=-\varepsilon h^{2+2s}\to 0. \] For $2\leq i\leq J-2$, \[ W_{i,i}=W_{0,0}=W_{J,J}=\frac{1}{2}+2\alpha \varepsilon \,h^{2+2s}-12\alpha ^2\varepsilon h^{2s}\to \frac{1}{2}. \] Similarly, \[ W_{1,1}=W_{J-1,J-1}=\frac{1}{2}+2\alpha \varepsilon h^{2+2s}-14\alpha ^2\varepsilon h^{2s}\to \frac{1}{2} \] and \[ W_{i,i-1}=W_{i,i+1}=\frac{W_{0,1}}{2} =\frac{W_{J,J-1}}{2}=-\alpha \varepsilon h^{2+2s}+8\alpha ^2\varepsilon h^{2s}\to 0 \] Finally, \[ W_{i,i-2}=W_{i,i+2}=\frac{W_{0,2}}{2}=\frac{ W_{J,J-2}}{2}=-2\alpha ^2\varepsilon \,h^{2s}\to 0. \] Recall that the technique assumption $l=o(h^{2+s})$ is a necessary requirement for the resolution of the present problem and may not be necessary in other PDEs. See for example \cite{Benmabrouk2,Benmabrouk1,Benmabrouk3} for NLS and Heat equations. Next, observing that for all $X$ in the space $\mathcal{M}_{(J+1)^2}(\mathbb{R})\times \mathcal{M}_{(J+1)^2}(\mathbb{R})$, \begin{align*} \| (\mathcal{L}_{W,A}-I)(X)\| & =\| (W-\frac{1}{2}I)X+X(A^{T}-\frac{1}{2}I)\| \\ & \leq [ \| W-\frac{1}{2}I\| +\| A^{T}- \frac{1}{2}I\| ] \| X\| , \end{align*} it results that \begin{equation} \| \mathcal{L}_{W,A}-I\| \leq \| W-\frac{1}{2}I\| +\| A^{T}-\frac{1}{2}I\| \leq \,C(\alpha )h^{2s}. \label{LWAtendsVersId} \end{equation} Consequently, the Lyapunov endomorphism $\mathcal{L}_{W,A}$ converges uniformly to the identity $I$ as $h$ goes towards $0$ and $l=o(h^{2+s})$ with $s>0$. Using Lemma \ref{LemmeInversion}, the operator $\mathcal{L}_{W,A} $ is invertible for $h$ small enough. \begin{remark} \label{rmk3.3} \rm The strict hypothesis $l=o(h^{2+s})$,$ s>0$ is theoretical and used to prove the invertibility (solvability) of the discrete system, but from the numerical point of view, we will see that even if this assumption is not satisfied, the algorithm converges faster the tri-diagonal classical methods,and with good error estimates. \end{remark} \section{Consistency, stability and convergence of the discrete method} The consistency of the proposed method is done by evaluating the local truncation error arising from the discretization of the system \begin{equation} \label{system1} \begin{gathered} u_{tt}-\Delta u-v_{xx}=0,\\ v=qu_{xx}+u^2. \end{gathered} \end{equation} The principal part of the first equation is \begin{equation} \label{consistency1} \begin{aligned} \mathcal{L}_{u,v}^1(t,x,y) & = \frac{l^2}{12}\frac{\partial^4u}{\partial t^4} -\frac{h^2}{12}\Big(\frac{\partial^4u}{\partial x^4}+\frac{ \partial^4u }{\partial y^4}\Big) -\alpha l^2\frac{\partial^2(\Delta u) }{\partial t^2}\\ &\quad - \frac{h^2}{12} \frac{\partial^2v}{\partial x^4} -\alpha l^2 \frac{\partial^4v }{\partial t^2\partial x^2} +O(l^2+h^2).\hfill \end{aligned} \end{equation} The principal part of the local error truncation due to the second part is \begin{equation} \label{consistency2} \mathcal{L}_{u,v}^2(t,x,y) = \frac{l^2}{2}\frac{\partial^2v}{\partial t^2} +\frac{l^4}{24} \frac{\partial^4v}{\partial t^4} -\frac{h^2}{12} \frac{\partial^4u}{\partial x^4} - \alpha l^2 \frac{\partial^4u}{\partial t^2\partial x^2}+O(l^2+h^2). \end{equation} It is clear that the two operators $\mathcal{L}_{u,v}^1$ and $\mathcal{L}_{u,v}^2$ tend toward 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. We now proceed by proving the stability of the method by applying the Lyapunov criterion. A linear system $\mathcal{L}(x_{n+1},x_n,x_{n-1},\dots)=0$ is stable in the sense of Lyapunov if for any bounded initial solution $x_{0}$ the solution $x_n$ remains bounded for all $n\geq0$. Here, we will precisely prove the following result. \begin{lemma} \label{LyapunovStabilityLemma} $\mathcal{P}_n$: The solution $(U^n,V^n) $ is bounded independently of $n$ whenever the initial solution $(U^0,V^0)$ is bounded. \end{lemma} We will proceed by recurrence on $n$. Assume firstly that $\|(U^0,V^0)\| \leq \eta $ for some $\eta $ positive. Using the system \eqref{matrixform1}, we obtain \begin{equation} \begin{gathered} \mathcal{L}_{W,A}(U^{n+1})=\mathcal{L}_{Z,B}(U^n)+b_2RV^n-\mathcal{L} _{W,A}(U^{n-1})-a_2R(F^{n-1}+F^n),\\ V^{n+1}=2c_2RU^{n+1}+2R(c_1U^n+c_2U^{n-1})-V^{n-1}+2\widehat{F^n}. \end{gathered} \label{LyapunovStability2} \end{equation} where $Z=B-2a_2c_1R^2$. Consequently, \begin{equation} \begin{aligned} &\| \mathcal{L}_{W,A}(U^{n+1})\| \\ &\leq \| \mathcal{L}_{Z,B}\|\, \| U^n\| +2|b_2|\,\| V^n\| +\| \mathcal{L}_{W,A}\| \,\| U^{n-1}\| +2|a_2|(\| F^{n-1}\| +\| F^n\| ) \end{aligned} \label{LyapunovStability3} \end{equation} and \begin{equation} \begin{aligned} \| V^{n+1}\| &\leq 4|c_2|\,\| U^{n+1}\| +4(|c_1|\, \|U^n\| +|c_2|\,\| U^{n-1}\| ) \\ &\quad +\| V^{n-1}\| +\| F^{n-1}\| +\| F^n\|. \end{aligned} \label{LyapunovStability4} \end{equation} Next, recall that for $l=o(h^{s+2})$ small enough and $s>0$, we have \begin{gather*} a_1=\frac{1}{2}+2\alpha h^{2s+2}\to \frac{1}{2},\quad a_2=-\alpha \,h^{2s+2}\to 0, \\ b_1=1-2(1-2\alpha )h^{2s+2}\to 1,\quad b_2=(1-2\alpha)h^{2s+2}\to 0, \\ c_1=(1-2\alpha )h^{-2}\to \infty, \quad c_2=\alpha \,h^{2s+2}\to 0, \\ a_2c_1=-\alpha (1-2\alpha )h^{2s}\to 0. \end{gather*} As a consequence, for $h$ small enough, \begin{equation} \| \mathcal{L}_{Z,B}\| \leq 2\| B\| +2|a_2c_1|\| R\|^2\leq 2\max (|b_1|,2|b_2|)+4|a_2c_1| \leq 2+4=6, \label{LtildeBB} \end{equation} and the following lemma is deduced from \eqref{LWAtendsVersId}. \begin{lemma}\label{LWABounded} For $h$ small enough, it holds for all $X\in \mathcal{M}_{(J+1)^2}(\mathbb{R})$ that \[ \frac{1}{2}\|X\|\leq(1-C(\alpha)h^{2s})\|X\| \leq\|\mathcal{L}_{W,A}(X)\| \leq(1+C(\alpha)h^{2s})\|X\| \leq\frac{3}{2}\|X\|. \] \end{lemma} Indeed, recall that equation \eqref{LWAtendsVersId} affirms that $\| \mathcal{L}_{W,A}-I\| \leq C(\alpha )h^{2s}$ for some constant $C(\alpha )>0$. Consequently, for any $X$ we obtain \[ (1-C(\alpha )h^{2s})\| X\| \leq \| \mathcal{L}_{W,A}(X)\| \leq (1+C(\alpha )h^{2s})\| X\| . \] For $h\leq \frac{1}{(2C(\alpha ))^{1/2s}}$, we obtain \[ \frac{1}{2}\leq (1-C(\alpha )h^{2s})<(1+C(\alpha )h^{2s})\leq \frac{3}{2} \] and thus Lemma \ref{LWABounded}. As a result, \eqref{LyapunovStability3} yields \begin{equation} \frac{1}{2}\| U^{n+1}\| \leq 6\| U^n\| +2\| V^n\| +\frac{3}{2}\| U^{n-1}\| +2(\| F^{n-1}\|+\| F^n\| ). \label{LyapunovStability3-1} \end{equation} For $n=0$, this implies \begin{equation} \| U^1\| \leq 12\| U^0\| +4\| V^0\| +3\| U^{-1}\| +4(\| F^{-1}\| +\| F^0\| ). \label{LyapunovStability3-2} \end{equation} Using the discrete initial condition \[ U^0=U^{-1}+l\varphi . \] We identify the function $\varphi $ to the matrix whose coefficients are $\varphi _{j,m}=\varphi (x_{j},y_{m})$. We obtain \begin{equation} \| U^{-1}\| \leq \| U^0\| +l\| \varphi \| . \label{U-1Bounds} \end{equation} Observing that \[ F_{j,m}^{-1}=F(U_{j,m}^{-1})=(U_{j,m}^0-l\varphi _{j,m})^2, \] it follows that \[ |F_{j,m}^{-1}|\leq |U_{j,m}^0|^2+2l|\varphi _{j,m}|.|U_{j,m}^0|+l^2|\varphi _{j,m}|^2 \] and consequently, \begin{equation} \| F^{-1}\| \leq \| U^0\| ^2+2l\| \varphi \| \,\| U^0\| +l^2\| \varphi \| ^2. \label{F-1Bounds} \end{equation} Hence, equation \eqref{LyapunovStability3-2} yields \begin{equation} \| U^1\| \leq (15+8l\| \varphi \| )\| U^0\| +4\| V^0\| +8\| F^0\| +3l\| \varphi \| +4l^2\| \varphi \| ^2. \label{LyapunovStability3-3} \end{equation} Now, the Lyapunov criterion for stability states exactly that for each $\varepsilon >0$ thee exists $\eta >0$ such that \begin{equation} \|(U^0,V^0)\| \leq \eta \;\Rightarrow \; \| (U^n,V^n)\|\leq \varepsilon ,\quad \forall n\geq 0. \label{LyapunovStability1} \end{equation} For $n=1$ and $\| (U^1,V^1)\| \leq \varepsilon $, we seek an $\eta>0$ for which $\| (U^0,V^0)\| \leq \eta $. Indeed, using \eqref{LyapunovStability3-3}, this means that, it suffices to find $\eta $ such that \begin{equation} 8\eta ^2+(19+8l\| \varphi \| )\eta +3l\| \varphi \| +4l^2\| \varphi \| ^2-\varepsilon <0. \label{LyapunovStability3-4} \end{equation} The discriminant of this second order inequality is \begin{equation} \Delta (l,h)=(19+8l\| \varphi \| )^2-32(3l\| \varphi \| +4l^2\| \varphi \| ^2-\varepsilon ). \label{Delta1} \end{equation} For $h,l$ small enough, this is estimated as \[ \Delta (l,h)\sim 361+32\varepsilon >0. \] Thus there are two zeros of the second order equality above \[ \eta _1=\frac{\sqrt{\Delta (l,h)}-(19+8l\| \varphi \| ) }{16}>0 \] and a second zero $\eta _2<0$ rejected. Consequently, choosing $\eta \in]0,\eta _1[$ we obtain \eqref{LyapunovStability3-4}. Finally, \eqref{LyapunovStability3-3} yields $\| U^1\| \leq \varepsilon $. Next, equation \eqref{LyapunovStability4}, for $n=0$, implies that \begin{equation} \| V^1\| \leq \,A(l,h,\varphi )\| U^0\| ^2+B(l,h,\varphi )\| U^0\| +C(l,h,\varphi )+16|c_2|\| V^0\| , \label{LyapunovStability4-1} \end{equation} where \begin{gather*} A(l,h,\varphi )=3+32|c_2|, \\ B(l,h,\varphi )=4\Big(|c_1|+8|c_2|(2+l\| \varphi \| )+l\| \varphi \| +\frac{1}{h^2}\Big) , \\ C(l,h,\varphi )=2(1+8|c_2|)l^2\| \varphi \| ^2+4l(4|c_2|+ \frac{1}{h^2})\| \varphi \| . \end{gather*} Choosing $\| (U^0,V^0)\| \leq \eta $, it suffices to study the inequality \begin{equation} A(l,h,\varphi )\eta ^2+(B(l,h,\varphi )+16|c_2|) \eta +C(l,h,\varphi )-\varepsilon \leq 0. \label{V1Bound} \end{equation} Its discriminant satisfies for $h,l$ small enough, \begin{equation} \Delta (l,h)\sim \frac{16}{h^{4}}(1+20\alpha +|1-2\alpha|) ^2 +\frac{128\alpha |q|}{h^2}\varepsilon >0. \label{Delta2} \end{equation} Here also there are two zeros, $\eta _1'=\frac{\sqrt{\Delta (l,h)}-(B(l,h,\varphi )+16|c_2|)}{2A(l,h,\varphi )}>0$ and a second one $\eta '<0$ and thus rejected. As a consequence, for $\eta \in]0,\eta _1'[$ we obtain $\| V^1\| \leq \varepsilon $. Finally, for $\eta \in ]0,\eta _{0}[$ with $\eta _{0}=\min (\eta _1,\eta_1')$, we obtain $\| (U^1,V^1)\| \leq \varepsilon $ whenever $\| (U^0,V^0)\| \leq \eta $. Assume now that the $ (U^{k},V^{k})$ is bounded for $k=1,2,\dots ,n$ (by $\varepsilon _1$) whenever $(U^0,V^0)$ is bounded by $\eta $ and let $\varepsilon >0$. We shall prove that it is possible to choose $\eta $ satisfying $\| (U^{n+1},V^{n+1})\| \leq \varepsilon $. Indeed, from (\eqref {LyapunovStability3-1}, we have \begin{equation} \| U^{n+1}\| \leq 19\varepsilon _1+8\varepsilon _1^2. \label{LyapunovStability3Ordren-1} \end{equation} So, one seeks, $\varepsilon _1$ for which $8\varepsilon_1^2+19\varepsilon _1-\varepsilon \leq 0$. Its discriminant $\Delta =361+32\varepsilon $, with one positive zero $\varepsilon _1= \frac{\sqrt{361+32\varepsilon }-19}{16}$. Then $\| U^{n+1}\| \leq \varepsilon $ whenever $\| (U^{k},V^{k})\| \leq \varepsilon _1$, $k=1,2,\dots ,n$. Next, using \eqref{LyapunovStability4} and \eqref{LyapunovStability3Ordren-1}, we have \begin{equation} \| V^{n+1}\| \leq (4|c_1|+80|c_2|+1) \varepsilon_1+(32|c_2|+2) \varepsilon _1^2. \label{LyapunovStability3Ordren-2} \end{equation} So, it suffices as previously to choose $\varepsilon _1$ such that \[ (32|c_2|+2) \varepsilon _1^2+(4|c_1|+80|c_2|+1) \varepsilon _1-\varepsilon \leq 0. \] $\Delta =(4|c_1|+80|c_2|+1)^2+4(32|c_2|+2)\varepsilon $, with positive zero \[ \varepsilon _1'=\frac{\sqrt{\Delta } -(4|c_1|+80|c_2|+1)}{2(32|c_2|+2)}. \] Then $\| V^{n+1}\| \leq \varepsilon $ whenever $\|(U^{k},V^{k})\| \leq \varepsilon _1'$, $k=1,2,\dots ,n$. Next, it holds from the recurrence hypothesis for $\varepsilon _{0}=\min(\varepsilon _1,\varepsilon _1')$, that there exists $\eta >0$ for which $\| (U^0,V^0)\| \leq \eta $ implies that $\|(U^{k},V^{k})\| \leq \varepsilon _{0}$, for $k=1,2,\dots ,n$, which by the next induces that $\| (U^{n+1},V^{n+1})\| \leq \varepsilon $. \begin{lemma} \label{laxequivresult} As the numerical scheme is consistent and stable, it is then convergent. \end{lemma} This lemma is a consequence of the well known Lax-Richtmyer equivalence theorem, which states that for consistent numerical approximations, stability and convergence are equivalent. Recall here that we have already proved in \eqref{consistency1} and \eqref{consistency2} that the used scheme is consistent. Next, Lemma \ref{LyapunovStabilityLemma}, Lemma \ref{LWABounded} and equation \eqref{LyapunovStability1} yields the stability of the scheme. Consequently, the Lax equivalence Theorem guarantees the convergence. So as Lemma \ref{laxequivresult}. \section{Numerical implementation} 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 $L_2$ discrete norm will be estimated. The matrix norm used will be \[ \|X\|_2=\Big(\sum_{i,j=1}^{N}|X_{ij}|^2\Big)^{1/2} \] for a matrix $X=(X_{ij})\in \mathcal{M}_{N+2}\mathbb{C}$. Denote $u^n$ the net function $u(x,y,t^n)$ and $U^n$ the numerical solution. We propose to compute the discrete error \begin{equation} \mathrm{Er}=\max_n\|U^n-u^n\|_2 \label{Er} \end{equation} on the grid $(x_{i},y_{j})$, $0\leq\,i,j\leq J+1$ and the relative error between the exact solution and the numerical one as \begin{equation} \mathrm{Relative\,Er}=\max_n\frac{\|U^n-u^n\|_2}{ \|u^n\|_2} \label{Errelative} \end{equation} on the same grid. \subsection{A polynomial-exponential example} We develop in this part a classical example based on polynomial function with an exponential envelop. We consider the inhomogeneous problem \begin{equation} \begin{gathered} u_{tt}=\Delta u+v_{xx}+g(x,y,t),\quad (x,y,t)\in\Omega\times(t_{0},T), \\ v=qu_{xx}+u^2,\quad (x,y,t)\in\Omega\times(t_{0},T), \\ (u,v)(x,y,t_{0})=(u_{0},v_{0})(x,y),\quad(x,y,t)\in \overline{\Omega}\times(t_{0},T), \\ \frac{\partial u}{\partial t}(x,y,t_{0})=\varphi(x,y),\quad(x,y)\in \overline{\Omega}, \\ \overrightarrow{\nabla}(u,v)=0,\quad(x,y,t)\in\partial\Omega\times(t_{0},T) \end{gathered} \end{equation} where $\Omega=[-1,1]^2$ and where the right hand term is \begin{align*} g(x,y,t) &=[ (x^2-1) ^2(x^{4}-58x^2+9)-48(35x^{4}-30x^2+3)] +[ y^{4}-14y^2+5] e^{-t}\\ &\quad -16(x^2-1)^2\big[ (x^2-1) ^{4}(15x^2-1) +(y^2-1) ^2(7x^2-1) \big] e^{-2t} \end{align*} The exact solution is \begin{equation} u(x,y,t)=[(x^2-1)^{4}+(y^2-1)^2]e^{-t}. \end{equation} In the following tables, numerical results are provided. We computed for different space and time steps the discrete $L_2$-error estimates defined by \eqref{Er}. The time interval is $[0,1]$ for a choice $t_{0}=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 $q$ and $\alpha $ are fixed to $q=0.01$ and $\alpha =0.25$. We just notice that some variations done on these latter parameters have induced an important variation in the error estimates which explains the effect of the parameter $q$ which has the role of a viscosity-type coefficient and the barycenter parameter $\alpha $ which calibrates the position of the approximated solution around the exact one. Finally, some comparison with our work in \cite{Benmabrouk1} has proved that Lyapunov type operators already result in fast convergent algorithms with a maximum time of execution of 2.014 sd for the present one. The classical tri-diagonal algorithms associated to the same problem with the same discrete scheme and the same parameters yielded a maximum time of 552.012 sd, so a performance of $23.10^{-4}$ faster algorithm for the present one. We recall that the tests are done on a Pentium Dual Core CPU 2.10 GHz processor and 250 Mo RAM. \begin{table}[th] \caption{} \label{table1} \begin{center} \begin{tabular}{|c|c|c|c|} \hline J & $l$ & Er & Relative Er \\ \hline 10 & 1/100 & $4,0.10^{-3}$ & 0,1317 \\ \hline 16 & 1/120 & $3,3.10^{-3}$ & 0,1323 \\ \hline 20 & 1/200 & $2,0.10^{-3}$ & 0,1335 \\ \hline 24 & 1/220 & $1,8.10^{-3}$ & 0,1337 \\ \hline 30 & 1/280 & $1,4.10^{-3}$ & 0,1340 \\ \hline 40 & 1/400 & $9,8.10^{-4}$ & 0.1344 \\ \hline 50 & 1/500 & $7,8.10^{-4}$ & 0,1346 \\ \hline \end{tabular} \end{center} \end{table} \subsection{A 2-particle interaction example} The example developed hereafter is a model of interaction of two particles or two waves. We consider the inhomogeneous problem \begin{equation} \begin{gathered} u_{tt}=\Delta u+v_{xx}+g(x,y,t),\quad (x,y,t)\in\Omega\times(t_{0},T), \\ v=qu_{xx}+u^2,\quad (x,y,t)\in\Omega\times(t_{0},T), \\ (u,v)(x,y,t_{0})=(u_{0},v_{0})(x,y),\quad (x,y,t)\in \overline{\Omega}\times(t_{0},T), \\ \frac{\partial u}{\partial t}(x,y,t_{0})=\varphi(x,y),\quad (x,y)\in \overline{\Omega}, \\ \overrightarrow{\nabla}(u,v)=0,\quad (x,y,t)\in\partial\Omega\times(t_{0},T) \end{gathered} \end{equation} where \[ g(x,y,t)=(4-6\psi^2(y))u^2-\psi^2(x)u. \] and $u$ is the exact solution given by \[ u(x,y,t)=2\psi^2(x)\psi^2(y)\theta(t) \] with \begin{gather*} \psi(x)=\cos(\frac{x}{2}),\quad\theta(t)=e^{-it}, \varphi(x,y)=-2i\psi^2(x)\psi^2(y) \end{gather*} As for the previous example, the following tables shows the numerical computations for different space and time steps the discrete $L_2$-error estimates defined by \eqref{Er} and the relative error \eqref{Errelative}. The time interval is $[-2\pi,+2\pi]$ for a choice $t_{0}=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 $q$ and $\alpha $ are fixed here-also the same as previously, $q=0.01$ and $\alpha =0.25$. Compared to the tri-diagonal scheme the present one leads a faster convergent algorithms \begin{table}[th] \caption{} \label{table2} \begin{center} \begin{tabular}{|c|c|c|c|} \hline J & $l$ & Er & Relative Er \\ \hline 10 & 1/100 & $4,6.10^{-3}$ & 0,2311 \\ \hline 16 & 1/120 & $4,4.10^{-3}$ & 0,2372 \\ \hline 20 & 1/200 & $2,4.10^{-3}$ & 0,2506 \\ \hline 24 & 1/220 & $2,3.10^{-3}$ & 0,2671 \\ \hline 30 & 1/280 & $2,0.10^{-3}$ & 0,3074 \\ \hline 40 & 1/400 & $1,4.10^{-3}$ & 0,3592 \\ \hline 50 & 1/500 & $7,6.10^{-4}$ & 0,2355 \\ \hline \end{tabular} \end{center} \end{table} \begin{remark} \label{rmk5.1} \rm For the convenience of the paper, we give here some computations of the determinants $\Delta (l,h)$ for different values of the parameters of the discrete scheme Firstly, for both examples above, we can easily see that $\| \varphi \| =2$ and thus, equation \eqref{Delta1} yields that \[ \Delta (l,h)=361+32\varepsilon +416l-256l^2. \] For the different values of $l$ as in the tables $1$ and $2$, we obtain a positive discriminant leading two zeros with a rejected one. For the discriminant of equation \eqref{Delta2} we obtain \[ \Delta (l,h)=\frac{676}{h^{4}}+\frac{8\varepsilon}{h^2}. \] Hence, the results explained previously hold. \end{remark} \section{Conclusion} This paper investigated the solution of the well-known Boussinesq equation in two-dimensional case by applying a two-dimensional finite difference discretization. The Boussinesq equation in its original form 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 using a reduction order idea. Next, the system has been transformed into an algebraic discrete system involving Lyapunov-Syslvester matrix terms by using a full time-space discretization. Solvability, consistency, stability and convergence are then established by applying well-known methods such as Lax-Richtmyer equivalence theorem and Lyapunov Stability and by examining the Lyapunov-Sylvester operators. The method was finally improved by developing numerical examples. It was shown to be efficient by means of error estimates as well as time execution algorithms compared to classical ones. \section{Appendix} \subsection{Main steps of the algorithm applied} \begin{itemize} \item Compute the matrices of the system \item Initialization: Compute the matrices $U^0$, $U^1$, $V^0$ and $V^1$ \item for $n\geq 2$, \[ U^n=\mathrm{lyap}(W,A,\mathcal{L}_{Z,B}(U^{n-1})+b_2RV^{n-1}-\mathcal{L} _{W,A}(U^{n-2})-a_2R(F^{n-2}+F^{n-1})), \] and \[ V^n=2c_2RU^n+2R(c_1U^{n-1}+c_2U^{n-2})-V^{n-2}+2\widehat{F^{n-1}}. \] \end{itemize} \subsection{The tridiagonal associated system} Consider the lexicographic mesh $k=j(J+1)+m$ for $0\leq j,m\leq J$, and denote $N=J(J+2)$, and \[ \Lambda _{N}=\{nJ+n-1: n\in \mathbb{N}\},\quad \widetilde{\Lambda }_{N}=\{n(J+1): n\in \mathbb{N}\},\quad \Theta _{N}=\Lambda_{N}\cup \widetilde{\Lambda }_{N}. \] Using the Kroncker product we obtain a tri-diagonal block system on the form \begin{equation} \begin{gathered} \widetilde{A}U^{n+1}+a_2\widetilde{R}V^{n+1} =\widetilde{B}U^n- \widetilde{A}U^{n-1} +b_2\widetilde{R}V^n-a_2\widetilde{R}V^{n-1}\\ V^{n+1}-2c_2\widetilde{R}U^{n+1}=2c_1\widetilde{R} U^n+2c_2\widetilde{R}U^{n-1}-V^{n-1}+2\widehat{F^n}. \end{gathered} \label{tridiagonalsystem1} \end{equation} The numerical solutions' matrices $U^n$ and $V^n$ are identified here as one-column $(N+1)$-vectors and the matrices $\widetilde{A}$, $\widetilde{B}$ and $\widetilde{R}$ are evaluated as follows. The matrix $\widetilde{A}$ \begin{gather*} \widetilde{A}_{j,j}=2a_1 \quad \forall j, 0\leq j\leq N, \\ \widetilde{A}_{j,j+1}=\frac{1}{2}\widetilde{A}_{0,1}=a_2\quad \forall j, 1\leq j\leq N,\; j\notin \Theta_{N},\text{ and $0$ on }\Lambda _{N}, \\ \widetilde{A}_{j-1,j}=\frac{1}{2}\widetilde{A}_{N,N-1}=a_2\quad \forall j, \; 1\leq j\leq N, \; j\notin \Theta _{N}, \text{ and $0$ on }\widetilde{\Lambda }_{N}, \\ \widetilde{A}_{j,j+J+1}=2a_2,\quad \forall j,\; 0\leq j\leq J, \\ \widetilde{A}_{j,j+J+1}=a_2,\quad \forall j,\;J+1\leq j\leq N-J-1 , \\ \widetilde{A}_{j-J-1,j}=a_2,\quad \forall j,\; J+1\leq j\leq N-J-1, \\ \widetilde{A}_{j-J-1,j}=2a_2,\quad \forall j,\; N-J\leq j\leq N. \end{gather*} The matrix $\widetilde{B}$ \begin{gather*} \widetilde{B}_{j,j}=2b_1\quad \forall j,\; 0\leq j\leq N, \\ \widetilde{B}_{j,j+1}=\frac{1}{2}\widetilde{B}_{0,1}=b_2,\quad \forall j, \; 1\leq j\leq N,;j\notin \Theta _{N}, \text{ and $0$ on } \Lambda _{N}, \\ \widetilde{B}_{j-1,j}=\frac{1}{2}\widetilde{B}_{N,N-1}=b_2,\quad \forall j,\; 1\leq j\leq N,\;j\notin \Theta_{N}, \text{ and $0$ on }\widetilde{\Lambda }_{N}, \\ \widetilde{B}_{j,j+J+1}=2b_2\quad \forall j,\;0\leq j\leq J, \\ \widetilde{B}_{j,j+J+1}=b_2\quad \forall j,\; J+1\leq j\leq N-J-1, \\ \widetilde{B}_{j-J-1,j}=b_2\quad \forall j,\;J+1\leq j\leq N-J-1, \\ \widetilde{B}_{j-J-1,j}=2b_2\quad \forall j,\; N-J\leq j\leq N. \end{gather*} The matrix $\widetilde{R}$ \begin{gather*} \widetilde{R}_{j,j}=-2\quad \forall j,\; 0\leq j\leq N, \\ \widetilde{R}_{j,j+J+1}=2\quad \forall j,\; 0\leq j\leq J, \\ \widetilde{R}_{j,j-J-1}=2\quad \forall j,\; N-J\leq j\leq N, \\ \widetilde{R}_{j,j+J+1}=\widetilde{R}_{j-J-1,j}=1\quad \forall j,\;J+1\leq j\leq N-J-1. \end{gather*} System \eqref{tridiagonalsystem1} can be written as a linear standard form \begin{equation} \begin{gathered} \widetilde{W}U^{n+1}=\widetilde{Z}U^n-\widetilde{W}U^{n-1} +b_2\widetilde{R}V^n-2a_2\widetilde{R}\widehat{F^n} \\ V^{n+1}=2\widetilde{R}(c_1I+c_2\widetilde{Z}) U^n+2(c_2-c_1) \widetilde{R}U^{n-1}+2b_2c_2\widetilde{R} ^2V^n-V^{n-1}+2\widehat{F^n}. \end{gathered} \label{tridiagonalsystem2} \end{equation} where $\widetilde{W}$ and $\widetilde{Z}$ are given by $\widetilde{W}=\widetilde{A}+2c_2a_2\widetilde{R}^2$ and $\widetilde{Z}=\widetilde{B}-2c_1a_2\widetilde{R}^2$. \subsection{Some facts on the convergence of solutions and associated spaces} Usually the problem of convergence depends on different quantities in the model and on the geometry of the domain. Denote $$ \Omega_h=\{(x_j,y_m)\in\mathbb{R}^2: 0\leq j,m\leq J\},\quad \Omega_t=\{t_n:n\in\mathbb{N}\} $$ and define the space of grid functions on $\Omega_h$ as $$ \mathcal{V}_h=\{U=(U_{j,m})_{j,m\in\mathbb{Z}} \text{ satisfying \eqref{CBD1}-\eqref{CBD2}}\}. $$ On the grid functions space, we usually define some appropriate norms to compute the error estimate between exact solutions of the continuous inhomogenous problem associated to \eqref{eqn1-1} and its discrete variant obtained through the discrete scheme. For $U\in \mathcal{V}_h$ and $V\in \mathcal{V}_h$ define the inner product $$ (U,V)_h=h^2\sum_{j,m=0}^{J}U_{j,m}.V_{j,m}. $$ This leads to Sobolev norms (or semi-norms) such as \begin{gather*} \|V\|_h=(V,V)^{1/2}_h,\quad \|V\|_{\infty,h}=\max_{0\leq j,m\leq J}|V_{j,m}|, \\ |V|_{1,h}=\Big[h^2\sum_{j,m=0}^{J}(|(U_x)_{j,m}|^2+|(U_y)_{j,m}|^2)\Big]^{1/2},\\ \|V\|_{2,h}=\Big[h^2\sum_{j,m=0}^{J}(|\Delta_{h}U_{j,m}|^2\Big]^{1/2}. \end{gather*} Next, as it appears in the continuous problem derivatives of order 4 of the unknown function $u$, we generally restrict on suitable regularity spaces. It is not sometimes necessary to go to higher derivatives. In the present case for example, we may consider functions that are of class $C^4$ with respect to $x$, $C^2$ with respect to $y$ and class $C^2$ with respect to $t$. We get using summation by parts $$ (\Delta_hV,U)_h= (V,\Delta_hU)_h ,\quad -(\Delta_hV,V)_h=|V|_{1,h}^2, \quad (\Delta_h^2V,V)_h=|V|_{2,h}^2. $$ As we work on a finite grid and thus a finite space of grid functions all these norms (semi-norms) are equivalent, and thus there is no essential difference between them. The norms $\|.\|_h$ and $\|\cdot\|_{\infty,h}$ reflects the $L^2$ convergence, while the semi-norms $|\cdot|_{1,h}$ and $|\cdot|_{1,h}$ reflects somehow the convergence of the discrete derivatives and thus the convergence in the discrete Sobolev space. For more details and backgrounds on these facts we refer to \cite{Benmabrouk1,Benmabrouk2,Benmabrouk3,Bratsos1,Dehghan2006,Dehghan2008, Dehghan-Salehi,Diaz1,Diaz2,Goncalves,Jafarizadeh1,Kano1,Song1,Wazwaz1}. \begin{thebibliography}{99} \bibitem{BaoDan1} T. Bao-Dan, Q. Yan-Hong, Chen-Ning; \emph{Exact Solutions for a Class of Boussinesq Equation}, Applied Mathematical Sciences, 3(6) (2009), 257-265. \bibitem{Bartels} R. H. Bartels, G. W. Stewart; \emph{Algorithm 432: solution of the matrix equation $AX+XB=C$}, Comm. ACM 15 (9) (1972) 820--826. \bibitem{Benmabrouk2} A. Ben Mabrouk, M. Ayadi; \emph{A linearized finite-difference method for the solution of some mixed concave and convex non-linear problems}, Appl. Math. Comput. 197 (2008), 1--10. \bibitem{Benmabrouk1} A. Ben Mabrouk, M. Ayadi; \emph{Lyapunov type operators for numerical solutions of PDEs}, Appl. Math. Comput. 204 (2008), 395--407. \bibitem{Benmabrouk3} A. Ben Mabrouk, M. L. Ben Mohamed, K. Omrani; \emph{Finite-difference approximate solutions for a mixed sub-superlinear equation}, Appl. Math. Comput. 187 (2007) 1007-1016. \bibitem{Bratsos1} A. G. Bratsos, Ch. Tsituras, D. G. Natsis; \emph{Linearized numerical schemes for the Boussinesq equation}. Appl. Num. Anal. Comp. Math. 2(1) (2005), 34-53. \bibitem{Clarkson1} P. A. Clarkson; \emph{Nonclassical symmetry reductions for the Boussinesq equation}, Chaos, Solitons and Fractals, 5 (1995), 2261--2301. \bibitem{Dehghan2006} M. Dehghan; \emph{Finite difference procedures for solving a problem arising in modeling and design of certain optoelectronic devices}, Mathematics and Computers in Simulation 71 (2006), 16--30. \bibitem{Dehghan2008} M. Dehghan; \emph{Use of Hes homotopy perturbation method for solving a partial differential equation arising in modeling of flow in porous media}, Journal of Porous Media, 11(8) (2008), 765--778. \bibitem{Dehghan-Salehi} M. Dehghan, R. Salehi; \emph{A meshless based numerical technique for traveling solitary wave solution of Boussinesq equation}, Applied Mathematical Modelling 36 (2012), 1939--1956. \bibitem{Diaz1} J. I. D\^iaz, G. Galiano; \emph{Existence and uniqueness of solutions of the Boussinesq system with nonlinear thermal di¤usion}, Topological Methods in Nonlinear Analysis, (volume dedicated to O.A. Ladyzhenskaya), 11, 58-92, 1998. \bibitem{Diaz2} J. I. D\^iaz, J. M. Rakotoson, P. G. Schmidt; \emph{Local strong solutions of a parabolic system related to the Boussinesq approximation for buoyancy-driven ow with viscous heating} Advances in Di¤erential Equations, 13(9-10) (2008), 977-1000. \bibitem{ElMikkawy1} M. El-Mikkawy; \emph{A note on a three-term recurrence for a tridiagonal matrix}, Appl. Math. Computa., 139 (2003), 503-511. \bibitem{ElMikkawy2} M. El-Mikkawy; \emph{A fast algorithm for evaluating $n$th order tri-diagonal determinants}, J. Computa. \& Appl. Math., 166 (2004), 581-584. \bibitem{ElMikkawy3} M. El-Mikkawy; \emph{On the inverse of a general tridiagonal matrix}, J. Computa. \& Appl. Math., 150 (2004), 669-679. \bibitem{ElMikkawy-Karawia} M. El-Mikkawy, A. Karawia; \emph{Inversion of general tridiagonal matrices}, Appl. Math. Letters 19 (2006), 712-720. \bibitem{El-Mikkawy-Atlan} M. El-Mikkawy, F. Atlan; \emph{A new recursive algorithm for inverting general image-tridiagonal matrices}. Applied Mathematics Letters, 44 (2015), 34-39. \bibitem{Goloub} G. H. Golub, S. Nash, C. Van Loan; \emph{A Hessenberg--Schur method for the matrix problem AX + XB = C}, IEEE Trans. Automat. Control AC-24 (6) (1979) 909--913. \bibitem{Goncalves} E. Gon\c{c}alv\`{e}s; \emph{Resolution numerique, discretisation des EDP et EDO. Cours}, Institut National Polytechnique de Grenoble, 2005. \bibitem{HenriCartan} H. Cartan; \emph{Differential Calculus}, Kershaw Publishing Company, London 1971, Translated from the original French text Calcul differentiel, first published by Hermann in 1967. \bibitem{Jafarizadeh1} M. A. Jafarizadeh, A. R. Esfandyari, M. Moslehi-fard; \emph{Exact solutions of Boussinesq equation}. Third International Conference on Geometry, Integrability and Quantization June 14-23, 2001, Varna, Bulgaria Ivailo M. Mladenov and Gregory L. Naber, Editors Coral Press, Sofia 2001, pp. 304-314. \bibitem{Jameson} A. Jameson; \emph{Solution of equation $AX+XB=C$ by inversion of an $M\times M$ or $N\times N$ matrix}. SIAM J. Appl Math. 16(5) (1968), 1020-1023. \bibitem{Jia-Li} J. Jia, S. Li; \emph{On the inverse and determinant of general bordered tridiagonal matrices}. Computers \& Mathematics with Applications, 69(6) (2015), 503-509. \bibitem{Kano1} T. Kano, T. Nishida; \emph{A mathematical justification for Korteweg-De Vries equation and Boussinesq equation of water surface waves}. Osaka J. Math. 23 (1986), 389-413. \bibitem{Kaya1} D. Kaya; \emph{Explicit solutions of generalized nonlinear Boussinesq equations}. Journal of Applied Mathematics 1(1) (2001), 29--37. \bibitem{Lai1} Sh. Lai, Y. Wub, Y. Zhou; \emph{Some physical structures for the (2+1)-dimensional Boussinesq water equation with positive and negative exponents}. Computers and Mathematics with Applications 56 (2008), 339--345. \bibitem{Liu1} Ch. Liu, Z. Dai; \emph{Exact periodic solitary wave solutions for the (2+1)-dimensional Boussinesq equation}. J. Math. Anal. Appl. 367 (2010), 444--450. \bibitem{Parlange1} J.-Y. Parlange, W. L. Hogarth, R. S. Govindaraju, M. B. Parlange, D. Lockington; \emph{On an Exact Analytical Solution of the Boussinesq Equation}, Transport in Porous Media 39 (2000), 339--345. \bibitem{Kirrini} P. Kirrinni; \emph{Fast algorithms for the Sylvester equation $AX-XB^{T}=C$}. Theoretical Computer Science. Volume 259, Issues 1--2, 28 May 2001, Pages 623-638. Elsevier. \bibitem{Rionero} S. Rionero; \emph{Stability results for hyperbolic and parabolic equations}. Transport Theory \& Statistical Physics, 25(3-5) (1996), 323-337. \bibitem{Shokri-Dehghan} A. Shokri, M. Dehghan; \emph{A Not-a-Knot meshless method using radial basis functions and predictor-corrector scheme to the numerical solution of improved Boussinesq equation}, Computer Physics Communications 181 (2010), 1990--2000. \bibitem{Simoncini} V. Simoncini; \emph{Computatioanl methods for linear matrix equations}, Course in Dipartimento di Matematica, Universita di Bologna, Piazza di Porta San Donato 5, I-40127 Bologna, Italia, March 12, (2013). \bibitem{Song1} M. Song, Sh. Shao; \emph{Exact solitary wave solutions of the generalized (2+1)-dimensional Boussinesq equation}. Applied Mathematics and Computation 217 (2010), 3557--3563. \bibitem{Varlamov1} V. Varlamov; \emph{Two-dimensional Boussinesq equation in a disc and anisotropic Sobolev spaces}. C. R. Mecanique 335 (2007), 548--558. \bibitem{Wazwaz1} A.-M. Wazwaz; \emph{Variants of the Two-Dimensional Boussinesq Equation with Compactons}, Solitons, and Periodic Solutions. Computers and Mathematics with Applications 49 (2005), 295-301. \bibitem{Yi1} Z. Yi, Y. Ling-ya, L. Yi-neng, Z. Hai-qiong; \emph{Periodic wave solutions of the Boussinesq equation}. J. Phys. A: Math. Theor. 40(21) (2007), 5539--5549. \bibitem{Yi2} Z. Yi, Y. Ling-Ya; \emph{Rational and Periodic Wave Solutions of Two-Dimensional Boussinesq Equation}. Commun. Theor. Phys. 49(4) (2008), 815-824. \end{thebibliography} \end{document}