High precision compact numerical approximation in exponential form for the system of 2D quasilinear elliptic BVPs on a discrete irrational region

This article presents a new approximation of order four in exponential form for two-dimensional (2D) quasilinear partial differential equation (PDE) of elliptic form with solution domain being irrational. It is further extended for application to a system of quasilinear elliptic PDEs with Dirichlet boundary conditions (DBCs). The main highlights of the method framed in this article are as under:• It uses a 9-point stencil with unequal mesh to approach the solution. The error analysis is discussed to authenticate the order of convergence of the proposed numerical approximation.• Various validating problems, for instance the Burgers’ equation, Poisson equation in cylindrical coordinates, Navier-Stokes (NS) equations in rectangular and cylindrical coordinates are solved using the proposed techniques to depict their stability. The proposed approximation produces solution free of oscillations for large values of Reynolds Number in the vicinity of a singularity.• The results of the proposed method are superior in comparison to the existing methods of [49] and [56].

The condition (I) is required to satisfy the elliptic condition, the condition (II) is required for a valid local truncation error. Conditions (III)-(IV) are the sufficient conditions for the convergence of the numerical solution (Jain et al. [1] ). Further, we presume that the solution of the boundary value problem (BVP) (1.1) -(1.2) is unique and exists in the domain h .
In the past, many high order compact numerical approximations for the nonlinear EPDEs have aroused renewed interest amongst the researchers and different type of techniques have been developed. A compact finite difference method (FDM) is restricted to cells surrounding any given grid and does not extend further, so it is convenient for computation since no special techniques are required for nodes near the boundary (Chawla [64] ). Many problems of physical significance like the diffusion equation with first derivative convection terms, viscous Burgers', and viscous Navier-Stokes (NS) equations in rectangular, cylindrical and spherical coordinates, are modeled by the EPDEs. Most of the numerical methods to solve EPDEs are iterative methods. We have attempted to solve numerically the coupled nonlinear equations like NS equations which show up in diverse areas of physics and mathematics. Li et al. [2] have expressed the NS equations in stream function-vorticity form and solved the system using nonlinear SOR iteration method. Yanwen et al. [3] have discretized the convection terms with high order upwind compact difference approximation and solved the incompressible NS equations in a rectangular domain. In 2009, Shah and Yuan [4] have proposed upwind compact approximations, of order three and five, for the solution of NS equations with artificial compressibility. These methods involved a flux-difference splitting technique. On the basis of the projection method, a high order FDM has been developed by Tian et al. [5] . This method is applied over the 2D incompressible NS equations in primitive variables. Using a fine grid mesh and a tri-diagonal solver, a new compact fourth order approximation has been employed by Erturk and Gökcöl [7] to solve the NS equations. Liu and Wang [8] have proposed a new high order FDM in terms of mean-vorticity formulation. The solutions obtained through this FDM were analyzed for the primitive equations of significant oceanic and atmospheric flow. A compact MAC FDM of high accuracy involving a staggered mesh has been developed by Ito and Qiao [9] for the solution of Stokes equations with Dirichlet boundary conditions (DBCs) on the velocity. A fourth-order compact nine-point 2D stencil has been formulated in [10][11][12] for the NS equations in steady stream functionvorticity form. A spatial approximation based on orthogonal spline collocation has been proposed by Fairweather et al. [13] for the 2D NS equations in fully coupled stream function-vorticity form. Various numerical approximations for 2D EPDEs with significant linear first derivative terms have been studied in [14][15][16] . In practical, numerical approximations for nonlinear EPDEs are of extreme interest in applied physics and applied mathematics. Different approximations have been developed for solving EPDEs, in particular finite element methods (FEMs), finite volume methods (FVMs) and FDMs [17][18][19][20][21][22] . Ananthakrishnaiah et al. [23] have analyzed a fourth order approximation for 2D mildly nonlinear EPDEs. Different types of fourth order approximations for the solution for 2D quasilinear elliptic BVPs and their normal derivatives have been studied in [24][25][26] . Zhang [27] and Saldanha [28] have derived compact numerical approximations for 2D elliptic BVPs and discussed convergence and performances of iterative methods. Arabshahi and Dehghan [29] have adopted preconditioning techniques to solve linear EPDEs.
Let us now document some of the related research works done in the past one decade. Several meshless techniques for solving the various counterparts of Eq. (1.1) have been proposed by various authors. In 2013, Aziz et al. [59] have proposed two new efficient collocation methods to obtain the solution of a linear EPDE with a variety of boundary conditions. Their methods were based on Haar and Legendre Wavelets. This work has been further extended, in the year 2017, for applicability to nonlinear EPDEs with improved efficiency by Aziz and Islam [61] . In 2015, Oberman and Zwiers [30] have introduced a combination of adaptive and monotone FDMs to non-linear elliptic and parabolic PDEs with free boundaries. Subsequently, in 2017, a higher order method using the similar approach has been developed by Hamfeldt and Salvador [31] . A meshless technique for solving a 3D linear EPDE has been designed by Gavete et al. [32] in 2016, and for a 2D non-linear EPDE by Gavete et al. [33] in 2017. In 2019, Oruç [34] have solved the Poisson equation with mixed boundary conditions and irregular domain using a generalized FDM. In the same year, two new numerical techniques for two and three dimensional EPDEs with regular interfaces have been proposed by Haider et al. [63] . These techniques were based on meshless collocation and Haar wavelet collocation. In 2021, Milewski [35] have proposed a meshless technique for solving a 2D linear EPDE with Dirichlet and Neumann conditions. Apart from the meshless approaches, Mohanty and Setia [36] have designed a highly accurate compact half-step discretization for a system of general non-linear 2D EPDEs in the year 2013. A compact FDM of high order accuracy for a 2D Poisson equation has been given by Zhai et al. [37] in the following year. Further, in 2015, Papanikos and Gousidou-Koutita [39] have studied and compared the numerical results obtained upon application of the central difference scheme and an FEM algorithm over a general linear second order EPDE. The same year, Islam et al. [60] reported a new numerical technique based on Haar Wavelet collocation and a meshless method involving different types of radial basis functions for a 2D Poisson equation with a two-point non-local boundary condition and an integral boundary condition. In 2016, Mittal et al. [41] have given a class of FDMs for solving 1D and 2D elliptic and parabolic equations. The accuracy of this scheme is atleast two. It is successfully applicable to time-dependent NS equations in cylindrical polar coordinates. RBF generated FDMs, successfully applicable to the NS equations in cylindrical polar coordinates, have also been proposed by Bayona et al. [40] in the year 2017. This is. Afterwards, in 2018, Mittal and Ray [42] have proposed a new second order FDM for 2D and 3D non-linear EPDEs. The EPDEs involved in [41] and [42] have discontinuous coefficients and singular source terms. Raeli et al. [45] have developed an FDM for the variable coefficient Poisson equation. Aziz et al . [62] have designed a new collocation method for three dimensional nonlinear EPDEs with DBCs. This method is based on Haar Wavelets. A second order FDM for a system of EPDEs in a square domain has been given by Pandey and Jaboob [48] . Zhang et al. [6] have proposed a highly accurate discontinuous Galerkin method for solving the 2D NS equations in conservative form on arbitrary grids. In 2020, Li and Zhang [50] have reduced the classical continuous FEM to a high order FDM applied to a 2D elliptic equation with smooth coefficients on a rectangular domain. A new bi-cubic spline collocation method of fourth order accuracy has been designed by Singh and Singh [57] for a linear EPDE with DBCs. Most recently, new high accuracy approximations in exponential form for solving two-point BVPs on a variable mesh, and 2D nonlinear EPDEs on an unequal uniform mesh have been reported in [51][52][53][54][55] . Mohanty and Kumar [56] , and Priyadarshini and Mohanty [ 47 , 49 ] have proposed new high accuracy numerical algorithms for 2D quasilinear EPDEs using unequal mesh. Singh et al. [58] have developed a collocation methodbased technique for second order linear EPDEs with irregularities in one and two dimensions.
In the proposed domain h , the grid sizes in both x -and y-directions cannot be equal. It is necessary to develop an appropriate stable compact highly accurate computational method for the EPDE (1.1) using an unequal mesh. In this piece of work, we design a new compact numerical method in exponential form of order four on the given irrational domain for the solution of quasilinear elliptic PDE (1.1) . We use nine-point stencil in the formulation of the method ( Fig. 1 ). The outline of the rest of the paper is as follows: In Section 1.2 , we frame the FDM in the exponential form on an irrational domain using unequal grid sizes in x-and y-directions. Further, in Section 1.3 , we deduce the compact numerical approximation. The obtained discretization is extended for applicability to a system of 2D EPDEs with DBCs. A proof of convergence of the numerical method is given in Section 1.4 . The error is proved to be of fourth-order accuracy. Section 2 comprises of some numerical illustrations for method validation. Different kinds of benchmark elliptic BVPs are included for numerical illustrations to demonstrate the fourth order accuracy and application of the method discussed.

