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

Specifications table

Method details

Scientific mathematical background

Consider the following equation

Eq. (1.1) is a two-dimensional (2D) elliptic PDE (EPDE) in quasilinear form, with solution domain and boundary , and being two real numbers on the positive x- and y-axis respectively, such that the fraction is irrational.

The Dirichlet values on the boundary are

The following assumptions for EPDE (1.1) with respect to the boundary conditions (1.2) are made:

in .

and .

and .

a positive constant, a positive constant.

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 .

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 function-vorticity 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 method-based technique for second order linear EPDEs with irregularities in one and two dimensions.

In the proposed domain , 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.

Fig. 1

Nine-point stencil with unequal mesh sizes in x- and y- directions.

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 . Let and be the mesh sizes in x- and y- directions, respectively, where The coordinate points of the mesh so generated are , where and , for and , such that ( and At each mesh point , the exact solution value is denoted by , and approximate solution value is denoted by . Further, let denote and denote .

For every mesh point and for S = a and b, we defineLet and .

Further, for each , 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 of the irrational domain , the proposed discretization of the EPDE (2) is given bywhere is the local truncation error (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 , the EPDE (2) can be expressed asFrom (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 getFor higher order approximations of the derivatives, let us definewhere 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 getwhereFrom (12.1)–(12.2), it is evident that for =0, that is, forwe haveThus, with the aid of (13.1)–(13.2), (5.26) simplifies toFurther, in order to enhance the accuracy of the designed FDM, letwhere j = 1,2,3,4, are arbitrary parameters to be estimated. With the aid of (9.7)–(9.10), (10.1) and (10.2), in the approximations (15.1) and (15.2), and simplifying, we getwhere

As a result of (16.1)–(16.2), from (5.27), we haveFinally, with the support of (10.1), (10.2), (14) and (17), from (6) and (8), the LTE can be computed asSince the LTE is required to be of , the coefficients of and in (18) should vanish, hence the parameters attain the following valuesand the LTE (18) is reduced to .

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 planeLet , and be the mesh sizes in radial and z-axis directions, respectively. Let > 0 be the mesh ratio aspect.

Replacing the variables by , and employing the approximation (6) to the PDE (19), we getwhere

The computational technique (20) is of . However, it is to be noticed that the coefficients and the function in (20) are not defined at . This means that the scheme (20) is unsuccessful for computation at and . The approximation usually declines in the neighborhood of the point (0,0) as and . We conquer this circumstance by changing partially (20) by using the following approximationswhere etc.

Now, with the assistance of (21)–(23), from (20), we have

The modified scheme (24) is of . This modified scheme is suitable to solve EPDE (19) in r-z plane, where (0,0) is a boundary point. This has been illustrated through Example 4 in this article.

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. subject to the condition (1.2), where are constants. We assume that constant, and hence show that the discretization (6) applied to (25) is convergent.

Applying (6) over the Eq. (25), we obtainfor each internal grid point ), with .

Let for ,

+boundary values, and . Further, with t denoting the transpose and for and , we defineThen, as are varied such that , let us re-write Eq. (26) asfor

(Tri-diagonal matrix),

(Tri-diagonal matrix),

(Block Tri-diagonal matrix).

Here, we have made the following assumptions

It is easy to see that the above two conditions yield and This renders strictly positive diagonal entries and strictly negative off-diagonal entries for the matrix .

Since is assumed to be the analytical solution calculated at each grid point ,It is to be noted that for every pair of such that , .

Now, let us defineThen, we can easily obtain the followingfor suitable , and , where , and Now, taking and , we obtainandThus, Eqs. (31.1)–(33.2) fetch us the followingfor , [ 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.

Let , and

Min on and Max on thenAlso, for and , letand

and where , and are some positive constants.

For sufficiently small , we obtainFurther, it is easy to see that the directed graph of the matrix is strongly connected. Therefore, by Varga [43], we conclude that is an irreducible matrix.

Let the sum of elements in the row of the matrix be denoted by Then, for , we obtain the followingwherewhereFor whereFor whereFinally, for and From the Eqs. (36)–(40), we getfor for and , where .for and , where .

For sufficiently small , we obtain

[for ,

[for and , where ,

[for and , where ],

[for and .

Thus, is monotone for sufficiently small value of . Hence, its inverse exists (Varga [43]). Let , where for and , .

Since

, we obtainNow, we re-write Eq. (35) aswhere

Using Eqs. (41)–(44) and (46) in Eq. (45), for appropriately small value of , we get

Thus, the proposed technique applied over Eq. (25) is convergent.

Method validation

Five standard problems of physical significance are numerically solved in this section. The closed-form 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 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.

(Poisson equation in plane)

The analytical solution is 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. 2a and b depict the plots of the analytical and numerical solutions, respectively, for =31.

Table 1

(Example 1): The MAEs.

N_x=N_y	Proposed FDM	Method [49]	Method [56]
15	1.2098(−08)	2.1503(−08)	1.0607(−03)
CPU time (sec)	(0.1518)
31	7.7217(−10)	1.5357(−09)	2.6858(−04)
CPU time (sec)	(0.8118)
63	4.8001(−11)	9.6354(−11)	6.8253(−05)
CPU time (sec)	(6.4516)

Fig. 2

a: The graph of analytical solution with N = N = 31. Fig. 2b: The graph of numerical solution with .

(Burgers' equation)For Eq. (49), the closed-form solution is . The MAEs for the solution are presented in Table 2. The plots of the closed-form and numerical solutions are depicted in Figs. 3a and b, respectively, for and =31.

Table 2

(Example 2): The MAEs.

N_x=N_y	Proposed FDM			Method [49]			Method [56]
	ε=0.1	ε=0.01	ε=0.001	ε=0.1	ε=0.01	ε=0.001	ε=0.1	ε=0.01,0.001
15 CPU time (sec)	2.3721(−07) (0.0795)	5.3311(−06) (1.1876)	2.2039(−05) (1.2340)	7.6613(−07)	1.2494(−05)	3.9373(−05)	1.6176(−03)	Unstable
31 CPU time (sec)	1.4791(−08) (0.3136)	3.3019(−07) (2.2133)	1.3712(−06) (2.3011)	4.2417(−08)	7.7677(−07)	2.2758(−06)	3.8695(−04)	Unstable
63 CPU time (sec)	9.2983(−10) (2.5563)	2.0867(−08) (8.2431)	8.6777(−08) (9.1661)	2.6162(−09)	4.7816(−08)	1.4042(−07)	9.4937(−05)	Unstable

Fig. 3

a: The graph of closed-form solution with and N = N = 31. Fig. 3b: The graph of numerical solution with and N = N = 31.

(NS equations in cartesian coordinates)

The analytical solutions are For varying values of R, the MAEs are reported in Table 3. The analytical and numerical solutions are plotted in Figs. 4a-d for =31 and

Table 3

(Example 3): The MAEs.

N_x=N_y	Proposed FDM			Method [49]			Method [56]
	R_e=10	R_e=10^2	R_e=10^4	R_e=10	R_e=10^2	R_e=10^4	R_e=10	R_e=10^2,10^4
15 u	2.3759(−06)	2.2705(−05)	1.6686(−03)	3.0563(−06)	7.1474(−05)	4.6314(−03)	1.2141(−03)	Unstable
v	6.5170(−07)	7.1474(−06)	1.4175(−03)	2.2897(−06)	2.3107(−05)	3.0187(−03)	6.5296(−04)
CPU time (sec)	(2.2386)	(2.8835)	(3.2231)
31 u	1.5405(−07)	1.3628(−06)	1.1188(−04)	1.9065(−07)	4.3910(−06)	2.8929(−04)	2.9334(−04)	Unstable
v	4.0903(−08)	4.3910(−07)	8.9654(−05)	1.4276(−07)	1.4416(−06)	1.8831(−04)	1.5696(−04)
CPU time (sec)	(6.1298)	(7.4762)	(8.2342)
63 u	9.7192(−09)	8.5442(−08)	7.0458(−06)	1.1845(−08)	1.7156(−07)	1.8045(−05)	7.2146(−05)	Unstable
v	2.5694(−09)	2.7756(−08)	5.6432(−06)	8.9755(−09)	9.0650(−08)	1.1785(−05)	3.8553(−05)
CPU time (sec)	(25.3462)	(28.1260)	(31.7743)

Fig. 4

a: The graph of analytical solution, u, with R=104 and N = N = 31. Fig. 4b: The graph of numerical solution, u, with R = 104 and N = N = 31. Fig. 4c: The graph of analytical solution, v, with R = 104 and  = N = 31. Fig. 4d: The graph of numerical solution, v, with R = 104 and N = N = 31.

(NS equations in cylindrical polar coordinates)

The analytical solutions are The MAEs for the solutions u and v are given in Table 4, taking various values of the Reynold constant R. The analytical and numerical plots for u and v are graphed in Figs. 5a-d for =31 and

Table 4

(Example 4): The MAEs.

N_x=N_y	Proposed FDM			Method [49]			Method [56]
	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
15 u	1.8317(−06)	2.0204(−05)	1.7178(−03)	7.6512(−06)	4.8900(−05)	3.0319(−03)	3.7991(−04)	Unstable
v	1.1949(−06)	1.6886(−05)	1.4376(−03)	4.7072(−06)	3.0462(−05)	2.9142(−03)	1.4492(−03)
CPU time (sec)	(3.1818)	(4.4332)	(5.1043)
31 u	1.1859(−07)	1.2679(−06)	1.3303(−04)	4.2316(−07)	2.9424(−06)	1.9105(−04)	9.1578(−05)	Unstable
v	8.2138(−08)	1.1302(−06)	8.9855(−05)	2.9392(−07)	1.9064(−06)	1.8073(−04)	3.5201(−04)
CPU time (sec)	(8.4562)	(12.2167)	(15.1543)
63 u	7.4638(−09)	7.9503(−08)	8.3304(−06)	2.6061(−08)	1.8196(−07)	1.1954(−05)	2.2508(−05)	Unstable
v	5.1803(−09)	7.0859(−08)	5.6823(−06)	1.8333(−08)	1.1844(−07)	1.1244(−05)	8.6797(−05)
CPU time (sec)	(30.3634)	(41.4245)	(48.3198)

Fig. 5

a: The graph of analytical solution, u, with R = 103 and N = N = 31. Fig. 5. b: The graph of numerical solution, u, with R = 103 and N = N = 31. Fig. 5. c: The graph of analytical solution, v, with R = 103 and N = N = 31. Fig. 5. d: The graph of numerical solution, v, with R = 103 and N = N = 31.

Clearly, this is a quasi-linear equation. The closed form solution is The MAEs are given in Table 5. The closed-form and numerical solutions are graphed in Figs. 6a and b, respectively, for =31.

Table 5

(Example 5): The MAEs.

N_x=N_y	Proposed Method	Method [49]	Method [56]
15	1.3913(−05)	3.9473(−05)	3.5901(−04)
CPU time (sec)	(0.1818)
31	8.6978(−07)	2.2858(−06)	9.1000(−05)
CPU time in (sec)	(0.8879)
63	5.4196(−08)	1.4142(−07)	2.2652(−05)
CPU time in (sec)	(7.8873)

Fig. 6

a: The graph of closed-form solution with N = N = 31. Fig. 6b: The graph of numerical solution with N = N = 31.

The order of the proposed method is . When or, the order of accuracy transforms as either , or of . In the given domain , The order of convergence can be computed using the formulawhere and refer to the MAEs of grid sizes and respectively, with . We have chosen, and The order of convergence for all problems are computed in Table 6.

Table 6

(The rate of convergence): .

Example	Parameters, if any		Order of convergence
1			4.00
2	ε=0.1		3.99
	ε=0.01		3.98
	ε=0.001		3.98
3	R_e=10	foru	3.98
		forv	3.99
	R_e=10^2	foru	3.99
		forv	3.98
	R_e=10^4	foru	3.98
		forv	3.98
4	R_e=10	foru	3.98
		forv	3.98
	R_e=10^2	foru	3.99
		forv	3.99
	R_e=10^3	foru	3.99
		forv	3.98
5			4.00

CRediT authorship contribution statement

R.K. Mohanty: Conceptualization, Funding acquisition, Investigation, Methodology, Project administration, Supervision, Writing – review & editing. Nikita Setia: Conceptualization, Data curation, Formal analysis, Funding acquisition, Investigation, Resources, Software, Validation, Visualization, Writing – original draft, Writing – review & editing. Gunjan Khurana: Conceptualization, Formal analysis, Funding acquisition, Resources, Software, Validation, Visualization, Writing – review & editing. Geetan Manchanda: Conceptualization, Formal analysis, Funding acquisition, Validation, Writing – review & editing.

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.

Subject area	Mathematics
More specific subject area	Partial Differential Equations
Name of your method	Finite Difference Method
Name and reference of original method	• M.M. Chawla, A fourth-order tridiagonal finite difference method for general non-linear two-point boundary value problems with mixed boundary conditions, IMA Journal of Applied Mathematics. 21(1) (1978) 83–93. https://doi.org/10.1093/imamat/21.1.83
	• M.K. Jain, R.K. Jain and R.K. Mohanty, A fourth order difference method for elliptic equations with non-linear first derivative terms, Numerical Methods for Partial Differential Equation., 5 (1989) 87-95. https://doi.org/10.1002/num.1690050203
Resource availability	MATLAB