Lyapunov-Sylvester operators for Generalized Nonlinear Euler-Poisson-Darboux System

TL;DR

This paper develops a new numerical method using Lyapunov-Sylvester operators for solving a nonlinear Euler-Poisson-Darboux system, demonstrating improved accuracy and efficiency over existing methods.

Contribution

It introduces a novel finite difference scheme transforming the system into a Lyapunov-Sylvester operator framework, proving its stability, convergence, and efficiency.

Findings

The scheme is uniquely solvable, stable, and convergent.

Numerical results show higher accuracy than existing methods.

The method leads to faster computational algorithms.

Abstract

In this paper a nonlinear Euler-Poisson-Darboux system is considered. In a first part, we proved the genericity of the hypergeometric functions in the development of exact solutions for such a system in some special cases leading to Bessel type differential equations. Next, a finite difference schemle in two-dimensional case has been developed. The continuous system is transformed into an algebraic quasi linear discrete one leading to generalized Lyapunov-Sylvester operators. The discrete algebraic system is proved to be uniquely solvable, stable and convergent based on Lyapunov criterion of stability and Lax-Richtmyer equivalence theorem for the convergence. A numerical example has been provided at the end to illustrate the efficiency of the theoretical results. The present method is thus proved to be more accurate than existing ones and lead to faster algorithms.

Tables

Table 1. Table 1:

J	ErI	ErII	Relative ErI	Relative ErII	TimeI	TimeII
24	4.01E-3	3.31E-3	0.1415	0.1382	20 s	05 s
49	4.01E-3	3.10E-4	0.1423	0.1380	49 s	7.1 s
99	3.80E-3	3.20E-5	0.1431	0.1381	466 s	20 s
149	3.10E-4	3.03E-5	0.1433	0.1383	1002 s	31.2 s
199	2.80E-4	2.19E-5	0.1436	0.1384	2114 s	35.6 s
399	2.30E-4	1.29E-6	0.1441	0.1384	3215 s	41.2 s
499	2.10E-4	1.11E-6	0.1445	0.1385	4001 s	46 s
999	2.01E-4	1.02E-6	0.1446	0.1386	4516 s	49.3 s

Equations