Discretization procedure
The 2D nonlinear EPDE is given by The DBCs for (2) are given by the Eq. (1.2) . The discretization is implemented on the solution domain h . Let h x > 0 and h y > 0 be the mesh sizes in x -and y-directions, respectively, where h x = h y . The coordinate points of the mesh so generated are ( x l , y m ) , where x l = lh x and y m = mh y , for l = 0 , 1 , . . . , (N x + 1) and m = 0 , 1 , . . . , (N y + 1) , such that ( N x + 1) h x = x a and (N y + 1) h y = y b . At each mesh point ( x l , y m ) , the exact solution value is denoted by U l,m , and approximate solution value is denoted by u l,m . Further, let a l,m denote a ( x l , y m ) and b l,m denote b( x l , y m ) .
For every mesh point ( x l , y m ) and for S = a and b , we define Let 10 6 a 00 and Further, for each ( x l , y m ) , we denote Using (4) , we generate an approximation using unequal grid sizes in x-and y-directions by defining the following Then, at every mesh point ( x l , y m ) of the irrational domain h , the proposed discretization of the EPDE (2) is given by where ˆ is the local truncation err or (LTE).
Incorporating the prescribed Dirichlet values from the boundary, the approximation (6) may be expressed as a block tri-diagonal matrix, that may be iteratively dealt with for a solution by choosing an appropriate block-iterative method ( [ 38 , 43 , 44 , 46 ]).

