RIADH Chteoui ,SABRINE Arfaoui and ANOUAR Ben Mabrouk

1 Laboratory of Algebra,Number Theory and Nonlinear Analysis LR15ES18,Department of Mathematics,Faculty of Sciences,5019 Monastir. Tunisia.

2 Department of Mathematics,Higher Institute of Applied Mathematics and Computer Science,University of Kairouan,Street of Assad Ibn Al-Fourat,Kairouan 3100,Tunisia.

3 Department of Mathematics,Faculty of Sciences,University of Tabuk,KSA.

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 scheme 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 numerical scheme developed in section 3. The present method is thus proved to be more accurate than existing ones and lead to faster algorithms.

Key Words: Finite difference method;Lyapunov-Sylvester operators;generalized Euler-Poisson-Darboux equation;hyperbolic equation;Lauricella hypergeometric functions.

1 Introduction

In the present paper a nonlinear Euler-Poisson-Darboux system is studied for two folds.In a first part,exact solutions based on general hypergeometric series and on Bessel functions in some special cases have been developed. The application of such forms showed that the system may be transformed in some cases to Bessel type differential equations. The general system of coupled equations is characterized by the the presence of some cross-correlated nonlinearities characterized by the presence of simultaneous superlinear two power laws which interchange the role according to the dynamical system studied.This type of problems is very important in plasma physics and engineering as it is used to model telegraphic phenomena,turbulense especially for plasma and for accelerating electrons.

In a second part, Lyapunov-Sylvester algebraic operators have been applied to approximate the numerical solutions of the 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 in2in the presence of mixed power laws nonlinearities. Denote fora∈Γa(x)=2a/xand forλ,γin,Fλ,γ(x)=(Γλ(x),Γγ(x)). We consider the evolutive system

with initial conditions

and boundary conditions

on a rectangular domain ?=[L0,L1]×[L0,L1]in2.t0≥0 is a real parameter fixed as the initial time,utis the first order partial derivative in time,uttis the second order partial derivative in time,is the Laplace operator onis the outward normal derivative operator along the boundary??.Finally,u,u0andu1are real valued functions withu0andu1are C2onuandvare the unknown candidates assumed to be C4onpandqare real parameters such thatp,q＞1.