$$\left\{\begin{array}[]{lll}u_{tt}+\Gamma_{a}(t)v_{t}=\Delta u+<F_{\lambda,\gamma},\nabla v>+|u|^{p-1}v,\\ v_{tt}+\Gamma_{a}(t)u_{t}=\Delta v+<F_{\lambda,\gamma},\nabla u>+|v|^{q-1}u\end{array}\right.$$

$$u(x,y,t_0)=u_0(x,y)and\frac{\partial u}{\partial t}(x,y,t_0)=u_1(x,y),(x,y)\in \Omega$$

$$\frac{\partial u}{\partial \eta }(x,y,t)=0,((x,y),t)\in \partial \Omega \times (t_0,+\infty ).$$

$$q{\psi }^{\prime \prime }(x)+\psi (x)+\alpha {\psi }^2(x)=ax+b,$$

$$\left\{\begin{array}[]{lll}u_{tt}+\Gamma_{a}(t)v_{t}=\Delta u+<F_{\lambda,\gamma},\nabla v>,\\ v_{tt}+\Gamma_{a}(t)u_{t}=\Delta v+<F_{\lambda,\gamma},\nabla u>.\end{array}\right.$$

$$u(x,y,t)=v(x,y,t)=\phi (x,y)=f(x)+g(y).$$

$$\Delta u+<F_{\lambda ,\gamma },\nabla v>=\Delta v+<F_{\lambda ,\gamma },\nabla u>=0.$$

$$f^{\prime \prime }(x)+\frac{2\lambda }{x}{f}^{\prime }(x)=g^{\prime \prime }(y)+\frac{2\gamma }{y}{g}^{\prime }(y).$$

$$\left\{\begin{array}[]{lll}f^{\prime\prime}(x)+\frac{2\lambda}{x}f^{\prime}(x)=K\\ g^{\prime\prime}(y)+\frac{2\gamma}{y}g^{\prime}(y)=-K.\end{array}\right.$$

$$f(x)=\frac{K_1}{1-2\lambda }|x{|}^{1-2\lambda }sign(x)+\frac{K}{1+2\lambda }\frac{x^2}{2},x\ne 0$$

$$g(y)=\frac{K_2}{1-2\gamma }|y{|}^{1-2\gamma }sign(y)-\frac{K}{1+2\gamma }\frac{y^2}{2},y\ne 0.$$

$$f(x)=K_1\log |x|+\frac{K}{6}{x}^3,x\ne 0$$

$$g(y)=K_2\log |y|+\frac{K}{6}{y}^3,y\ne 0.$$

$$f(x)=\frac{x^2}{2}\left(K\log |x|+K_1-\frac{1}{2}\right),x\ne 0$$

$$g(y)=\frac{y^2}{2}\left(-K\log |y|+K_2-\frac{1}{2}\right),y\ne 0.$$

$$u(x,y,t)=v(x,y,t)=\phi (x,y)=f(x)g(y).$$

$$g(y)\left(f^{\prime \prime }(x)+\frac{2\lambda }{x}{f}^{\prime }(x)\right)=f(x)\left(g^{\prime \prime }(y)+\frac{2\gamma }{y}{g}^{\prime }(y)\right).$$

$$\left\{\begin{array}[]{lll}f^{\prime\prime}(x)+\frac{2\lambda}{x}f^{\prime}(x)=Kf(x)\\ g^{\prime\prime}(y)+\frac{2\gamma}{y}g^{\prime}(y)=-Kg(y).\end{array}\right.$$

$$f(x)=|x{|}^{-\lambda }\left(a_0\cos (\sqrt{K}x)+\frac{a_1}{\sqrt{K}}\sin (\sqrt{K}x)\right)$$

$$g(y)=|y{|}^{-\gamma }(b_0\cos (\sqrt{K}y)+\frac{b_1}{\sqrt{K}}\sin (\sqrt{K}y)$$

$$y(x)=|x{|}^{\nu }\sum _{n=0}^{+\infty }{a}_n{x}^n.$$

$$\left\{\begin{array}[]{lll}a_{0}(1+2\lambda\nu)=0,\\ a_{1}(1+2(\lambda+\nu)+2\lambda\nu)=0,\\ (1+2\lambda\nu+2(\lambda+\nu)n+n(n-1))a_{n}=Ka_{n-2},\;\forall\,n\geq 2.\end{array}\right.$$

$$a_n=\frac{K}{n(n-1)}{a}_{n-2},n\ge 2$$

$$(1+2\lambda \nu +n(n-1))a_n=K{a}_{n-2},\forall n\ge 2.$$

$$(n-p)(n+p-1)a_n=K{a}_{n-2},\forall n\ge 2$$

$$a_{2n}=\frac{1}{(p-1)!}\frac{(n+k)!K^{n-k}}{(n-k)!(2m+p)!}{a}_p$$

$$f(x)=\frac{a_p|x{|}^{-\lambda }}{(p-1)!}\sum _{n=k}^{+\infty }\frac{(n+k)!K^{n-k}}{(n-k)!(2m+p)!}{x}^{2n}.$$

$$a_{2n+1}=(2p-1)!\frac{(n+k)!K^{n-k}}{(n-k)!(2m+p)!}{a}_p$$

$$f(x)=a_p(2p-1)!|x{|}^{-\lambda }\sum _{n=k}^{+\infty }\frac{(n+k)!K^{n-k}}{(n-k)!(2m+p)!}{x}^{2n+1}.$$

$$n(n+s)a_n=K{a}_{n-2},\forall n\ge 2.$$

Full text

1 Introduction

In this work we use Lyapunov-Sylvester algebraic operators to approximate the solutions of a generalized Euler-Poisson-Darboux (EPD) system in two-dimensional case. The present article is precisely devoted to the development of a numerical method based on two-dimensional finite difference scheme to approximate the solution of the generalized EPD system in $\mathbb{R}^{2}$ in the presence of mixed power laws nonlinearities. Denote for $a\in\mathbb{R}$, $\Gamma_{a}(x)=\displaystyle\frac{2a}{x}$ and for $\lambda,\gamma$ in $\mathbb{R}$, $F_{\lambda,\gamma}(x)=(\Gamma_{\lambda}(x),\Gamma_{\gamma}(x))$. We consider the following evolutive system.

[TABLE]

with initial conditions

[TABLE]

and boundary conditions

[TABLE]

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_{t}$ is the first order partial derivative in 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 on $\mathbb{R}^{2}$. $\frac{\partial}{\partial\eta}$ is the outward normal derivative operator along the boundary $\partial\Omega$. Finally, $u$, $u_{0}$ and $u_{1}$ are real valued functions with $u_{0}$ and $u_{1}$ are $\mathcal{C}^{2}$ on $\overline{\Omega}$. $u$ and $v$ are the unknown candidates supposed to be $\mathcal{C}^{4}$ on $\overline{\Omega}$. $p$ and $q$ are reel parameters such that $p,q>1$.

Several papers have been devoted to the study of existence and uniqueness of solutions of the linear problem 1 (without the nonlinear parts) for the case where $u\equiv v$ and where the system is in fact the well known EPD equation. In some of these studies, 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

[TABLE]

for some constants $a$ and $b$ depending on the initial-boundary conditions. The existence and uniqueness problems are overcame using ODEs.

In [25], a generalized form of EPD equation is considered and particular solutions are constructed in an explicit form expressed by the Lauricella hypergeometric function of three variables. Properties of each constructed solutions have been investigated in sections of surfaces of the characteristic cone. The authors proved that found solutions have singularity $\frac{1}{r}$ as $r\rightarrow 0$, where $r^{2}=(x-x_{0})^{2}+(y-y_{0})^{2}-(t-t_{0})^{2}$.

In [27], a singular Cauchy problem for the multi-dimensional EPD equation with spectral parameter has been investigated with the help of the generalized Erdelyi-Kober fractional operator. Solution of the considered problem is found in explicit form for various values of the parameter $p$ of the equation.

In the present work, we intend to apply some algebraic operators to develop numerical solutions for EPD system. It consists of the well known Lyapunov-Sylvester operators. 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 classical idea to obtain 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,[18]. 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 [26].

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 (See [2] and [16]).

In [20] the author described 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}\log N)$ 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.

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 quasi linear Lyapunov-Sylvester system. The motivation behind the application of Lyapunov-Syslvester operators was already evoked in [3]. 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.

• •

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.

• •

The method developed is not just a resolution of a PDE. We recall already that the resolution itself is not a negligible aim. Further, it proves the efficiency of algebraic operators other than classical tri-diagonal ones.

• •

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.

• •

The fact of 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.

In the present work, existence and multiplicity of the solutions of problem (1)-(3) are developed in some special cases. We showed that special functions such as hypergeometric series and Bessel function are generic for developing such solutions. Section 3 is devoted to the development of a 2-dimensional discrete scheme to transform the continuous problem (1)-(3) to a discrete one. A system of generalized Lyapunov-Sylvester equations is obtained. The solvability of such a discrete system is proved next in section 4. Section 5 is concerned with the consistency, stability and the convergence of the discrete Lyapunov-Sylvester problem obtained in section 3. The crucial idea is the application of the truncation error for consistency, Lyapunov cretirion for stability and the Lax equivalence theorem for the convergence. Section 6 is devoted to the development of a numerical example. The performance of the discrete scheme is proved by means of error estimates as well as fast algorithms. The conclusion is finally subject of section 7.

2 On the existence and multiplicity of the solutions of the continuous problem

In this part, we review the generalized linear EPD system. We will show the genericity of hypergeometric functions in the development of solutions. The generalized linear EPD system is

[TABLE]

2.1 A first class of time-independent solutions

In this section we propose to develop a first class of solutions of problem (4) on the time-independent additive form

[TABLE]

The stationary problem associated to the system (4) is

[TABLE]

Substituting $f$ and $g$ in the last equations yield that

[TABLE]

Hence, there exists a constant $K\in\mathbb{C}$ for which

[TABLE]

Therefore, whenever $\lambda,\gamma\notin\{\pm\frac{1}{2}\}$, by applying classical resolution of ODEs we get the general solutions

[TABLE]

and

[TABLE]

For the case when $\lambda=\gamma=\displaystyle\frac{1}{2}$ we get

[TABLE]

and

[TABLE]

Similarly, for the case $\lambda=\gamma=-\displaystyle\frac{1}{2}$ we get

[TABLE]

and

[TABLE]

The remaining solutions (($\lambda\not=\pm\frac{1}{2}$ and $\gamma=\pm\frac{1}{2}$) and ($\lambda=\pm\frac{1}{2}$ and $\gamma\not=\pm\frac{1}{2}$)) are composed of linear combinaisions of the developed cases above.

2.2 Second class of time-independent solutions

In this section, we continue to develop a second class of time-independent but multiplicative solutions. We consider solutions of the form

[TABLE]

Substituting as for the previous section $\varphi$ in the stationary system (4) we get

[TABLE]

Hence, there exists a constant $K\in\mathbb{C}$ for which

[TABLE]

For example, when $\lambda=\gamma=\pm\displaystyle\frac{1}{\sqrt{2}}$, an explicit solution may be obtained by

[TABLE]

and

[TABLE]

As previously, the remaining solutions (($\lambda\not=\pm\frac{1}{\sqrt{2}}$ and $\gamma=\pm\frac{1}{\sqrt{2}}$) and ($\lambda=\pm\frac{1}{\sqrt{2}}$ and $\gamma\not=\pm\frac{1}{\sqrt{2}}$)) are composed of linear combinaisions of the solutions developed above.

2.3 A hypergeometric/Bessel type solution

In this section we will prove that the EPD time-independent system has solutions that may be expressed by means of the famous special functions such as the hypergeometric one and its general variants especially Bessel function. For this we assume that the solutions $u$ and $v$ are on the form (7) and that $f$ and $g$ are of the form

[TABLE]

Substituting in the first equation of system (6), we obtain the following recurrence system

[TABLE]

So, for $\lambda=-\nu=\pm\displaystyle\frac{1}{\sqrt{2}}$, we get

[TABLE]

which yields the same oscillating singular solutions provided in (9) and (10).

For $\lambda=-\nu$ but $1+2\lambda\nu\not=0$, we get firstly $a_{0}=a_{1}=0$. It remains in (11) just one equation to handle,

[TABLE]

Assume now that there exists $p=2k\in\mathbb{N}$ for which $1+2\lambda\nu=-p(p-1)$. Equation (12) becomes

[TABLE]

which yields that $a_{2n+1}=0$, for all $n$, and that

[TABLE]

and thus $f$ (and similarly $g$) may be expressed by means of a hypergeometric series

[TABLE]

For $p=(2k+1)\in\mathbb{N}$ satisfying the same hypothesis $1+2\lambda\nu=-p(p-1)$, we get $a_{2n}=0$, for all $n$, and

[TABLE]

and thus $f$ (and similarly $g$) may be expressed by

[TABLE]

Assume now that $\lambda\not=\nu$ and that already $1+2\lambda\nu=0$, and denote $r=\lambda+\nu$ and $s=2r-1$. From equation (11) we obtain immediately $a_{1}=0$ and that

[TABLE]

So, when $s\in-\mathbb{N}$, this yields that $a_{2n+}=0$, for all $n$ and that

[TABLE]

As a result,

[TABLE]

Now, when $\lambda\not=\nu$, $1+2\lambda\nu=0$ and $s=2(\lambda+\nu)-1=2r-1=-2k\in-\mathbb{N}$, the recurence relation (11) permits to obtain $a_{1}=0$ and

[TABLE]

which in turn yields that $a_{2n+1}=0$, for all $n$ and that

[TABLE]

As a result,

[TABLE]

Next, for $\lambda\not=-\nu$, $1+2\lambda\nu\not=0$ but $1+2\lambda\nu+2(\lambda+\nu)=0$ and $\lambda+\nu=-k\in-\mathbb{N}$, the recurrence (11) becomes

[TABLE]

Therefore, $a_{2n}=0$ for $n\leq k-1$ and for $n\geq k$, we get

[TABLE]

Similarly, for $n\geq k$ we get

[TABLE]

and for $n\leq k-1$, we have

[TABLE]

Thus we get for $\lambda\not=-1$,

[TABLE]

For $\lambda=-1$ and $\nu=\displaystyle\frac{1}{2}$, we get $a_{0}=0$ and

[TABLE]

and

[TABLE]

Hence,

[TABLE]

Now, for $\lambda=-1$ and $\nu=\displaystyle\frac{3k-k^{2}-1}{2(k-1)}$ for some $k\in\mathbb{N}$, $k\geq 3$ (to guaranty that $\nu\not=\displaystyle\frac{1}{2}$, we get $a_{0}=a_{1}=0$ and

[TABLE]

where $2\eta=1-2\nu-\sqrt{(1-2\nu)^{2}+4}$. Of course, in the present case $\eta\notin\mathbb{N}$. More precisely, for $\nu<\displaystyle\frac{1}{2}$, we have $0<\eta<1$ and for $\nu>\displaystyle\frac{1}{2}$, it holds that $\eta<0$. Hence, whenever $k=2p$, we obtain $a_{2n+1}=0$ for all $n$ and $a_{2n}=0$ just for $n\leq p-1$. The recurrence relation (23) yields that

[TABLE]

where $\zeta=\displaystyle\frac{\eta}{2}$. Consequently,

[TABLE]

where $\nu_{k}=\displaystyle\frac{(3k-k^{2}-1)}{2(k-1)}$. Similarly, whenever whenever $k=2p+1$, we obtain $a_{2n}=0$ for all $n$ and $a_{2n+1}=0$ just for $n\leq p-1$ and from the recurrence relation (23) we get

[TABLE]

Consequently,

[TABLE]

It remains now to study the case when $1+2\lambda\nu\not=0$ and $1+2(\lambda+\nu)+2\lambda\nu\not=0$. It holds so that $a_{0}=a_{1}=0$ and remains in (11) that

[TABLE]

So assume now that $(\lambda,\nu)\in D(\Omega,\displaystyle\frac{\sqrt{5}}{2})$ the open disc of center $\Omega(\displaystyle\frac{1}{2},\displaystyle\frac{1}{2})$ and radius $R=\displaystyle\frac{\sqrt{5}}{2}$. Then the second order equation $n^{2}+(2(\lambda+\nu)-1)n+(1+2\lambda\nu))>0$ for all $n$ which yields with the fact that $a_{0}=a_{1}=0$ the null solution.

Next, assume that $(\lambda,\nu)\in C(\Omega,\displaystyle\frac{\sqrt{5}}{2})$ the circle of center $\Omega(\displaystyle\frac{1}{2},\displaystyle\frac{1}{2})$ and radius $R=\displaystyle\frac{\sqrt{5}}{2}$. We get

[TABLE]

with $\eta=\displaystyle\frac{1-2(\lambda+\nu)}{2}$. So, whenever there exists $k\in\mathbb{N}$ for which $\eta=2k$, we obtain $a_{2n+1}=0$ for all $n$, and

[TABLE]

and $a_{2n}=0$ otherwise. Hence, the generated solution will be on the form

[TABLE]

where $\nu=\displaystyle\frac{1}{2}\pm\sqrt{\displaystyle\frac{5}{4}-(\lambda-\displaystyle\frac{1}{2})^{2}}$.

Now, whenever there exists $k\in\mathbb{N}$ for which $\eta=2k+1$, we obtain $a_{2n}=0$ for all $n$, and

[TABLE]

and $a_{2n}=0$ otherwise. Hence, the generated solution will be on the form

[TABLE]

for the same expression of $\nu=\displaystyle\frac{1}{2}\pm\sqrt{\displaystyle\frac{5}{4}-(\lambda-\displaystyle\frac{1}{2})^{2}}$ as in the last previous case.

Next, assume that $(\lambda,\nu)\notin\overline{D}(\Omega,\displaystyle\frac{\sqrt{5}}{2})$ the complementary of the closed disc of center $\Omega(\displaystyle\frac{1}{2},\displaystyle\frac{1}{2})$ and radius $R=\displaystyle\frac{\sqrt{5}}{2}$ and denote

[TABLE]

and

[TABLE]

so that the recurrence (11) becomes

[TABLE]

Next, as previously we get one of the following solutions.

[TABLE]

and $a_{2n}=0$ otherwise, where $k=\displaystyle\frac{\eta_{1}}{2}\in\mathbb{N}$ and $\zeta_{2}=\displaystyle\frac{\eta_{2}}{2}$. In this case, of course, $a_{2n+1}=0$ for all $n$. Next, whenever $k=\displaystyle\frac{\eta_{1}-1}{2}\in\mathbb{N}$, we get

[TABLE]

Hence, the generated solution are on the respective form

[TABLE]

and

[TABLE]

2.4 A last hypergeometric in time class of solutions

In this section we will prove that the EPD system has solutions of the form

[TABLE]

where $\psi$ is an hypergeometric type function. Indeed, substituting in (4) we obtain

[TABLE]

where $K$ is a real constant.

The function $\psi$ takes the form of the solution of problem (8) developed in section 2.2. Next, to develop a solution $\varphi$ of the problem

[TABLE]

we may proceed as in the previous sections by assuming either

[TABLE]

or

[TABLE]

In both cases, we get

[TABLE]

where $f_{H}$ is the solution of the homogeneous problem (8) and $C$ is a suitable constant. Indeed, for the first choice we get the following system

[TABLE]

For $K=0$, we get analogous problem as (6) by replacing $K$ by $\widetilde{K}$ and thus the solutions may be obtained from section 2.1. Next, whenever $K\not=0$, a particular solution is

[TABLE]

and the homogeneous problem in this case is the same as (8).

3 Discrete two-dimensional Nonlinear EPD system

The object of this section is to explain the discretization scheme proposed to transform problem (1)-(3) into a discrete quasi-linear one. Let $n\in\mathbb{N}$ and consider a time step $l=\Delta t$ and a space one $h=\displaystyle\frac{L_{1}-L_{0}}{J+1}$. Next, denote for $k\in\mathbb{N}$ and $j,m\in\{0,...,J+1\}$

[TABLE]

so that the cube $[L_{0},L_{1}]\times[L_{0},L_{1}]$ is subdivided into cubes $C_{j,m}=[x_{j},x_{j+1}]\times[y_{m},y_{m+1}]$. For a function $z$ defined on the cube $[L_{0},L_{1}]\times[L_{0},L_{1}]$, we denote by small $z_{j,m}^{k}$ the net function $z(x_{j},y_{m},t^{k})$ and capital $Z_{j,m}^{k}$ the numerical approximation. Consider next the discrete finite difference operators

[TABLE]

[TABLE]

[TABLE]

[TABLE]

Denote also $G(u,v)=|u|^{p-1}v$ and $H(u)=|v|^{p-1}u$.

The discretization of the first equation of problem (1) yields that

[TABLE]

Denote next

[TABLE]

Par suite

[TABLE]

or equivalently,

[TABLE]

This may be written in a matrix-vector form

[TABLE]

where $A$, $\Theta$ and $\Lambda$ are the matrices given by

[TABLE]

and for $1\leq j,m\leq J$,

[TABLE]

[TABLE]

and

[TABLE]

Denote next,

[TABLE]

[TABLE]

[TABLE]

We get

[TABLE]

Similarly, for the second equation we get

[TABLE]

4 Solvability of the discrete problem

In [3], 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 [3], but we develop different arguments. The first main result is stated as follows.

Theorem 4.1

The system (34)-(35) is uniquely solvable whenever $U^{0}$ and $U^{1}$ are known.

Proof. It reposes on the inverse of Lyapunov-Syslvester operators. Consider the endomorphism $\Phi$ defined by

[TABLE]

where for matrices $A$ and $B$, $\mathcal{L}_{A,B}$ is the Lyapunov-Sylvester operator defined by

[TABLE]

and where $\mathcal{L}_{A}\equiv\mathcal{L}_{A,A}$. To prove Theorem 4.1, it suffices to show that $ker\Phi$ is reduced to [math]. Indeed, whenever $l=o(h)$ and $l,h\rightarrow 0$, we get

[TABLE]

Next, whenever $\Phi(X,Y)=0$, we get

[TABLE]

Which means that $X=Y=0$. Next, we apply the following result.

Lemma 4.1

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.

Observing that the operator $\Phi$ obtained above is invertible, we get that $\Phi_{l,h}$ is invertible for $l,h$ small enough.

5 Consistency, stability and convergence of the discrete method

Recall firstly that the consistency of the numerical scheme is always done by evaluating the local truncation error arising from the discrete and the continuous problem. Applying Taylor’s expansion in the discrete equations raised in section 2, we get the following truncation principal part for the first equation in system (1)

[TABLE]

and for the second equation, we get

[TABLE]

where $\Delta_{2}=\displaystyle\frac{\partial^{4}}{\partial x^{4}}+\frac{\partial^{4}}{\partial y^{4}}$ and $\nabla_{3}=(\displaystyle\frac{\partial^{3}}{\partial x^{3}},\frac{\partial^{3}}{\partial y^{3}})$. WE then have the following lemma.

Lemma 5.1

The discrete scheme is consistent with order $(l^{2}+h^{2})$.

Now, we will examine the stability of the scheme. We will apply the same method as in [3], [5], [6] and [9] based on the Lyapunov criterion of stability. A linear system $\mathcal{L}(u_{n+1},u_{n},u_{n-1},\dots)=0$ is stable in the sense of Lyapunov iff for any bounded initial value $u_{0}$, the solution $u_{n}$ ramains bounded for all $n\geq 0$.

Lemma 5.2

$\mathcal{P}_{n}$: The solution $(U^{n},V^{n})$ is bounded independently of $n$ whenever the initial solution $(U^{0},V^{0})$ is bounded.

Before going on proving this result, we stress on the fact that contrarily to previous studies such as [3], [5], [6] and [9], we are confronted in the present study to Lyapunov-Sylvester operators that are non commutative, which raises some new difficulties in the proof of the stability and thus imposes different ideas.

Proof. Recall firstly that the discrete scheme in the matrix form may be written as

[TABLE]

where $\Phi_{l,h}$ is the Lyapunov-Sylvester operator defined in (36) and where the linear operator $F=(F_{1},F_{2})$ (which is a melange of Lyapunov-Sylvester operators) is defined by

[TABLE]

and

[TABLE]

where $G(X,Y,Z,T)$ is the matrix with coefficients

[TABLE]

and $H(X,Y,Z,T)$ is the matrix with coefficients

[TABLE]

We remark immediately that

[TABLE]

and

[TABLE]

Consequently,

[TABLE]

and

[TABLE]

Recall also that

[TABLE]

Observe next that for $l,h$ small enough

[TABLE]

So, for $n$ large enough, we may obtain

[TABLE]

On the other hand, we have

[TABLE]

Next, observe that

[TABLE]

and consequently, by denoting

[TABLE]

we get

[TABLE]

Conseqyently,

[TABLE]

Next, by choosing the time step $l$ and the space step $h$ so that $4\sigma C_{\alpha}<1$, we obtain

[TABLE]

Finally, the lemma follows by recurrence on $n$ by observing for the choice above that

[TABLE]

and the continuity of $F$.

Now, it remains finally to check the convergence of the discrete scheme. This is done by a direct application of the following well-known result [22].

Theorem 5.1

(Lax Equivalence Theorem). For a consistent finite difference scheme, stability is equivalent to convergence.

Lemma 5.3

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

6 Numerical implementation

In this section we present a numerical example that links between the EPD system studied here and the classical EPD equation in some parts so that we point out two aims; we firstly validate the theoretical results developed previously and secondly to show that such a system may be considered as a perturbation of the classical EPD equation.

To measure the closeness of the numerical solution and the exact one, the error is evaluated via an $L_{2}$ matrix norm

[TABLE]

for a matrix $X=(X_{ij})\in\mathcal{M}_{J+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

[TABLE]

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

[TABLE]

on the same grid. Next, we will assume that

[TABLE]

We denote also

[TABLE]

We consider the inhomogeneous system

[TABLE]

where $G_{1}$ and $G_{2}$ are explicited respectively by

[TABLE]

and

[TABLE]

The exact solution of such a system is given by

[TABLE]

We fix the domain bounds as follows:

[TABLE]

Finally tu fulful the assumption $l=o(h)$ we assume that $l=h\sqrt{h}$. The following table resumes the error estimates, the relative error as well as the time of execusion for the corresponding algorithm (denoted here II) compared to the time execusion of the classical tri-diagonal one (denoted I) for different values of the space step.

In numerical studies of PDEs one important task is the convergence of algorithms and especially the rate of convergence. Different methods have been developed to get fast algorithms. In the present work, Table 1 above includes a comparison between existing method (denote Method I) and the present scheme based on the generalized Lyapunov-Syslvester form (denote method II) for different values of the parameter $J$. It is remarkable that our method is faster than the classical one. This is perfect as nowadays focuses are on big and/or cloud data and thus may seek fast and accurate algirithms. In Table 1, we noticed easily that an accelerated procedure is pointed out. For $J=24$ for example, a 4-times faster algorithm is obtained. For $J=99$ the running time is reduced to be more than 20-times. Increasing more the mesh size ($J$) results in more and more best faster algorithms with a rate of running time proportion crossing 90-times for $J=999$.

7 Conclusion

In this work, computational method has been developed for numerical solutions of Poisson-Darbous-Euler system in 2-D case. The method has yielded algebraic systems based on general Lyapunov-Sylvester equations. Fast and accurate algorithms have been obtained compared with the associated tri-diagonal classical algorithms always applied in such problems.