Derivation of the algorithms in the exponential form
At any mesh point ( x l , y m ) , the EPDE (2) can be expressed as From (7) , it is simple to observe that Simplifying the approximations (5.1) -(5.19) , we obtain Employing the approximations (9.2) , (9.3) , (9.5) and (9.6) in (5.20) and (5.21) , and simplifying, we get For higher order approximations of the derivatives, let us define where α 1 j , β 1 j , j = 1,2,3,4, are the parameters to be estimated. Using the approximations (9.7) -(9.10) , (10.1) and (10.2) , in the expressions (11.1) and (11.2) , and simplifying, we get where Further, in order to enhance the accuracy of the designed FDM, let Finally, with the support of (10.1) , (10.2) , (14) and (17) , from (6) and (8) , the LTE can be computed as   a 00 , a 01 , b 01 ) , The technique discussed in [49] can be used to obtain numerical approximation of order four for the quasilinear EPDE (1.1) . Further, the proposed method (6) for the single equation can be stretched out to the system of quasilinear EPDE in vector-matrix form, as discussed in [36] .

Note: Application to elliptic BVP in r-z plane
Let us now discuss the implementation of the proposed discretization to the following EPDE in r-z plane ) , and h r > 0 , h z > 0 be the mesh sizes in radial and zaxis directions, respectively. Let γ = h r /h z > 0 be the mesh ratio aspect.
However, it is to be noticed that the coefficients 1 r and the function g( r, z ) in (20) are not defined at ( r, z ) = ( 0 , 0 ) . This means that the scheme (20) is unsuccessful for computation at l = 1 and m = 1 . The approximation usually declines in the neighborhood of the point (0,0) as h r → 0 and h z → 0 . We conquer this circumstance by changing partially (20) by using the following approximations where g l,m = g( r l , z m ) etc.
The modified scheme (24)  In a similar manner, we can derive modified method for nonlinear elliptic BVPs in r-z plane.