In the present work,existence and multiplicity of the solutions of problem(1.1)-(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 two-dimensional discrete scheme to transform the continuous problem (1.1)-(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 Review and motivations

EPD equation as well as systems are widely applied and related to many phenomena such as telegraphic models, acceleration models of electrons by the plasma turbulence and general dynamical systems.See[1-11].

Remark that problem(1.1)may be seen as perturbation of the linear problem

Wheneveru≡v,this becomes the famous linear EPD equation

Several papers have been devoted to the study of existence and uniqueness of solutions of problem(2.2). In some cases,exact solutions are developed such as solitary,stationary,time-independent, one-dimensional, ..etc. For example, in the case of a one-direction viscous fluid we may seek solutions of the formu(x,y,t)=αψ(x).In this case,the problem is transformed into a one variable ordinary differential equation

for some constantsα,q,aandbdepending on the initial-boundary conditions. The existence and uniqueness problems are studied using ODEs. In[12]and[13],weak solutions of an EPD initial value problem are constructed using distributional methods. Firstly the Fourier transform with respect to the space variables is applied leading to a Bessel type differential equation which is easily solved. Next,the inverse Fourier transform is applied to obtain the solutions of the original EPD initial value problem. Conditions for the regularity of the solutions are also provided depending on the problem parameters.In[14],a Cauchy problem with modified conditions for the EPD equation is studied.Explicit solutions are investigated in terms of Gauss and Appell hypergeometric functions.In [15], a time Fuchsian type singular ramified Cauchy data problem for EPD equation is considered. An expansion of the solutions in a series of hypergeometric functions has been established and the nature of the singularities of the solutions has been investigated. In[16], the wave equation in space of any odd number of dimensions is considered without need for evaluating divergent integrals leading to a simple method to solve a regular Cauchy problem for EPD equation in the same space. In[17],criteria for blow-up solutions in finite time are given for an hyperbolic initial nonlocal boundary-value EPD problem.Under suitable assumptions on the problem parameters blow-up solution have been proved to exist with finite blow-up time.In[18],random motions on the line and in the space Rdin higher dimensions governed by a non-homogeneous Poisson process has been analysed with explicit distributions of the position of the randomly moving particles seen as solutions of initial-value problems governed by the EPD equation. In [19],a singular Cauchy EPD problem is studied by applying Riemann function and Appell?s hypergeometric one to evaluate the solutions.In[20],a general analytical solution of one dimensional regular EPD Cauchy problem has been investigated using similarity transformation to get an ordinary differential equation. In[21],solutions of a Cauchy problem for a general form of the EPD equation have been represented in the compact integral form using Bessel operators via generalized translation and spherical mean operators.Under suitable conditions a distributional solution has been proved to be a classical one.In [22], 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 solution have been investigated in sections of surfaces of the characteristic cone. The authors proved that the solutions have singularitywherer2=(x?x0)2+(y?y0)2?(t?t0)2. In[23],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 parameters of the equation.

For the case of systems of the form(2.1),in the best knowlege of the author there is few works dealing with. In[24],an explicit solution of a mixed problem for the EPD equation is obtained in the quarter spacex＞0,?∞＜y＜∞,t＞0. Uniqueness of solution has been proved to hold for some special values of the problem parameters and is lost for some others. Even though existence of nontrivial solutions is already holding. The method applied is based on analytic continuation of a generalized Riemann-Liouville integral developed by Riesz. In[25],a derivation of the EPD equation from the Laplace equation in R3is presented.Using symmetries of the Laplace equation,linear differential relations of the first order between solutions of the EPD equation are obtained leading to systems of EPD types.The method is based on the solution of the over-determined system of partial differential equations for unknown coefficients and a second one is based on the group analysis approach. In[26], the well known Whitham-type equations have been considered and transformed to obtain in a single-phase case a linear over-determined system of EPD type. It is proved that such a system has a unique solution,and that the solution can be explicitly written. In [27], a system of EPD equations in matrix form has been investigated. Cauchy-Goursat and Darboux problems has been formulated for the case when the eigenvalues of the coefficient matrix are in(0,1/2). The Jordan form is applied allowing to transform the original EPD system into independent systems associated to each Jordan cell. Next,solutions of such independent sub-systems are constructed. The solutions of the original problems are obtained as direct sum of the solutions of Jordan cells systems. In[28], Dirichlet problem for a system of EPD equations is considered in a half-plane and strip perpendicular to the line of degeneracy.Under suitable conditions two classes of solution have been investigated. In[29], existence and uniqueness of the solutions for a singular mixed EPD problems and modified versions has been studied.In[30]and[31]systems of mixed equations issued from different PDEs are studied such as degenerate hyperbolic systems of conservation laws and the link with such laws and EPD equation is investigated through the entropy and Young measure. Such a connection is next applied for solving hyperbolic systems of conservation laws with parabolic or hyperbolic degeneracy. In [32], the classical EPD equation in multi variables space has been investigated. The original EPD equation is transformed into an equivalent systems of equations. Next,by using explicit form of Appell and Lauricella hypergeometric functions the solutions of such system has been explicitly written to yield the solution of the original EPD equation. In [33], a general exact solution for an-dimensional regular Cauchy problem of EPD equation has been studied using quitely the same transformations but with necessary modifications to get a system of ordinary differential equations and next to left to the original problem.

In the present work, we intend in a first part to investigated exact solutions based on general hypergeometric series and on Bessel functions in some special caseswhere the system may be transformed to Bessel type differential equations. In the second part,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 matricesA∈m×m,B∈n×n, andC∈m×n, the Sylvester equation is given by the formAX+XB=C. A classical idea to obtain the solutionXis to rewrite the Sylvester equation in standardmn×mnlinear systemusing the Kronecker product, [34]. The Sylvester equation can be solved by Gaussian elimination withO(m3n3) flops. This approach dramatically increases the complexity of the computation, and also cannot preserve the intrinsic properties of the problem in practice[35].

In numerical analysis,solving the Sylvester equation using the Bratels-Stewart and/or the Golub-Nash-Van Loan algorithms usesO(m3+n3)floating point operations if one assumes that anm×mmatrix can be reduced to Schur form (See [36] and [37]). In [38],the author described an algorithm that computes the solutionXover an arbitrary field F. The complexity of the algorithm forA∈m×m,B∈n×nandm,n≤NisO(NβlogN)arithmetic operations in F, whereβ＞2 is such thatm×mmatrices can be multiplied withO(mβ) arithmetic operations. This algorithm is competitive in terms of arithmetic operation and faster than the classical algorithms.

Recently,in[39],Lyapunov-Sylvester algebraic operators are used to approximate the solutions of the continuous problem due to Boussinesq equation in higher dimensions without adapting classical developments based on separation of variables, radial solutions, etc. The continuous problem is transformed into a system of Lyapunov-Sylvester equations in the matrix space. Next, investigating solvability of the numerical algebraic Lyapunov Sylvester problem has been itself done by using a system of Lyapunov-Sylvester equations in matrix form. Compared to the classical numerical methods such as tri-diagonal transformations,the new method leads to best algorithms when regarded for convergence rates, time execution and error estimates. Besides,there is no need for the present method to compute eigenvalues and precisely bounds/estimates of eigenvalues or direct inverses which remain complicated problems in general linear algebra and especially for generalized Lyapunov-Sylvester operators. Lyapunov-Sylvester operators has been proved to be good candidates as solvers compared to tri-diagonal and/or block tri-diagonal ones. The later methods are unadvised because of the fact that they are costing methods from both the machine memory and time. They also need for higher dimensional cases 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.This may induce time and accuracy losing and error degradation.

In the present paper,time and space partial derivatives are replaced by finite-difference approximations in order to transform the continuous problem (1.1)-(1.3) into quasi linear Lyapunov-Sylvester system. The motivation behind the application of Lyapunov-Sylvester operators are various.We recall 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. The following are some reasons to apply the method developed here(See[40]).

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 is the 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.

The use of Lyapunov-Sylvester operators is anymore standard. It may be related to Sturm operators.More precisely,consider differential operators of the form

whereαis a real parameter. The EPD equation may be seen on the form

From an algebraic point of view, the study of Lyapunov-Sylvster operators is still fascinating and is not standard.It is sometimes related to the notion of derivation on algebras.More precisely,when the operators LA,Bare defined forA,Bin some ringRby

we often call them generalized derivations on the ringR. These have been studied in the context of algebras on some normed spaces. Indeed,a generalized derivation on an algebraRis any map of the form LA,BwhereAandBare fixed elements inR. In the theory of operator algebras,they are considered as an important class of the elementary operatorswhere as previously theAi’s and theBi’s are elements of the algebraR. Remark that the operator LA,Bsatisfies

where DB(Y)=YB?BYis the so-called inner derivation. Remark that

which looks like to the product derivation rule for functions(f g)′=f′g+f g′. Hence,the nomination of derivation on algebraic structures.So,it is questionable that an integration operator IA,Bcould or not be defined so that one has for anyXandYinR,

For example,denoting

the associated Lyapunov-Sylvester operator or generalized derivation operator LA,Bis defined by

OnM2()we get an isomorphism.However,onM2(Z/3Z),it is not injective. Consider next the operator

where

We get an isomorphism for bothandFurthermore,we have

Surprisingly, this latter is not a Lyaponov-Sylvester neither a generalised derivation in the algebraic sense above. More precisely,no one of the following forms can occur:

So,what could be the reciprocal operatorIn fact,we know the general form

Here also,a second problem appears;what could be the values of the parametersmandnand next the matrices appearing there?

More details on derivation on algebras and especially on rings may be found in[41-49]. More details on LS operators may be found in[3-6,8-10,34,37-39,50-60,62-67].

3 Existence and multiplicity for 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

3.1 A first class of time-independent solutions

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

The stationary problem associated to the system(3.1)is

Substitutingfandgin the last equations yield that

Hence,there exists a constantfor which

Therefore, wheneverby applying classical resolution of ODEs we get the general solutions

and

For the case whenwe get

and

Similarly,for the casewe get

and

3.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

Substituting as for the previous section?in the system(3.1)we get

Hence,there exists a constantK∈C for which

For example,whenan explicit solution may be obtained by

and

3.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. To do this we assume that the solutionsuandvare on the form(3.4)and thatfandgare of the form

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

So,forwe get

which yields the same oscillating singular solutions provided in(3.6)and(3.7).

Forλ=?νand 1+2λν/=0,we get firstlya0=a1=0.It remains in(3.8)just one equation to handle,

Assume now that there existsp=2k∈N for which 1+2λν=?p(p?1). Equation (3.9)becomes

which yields thata2n+1=0,for alln,and that

and thusf(and similarlyg)may be expressed by means of a hypergeometric series

Forp=(2k+1)∈N satisfying the same hypothesis 1+2λν=?p(p?1),we geta2n=0,for alln,and

and thusf(and similarlyg)may be expressed by

Assume now thatand that already 1+2λν=0,and denoter=λ+νands=2r?1.From equation(3.8)we obtain immediatelya1=0 and thus

So,whens∈?,this yields thata2n+1=0,for allnand that

As a result,

Now,when1+2λν=0 ands=2(λ+ν)?1=2r?1=?2k∈?,the recurence relation(3.8)permits to obtaina1=0 and

which in turn yields thata2n+1=0,for allnand that

As a result,

Next,forλ/=?ν,1+2but 1+2λν+2(λ+ν)=0 andλ+ν=?k∈?,the recurrence(3.8)becomes

Therefore,a2n=0 forn≤k?1 and forn≥k,we get

Similarly,forn≥kwe get

and forn≤k?1,

Thus we get for

Forλ=?1 andwe geta0=0

and

Hence,

Now, forλ=?1 andfor somek∈,k≥3 (to guaranty thatwe geta0=a1=0 and

whereConsequently,

Consequently,

It remains now to study the case whenandIt holds thata0=a1=0 and it remains in(3.8)that

withSo,whenever there existsfor whichη=2k, we obtaina2n+1=0 for alln,and

anda2n=0 otherwise.Hence,the generated solution will be on the form

Now,whenever there existsk∈N for whichη=2k+1,we obtaina2n=0 for alln,and

anda2n=0 otherwise.Hence,the generated solution will be on the form

for the same expression ofas in the last previous case.

and

so that the recurrence(3.8)becomes

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

anda2n=0 otherwise,whereIn this case,of course,a2n+1=0 for alln. Next,wheneverwe get

Hence,the generated solutions are on the respective forms

and

3.4 A last hypergeometric in time class of solutions

In this section we will develop solutions of problem(2.1)of the form

whereψis an hypergeometric type function. Indeed,substituting in(3.1)we obtain

whereKis a real constant.

The functionψtakes the form of the solution of problem (3.5) developed in Section 2.2. Next,to develop a solution?of the problem

we may proceed as in the previous sections by assuming either

or

In both cases,we get

wherefHis the solution of the homogeneous problem(3.5)andCis a suitable constant.Indeed,for the first choice we get the following system

ForK=0,we get analogous problem as(3.3)by replacingKbyand thus the solutions may be obtained from Section 2.1. Next,whenevera particular solution is

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

4 Discrete two-dimensional nonlinear EPD system

The object of this section is to explain the discretization scheme proposed to transform problem (1.1)-(1.3) into a discrete quasi-linear one. Letand consider a time stepl=?tand a space oneNext,denote forandj,m∈{0,...,J+1}

so that the cube[L0,L1]×[L0,L1]is subdivided into cubesCj,m=[xj,xj+1]×[ym,ym+1]. For a functionzdefined on the cube[L0,L1]×[L0,L1],we denote by smallthe net functionz(xj,ym,tk)and capitalthe numerical approximation.Consider next the discrete finite difference operators

whereis a barycentric parameter of calibration. Denote alsoG(u,v)=|u|p?1vandH(u)=|v|p?1u.

The discretization of the first equation of problem(1.1)will be

Denote next

We have

or equivalently,

This may be written in a matrix-vector form

whereA,Θ and Λ are the matrices given by

and for 1≤j,m≤J,

Denote next

and

We get

Similarly,for the second equation we get

5 Solvability of the discrete problem

In [40], 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 matrixXhorizontally which leads to an equivalent linear operator characterized by a fringe-tridiagonal matrix. Standard computations have been applied to prove the invertibility of such an operator. Here, we do not apply the same computations as in [40], but we develop different arguments. The first main result is stated as follows.

Theorem 5.1.The system(4.1)-(4.2)is uniquely solvable whenever U0and U1are known.

Proof.It reposes on the inverse of Lyapunov-Sylvester operators.Consider the endomorphism Φ defined by

where for matricesAandB,LA,Bis the Lyapunov-Sylvester operator defined by

and where LA≡LA,A. To prove Theorem 5.1, it suffices to show thatkerΦ is reduced to 0. Indeed,wheneverl=o(h)andl,h→0,we get

Next,whenever Φ(X,Y)=0,we get

Which means thatX=Y=0.

Next,we apply the following result.

Lemma 5.1.Let E be a finite dimensional () vector space and(Φn)n be a sequence of endomorphisms converging uniformly to an invertible endomorphismΦ. Then, there exists n0such that,for any n≥n0,the endomorphismΦn is invertible.

Observing that the operator Φ obtained above is invertible,we get that Φl,his invertible forl,hsmall enough.

6 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 for the discrete equations raised in Section 2, we get the following truncation principal part for the first equation in system(1.1)

and for the second equation,we get

Lemma 6.1.The discrete scheme is consistent with order O(l2+h2).

Now,we will examine the stability of the scheme. We will apply the same method as in[39,40,50,68]based on the Lyapunov criterion of stability. A linear system

is stable in the sense of Lyapunov iff for any bounded initial valueu0, the solutionunramains bounded for alln≥0.

Lemma 6.2.Pn: The solution(Un,Vn)is bounded independently of n whenever the initial solution(U0,V0)is bounded.

Before going on proving this result, we stress on the fact that contrarily to previous studies such as[39,50,68]and[40],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

where Φl,his the Lyapunov-Sylvester operator defined in(5.1) and where the linear operatorF=(F1,F2)(which is a mixture of Lyapunov-Sylvester operators)is defined by

and

whereG(X,Y,Z,T)is the matrix with coefficients

andH(X,Y,Z,T)is the matrix with coefficients

We remark immediately that

and

Consequently,

and

Recall also that

Observe next that forl,hsmall enough

So,fornlarge enough,we obtain

On the other hand,we have

Next,observe that

and consequently,by denoting

we get

Consequently,

Next,by choosing the time stepland the space stephso that 4σCα＜1,we obtain

Finally,the lemma follows by recurrence onnby observing for the choice above that

and the continuity ofF.

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[79].

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

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

7 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 findings in Section 5 by measuring the closeness of the numerical solution obtained by using Lyapunov-Sylvester method to the exact one. The error is evaluated via anL2matrix norm

for a matrixX=(Xij)∈MJ+2C.Denoteunthe net functionu(x,y,tn)andUnthe numerical solution. We propose to compute the discrete error

on the grid(xi,yj),0≤i,j≤J+1 and the relative error between the exact solution and the numerical one as

on the same grid. To do this we fix the parameters

and denote also

The second aim is to show that effectively the system(1.1) may be considered as a perturbation of the classical EPD equation. Indeed,consider the inhomogeneous system

whereG1andG2are explicited respectively by

and

The exact solution of such a system is given by

We fix the domain bounds as follows:

Probelm(7.3)is clearly a perturbation of(1.1). Foru≡vandp=q,we get the inhomogeneous EDP equation

To fulful the assumptionl=o(h) we assume thatTable 1 below 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 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 parameterJ. 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.ForJ=24 for example,a 4-times faster algorithm is obtained. ForJ=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 forJ=999.

Table 1: Error estimates and running time

8 Conclusion

In this work,firstly we developed some exact solutions of an EPD system based on special functions proving especially the genericity of the hypergeometric functions in the development of these solutions. Next, computational method has been developed for numerical solutions of Poisson-Darboux-Euler system in 2-D case based on standard two-dimensional finite difference scheme. The method has yielded non standard algebraic systems generated by general Lyapunov-Sylvester operators where the invertibility needs more effort than classical systems. Fast and accurate algorithms have been obtained compared with the associated tri-diagonal classical algorithms always applied in such problems.