Error analysis
In this subsection, the error analysis of the proposed discretization is discussed, when applied over the constant coefficient counterpart of the Eq. (1.1) viz.
Then, as ( l, m ) are varied such that 1 ≤ l, m ≤ N − 1 , let us re-write Eq. (26) as Here, we have made the following assumptions It is easy to see that the above two conditions yield a > 0 and b > 0 . This renders strictly positive diagonal entries and strictly negative off-diagonal entries for the matrix D .
Since U l,m is assumed to be the analytical solution calculated at each grid point ( x l , y m ) , It is to be noted that for every pair of ( l, m ) such that 1 ≤ l, m ≤ N, Then, we can easily obtain the following l,m , where Q = ρ, α and β. Now, taking Q = α and β, we obtain being the block tri-diagonal matrix with the following entries Ignoring the existence of the round-off errors, relations (27) , (30) and (34) give us the following equation called the error equation.
For sufficiently small h x , we obtain Further, it is easy to see that the directed graph of the matrix D + P is strongly connected. Therefore, by Varga [43] , we conclude that D + P is an irreducible matrix. Let the sum of elements in the k th row of the matrix D + P be denoted by S k . Then, for k = 1 , N − 1 , we obtain the following Finally, for 2 ≤ q ≤ N − 2 and 2 ≤ r ≤ N − 2 From the Eqs. (36) -(40) , we get For sufficiently small h x , we obtain S k > 5 h 2 Thus, D + P is monotone for sufficiently small value of h x . Hence, its inverse exists (Varga [43] ). Let Now, we re-write Eq. (35) as where Using Eqs. (41) - (44) and (46) in Eq. (45) , for appropriately small value of h x , we get Thus, the proposed technique applied over Eq. (25) is

Method validation
Five standard problems of physical significance are numerically solved in this section. The closedform solutions of these problems are available. The functions on the right-hand side and the DBCs can be attained from the analytical solution as a test procedure. Both nonlinear and linear difference equations in block forms are solved with the aid of suitable iterative methods ( [ 38 , 43 , 44 , 46 ]).
During computation, the iterations are terminated when | u ( k +1 ) − u (k ) | ≤ 10 −12 is achieved. The initial solution is chosen as vector 0 for all iterative methods and the Maximum Absolute Errors (MAEs) are computed. All the results were computed using MATLAB codes.

Example 1. (Poisson equation in r − z plane)
The analytical solution is u ( r, z ) = coshr. coshz. By the help of the proposed method (6) and using the technique discussed in [24] , we solve the Eq. (48) . The MAEs are presented in Table 1 . Figs. 2 a and b depict the plots of the analytical and numerical solutions, respectively, for N x = N y = 31.

Example 2. (Burgers' equation)
For Eq. (49) , the closed-form solution is u = e x . sin ( π y ) . The MAEs for the solution u are presented in Table 2 . The plots of the closed-form and numerical solutions are depicted in Figs. 3 a and b, respectively, for ε = 0 . 01 and N x = N y = 31.

Example 3. (NS equations in cartesian coordinates)
1       The analytical solutions are u = sin ( π x ) sin ( π y ) , v = cos ( π x ) cos ( π y ) . For varying values of R e , the MAEs are reported in Table 3 . 1 R e u rr + 1 r u r + u zz − 1 1 R e v rr + 1 r v r + v zz = u v r + vv z + g ( r, z ) , 0 < r < 1 , 0 < z < √ 2 . (51. 2) The analytical solutions are u = r 3 sinhz , v = − 4 r 2 coshz. The MAEs for the solutions u and v are given in Table 4 , taking various values of the Reynold constant R e . The analytical and numerical plots for u and v are graphed in Figs. 5 a-d for N x = N y = 31 and R e = 10 3 .
Clearly, this is a quasi-linear equation. The closed form solution is u ( x, y ) = e x siny . The MAEs are given in where e h x 1 and e h x 2 refer to the MAEs of grid sizes h x 1 =  Table 6 . R e = 10 R e = 10 2 R e = 10 3 R e = 10 R e = 10 2 R e = 10 3 R e = 10 R e = 10 2 , 10 3

Declaration of Competing Interest
The authors declare that they have no known competing financial interests or personal relationships that could have appeared to influence the work reported in this paper.

Data Availability
No data was used for the research described in the article.