Skip to main content

Theory and Modern Applications

Operator compact exponential approximation for the solution of the system of 2D second order quasilinear elliptic partial differential equations

Abstract

In this paper, we suggest a new exponential implicit method based on full step discretization of order four for the solution of quasilinear elliptic partial differential equation of the form \(A ( x,y,z ) z_{xx} +C ( x,y,z ) z_{yy} =k ( x,y,z, z_{x}, z_{y} )\), \(0< x,y<1\). In this method a single compact cell consisting of nine nodal points is used. Convergence analysis of the said method is discussed in detail. The developed method is successfully applied to solving problems in polar coordinates. The method for scalar equation is eventually applied to solving the system of quasilinear elliptic equations. To measure the rationality and precision, the method is applied to solving several noteworthy problems and numerical results are provided to exhibit the effectiveness of the method.

1 Introduction

We examine a quasilinear elliptic equation in two space dimensions of the form

$$ A(x,y,z)z_{xx} + C(x,y,z)z_{yy} = k(x, y,z, z_{x}, z_{y}), $$
(1.1)

where \((x,y)\in \varGamma = (0,1) \times (0,1)\).

$$ z(x, y) = z_{0}(x, y),\quad (x,y) \in \partial \varGamma. $$
(1.2)

The partial differential equations (PDEs) of the type (1.1) with non-constant coefficients illustrate several real world problems of phenomenological importance like Poisson’s equation, convection–diffusion equation, Burgers’ equation and the nonlinear steady-state Navier–Stokes (NS) equations of motion.

We presuppose the following about the boundary value problem (1.1)–(1.2):

  1. 1.

    \(A.C>0\) in Γ,

  2. 2.

    \(z(x,y)\in C^{6}\),

  3. 3.

    \(A, C\in C^{4}\),

  4. 4.

    k is continuous,

  5. 5.

    \(\frac{\partial k}{\partial z} \geq0\),

  6. 6.

    \(\vert \frac{\partial k}{\partial z_{x}} \vert \leq H\), \(\vert \frac{\partial k}{\partial z_{y}} \vert \leq I\),

where H and Iare finite positive real bounds and \(C ^{p}\) is the class of functions with smooth partial derivatives up to order p (Jain et al. [1]).

The ellipticity condition of Eq. (1.1) is ensured by (1). Assumptions (2)–(4) are required to enable the Taylor series expansions, and (5)–(6) are the adequate conditions for the existence and uniqueness of the solution of the boundary value problem (1.1)–(1.2) (Jain et al. [2]).

Elliptic equations, a type of partial differential equations illustrate behavior that remains static with time. In other words processes which are in equilibrium, like heat flow or fluid flow through a medium without any accumulations. The elliptic equation satisfies a differential equation within a domain along with the values near the boundary of the region (boundary values), representing the effect from outside the domain. These conditions fall in two categories, one representing fixed temperature distribution at boundary points (Dirichlet problem) another, where heat is added or removed across the boundary in such a manner so that constant temperature is maintained throughout (Neumann problem).

The elliptic equations are used to model many natural phenomena like heat dissipation in a metal sheet. They are used in aircraft design as well as in weather prediction (weather, flow over wing, turbulence etc.). Many numerical strategies to solve a specific PDE are iterative methods based on finite differencing whereby a mesh is generated to represent a physical domain, the points on the mesh or grid are initialized usually with an approximate solution and then repeatedly updated to obtain an increasingly accurate solution. This process may be repeated for a fixed number of iterations or until the solution has reached the desired level of accuracy. Each iteration stores the newly calculated values in another array and swaps the arrays at the end of the iteration.

The Laplace equation, \(z _{xx} + z _{yy} = 0\), is one of the most fundamental 2D elliptic equation which represents the steady-state condition. It arises in electrostatics, gravitation, steady-state flow of inviscid fluids and steady-state heat conduction. We have numerically solved substantially important equations like Navier–Stokes equations which are beneficial in explaining multitude of phenomena arising in science and engineering-modeling weather, ocean currents, water flowing in a pipe and air flowing around a wing. The Navier–Stokes equations find their way in a host of applications viz. creating blueprints of aircraft and cars, flow of blood in animals, designing power stations, analyzing soil, air and water pollution and a plethora of real world phenomena. These equations govern the motion of fluids and can be seen as Newton’s second law of motion for fluids. These are always solved together with the continuity equation. The Navier–Stokes equations symbolize the conservation of momentum while the continuity equation describes the conservation of mass. They form the core of fluid dynamics. Due to their complex behavior, these equations allow for a few restricted analytical solutions. In 1995, Li, Tang, Fornberg [3] discussed a compact finite difference method for the equilibrium state of non-compressible Navier–Stokes equations.

In the year 2006, Erturk, Gökcöl [4] developed a fourth order compact method to solve Navier–Stokes equations with high Reynolds numbers. Then, in the year 2008, Liu, Wang proposed a fourth order numerical scheme for the primitive equations formulated in mean vortices [5]. In the same year Ito and Qiao [6] discussed compact MAC finite difference scheme of high order for the Stokes equations.

We have also numerically solved the Poisson equation, which arises in electrostatics and elasticity theory. The solution of a Poisson equation describes the steady state of a system. For example, the stabilized temperature of a steel rod with one end held in your hand and the other end in the air is a solution of a certain Poisson equation. To approximate the solution of a Poisson equation numerically, one needs to solve a diagonally dominant linear system.

For linear elliptic equations many numerical schemes have been discussed which date back to the year 1984 (see [7,8,9,10,11]). Solving fully nonlinear elliptic partial differential equations numerically finds strong interest among the research community. The applicability of these equations in innumerable areas of science, such as transportation theory, optimization, fluid dynamics and differential geometry, provides strong impetus to pursue deeper research into this field. Among various numerical methods in the literature for solving fully nonlinear equations; finite difference and finite element type methods are quite popular (see [1, 12,13,14,15]). In 1994, Jain et al. [16] developed fourth order difference method for quasilinear Poisson equation in cylindrical symmetry. The following year Ananthakrishnaiah, Saldanha [17] discussed a fourth order finite difference scheme for two-dimensional nonlinear elliptic partial differential equations. Thereafter, Mohanty et al. [18,19,20,21], Zhang [22], Saldanha [23] discussed order \(h ^{4}\) difference methods for a class of elliptic boundary value problem. Dehghan et al. [24] proposed preconditioning techniques to obtain faster convergence of the higher order methods applied to linear elliptic PDEs. Using the splitting technique, Mohanty et al. [25,26,27,28], Khattar et al. [29] and Singh et al. [30, 31] have proposed high accuracy numerical methods for the solution of nonlinear bi- and triharmonic elliptic boundary value problems.

This paper is devoted to the construction and analysis of the fourth order exponentially fitted discretization of a second order quasilinear elliptic PDE using full step grid points on a uniform mesh with Dirichlet boundary conditions. The exponentially fitted scheme is one of the upcoming classes of robust difference schemes. Such schemes exhibit good convergence and stability. Furthermore, they do not produce spurious oscillations as the previously known finite difference schemes. The supremacy of the exponentially fitted scheme is reflected by the numerical results in terms of maximum absolute errors.

A good first theoretical foundation of the technique of exponential fitting for multistep methods was given by Gautschi [32] and Lyche [33]. Since then, a lot of exponentially fitted linear multistep methods have been constructed; most of them were developed for second order differential equations where the first derivative is absent and applied to solving equations of the Schrödinger type.

This paper is structured in the following manner: In Sect. 2, we formulate the full step fourth order compact discretization scheme for the solution of a nonlinear elliptic equation with non-constant coefficients. Section 3 is dedicated to a detailed derivation of the numerical scheme under consideration. Thereafter, in Sect. 4, we establish the fourth order convergence of the method for a scalar equation under suitable conditions. Further, we extend our method to the system of quasilinear elliptic PDEs with non-constant coefficients, subject to Dirichlet boundary conditions in Sect. 5. In Sect. 6, we discuss our method to solve the elliptic equation in polar coordinates. The difficulties like the deterioration of the solution in the vicinity of singularity which we encountered in the past for obtaining the high accuracy numerical solution for the singular elliptic problems are solved by modifying our method and thus the method becomes valid to compute the singular problem in the entire solution domain. In Sect. 7, we have applied our method to solving nonlinear bi- and triharmonic problems. In Sect. 8, we solve a set of linear and nonlinear elliptic problems of physical importance to present and investigate the precision of the proposed method. The last section is devoted to concluding remarks.

2 Formulation of the numerical algorithm

We first consider the following two-dimensional elliptic PDE:

$$ A(x,y)z_{xx} + C(x,y)z_{yy} = k(x, y,z, z_{x}, z_{y}),\quad 0 < x,y < 1, $$
(2.1)

where \((x,y)\in \varGamma = (0,1) \times (0,1)\), satisfying Dirichlet boundary conditions given by (1.2).

We overlay on the solution domain Γ a square grid with spacing \(h>0\) in both x- and y-directions. Each nodal point is given by (\(x _{p} ,y _{q}\)), where \(x _{p}=ph\) and \(y _{q}=qh\), \(0 \leq p\), \(q \leq N\)+1 and \((N+1)h=1\).

Further, at each grid point (\(x _{p} ,y _{q}\)), let

  1. (a)

    \(Z_{p,q}\) and \(z_{p,q}\) be the exact and approximate values of \(z(x _{p} ,y _{q})\), respectively,

  2. (b)

    \(K_{p,q} = k(x_{p}, y_{q},Z_{p,q}, Z_{x_{p,q}}, Z_{y_{p,q}})\),

  3. (c)

    \(k_{p,q} = k(x_{p}, y_{q},z_{p,q}, z_{x_{p,q}}, z_{y_{p,q}})\),

  4. (d)

    \(S_{lm} = \frac{\partial^{l+m} S}{\partial x^{l} \partial y^{m}}\), for \(l,m=0,1,2,\dots\) and for \(S=A, C,Z\).

At each nodal point (\(x _{p} ,y _{q}\)), the differential equation (2.1) can be written as

$$ A_{00} Z_{20} + C_{00} Z_{02} = K_{p,q}. $$
(2.2)

Let \(\delta_{x} Z_{l} = ( Z_{l+ {\frac{1}{2}}} - Z_{l- {\frac{1}{2}}} )\) and \(\mu_{x} Z_{l} = {\frac{1}{2}} ( Z_{l+ {\frac{1}{2}}} + Z_{l- {\frac{1}{2}}} )\).

For the fourth order discretization of PDE (2.1), we use the following approximations:

$$\begin{aligned} &\overline{Z}_{x_{p,q}} = \frac{1}{2h} [ 2\mu_{x} \delta_{x} ]Z_{p,q}= Z_{x_{p,q}} + \frac{h^{2}}{6}Z_{30} + O\bigl(h^{4}\bigr), \end{aligned}$$
(2.3a)
$$\begin{aligned} &\overline{Z}_{x_{p \pm 1,q}} = \frac{1}{2h} \bigl[ 2\mu_{x} \delta_{x} \pm 2\delta_{x}^{2} \bigr]Z_{p,q} = Z_{x_{p \pm 1,q}} - \frac{h^{2}}{3}Z_{30} \pm O\bigl(h^{3}\bigr), \end{aligned}$$
(2.3b)
$$\begin{aligned} &\overline{Z}_{x_{p,q \pm 1}} = \frac{1}{2h} [ 2\mu_{x} \delta_{x} ]Z_{p,q \pm 1} = Z_{x_{p,q \pm 1}} + \frac{h^{2}}{6}Z_{30} \pm O\bigl(h^{3}\bigr), \end{aligned}$$
(2.3c)
$$\begin{aligned} &\overline{Z}_{y_{p,q}} = \frac{1}{2h} [ 2\mu_{y} \delta_{y} ]Z_{p,q}= Z_{y_{p,q}} + \frac{h^{2}}{6}Z_{03} + O\bigl(h^{4}\bigr), \end{aligned}$$
(2.4a)
$$\begin{aligned} &\overline{Z}_{y_{p \pm 1,q}} = \frac{1}{2h} [ 2\mu_{y} \delta_{y} ]Z_{p \pm 1,q}= Z_{y_{p \pm 1,q}} + \frac{h^{2}}{6}Z_{03} \pm O\bigl(h^{3}\bigr), \end{aligned}$$
(2.4b)
$$\begin{aligned} &\overline{Z}_{y_{p,q \pm 1}} = \frac{1}{2h} \bigl[ 2\mu_{y} \delta_{y} \pm 2\delta_{y}^{2} \bigr]Z_{p,q} = Z_{y_{p,q \pm 1}} - \frac{h^{2}}{3}Z_{03} \pm O\bigl(h^{3}\bigr), \end{aligned}$$
(2.4c)
$$\begin{aligned} &\overline{Z}_{xx_{p,q}} = \frac{1}{h^{2}} \bigl[ \delta_{x}^{2} \bigr]Z_{p,q}= Z_{xx_{p,q}} + \frac{h^{2}}{12}Z_{40} + O \bigl(h^{4}\bigr), \end{aligned}$$
(2.5a)
$$\begin{aligned} &\overline{Z}_{xx_{p,q \pm 1}} = \frac{1}{h^{2}} \bigl[ \delta_{x}^{2} \bigr]Z_{p,q \pm 1}= Z_{xx_{p,q \pm 1}} + \frac{h^{2}}{12}Z_{40} \pm O\bigl(h^{3}\bigr), \end{aligned}$$
(2.5b)
$$\begin{aligned} &\overline{Z}_{yy_{p,q}} = \frac{1}{h^{2}} \bigl[ \delta_{y}^{2} \bigr]Z_{p,q}= Z_{yy_{p,q}} + \frac{h^{2}}{12}Z_{04} + O \bigl(h^{4}\bigr), \end{aligned}$$
(2.6a)
$$\begin{aligned} &\overline{Z}_{yy_{p \pm 1,q}} = \frac{1}{h^{2}} \bigl[ \delta_{y}^{2} \bigr]Z_{p \pm 1,q}= Z_{yy_{p \pm 1,q}} + \frac{h^{2}}{12}Z_{04} \pm O\bigl(h^{3}\bigr). \end{aligned}$$
(2.6b)

Define

$$\begin{aligned} &\overline{K}_{p \pm 1,q} = k(x_{p \pm 1},y_{q},Z_{p \pm 1,q}, \overline{Z}_{x_{p \pm 1,q}},\overline{Z}_{y_{p \pm 1,q}}), \end{aligned}$$
(2.7a)
$$\begin{aligned} &\overline{K}_{p,q\pm 1} = k(x_{p},y_{q \pm 1},Z_{p,q \pm 1}, \overline{Z}_{x_{p,q\pm 1}},\overline{Z}_{y_{p,q \pm1}}). \end{aligned}$$
(2.7b)

Let

$$\begin{aligned} &\begin{aligned}[b] \overline{\overline{Z}}_{x_{p,q}} ={}& \overline{Z}_{x_{p,q}} - \frac{h}{12A_{00}}(\overline{K}_{p + 1,q} - \overline{K}_{p - 1,q}) + \frac{hC_{00}}{12A_{00}}(\overline{Z}_{yy_{p + 1,q}} - \overline{Z}_{yy_{p - 1,q}})\\ &{} + \frac{h^{2}}{6A_{00}}(A_{10}\overline{Z}_{xx_{p,q}} + C_{10} \overline{Z}_{yy_{p,q}}), \end{aligned} \end{aligned}$$
(2.8a)
$$\begin{aligned} & \begin{aligned}[b]\overline{\overline{Z}}_{y_{p,q}} ={}& \overline{Z}_{y_{p,q}} - \frac{h}{12C_{00}}(\overline{K}_{p,q + 1} - \overline{K}_{p,q - 1}) + \frac{hA_{00}}{12C_{00}}(\overline{Z}_{xx_{p,q + 1}} - \overline{Z}_{xx_{p,q - 1}})\\ &{} + \frac{h^{2}}{6C_{00}}(A_{01}\overline{Z}_{xx_{p,q}} + C_{01} \overline{Z}_{yy_{p,q}}), \end{aligned} \end{aligned}$$
(2.8b)

and let

$$ \overline{\overline{K}}_{p,q} = k(x_{p,q},y_{p,q},Z_{p,q}, \overline{\overline{Z}}_{x_{p,q}},\overline{\overline{Z}}_{y_{p,q}}). $$
(2.9)

Further, let

$$\begin{aligned} & \begin{aligned}[b]\hat{\hat{Z}}_{x_{p,q}} ={}& \overline{Z}_{x_{p,q}} - \frac{h}{8A_{00}}(\overline{K}_{p + 1,q} - \overline{K}_{p - 1,q}) + \frac{hC_{00}}{8A_{00}}(\overline{Z}_{yy_{p + 1,q}} - \overline{Z}_{yy_{p - 1,q}})\\ &{} + \frac{h^{2}}{4A_{00}}(A_{10}\overline{Z}_{xx_{p,q}} + C_{10} \overline{Z}_{yy_{p,q}}), \end{aligned} \end{aligned}$$
(2.10a)
$$\begin{aligned} &\begin{aligned}[b] \hat{\hat{Z}}_{y_{p,q}} ={}& \overline{Z}_{y_{p,q}} - \frac{h}{8C_{00}}(\overline{K}_{p,q + 1} - \overline{K}_{p,q - 1}) + \frac{hA_{00}}{8C_{00}}(\overline{Z}_{xx_{p,q + 1}} - \overline{Z}_{xx_{p,q - 1}}) \\ &{}+ \frac{h^{2}}{4C_{00}}(A_{01}\overline{Z}_{xx_{p,q}} + C_{01} \overline{Z}_{yy_{p,q}}), \end{aligned} \end{aligned}$$
(2.10b)

and let

$$ \hat{\hat{K}}_{p,q} = k(x_{p,q},y_{p,q},Z_{p,q}, \hat{\hat{Z}}_{x_{p,q}},\hat{\hat{Z}}_{y_{p,q}}). $$
(2.11)

We denote

$$\begin{aligned} &P_{1} = A_{10}/A_{00},\qquad P_{2} = C_{01}/C_{00}, \\ &I_{1} = A_{00} + \frac{1}{12}h^{2}(A_{20} + A_{02} - 2P_{1}A_{10} - 2P_{2}A_{01}), \\ & I_{2} = C_{00} + \frac{1}{12}h^{2}(C_{20} + C_{02} - 2P_{1}C_{10} - 2P_{2}C_{01}), \\ &I_{3} = \frac{1}{12}h(A_{01} - P_{2}A_{00}), \qquad I_{4} = \frac{1}{12}h(C_{10} - P_{1}C_{00}),\qquad I_{5} = \frac{1}{12}(A_{00} + C_{00}), \\ &J_{1} = 1 - hP_{1},\qquad J_{2} = 1 + hP_{1},\qquad J_{3} = 1 - hP_{2},\qquad J_{4} = 1 + hP_{2}. \end{aligned}$$

Assume that \(K_{p,q} \neq0\), then at each interior nodal point \(( x_{p}, y_{q} )\), the differential equation (2.1) is discretized as

$$\begin{aligned} &\bigl[I_{1}\delta_{x}^{2} + I_{2} \delta_{y}^{2} + I_{3}\bigl(2\delta_{x}^{2} \mu_{y}\delta_{y}\bigr) + I_{4}\bigl(2 \delta_{y}^{2}\mu_{x}\delta_{x}\bigr) + I_{5}\bigl(\delta_{x}^{2}\delta_{y}^{2} \bigr)\bigr]Z_{p,q} \\ &\quad = h^{2}\overline{\overline{K}}_{p,q}\exp \biggl( \frac{J_{1}\overline{K}_{p + 1,q} + J_{2}\overline{K}_{p - 1,q} + J_{3}\overline{K}_{p,q + 1} + J_{4}\overline{K}_{{p,q - 1}} - 4\hat{\hat{K}}_{{p,q}}}{12\overline{\overline{K}}_{p,q}} \biggr) + \overline{E}_{p,q}, \\ &\quad 1\leq p,q\leq N, \end{aligned}$$
(2.12)

where \(\overline{E}_{p,q} = O(h^{6})\).

3 Deriving the numerical scheme

For the derivation of the scheme (2.12), at the nodal point (\(x _{p,} y _{q}\)), we denote

$$ \alpha_{p,q} = \biggl( \frac{\partial k}{\partial z_{x}} \biggr)_{p,q},\qquad \beta_{p,q} = \biggl( \frac{\partial k}{\partial z_{y}} \biggr)_{p,q}. $$
(3.1)

With the help of (2.3b), (2.4b) and (3.1) from (2.7a), we get

$$\begin{aligned} \overline{K}_{p \pm 1,q} =& k\biggl(x_{p \pm 1},y_{q,}Z_{p \pm 1,q},Z_{x_{p \pm 1,q}} - \frac{h^{2}}{3}Z_{30} \pm O\bigl(h^{3} \bigr),Z_{y_{p \pm 1,q}} + \frac{h^{2}}{6}Z_{30} \pm O \bigl(h^{3}\bigr)\biggr) \\ =& k(x_{p \pm 1},y_{q,}Z_{p \pm 1,q},Z_{x_{p \pm 1,q}},Z_{y_{p \pm 1,q}}) - \frac{h^{2}}{3}Z_{30}\alpha_{p \pm 1,q} + \frac{h^{2}}{6}Z_{03} \beta_{p \pm 1,q} + O\bigl( \pm h^{3} + h^{4}\bigr)) \\ =& K_{p \pm 1,q} - \frac{h^{2}}{3}Z_{30}\bigl( \alpha_{p,q} \pm O(h)\bigr) + \frac{h^{2}}{6}Z_{03}\bigl( \beta_{p,q} \pm O(h)\bigr) + O\bigl( \pm h^{3} + h^{4}\bigr)) \\ =& K_{p \pm 1,q} + \frac{h^{2}}{6}E_{1} + O\bigl( \pm h^{3} + h^{4}\bigr). \end{aligned}$$
(3.2)

Similarly, using (2.3c), (2.4c) and (3.1) from (2.7b), we get

$$ \overline{K}_{p,q \pm 1} = K_{p,q \pm 1} + \frac{h^{2}}{6}E_{2} + O\bigl( \pm h^{3} + h^{4}\bigr), $$
(3.3)

where

$$\begin{aligned} &E_{1} = - 2\alpha_{p,q}Z_{30} + \beta_{p,q}Z_{03}, \end{aligned}$$
(3.4)
$$\begin{aligned} &E_{2} = \alpha_{p,q}Z_{30} - 2 \beta_{p,q}Z_{03}. \end{aligned}$$
(3.5)

Let us consider the following linear combinations:

$$\begin{aligned} & \overline{\overline{Z}}_{x_{p,q}} = \overline{Z}_{x_{p,q}} + ha_{1}(\overline{K}_{p + 1,q} - \overline{K}_{p - 1,q}) + ha_{2}(\overline{Z}_{yy_{p + 1,q}} - \overline{Z}_{yy_{p - 1,q}}) \\ &\hphantom{\overline{\overline{Z}}_{x_{p,q}} =}{}+ 2h^{2}(a_{3}\overline{Z}_{xx_{p,q}} + a_{4} \overline{Z}_{yy_{p,q}}), \end{aligned}$$
(3.6a)
$$\begin{aligned} &\overline{\overline{Z}}_{y_{p,q}} = \overline{Z}_{y_{p,q}} + hb_{1}(\overline{K}_{p,q + 1} - \overline{K}_{p,q - 1}) + hb_{2}(\overline{Z}_{xx_{p,q + 1}} - \overline{Z}_{xx_{p,q - 1}}) \\ &\hphantom{\overline{\overline{Z}}_{y_{p,q}} = }{} + 2h^{2}(b_{3}\overline{Z}_{xx_{p,q}} + b_{4} \overline{Z}_{yy_{p,q}}), \end{aligned}$$
(3.6b)
$$\begin{aligned} &\hat{\hat{Z}}_{x_{p,q}} = \overline{Z}_{x_{p,q}} + ha'_{1}(\overline{K}_{p + 1,q} - \overline{K}_{p - 1,q}) + ha'_{2}( \overline{Z}_{yy_{p + 1,q}} - \overline{Z}_{yy_{p - 1,q}}) \\ &\hphantom{\hat{\hat{Z}}_{x_{p,q}}=}{}+ 2h^{2} \bigl(a'_{3}\overline{Z}_{xx_{p,q}} + a'_{4}\overline{Z}_{yy_{p,q}}\bigr), \end{aligned}$$
(3.7a)
$$\begin{aligned} &\hat{\hat{Z}}_{y_{p,q}} = \overline{Z}_{y_{p,q}} + hb'_{1}(\overline{K}_{p,q + 1} - \overline{K}_{p,q - 1}) + hb'_{2}( \overline{Z}_{xx_{p,q + 1}} - \overline{Z}_{xx_{p,q - 1}}) \\ &\hphantom{\hat{\hat{Z}}_{y_{p,q}} =}{}+ 2h^{2} \bigl(b'_{3}\overline{Z}_{xx_{p,q}} + b'_{4}\overline{Z}_{yy_{p,q}}\bigr), \end{aligned}$$
(3.7b)

where \(a_{i}\), \(b_{i}\), \(a'_{i}\), \(b'_{i}\), \(i = 1(1)4\) are parameters to be determined.

Using (2.3a), (2.5a), (2.6a), (2.6b), (3.2) in (3.6a) and (2.4a), (2.5a), (2.5b), (2.6a), (3.3) in (3.6b), we get

$$\begin{aligned} \overline{\overline{Z}}_{x_{p,q}} =& Z_{x_{p,q}} + \frac{h^{2}}{6}(1 + 12a_{1}A_{00})Z_{30} + 2h^{2}(a_{1}C_{00} + a_{2})Z_{12} + 2h^{2}(a_{1}A_{10} + a_{3})Z_{20} \\ &{}+ 2h^{2}(a_{1}C_{10} + a_{4})Z_{02} + O\bigl(h^{4}\bigr), \end{aligned}$$
(3.8)
$$\begin{aligned} \overline{\overline{Z}}_{y_{p,q}} =& Z_{y_{p,q}} + \frac{h^{2}}{6}(1 + 12b_{1}C_{00})Z_{03} + 2h^{2}(b_{1}A_{00} + b_{2})Z_{21} + 2h^{2}(b_{1}A_{01} + b_{3})Z_{20} \\ &{}+ 2h^{2}(b_{1}C_{01} + b_{4})Z_{02} + O\bigl(h^{4}\bigr). \end{aligned}$$
(3.9)

Note that

$$\overline{\overline{Z}}_{x_{p,q}} = Z_{x_{p,q}} + O \bigl(h^{4}\bigr), $$

and

$$\overline{\overline{Z}}_{y_{p,q}} = Z_{y_{p,q}} + O \bigl(h^{4}\bigr), $$

if

$$\begin{aligned} &a_{1} = \frac{ - 1}{12A_{00}},\qquad a_{2} = \frac{C_{00}}{12A_{00}},\qquad a_{3} = \frac{A_{10}}{12A_{00}},\qquad a_{4} = \frac{C_{10}}{12A_{00}}, \\ & b_{1} = \frac{ - 1}{12C_{00}},\qquad b_{2} = \frac{A_{00}}{12C_{00}},\qquad b_{3} = \frac{A_{01}}{12C_{00}},\qquad b_{4} = \frac{C_{01}}{12C_{00}}. \end{aligned}$$

For this choice of parameters, from (2.9), it is easy to verify that

$$ \overline{\overline{K}}_{p,q} = K_{p,q} + O \bigl(h^{4}\bigr). $$
(3.10)

Using (2.3a), (2.5a), (2.6a), (2.6b), (3.1) in (3.6a) and (2.4a), (2.5a), (2.5b), (2.6a), (3.2) in (3.6b) we get

$$\begin{aligned} &\hat{\hat{Z}}_{x_{p,q}} = Z_{x_{p,q}} + \frac{h^{2}}{6}E_{3} + O\bigl(h^{4}\bigr). \end{aligned}$$
(3.11)
$$\begin{aligned} &\hat{\hat{Z}}_{y_{p,q}} = Z_{y_{p,q}} + \frac{h^{2}}{6}E_{4} + O\bigl(h^{4}\bigr). \end{aligned}$$
(3.12)

where

$$\begin{aligned} & \begin{aligned}[b] E_{3} = {}&\bigl(1 + 12a'_{1}A_{00} \bigr)Z_{30} + \bigl(12a'_{1}C_{00} + 12a'_{2}\bigr)Z_{12} + \bigl(12a'_{1}A_{10} + 12a'_{3}\bigr)Z_{20} \\ &{}+ \bigl(12a'_{1}C_{10} + 12a'_{4}\bigr)Z_{02}, \end{aligned} \end{aligned}$$
(3.13)
$$\begin{aligned} &\begin{aligned}[b] E_{4} = {}&\bigl(1 + 12b'_{1}C_{00} \bigr)Z_{03} + \bigl(12b'_{1}A_{00} + 12b'_{2}\bigr)Z_{21} + \bigl(12b'_{1}A_{01} + 12b'_{3}\bigr)Z_{20}\\ &{} + \bigl(12b'_{1}C_{01} + 12b'_{4}\bigr)Z_{02}. \end{aligned} \end{aligned}$$
(3.14)

Using (3.10) and (3.11) in (2.11) we get

$$ \hat{\hat{K}}_{p,q} = K_{p,q} + \frac{h^{2}}{6}E_{3} \alpha_{p,q} + \frac{h^{2}}{6}E_{4}\beta_{p,q} + O \bigl(h^{4}\bigr). $$
(3.15)

With the help of a Taylor expansion, it is easy to verify that

$$\begin{aligned} &\bigl[I_{1}\delta_{x}^{2} + I_{2} \delta_{y}^{2} + I_{3}\bigl(2\delta_{x}^{2} \mu_{y}\delta_{y}\bigr) + I_{4}\bigl(2 \delta_{y}^{2}\mu_{x}\delta_{x}\bigr) + I_{5}\bigl(\delta_{x}^{2}\delta_{y}^{2} \bigr)\bigr]Z_{p,q} \\ &\quad = h^{2}K_{p,q}\exp \biggl( \frac{J_{1}K_{{p + 1,q}} + J_{2}K_{{p - 1,q}} + J_{3}K_{{p,q + 1}} + J_{4}K_{{p,q - 1}} - 4K_{{p,q}}}{12K_{p,q}} \biggr) + O\bigl(h^{6}\bigr), \\ &\quad K_{p,q} \neq0. \end{aligned}$$
(3.16)

With the help of approximations (3.2), (3.3), (3.10), (3.15), from (2.12) and (3.16), we get the local truncation error

$$ \overline{E}_{p,q} = \frac{ - h^{4}}{36} [ E_{1} + E_{2} - 2E_{3}\alpha_{p,q} - 2E_{4} \beta_{p,q} ] + O\bigl(h^{6}\bigr). $$
(3.17)

For the proposed method (2.12) to be of O(\(h ^{4}\)), the coefficient of \(h ^{4}\) in (3.17) must be zero. Thus we obtain the values of the coefficients,

$$\begin{aligned} & a'_{1} = \frac{ - 1}{8A_{00}},\qquad a'_{2} = \frac{C_{00}}{8A_{00}},\qquad a'_{3} = \frac{A_{10}}{8A_{00}}, \qquad a'_{4} = \frac{C_{10}}{8A_{00}}, \\ &b'_{1} = \frac{ - 1}{8C_{00}},\qquad b'_{2} = \frac{A_{00}}{8C_{00}},\qquad b'_{3} = \frac{A_{01}}{8C_{00}}, \qquad b'_{4} = \frac{C_{01}}{8C_{00}}. \end{aligned}$$

We now consider the numerical method of \(O(h^{4})\) for the solution of the 2D quasilinear elliptic equation (1.1). Using the technique discussed in [19], we can get the fourth order method for the quasilinear equation (1.1).

4 Study of convergence

Consider the following 2D nonlinear elliptic partial differential equation:

$$ Az_{xx} + Cz_{yy} = k(x,y,z,z_{x},z_{y}), $$
(4.1)

defined in the region Γ and subject to \(z(x,y) = z_{0}(x,y), (x,y) \in \partial \varGamma\), where A and Care positive constants.

Then the difference method (2.12) for Eq. (4.1) becomes

$$\begin{aligned} &\lambda_{1}(Z_{p + 1,q} + Z_{p - 1,q}) + \lambda_{2}(Z_{p,q + 1} + Z_{p,q - 1}) \\ &\qquad {} + \lambda_{3}(Z_{p + 1,q + 1} + Z_{p + 1,q - 1} + Z_{p - 1,q + 1} + Z_{p - 1,q - 1} - 20Z_{p,q}) \\ &\quad = 6 h^{2}\overline{\overline{K}}_{p,q}\exp \biggl( \frac{J_{1}\overline{K}_{{p + 1,q}} + J_{2}\overline{K}_{{p - 1,q}} + J_{3}\overline{K}_{{p,q + 1}} + J_{4}\overline{K}_{{p,q - 1}} - 4\hat{\hat{K}}_{{p,q}}}{12\overline{\overline{K}}_{p,q}} \biggr) + \overline{E}_{p,q}, \\ &\quad 1\leq p,q \leq N, \end{aligned}$$
(4.2)

where \(\overline{E}_{p,q} = O(h^{6})\), \(\lambda_{1} = 5A - C\), \(\lambda_{2} = 5C - A\) and \(\lambda_{3} = \frac{A + C}{2}\). The conditions which are usually imposed on (4.2) are \(\lambda _{1}\)>0 and \(\lambda _{2}\)>0.

The scheme (4.2) can easily be written in matrix form.

Let \(\boldsymbol{S} = [\boldsymbol{S}_{1}, \boldsymbol{S}_{2}, \boldsymbol{S}_{1}]_{N^{2} \times N^{2}}\) be a triblock-diagonal matrix, where \(S_{1} = [ - \lambda_{3}, - \lambda_{2}, - \lambda_{3}\ ]_{N \times N}\) and \(S_{2} = [ - \lambda_{1}, 20\lambda_{3}, - \lambda_{1}\ ]_{N \times N}\) are tri-diagonal matrices.

Let

$$\phi_{p,q} = 6h^{2}\overline{\overline{K}}_{p,q}\exp \biggl( \frac{J_{1}\overline{K}_{{p + 1,q}} + J_{2}\overline{K}_{{p - 1,q}} + J_{3}\overline{K}_{{p,q + 1}} + J_{4}\overline{K}_{{p,q - 1}} - 4\hat{\hat{K}}_{{p,q}}}{12\overline{\overline{K}}_{p,q}} \biggr). $$

Then the method (4.2) in matrix form may be written as

$$ \boldsymbol{SZ} + \boldsymbol{\phi Z} + \boldsymbol{E} = \boldsymbol{0}, $$
(4.3)

where E is the local truncation error vector.

Thus the method involves computing a numerical value z for the exact value Z by solving a system of \(( N^{2} \times N^{2} )\) equations:

$$ \boldsymbol{Sz} + \boldsymbol{\phi z} = \boldsymbol{0}. $$
(4.4)

Let

$$ \varepsilon_{p,q} = z_{p,q} - Z_{p,q} \quad \bigl(p=1(1)N, q=1(1)N\bigr), $$
(4.5)

and

$$\boldsymbol{T} = \boldsymbol{z} - Z $$

Let

$$\begin{aligned} &\overline{k}_{p \pm 1,q} = k(x_{p \pm 1},y_{q},z_{p \pm 1,q}, \bar{z}_{x_{p \pm 1,q}},\bar{z}_{y_{p \pm 1,q}}) \approx \overline{K}_{p \pm 1,q}, \end{aligned}$$
(4.6a)
$$\begin{aligned} &\overline{k}_{p,q \pm 1} = k(x_{p},y_{q \pm 1},z_{p,q \pm 1}, \bar{z}_{x_{p,q \pm 1}},\bar{z}_{y_{p,q \pm 1}}) \approx \overline{K}_{p,q \pm 1}, \end{aligned}$$
(4.6b)
$$\begin{aligned} &\overline{\overline{k}}_{p,q} = k(x_{p},y_{q},z_{p,q}, \bar{\bar{z}}_{x_{p,q}},\bar{\bar{z}}_{y_{p,q}}) \approx \overline{ \overline{K}}_{p,q}, \end{aligned}$$
(4.6c)
$$\begin{aligned} &\hat{\hat{k}}_{p,q} = k(x_{p},y_{q},z_{p,q}, \hat{\hat{z}}_{x_{p,q}},\hat{\hat{z}}_{y_{p,q}}) \approx \hat{ \hat{K}}_{p,q}. \end{aligned}$$
(4.6d)

We may write

$$\begin{aligned} & \begin{aligned}[b] \overline{k}_{p \pm 1,q} - \overline{K}_{p \pm 1,q} ={}& \varepsilon_{p \pm 1,q}V_{p \pm 1,q}^{(1)} \\ &{}+ (\overline{z}_{x_{p \pm 1,q}} - \overline{Z}_{x_{p \pm 1,q}})G_{p \pm 1,q}^{(1)} + ( \overline{z}_{y_{p \pm 1,q}} - \overline{Z}_{y_{p \pm 1,q}})H_{p \pm 1,q}^{(1)}, \end{aligned} \end{aligned}$$
(4.7a)
$$\begin{aligned} & \begin{aligned}[b]\overline{k}_{p,q \pm 1} - \overline{K}_{p,q \pm 1} ={}& \varepsilon_{p,q \pm 1}V_{p,q \pm 1}^{(2)} + (\overline{z}_{x_{p,q \pm 1}} - \overline{Z}_{x_{p,q \pm 1}})G_{p,q \pm 1}^{(2)}\\ &{} + ( \overline{z}_{y_{p,q \pm 1}} - \overline{Z}_{y_{p,q \pm 1}})H_{p,q \pm 1}^{(2)}, \end{aligned} \end{aligned}$$
(4.7b)
$$\begin{aligned} &\overline{\overline{k}}_{p,q} - \overline{\overline{K}}_{p,q} = \varepsilon_{p,q}V_{p,q}^{(3)} + (\bar{ \bar{z}}_{x_{p,q}} - \overline{\overline{Z}}_{x_{p,q}})G_{p,q}^{(3)} + (\bar{\bar{z}}_{y_{p,q}} - \overline{\overline{Z}}_{y_{p,q}})H_{p,q}^{(3)}, \end{aligned}$$
(4.7c)
$$\begin{aligned} &\hat{\hat{k}}_{p,q} - \hat{\hat{K}}_{p,q} = \varepsilon_{p,q}V_{p,q}^{(4)} + (\hat{ \hat{z}}_{x_{p,q}} - \hat{\hat{Z}}_{x_{p,q}})G_{p,q}^{(4)} + (\hat{\hat{z}}_{y_{p,q}} - \hat{\hat{Z}}_{y_{p,q}})H_{p,q}^{(4)}, \end{aligned}$$
(4.7d)

for suitable \(R_{p \pm 1,q}^{(1)}\), \(R_{p,q \pm 1}^{(2)}\), \(R_{p,q}^{(3)}\) and \(R_{p,q}^{(4)}\), where \(R = V, G\) and H.

Also, for \(R = G\) and H we may write

$$\begin{aligned} &R_{p \pm 1,q}^{(1)} = R_{p,q}^{(1)} \pm hR_{x_{p,q}}^{(1)} + O\bigl(h^{2}\bigr), \end{aligned}$$
(4.8a)
$$\begin{aligned} &R_{p,q \pm 1}^{(2)} = R_{p,q}^{(2)} \pm hR_{y_{p,q}}^{(2)} + O\bigl(h^{2}\bigr) \end{aligned}$$
(4.8b)

and

$$\begin{aligned} &V_{p \pm 1,q}^{(1)} = V_{p,q}^{(1)} \pm O(h), \end{aligned}$$
(4.8c)
$$\begin{aligned} &V_{p,q \pm 1}^{(2)} = V_{p,q}^{(2)} \pm O(h). \end{aligned}$$
(4.8d)

With the help of Eqs. (4.7a)–(4.7d) and (4.8a)–(4.8d), we get

$$ \boldsymbol{\phi z} - \boldsymbol{\phi Z} = \boldsymbol{QT}. $$
(4.9)

where \(\boldsymbol{Q} = (Q_{i,j})\), \(((i = 1(1)N^{2}), j = 1(1)N^{2})\) is the triblock-diagonal matrix with

$$\begin{aligned} &Q_{(q - 1)N + p,(q - 1)N + p} \\ &\quad = h^{2}\bigl[\alpha_{1}G_{p,q}^{(1)} \bigl(2G_{p,q}^{(3)} - G_{p,q}^{(4)}\bigr) + \beta_{1}H_{p,q}^{(2)}\bigl(2H_{p,q}^{(3)} - H_{p,q}^{(4)}\bigr) + 2\bigl(3V_{p,q}^{(3)} - V_{p,q}^{(4)}\bigr) \\ &\qquad {}- 2G_{x_{p,q}}^{(1)} - 2H_{y_{p,q}}^{(2)} \bigr] + O\bigl(h^{4}\bigr),\quad \bigl(\bigl(p = 1(1)N\bigr),q = 1(1)N \bigr); \\ &Q_{(q - 1)N + p,(q - 1)N + p \pm 1} \\ &\quad = \pm \frac{h}{2}\bigl[G_{p,q}^{(1)} + 2 \bigl(3G_{p,q}^{(3)} - G_{p,q}^{(4)}\bigr) - \alpha_{2}\bigl(2G_{p,q}^{(3)} - G_{p,q}^{(4)} \bigr)\bigr] \\ &\qquad {} + \frac{h^{2}}{2}\bigl[V_{p,q}^{(1)} + 2G_{x_{p,q}}^{(1)} - \alpha_{1}G_{p,q}^{(1)} \bigl(2G_{p,q}^{(3)} - G_{p,q}^{(4)}\bigr) \bigr] + O\bigl(h^{3}\bigr),\\ &\quad \bigl(\bigl(p = 1(1)N - 1,2(1)N \bigr),q = 1(1)N\bigr); \\ &Q_{(q - 1)N + p,(q - 1 \pm 1)N + p} \\ &\quad = \pm \frac{h}{2}\bigl[H_{p,q}^{(2)} + 2 \bigl(3H_{p,q}^{(3)} - H_{p,q}^{(4)}\bigr) - \beta_{2}\bigl(2H_{p,q}^{(3)} - H_{p,q}^{(4)} \bigr)\bigr] \\ &\qquad {} + \frac{h^{2}}{2}\bigl[V_{p,q}^{(2)} + 2H_{y_{p,q}}^{(2)} - \beta_{1}H_{p,q}^{(2)} \bigl(2H_{p,q}^{(3)} - H_{p,q}^{(4)}\bigr) \bigr] + O\bigl(h^{3}\bigr),\\ & \quad \bigl(\bigl(p = 1(1)N\bigr),q = 1(1)N - 1,2(1)N\bigr); \\ & Q_{(q - 1)N + p,qN + p \pm 1} \\ &\quad = \frac{h}{4}\bigl[ \pm G_{p,q}^{(2)} + H_{p,q}^{(1)} \pm \alpha_{2}\bigl(2G_{p,q}^{(3)} - G_{p,q}^{(4)}\bigr) + \beta_{2} \bigl(2H_{p,q}^{(3)} - H_{p,q}^{(4)}\bigr) \bigr] \\ &\qquad {} \pm \frac{h^{2}}{8}\bigl[2H_{x_{p,q}}^{(1)} + 2G_{y_{p,q}}^{(2)} - \alpha_{1}H_{p,q}^{(1)} \bigl(2G_{p,q}^{(3)} - G_{p,q}^{(4)}\bigr) - \beta_{1}G_{p,q}^{(2)}\bigl(2H_{p,q}^{(3)} - H_{p,q}^{(4)}\bigr)\bigr] + O\bigl(h^{3}\bigr), \\ &\quad \bigl(\bigl(p = 1(1)N -1,2(1)N\bigr),q = 1(1)N - 1\bigr); \\ &Q_{(q - 1)N + p,(q - 2)N + p \pm 1} \\ &\quad = \frac{h}{4}\bigl[ \pm G_{p,q}^{(2)} - H_{p,q}^{(1)} \pm \alpha_{2}\bigl(2G_{p,q}^{(3)} - G_{p,q}^{(4)}\bigr) - \beta_{2} \bigl(2H_{p,q}^{(3)} - H_{p,q}^{(4)}\bigr) \bigr] \\ &\qquad {}\pm \frac{h^{2}}{8}\bigl[ - 2H_{x_{p,q}}^{(1)} - 2G_{y_{p,q}}^{(2)} + \alpha_{1}H_{p,q}^{(1)} \bigl(2G_{p,q}^{(3)} - G_{p,q}^{(4)}\bigr) + \beta_{1}G_{p,q}^{(2)}\bigl(2H_{p,q}^{(3)} - H_{p,q}^{(4)}\bigr)\bigr] + O\bigl(h^{3}\bigr), \\ &\quad \bigl(\bigl(p = 1(1)N -1,2(1)N\bigr),q = 2(1)N\bigr); \end{aligned}$$

where

$$\alpha_{1} = \frac{1}{A_{00}},\qquad \alpha_{2} = \frac{C_{00}}{A_{00}},\qquad \beta_{1} = \frac{1}{C_{00}}, \qquad \beta_{2} = \frac{A_{00}}{C_{00}}. $$

With the help of (4.9), from (4.3) and (4.4), we get the following equation:

$$ (\boldsymbol{S} + \boldsymbol{Q})\boldsymbol{T} = \boldsymbol{E}. $$
(4.10)

Let

$$V_{*}^{(2)} = \mathop{\mathrm{Min}}\limits _{(x,y) \in \overline{\varGamma}} \frac{\partial k}{\partial Z}\quad \mbox{and} \quad V_{(2)}^{*} = \mathop{\mathrm{Max}}\limits _{(x,y) \in \overline{\varGamma}} \frac{\partial k}{\partial Z}, $$

where \(\overline{\varGamma} = \varGamma \cup \partial \varGamma\).

Then \(0 < V_{*}^{(2)} \le V_{p \pm 1,q}^{(1)}, V_{p,q \pm 1}^{(2)}, V_{p,q}^{(3)}, V_{p,q}^{(4)} \le V_{(2)}^{*}\), and for \(R=G,H\), let \(0 < \vert R_{p \pm 1,q}^{(1)} \vert , \vert R_{p,q \pm 1}^{(2)} \vert , \vert R_{p,q}^{(3)} \vert , \vert R_{p,q}^{(4)} \vert \le R\), and \(\vert R_{{x_{p,q}}}^{(1)} \vert \le R^{(1)}\), \(\vert R_{{y_{p,q}}}^{(2)} \vert \le R^{(2)}\) for some positive constants \(R^{(1)}\) and \(R^{(2)}\).

Further it is easy to verify that, for howsoever small h,

$$\begin{aligned} & \vert Q_{(q - 1)N + p,(q - 1)N + p \pm 1} \vert < \lambda_{1},\quad \bigl(\bigl(p = 1(1)N - 1,2(1)N\bigr),q = 1(1)N\bigr); \\ & \vert Q_{(q - 1)N + p,(q - 1 \pm 1)N + p} \vert < \lambda_{2},\quad \bigl(\bigl(p = 1(1)N\bigr),q = 1(1)N - 1,2(1)N\bigr); \\ & \vert Q_{(q - 1)N + p,qN + p \pm 1} \vert < \lambda_{3},\quad \bigl(\bigl(p = 1(1)N - 1,2(1)N\bigr),q = 1(1)N - 1\bigr); \\ & \vert Q_{(q - 1)N + p,(q - 2)N + p \pm 1} \vert < \lambda_{3},\quad \bigl(\bigl(p = 1(1)N - 1,2(1)N\bigr),q = 2(1)N\bigr); \end{aligned}$$

Further, the directed graph of (\(\boldsymbol{S} + \boldsymbol{Q}\)) shows that it is an irreducible matrix (see Fig. 1). The arrows indicate the paths \(i \to j\) for every nonzero entry of the matrix \((\boldsymbol{S} + \boldsymbol{Q})\). For any ordered pair of nodes i and j, there exists a direct path \((\overrightarrow{i,l_{1}} ), (\overrightarrow{l_{1},l_{2}} ),\ldots, (\overrightarrow{l_{k},j} )\) connecting i to j. Thus, the graph is strongly connected. So, the matrix (\(\boldsymbol{S} + \boldsymbol{Q}\)) is irreducible. (See Varga [34], Young [35].)

Figure 1
figure 1

Directed graph of (\(\boldsymbol{S} + \boldsymbol{Q}\))

Let \(M_{k}\) be the sum of the elements in the kth row of \((\boldsymbol{S} + \boldsymbol{Q})\), then for \(k = 1\) and N, we have

$$ M_{k} = 11\lambda_{3} + \frac{h}{4} \biggl(d_{k} + \frac{h}{2}p_{k}\biggr) + \frac{h^{2}}{2}\bigl[V_{k,1}^{(1)} + V_{k,1}^{(2)} + 4\bigl(3V_{k,1}^{(3)} - V_{k,1}^{(4)} \bigr)\bigr] + O\bigl(h^{3}\bigr), $$
(4.11a)

where

$$\begin{aligned} & \begin{aligned} d_{k} = {}&\pm 2G_{k,1}^{(1)} \pm G_{k,1}^{(2)} \pm 4\bigl(3G_{k,1}^{(3)} - G_{k,1}^{(4)}\bigr) \mp \alpha_{2} \bigl(2G_{k,1}^{(3)} - G_{k,1}^{(4)}\bigr) \\ &{} + H_{k,1}^{(1)} + 2H_{k,1}^{(2)} + 4 \bigl(3H_{k,1}^{(3)} - H_{k,1}^{(4)}\bigr) - \beta_{2}\bigl(2H_{k,1}^{(3)} - H_{k,1}^{(4)} \bigr) \end{aligned} \\ & \begin{aligned} p_{k} ={}& 4\alpha_{1}G_{k,1}^{(1)} \bigl(2G_{k,1}^{(3)} - G_{k,1}^{(4)}\bigr) + 4\beta_{1}H_{k,1}^{(2)}\bigl(2H_{k,1}^{(3)} - H_{k,1}^{(4)}\bigr) \mp \beta_{1}G_{k,1}^{(2)} \bigl(2H_{k,1}^{(3)} - H_{k,1}^{(4)}\bigr) \\ &{} \mp \alpha_{1}H_{k,1}^{(1)} \bigl(2G_{k,1}^{(3)} - G_{k,1}^{(4)}\bigr) - 8G_{x_{k,1}}^{(1)} - 8H_{y_{k,1}}^{(2)} \pm 2G_{y_{k,1}}^{(2)} \pm 2H_{x_{k,1}}^{(1)}, \end{aligned} \\ &\begin{aligned} M_{(N - 1)N + k} ={}& 11\lambda_{3} + \frac{h}{4} \biggl(d_{(N - 1)N + k}+ \frac{h}{2}p_{(N - 1)N + k}\biggr) \\ &{}+ \frac{h^{2}}{2}\bigl[V_{k,N}^{(1)} + V_{k,N}^{(2)} + 4\bigl(3V_{k,N}^{(3)} - V_{k,N}^{(4)} \bigr)\bigr] + O\bigl(h^{3}\bigr), \end{aligned} \end{aligned}$$
(4.11b)

where

$$\begin{aligned} d_{(N - 1)N + k} ={}& \pm 2G_{k,N}^{(1)} \pm G_{k,N}^{(2)} \pm 4\bigl(3G_{k,N}^{(3)} - G_{k,N}^{(4)}\bigr) \mp \alpha_{2} \bigl(2G_{k,N}^{(3)} - G_{k,N}^{(4)}\bigr) \\ &{} - H_{k,N}^{(1)} - 2H_{k,N}^{(2)} - 4\bigl(3H_{k,N}^{(3)} - H_{k,N}^{(4)} \bigr) + \beta_{2}\bigl(2H_{k,N}^{(3)} - H_{k,N}^{(4)}\bigr) \\ p_{(N - 1)N + k} = {}&4\alpha_{1}G_{k,N}^{(1)} \bigl(2G_{k,N}^{(3)} - G_{k,N}^{(4)}\bigr) + 4\beta_{1}H_{k,N}^{(2)}\bigl(2H_{k,N}^{(3)} - H_{k,N}^{(4)}\bigr) \pm \beta_{1}G_{k,N}^{(2)} \bigl(2H_{k,N}^{(3)} - H_{k,N}^{(4)}\bigr) \\ &{} \pm \alpha_{1}H_{k,N}^{(1)} \bigl(2G_{k,N}^{(3)} - G_{k,N}^{(4)}\bigr) - 8G_{x_{k,N}}^{(1)} - 8H_{y_{k,N}}^{(2)} \mp 2G_{y_{k,N}}^{(2)} \mp 2H_{x_{k,N}}^{(1)}. \end{aligned}$$

For \(( 2 \le j \le N - 1 )\):

$$\begin{aligned} M_{(j - 1)N + k} ={}& 6A + \frac{h}{2} [ d_{(j - 1)N + k} + hp_{(j - 1)N + k} ] \\ &{}+ \frac{h^{2}}{2}\bigl[V_{k,j}^{(1)} + 2V_{k,j}^{(2)} + 4\bigl(3V_{k,j}^{(3)} - V_{k,j}^{(4)}\bigr)\bigr] + O\bigl(h^{3} \bigr), \end{aligned}$$
(4.11c)

where

$$\begin{aligned} &d_{(j - 1)N + k} = \pm \bigl[G_{k,j}^{(1)} + G_{k,J}^{(2)} + 2\bigl(3G_{k,j}^{(3)} - G_{k,j}^{(4)}\bigr)\bigr], \\ &p_{(j - 1)N + k} = \alpha_{1}\bigl[G_{k,j}^{(1)} \bigl(2G_{k,j}^{(3)} - G_{k,j}^{(4)}\bigr) \bigr] - 2G_{x_{k,j}}^{(1)}. \end{aligned}$$

For \(( i = 2(1)N - 1 )\):

$$\begin{aligned} M_{(k - 1)N + i} ={}& 6B + \frac{h}{2} [ d_{(k - 1)N + i} + hp_{(k - 1)N + i} ] \\ &{}+ \frac{h^{2}}{2}\bigl[2V_{i,k}^{(1)} + V_{i,k}^{(2)} + 4\bigl(3V_{i,k}^{(3)} - V_{i,k}^{(4)}\bigr)\bigr] + O\bigl(h^{3} \bigr), \end{aligned}$$
(4.11d)

where

$$\begin{aligned} &d_{(k - 1)N + i} = \pm \bigl[H_{i,k}^{(1)} + H_{i,k}^{(2)} + 2\bigl(3H_{i,k}^{(3)} - H_{i,k}^{(4)}\bigr)\bigr], \\ &p_{(k - 1)N + i} = \beta_{1}\bigl[H_{i,k}^{(2)} \bigl(2H_{i,k}^{(3)} - H_{i,k}^{(4)}\bigr) - 2H_{y_{i,k}}^{(2)}\bigr], \end{aligned}$$

and, finally, for \(((2 \le i \le N - 1), 2 \le j \le N - 1)\):

$$ M_{(j - 1)N + i} = h^{2}\bigl[V_{i,j}^{(1)} + V_{i,j}^{(2)} + 2\bigl(3V_{i,j}^{(3)} - V_{i,j}^{(4)}\bigr)\bigr] + O\bigl(h^{4} \bigr). $$
(4.11e)

With the help of Eqs. (4.11a)–(4.11e), for \(k = 1, N, (N - 1)N + 1\) and \(N^{2}\)

$$\begin{aligned} & \vert d_{k} \vert \le (19 + 3\alpha_{2})G + (19 + 3 \beta_{2})H, \\ & \vert p_{k} \vert \le 12\bigl(\alpha_{1}G^{2} + \beta_{1}H^{2}\bigr) + 3(\alpha_{1} + \beta_{1})GH + 8\bigl(G^{(1)} + H^{(2)}\bigr) + 2 \bigl(G^{(2)} + H^{(1)}\bigr), \end{aligned}$$

and for \(k = i\) and \((N - 1)N + i\); \(i = 2(1)N - 1\)

$$\begin{aligned} & \vert d_{k} \vert \le 10H, \\ & \vert p_{k} \vert \le 3\beta_{1}H^{2} + 2H^{(2)}, \end{aligned}$$

and for \(k = (j - 1)N + 1\) and jN; \(j = 2(1)N - 1\)

$$\begin{aligned} & \vert d_{k} \vert \le 10G, \\ &\vert p_{k} \vert \le 3\alpha_{1}G^{2} + 2G^{(1)}. \end{aligned}$$

Hence for h (howsoever small) using (4.11) it is easy to see that

$$\begin{aligned} &M_{k} > 9h^{2}V_{*}^{(2)}; \quad k = 1, N, (N - 1)N + 1 \mbox{ and } N^{2}, \end{aligned}$$
(4.12a)
$$\begin{aligned} &M_{k} > \frac{19}{2}h^{2}V_{*}^{(2)};\quad k = i \mbox{ and } (N - 1)N + i; i = 2(1)N - 1, \end{aligned}$$
(4.12b)
$$\begin{aligned} &M_{k} > \frac{19}{2}h^{2}V_{*}^{(2)};\quad k = (j - 1)N + 1 \mbox{ and } jN; j = 2(1)N - 1, \end{aligned}$$
(4.12c)
$$\begin{aligned} &M_{(j - 1)N + i} \ge 10h^{2}V_{*}^{(2)};\quad \bigl( \bigl(i = 2(1)N - 1\bigr), j = 2(1)N - 1\bigr). \end{aligned}$$
(4.12d)

Thus, for h (howsoever small), (\(\boldsymbol{S} + \boldsymbol{Q}\)) is monotone. Hence, \((\boldsymbol{S} + \boldsymbol{Q})^{ - 1}\) exists and \((\boldsymbol{S} + \boldsymbol{Q})^{ - 1}= \boldsymbol{J} >\boldsymbol{0} \) where \(\boldsymbol{J} = [J_{i,j}]_{N^{2} \times N^{2}}\).

Since \(\sum_{j = 1}^{N^{2}} J_{l,j}M_{j} = 1\), \(l = 1(1)N^{2}\), using (4.12a)–(4.12d) with \(l = 1(1)N^{2}\), we obtain

$$\begin{aligned} &J_{l,k} \le \frac{1}{M_{k}} < \frac{1}{9h^{2}V_{*}^{(2)}},\quad k = 1,N,(N - 1)N + 1,N^{2}, \end{aligned}$$
(4.13a)
$$\begin{aligned} &\sum_{i = 2}^{N - 1} J_{l,k} \le \frac{1}{\min_{2 \le i \le N - 1}M_{k}} < \frac{2}{19h^{2}V_{*}^{(2)}};\quad k = i,(N - 1)N + i, \end{aligned}$$
(4.13b)
$$\begin{aligned} &\sum_{j = 2}^{N - 1} J_{l,k} \le \frac{1}{\min_{2 \le j \le N - 1}M_{k}} < \frac{2}{19h^{2}V_{*}^{(2)}};\quad k = (j - 1)N + 1,jN, \end{aligned}$$
(4.13c)
$$\begin{aligned} &\sum_{i = 2}^{N - 1} \sum _{j = 2}^{N - 1} J_{l,k} \le \frac{1}{\min_{{\scriptstyle 2 \le i \le N - 1 \atop \scriptstyle 2 \le j \le N - 1}}M_{k}} \le \frac{1}{10h^{2}V_{*}^{(2)}},\quad k = (j - 1)N + i. \end{aligned}$$
(4.13d)

Since E is the local truncation error vector, we define \(\Vert \boldsymbol{E} \Vert =\max_{p,q} \vert \bar{E}_{p,q} \vert \leq C ^{*} h ^{6}\), where \(C ^{*}\) is a positive constant.

Now, from Eq. (4.10), we obtain

$$ \Vert \boldsymbol{T} \Vert \le \Vert \boldsymbol{J} \Vert \Vert \boldsymbol{E} \Vert , $$
(4.14)

where

$$\begin{aligned} \Vert \boldsymbol{J} \Vert =& \max_{1 \le l \le N^{2}} \Biggl[ \Biggl( J_{l,1} + \sum_{i = 2}^{(N - 1)} J_{l,i} + J_{l,N} \Biggr) + \Biggl( \sum _{j = 2}^{N - 1} J_{l,(j - 1)N + 1} + \sum _{i = 2}^{N - 1} \sum_{j = 2}^{N - 1} J_{l,(j - 1)N + i} + \sum_{j = 2}^{N - 1} J_{l,jN} \Biggr) \\ &{} + \Biggl( J_{l,(N - 1)N + 1} + \sum_{i = 2}^{N - 1} J_{l,(N - 1)N + i} + J_{l,N^{2}} \Biggr) \Biggr]. \end{aligned}$$
(4.15)

Using inequalities (4.13a)–(4.13d) in Eq. (4.15), we get

$$ \Vert \boldsymbol{J} \Vert \le \frac{1651}{1710h^{2}V_{*}^{(2)}}. $$
(4.16)

Finally, with the help of (4.16) for sufficiently small h, from (4.14) we obtain

$$ \Vert \boldsymbol{T} \Vert \le O\bigl(h^{4}\bigr). $$
(4.17)

This proves the convergence of the fourth order of the method (4.2) for the elliptic equation (4.1).

5 Method for system of equations

In this section, we extend our method to the system of quasilinear PDEs of the form

$$\begin{aligned} &A^{(i)}z^{(i)}_{xx} + C^{(i)}z^{(i)}_{yy} = k^{(i)}\bigl(x, y,z^{(1)},z^{(2)}, \ldots,z^{(n)}, z_{x}^{(1)},z_{x}^{(2)}, \ldots,z_{x}^{(n)},z_{y}^{(1)},z_{y}^{(2)}, \ldots,z_{y}^{(n)} \bigr), \\ &\quad i = 1(1)n, \end{aligned}$$
(5.1)

where \((x,y)\in \varGamma = (0,1) \times (0,1)\), \(A^{(i)} = A^{(i)}(x, y)\) and \(C^{(i)} = C^{(i)}(x, y)\).

The boundary conditions of Dirichlet type are given by

$$ z^{(i)}(x,y) = z_{0}^{(i)}(x,y),\quad (x,y) \in \partial \varGamma. $$
(5.2)

We assume \(Z_{p,q}^{(i)}\) and \(z_{p,q}^{(i)}\) to be the exact and approximate values of \(z^{(i)}(x_{p},y_{q})\) respectively.

For each \(i=1(1)n\), let

$$K_{p,q}^{(i)} = k^{(i)}\bigl(x_{p}, y_{q},Z_{p,q}^{(1)},Z_{p,q}^{(2)}, \ldots,Z_{p,q}^{(n)}, Z_{x_{p,q}}^{(1)},Z_{x_{p,q}}^{(2)}, \ldots,Z_{x_{p,q}}^{(n)},Z_{y_{p,q}}^{(1)},Z_{y_{p,q}}^{(2)}, \ldots,Z_{y_{p,q}}^{(n)} \bigr). $$

We define the following approximations:

$$\begin{aligned} &\overline{Z}_{x_{p,q}}^{(i)} = \frac{1}{2h} [ 2 \mu_{x}\delta_{x} ]Z_{p,q}^{(i)}, \end{aligned}$$
(5.3a)
$$\begin{aligned} &\overline{Z}_{x_{p \pm 1,q}}^{(i)} = \frac{1}{2h} \bigl[ 2 \mu_{x}\delta_{x} \pm 2\delta_{x}^{2} \bigr]Z_{p,q}^{(i)}, \end{aligned}$$
(5.3b)
$$\begin{aligned} &\overline{Z}_{x_{p,q \pm 1}}^{(i)} = \frac{1}{2h} [ 2 \mu_{x}\delta_{x} ]Z_{p,q \pm 1}^{(i)}, \end{aligned}$$
(5.3c)
$$\begin{aligned} &\overline{Z}_{y_{p,q}}^{(i)} = \frac{1}{2h} [ 2 \mu_{y}\delta_{y} ]Z_{p,q}^{(i)}, \end{aligned}$$
(5.4a)
$$\begin{aligned} &\overline{Z}_{y_{p \pm 1,q}}^{(i)} = \frac{1}{2h} [ 2 \mu_{y}\delta_{y} ]Z_{p \pm 1,q}^{(i)}, \end{aligned}$$
(5.4b)
$$\begin{aligned} &\overline{Z}_{y_{p,q \pm 1}}^{(i)} = \frac{1}{2h} \bigl[ 2 \mu_{y}\delta_{y} \pm 2\delta_{y}^{2} \bigr]Z_{p,q}^{(i)}, \end{aligned}$$
(5.4c)
$$\begin{aligned} &\overline{Z}_{xx_{p,q}}^{(i)} = \frac{1}{h^{2}} \bigl[ \delta_{x}^{2} \bigr]Z_{p,q}^{(i)}, \end{aligned}$$
(5.5a)
$$\begin{aligned} &\overline{Z}_{xx_{p,q \pm 1}}^{(i)} = \frac{1}{h^{2}} \bigl[ \delta_{x}^{2} \bigr]Z_{p,q \pm 1}^{(i)}, \end{aligned}$$
(5.5b)
$$\begin{aligned} &\overline{Z}_{yy_{p,q}}^{(i)} = \frac{1}{h^{2}} \bigl[ \delta_{y}^{2} \bigr]Z_{p,q}^{(i)}, \end{aligned}$$
(5.6a)
$$\begin{aligned} &\overline{Z}_{yy_{p \pm 1,q}}^{(i)} = \frac{1}{h^{2}} \bigl[ \delta_{y}^{2} \bigr]Z_{p \pm 1,q}^{(i)}. \end{aligned}$$
(5.6b)

Define

$$\begin{aligned} &\overline{K}_{p \pm a,q \pm b}^{(i)} = k^{(i)} \bigl(x_{p \pm a},y_{q \pm b},Z_{p \pm a,q \pm b}^{(1)},Z_{p \pm a,q \pm b}^{(2)}, \ldots,Z_{p \pm a,q \pm b}^{(n)},\overline{Z}_{x_{p \pm a,q \pm b}}^{(1)}, \overline{Z}_{x_{p \pm a,q \pm b}}^{(2)},\ldots, \\ &\hphantom{\overline{K}_{p \pm a,q \pm b}^{(i)} =}\overline{Z}_{x_{p \pm a,q \pm b}}^{(n)}, \overline{Z}_{y_{p \pm a,q \pm b}}^{(1)},\overline{Z}_{y_{p \pm a,q \pm b}}^{(2)}, \ldots,\overline{Z}_{y_{p \pm a,q \pm b}}^{(n)}\bigr); \\ &\quad (a = 0, b = 1; a = 1,b = 0). \end{aligned}$$
(5.7)

Let

$$\begin{aligned} \overline{\overline{Z}}_{x_{p,q}}^{(i)} = {}&\overline{Z}_{x_{p,q}}^{(i)} - \frac{h}{12A_{00}^{(i)}}\bigl(\overline{K}_{p + 1,q}^{(i)} - \overline{K}_{p - 1,q}^{(i)}\bigr) + \frac{hC_{00}^{(i)}}{12A_{00}^{(i)}}\bigl( \overline{Z}_{yy_{p + 1,q}}^{(i)} - \overline{Z}_{yy_{p - 1,q}}^{(i)} \bigr) \\ &{}+ \frac{h^{2}}{6A_{00}^{(i)}}\bigl(A_{10}^{(i)} \overline{Z}_{xx_{p,q}}^{(i)} + C_{10}^{(i)} \overline{Z}_{yy_{p,q}}^{(i)}\bigr), \end{aligned}$$
(5.8a)
$$\begin{aligned} \overline{\overline{Z}}_{y_{p,q}}^{(i)} = {}&\overline{Z}_{y_{p,q}}^{(i)} - \frac{h}{12C_{00}^{(i)}}\bigl(\overline{K}_{p,q + 1}^{(i)} - \overline{K}_{p,q - 1}^{(i)}\bigr) + \frac{hA_{00}^{(i)}}{12C_{00}^{(i)}}\bigl( \overline{Z}_{xx_{p,q + 1}}^{(i)} - \overline{Z}_{xx_{p,q - 1}}^{(i)} \bigr) \\ &{}+ \frac{h^{2}}{6C_{00}^{(i)}}\bigl(A_{01}^{(i)} \overline{Z}_{xx_{p,q}}^{(i)} + C_{01}^{(i)} \overline{Z}_{yy_{p,q}}^{(i)}\bigr), \end{aligned}$$
(5.8b)
$$\begin{aligned} \hat{\hat{Z}}_{x_{p,q}}^{(i)} = {}&\overline{Z}_{x_{p,q}}^{(i)} - \frac{h}{8A_{00}^{(i)}}\bigl(\overline{K}_{p + 1,q}^{(i)} - \overline{K}_{p - 1,q}^{(i)}\bigr) + \frac{hC_{00}^{(i)}}{8A_{00}^{(i)}}\bigl( \overline{Z}_{yy_{p + 1,q}}^{(i)} - \overline{Z}_{yy_{p - 1,q}}^{(i)} \bigr) \\ &{}+ \frac{h^{2}}{4A_{00}^{(i)}}\bigl(A_{10}^{(i)} \overline{Z}_{xx_{p,q}}^{(i)} + C_{10}^{(i)} \overline{Z}_{yy_{p,q}}^{(i)}\bigr), \end{aligned}$$
(5.9a)
$$\begin{aligned} \hat{\hat{Z}}_{y_{p,q}}^{(i)} = {}&\overline{Z}_{y_{p,q}}^{(i)} - \frac{h}{8C_{00}^{(i)}}\bigl(\overline{K}_{p,q + 1}^{(i)} - \overline{K}_{p,q - 1}^{(i)}\bigr) + \frac{hA_{00}^{(i)}}{8C_{00}^{(i)}}\bigl( \overline{Z}_{xx_{p,q + 1}}^{(i)} - \overline{Z}_{xx_{p,q - 1}}^{(i)} \bigr) \\ &{}+ \frac{h^{2}}{4C_{00}^{(i)}}\bigl(A_{01}^{(i)} \overline{Z}_{xx_{p,q}}^{(i)} + C_{01}^{(i)} \overline{Z}_{yy_{p,q}}^{(i)}\bigr) \end{aligned}$$
(5.9b)

and let

$$\begin{aligned} &\overline{\overline{K}}_{p,q}^{(i)} = k^{(i)} \bigl(x_{p},y_{q},Z_{p,q}^{(1)},Z_{p,q}^{(2)},\ldots,Z_{p,q}^{(n)}, \overline{\overline{Z}}_{x_{p,q}}^{(1)},\overline{ \overline{Z}}_{x_{p,q}}^{(2)},\ldots,\overline{\overline{Z}}_{x_{p,q}}^{(n)} \overline{\overline{Z}}_{y_{p,q}}^{(1)},\overline{ \overline{Z}}_{y_{p,q}}^{(2)},\ldots,\overline{\overline{Z}}_{y_{p,q}}^{(n)} \bigr), \end{aligned}$$
(5.10)
$$\begin{aligned} &\hat{\hat{K}}_{p,q}^{(i)} = k^{(i)} \bigl(x_{p},y_{q},Z_{p,q}^{(1)},Z_{p,q}^{(2)}, \ldots,Z_{p,q}^{(n)},\hat{\hat{Z}}_{x_{p,q}}^{(1)}, \hat{\hat{Z}}_{x_{p,q}}^{(2)},\ldots,\hat{\hat{Z}}_{x_{p,q}}^{(n)}, \hat{\hat{Z}}_{y_{p,q}}^{(1)},\hat{\hat{Z}}_{y_{p,q}}^{(2)}, \ldots,\hat{\hat{Z}}_{y_{p,q}}^{(n)}\bigr). \end{aligned}$$
(5.11)

We denote

$$\begin{aligned} &P_{1}^{(i)} = A_{10}^{(i)}/A_{00}^{(i)}, \qquad P_{2}^{(i)} = C_{01}^{(i)}/C_{00}^{(i)}, \\ &I_{1}^{(i)} = A_{00}^{(i)} + \frac{1}{12}h^{2}\bigl(A_{20}^{(i)} + A_{02}^{(i)} - 2P_{1}^{(i)}A_{10}^{(i)} - 2P_{2}^{(i)}A_{01}^{(i)}\bigr), \\ &I_{2}^{(i)} = C_{00}^{(i)} + \frac{1}{12}h^{2}\bigl(C_{20}^{(i)} + C_{02}^{(i)} - 2P_{1}^{(i)}C_{10}^{(i)} - 2P_{2}^{(i)}C_{01}^{(i)}\bigr), \\ &I_{3}^{(i)} = \frac{1}{12}h\bigl(A_{01}^{(i)} - P_{2}^{(i)}A_{00}^{(i)}\bigr),\qquad I_{4}^{(i)} = \frac{1}{12}h\bigl(C_{10}^{(i)} - P_{1}^{(i)}C_{00}^{(i)}\bigr),\qquad I_{5}^{(i)} = \frac{1}{12}\bigl(A_{00}^{(i)} + C_{00}^{(i)}\bigr), \\ &J_{1}^{(i)} = 1 - hP_{1}^{(i)},\qquad J_{2}^{(i)} = 1 + hP_{1}^{(i)},\qquad J_{3}^{(i)} = 1 - hP_{2}^{(i)},\qquad J_{4}^{(i)} = 1 + hP_{2}^{(i)}. \end{aligned}$$

Assume that \(K_{p,q}^{(i)} \ne 0\), \(i = 1(1) N\). Then at each nodal point \(( x_{ p}, y_{ q})\) the given system of differential equations (5.1) is discretized by

$$\begin{aligned} & \bigl[I_{1}^{(i)}\delta_{x}^{2} + I_{2}^{(i)}\delta_{y}^{2} + I_{3}^{(i)}\bigl(2\delta_{x}^{2} \mu_{y}\delta_{y}\bigr) + I_{4}^{(i)} \bigl(2\delta_{y}^{2}\mu_{x}\delta_{x} \bigr) + I_{5}^{(i)}\bigl(\delta_{x}^{2} \delta_{y}^{2}\bigr)\bigr]Z_{p,.q}^{(i)} \\ & \quad =h^{2}\overline{\overline{K}}_{p,q}^{(i)} \exp \biggl(\frac{J_{1}^{(i)}\overline{K}_{p + 1,q}^{(i)} + J_{2}^{(i)}\overline{K}_{p - 1,q}^{(i)} + J_{3}^{(i)}\overline{K}_{p,q + 1}^{(i)} + J_{4}^{(i)}\overline{K}_{p,q - 1}^{(i)} - 4\hat{\hat{K}}_{p,q}^{(i)}}{12\overline{\overline{K}}_{p,q}^{(i)}}\biggr) \\ &\qquad {} + O\bigl(h^{6} \bigr). \end{aligned}$$
(5.12)

Using the technique discussed in [19], we can derive the fourth order schemes for the system of quasilinear elliptic PDEs.

6 Application to elliptic equations in polar coordinates

In this section, we aim to derive a stable finite difference scheme for a class of two-dimensional quasilinear elliptic equations and ensure that the numerical methods developed here retain their order and accuracy everywhere including the region in the vicinity of the singularity \(x= 0\).

Let us consider the elliptic equation of the form

$$ z_{xx} + C(x)z_{yy} = D(x)z_{x} + G(x, y),\quad 0 < x,y < 1, $$
(6.1)

where \(D(x) = \frac{1}{x}\).

In this section, we denote

$$ S_{lm} = \frac{\partial^{l+m} S}{\partial x^{l} \partial y^{m}}\quad \mbox{for } l, m=0,1,2,\dots\mbox{ and for } S=C, D, G. $$
(6.2)

With the help of approximations (2.3a)–(2.11) and using the method (2.12), we obtain the following difference scheme for Eq. (6.1):

$$\begin{aligned} &\bigl[L_{1}\delta_{x}^{2} + L_{2} \delta_{y}^{2} + L_{3}\bigl(2\delta_{y}^{2} \mu_{x}\delta_{x}\bigr) + L_{4}\bigl( \delta_{x}^{2}\delta_{y}^{2}\bigr) \bigr]z_{p,q} \\ &\quad =\frac{h^{2}}{12}\biggl[8D_{00}\overline{z}_{x_{p,q}} + \biggl(1 - \frac{hD_{00}}{2}\biggr)D_{l + 1}\overline{z}_{x_{p + 1,q}} + \biggl(1 + \frac{hD_{00}}{2}\biggr)D_{p - 1}\overline{z}_{x_{p - 1,q}} \\ &\qquad {}+ D_{00}(\overline{z}_{x_{p,q + 1}} + \overline{z}_{x_{p,q - 1}}) + \frac{hD_{00}C_{00}}{2}(\overline{z}_{yy_{p + 1,q}} - \overline{z}_{yy_{p - 1,q}}) + h^{2}C_{10}D_{00} \overline{z}_{yy_{p,q}}\biggr] \\ &\qquad {}+ \frac{h^{2}}{12}\biggl[8G_{p,q} + G_{p + 1,q} + G_{p - 1,q} + G_{p,q + 1} + G_{p,q - 1} - \frac{hD_{00}}{2}(G_{p + 1,q} - G_{p - 1,q})\biggr], \\ &\quad p,q= 1,2,\dots, N, \end{aligned}$$
(6.3)

where

$$L_{1} = 1,\qquad L_{2} = C_{00} + \frac{h^{2}}{12}C_{20},\qquad L_{3} = \frac{h}{12}C_{10}, \qquad L_{4} = \frac{1 + C_{00}}{12}. $$

Note that the scheme (6.3) is of \(O(h ^{4})\) for the difference solution of (6.1). However, the scheme fails when the solution is to be determined at \(p=1\). We overcome this difficulty by modifying the method in such a way that the solutions retain the order and accuracy even in the vicinity of the singularity \(x=0\).

We consider the following approximations:

$$\begin{aligned} &D_{p \pm 1} = D_{p} \pm hD_{10} + \frac{h^{2}}{2}D_{20} + O\bigl( \pm h^{3} + h^{4} \bigr), \end{aligned}$$
(6.4a)
$$\begin{aligned} &G_{p \pm 1,q} = G_{00} \pm hG_{10} + \frac{h^{2}}{2}G_{20} + O\bigl( \pm h^{3} + h^{4} \bigr), \end{aligned}$$
(6.4b)
$$\begin{aligned} &G_{p,q \pm 1} = G_{00} \pm hG_{01} + \frac{h^{2}}{2}G_{02} + O\bigl( \pm h^{3} + h^{4} \bigr). \end{aligned}$$
(6.4c)

Now using the approximations (6.4a)–(6.4c) and neglecting the higher order terms, we get

$$\begin{aligned} & \bigl[L_{1}\delta_{x}^{2} + L_{2} \delta_{y}^{2} + L_{3}\bigl(2\delta_{y}^{2} \mu_{x}\delta_{x}\bigr) + L_{4}\bigl( \delta_{x}^{2}\delta_{y}^{2}\bigr) \bigr]z_{p,q} \\ &\quad = \bigl[M_{1}\delta_{x}^{2} + M_{2}\delta_{y}^{2} + M_{3}\bigl(2 \delta_{y}^{2}\mu_{x}\delta_{x}\bigr) + M_{4}\bigl(\delta_{x}^{2}\delta_{y}^{2} \bigr)\bigr]z_{p,q} + \sum G, \end{aligned}$$
(6.5)

where

$$\begin{aligned} &M_{1} = \frac{h^{2}}{12}\bigl(2D_{10} - D_{00}^{2}\bigr),\quad M_{2} = \frac{h^{2}}{12}(D_{00}C_{10}), \\ &M_{3} = \frac{h^{2}}{24}\biggl(\frac{12}{h}D_{00} + hD_{20} - hD_{00}D_{10}\biggr),\qquad M_{4} = h\frac{D_{00}}{24}(1 + c_{00}), \\ &\sum G = \frac{h^{2}}{12}\bigl[12G_{00} + h^{2}(G_{20} + G_{02}) - h^{2}D_{00}G_{10} \bigr]. \end{aligned}$$

Now consider the two dimensional Poisson equation in the rθ plane,

$$ z_{rr} + \frac{1}{r}z_{r} + \frac{1}{r^{2}}z_{\theta \theta} = G(r, \theta ), \quad 0 < r,\theta < 1. $$
(6.6)

Replacing the variables x, y by r, θ, respectively, we get the required difference scheme for the solution of differential equation (6.6) and the coefficients are given by

$$\begin{aligned} & L_{1} = 1,\qquad L_{2} = \frac{1}{r^{2}}\biggl(1 + \frac{h^{2}}{2r^{2}}\biggr),\qquad L_{3} = \frac{ - h}{6r^{3}},\qquad L_{4} = \frac{1 + r^{2}}{12r^{2}}, \\ &M_{1} = \frac{h^{2}}{12r^{2}},\qquad M_{2} = \frac{h^{2}}{6r^{4}}, \\ &M_{3} = - \frac{h}{2r}\biggl(1 + \frac{h^{2}}{12r^{2}}\biggr), \qquad M_{4} = - \frac{h}{24r^{3}}\bigl(1 + r^{2}\bigr), \\ &\sum G = \frac{h^{2}}{12}\biggl[12G_{00} + h^{2}(G_{20} + G_{02}) + \frac{h^{2}}{r}G_{10} \biggr]. \end{aligned}$$

Next, we consider the two dimensional Poisson equation in the rw plane,

$$ z_{rr} + \frac{1}{r}z_{r} + z_{ww} = G(r,w),\quad 0 < r,w < 1. $$
(6.7)

Replacing the variables x, y by r, w, respectively, we get the required difference scheme for the solution of the differential equation (6.7) and the coefficients are given by

$$\begin{aligned} & L_{1} = 1,\qquad L_{2} = 1,\qquad L_{3} = 0, \qquad L_{4} = \frac{1}{6}, \\ &M_{1} = \frac{h^{2}}{12r^{2}},\qquad M_{2} = 0,\qquad M_{3} = - \frac{h}{2r}\biggl(1 + \frac{h^{2}}{12r^{2}}\biggr),\qquad M_{4} = - \frac{h}{12r}, \\ &\sum G = \frac{h^{2}}{12}\biggl[12G_{00} + h^{2}(G_{20} + G_{02}) + \frac{h^{2}}{r}G_{10} \biggr]. \end{aligned}$$

Note that the scheme (6.3) along with the approximations (6.4a)–(6.4c) are of \(O(h ^{4})\) and free from the terms \(1/(p \pm1)\) and \(1/(q \pm1)\), hence it is very easily solved for \(p,q= 1, 2,\ldots,N\) in the region \(0 < r,\theta < 1\) and \(0 < r,w < 1\). In a similar manner, we can discuss the numerical schemes for nonlinear elliptic equations in polar coordinates.

7 Application to nonlinear bi- and triharmonic problems

7.1 Nonlinear biharmonic equation

We consider the 2D nonlinear biharmonic elliptic partial differential equation with a forcing function of the form

$$ \nabla^{4}z(x,y) = f\bigl(x,y,z,\nabla^{2}z,z_{x}, \nabla^{2}z_{x},z_{y},\nabla^{2}z_{y} \bigr),\quad \varGamma:0 < x,y < 1, $$
(7.1)

where the Laplacian operator \(\nabla^{2}\) is defined by \(\nabla^{2}z \equiv \frac{\partial^{2}z}{\partial x^{2}} + \frac{\partial^{2}z}{\partial y^{2}}\).

The boundary conditions of second kind z and (\(\partial^{2}z/\partial n^{2}\)) are prescribed on the boundary.

The values of z, \(z _{yy}\) are prescribed on the boundary \(y = 0\), \(y = 1\); and the values of z, \(z _{xx}\) are prescribed on the boundary \(x = 0\), \(x = 1\). As the grid lines are parallel to the coordinate axes and the values of zare exactly known on the boundary, this implies that the successive tangential partial derivatives of z are known exactly on the boundary. For example, on the line \(y = 0\), the values of \(z(x,0)\) and \(z _{yy}(x,0)\) are known, i.e. the values of \(z _{x}(x,0)\), \(z _{xx}(x,0)\), etc. are known on the line \(y = 0\). This implies the values of \(z(x,0)\) and \(\nabla^{2}z(x,0) \equiv z_{xx}(x,0) + z_{yy}(x,0)\) are known on the line \(y = 0\). Similarly, the values of z and \(\nabla^{2}z\) are known on all sides of the square region.

Let us denote \(\nabla^{2}z=v\). Then we can rewrite Eq. (7.1) in a coupled manner as

$$\begin{aligned} &\nabla^{2}z = v(x,y),\quad (x,y) \in \varGamma, \end{aligned}$$
(7.2a)
$$\begin{aligned} &\nabla^{2}v = f(x,y,z,v,z_{x},v_{x},z_{y},v_{y}), \quad (x,y) \in \varGamma. \end{aligned}$$
(7.2b)

In this case, the values of z and vare exactly known on the boundary of Γ.

Applying the method (5.12) to the system of equations (7.2a)–(7.2b), a numerical method of order four for the solution of the biharmonic equation (7.1) can be written as

$$\begin{aligned} &\biggl[\delta_{x}^{2} + \delta_{y}^{2} + \frac{1}{6}\bigl(\delta_{x}^{2} \delta_{y}^{2}\bigr)\biggr]Z_{p,q} \\ &\quad = h^{2}V_{p,q}\exp \biggl( \frac{V_{{p + 1,q}} + V_{{p - 1,q}} + V_{{p,q + 1}} + V_{{p,q - 1}} - 4V_{{p,q}}}{12V_{p,q}} \biggr) + O\bigl(h^{6}\bigr), \\ &\quad 1\leq p,q \leq N, \end{aligned}$$
(7.3a)
$$\begin{aligned} &\biggl[\delta_{x}^{2} + \delta_{y}^{2} + \frac{1}{6}\bigl(\delta_{x}^{2} \delta_{y}^{2}\bigr)\biggr]V_{p,q} \\ &\quad =h^{2}\overline{\overline{F}}_{p,q}\exp \biggl( \frac{\overline{F}_{{p + 1,q}} + \overline{F}_{p - 1,q} + \overline{F}_{p,q + 1} + \overline{F}_{{p,q - 1}} - 4\hat{\hat{F}}_{p,q}}{12\overline{\overline{F}}_{p,q}} \biggr) + O\bigl(h^{6}\bigr), \\ &\quad 1\leq p,q \leq N, \end{aligned}$$
(7.3b)

where

$$\begin{aligned} &\overline{F}_{p \pm 1,q} = f(x_{p \pm 1},y_{q,}z_{p \pm 1,q},v_{p \pm 1,q}, \bar{z}_{x_{p \pm 1,q}},\bar{v}_{x_{p \pm 1,q}},\bar{z}_{y_{p \pm 1,q}}, \bar{v}_{y_{p \pm 1,q}}), \end{aligned}$$
(7.4a)
$$\begin{aligned} &\overline{F}_{p,q \pm 1} = f(x_{p},y_{q \pm 1,}z_{p,q \pm 1},v_{p,q \pm 1}, \bar{z}_{x_{p,q \pm 1}},\bar{v}_{x_{p,q \pm 1}},\bar{z}_{y_{p,q \pm 1}}, \bar{v}_{y_{p,q \pm 1}}), \end{aligned}$$
(7.4b)
$$\begin{aligned} &\bar{\bar{F}}_{p,q} = f(x_{p},y_{q,}z_{p,q},v_{p,q}, \bar{\bar{z}}_{x_{p,q}},\bar{\bar{v}}_{x_{p,q}},\bar{ \bar{z}}_{y_{p,q}},\bar{\bar{v}}_{y_{p,q}}), \end{aligned}$$
(7.4c)
$$\begin{aligned} &\hat{\hat{F}}_{p,q} = f(x_{p},y_{q,}z_{p,q},v_{p,q}, \hat{\hat{z}}_{x_{p,q}},\hat{\hat{v}}_{x_{p,q}},\hat{ \hat{z}}_{y_{p,q}},\hat{\hat{v}}_{y_{p,q}}). \end{aligned}$$
(7.4d)

The approximations associated with (7.4a)–(7.4d) are already defined in Sect. 5.

7.2 Nonlinear triharmonic equation

Next we consider the nonlinear triharmonic equation with a forcing function of the form

$$ \nabla^{6}z(x,y) = g\bigl(x,y,z,\nabla^{2}z, \nabla^{4}z,z_{x},\nabla^{2}z_{x}, \nabla^{4}z_{x},z_{y},\nabla^{2}z_{y}, \nabla^{4}z_{y}\bigr),\quad \varGamma:0 < x,y < 1. $$
(7.5)

The values of z, (\(\partial^{2}z/\partial n^{2}\)) and (\(\partial^{4}z/\partial n^{4}\)) are known on the boundary.

For Eq. (7.5), the boundary values of \(z, z _{yy}\), \(z _{\mathit{yyyy}}\) are prescribed on the line \(y = 0\), \(y = 1\); and the boundary values of \(z, z _{xx}\), \(z _{\mathit{xxxx}}\) are prescribed on the line \(x = 0\), \(x = 1\). As discussed in the biharmonic case, the values of z, \(\nabla^{2}z\) and \(\nabla^{4}z\) are known on all sides of the square region Γ.

Let \(\nabla^{2}z =v\) and \(\nabla^{2}v =w\). Then we rewrite Eq. (7.5) in a split form as

$$\begin{aligned} &\nabla^{2}z =v(x,y),\quad (x,y) \in \varGamma, \end{aligned}$$
(7.6a)
$$\begin{aligned} &\nabla^{2}v = w(x,y),\quad (x,y) \in \varGamma, \end{aligned}$$
(7.6b)
$$\begin{aligned} &\nabla^{2}w = g(x,y,z,v,w,z_{x},v_{x},w_{x},z_{y},v_{y},w_{y}), \quad (x,y) \in \varGamma. \end{aligned}$$
(7.6c)

Applying the method (5.12) to the system of equations (7.6a)–(7.6c), a numerical method of order four for the solution of triharmonic equation (7.5) can be written as

$$\begin{aligned} &\biggl[\delta_{x}^{2} + \delta_{y}^{2} + \frac{1}{6}\bigl(\delta_{x}^{2} \delta_{y}^{2}\bigr)\biggr]Z_{p,q} \\ &\quad = h^{2}V_{p,q}\exp \biggl( \frac{V_{{p + 1,q}} + V_{{p - 1,q}} + V_{{p,q + 1}} + V_{{p,q - 1}} - 4V_{{p,q}}}{12V_{p,q}} \biggr) + O\bigl(h^{6}\bigr), \\ &\quad 1\leq p,q \leq N, \end{aligned}$$
(7.7a)
$$\begin{aligned} &\biggl[\delta_{x}^{2} + \delta_{y}^{2} + \frac{1}{6}\bigl(\delta_{x}^{2} \delta_{y}^{2}\bigr)\biggr]V_{p,q} \\ &\quad =h^{2}W_{p,q}\exp \biggl( \frac{W_{{p + 1,q}} + W_{{p - 1,q}} + W_{{p,q + 1}} + W_{{p,q - 1}} - 4W_{{p,q}}}{12W_{p,q}} \biggr) + O\bigl(h^{6}\bigr), \\ &\quad 1\leq p,q \leq N, \end{aligned}$$
(7.7b)
$$\begin{aligned} &\biggl[\delta_{x}^{2} + \delta_{y}^{2} + \frac{1}{6}\bigl(\delta_{x}^{2} \delta_{y}^{2}\bigr)\biggr]W_{p,q} \\ &\quad = h^{2}\overline{\overline{G}}_{p,q}\exp \biggl( \frac{\overline{G}_{{p + 1,q}} + \overline{G}_{{p - 1,q}} + \overline{G}_{{p,q + 1}} + \overline{G}_{{p,q - 1}} - 4\widehat{\widehat{G}}_{p,q}}{12\overline{\overline{G}}_{p,q}} \biggr) + O\bigl(h^{6}\bigr), \\ &\quad 1\leq p,q \leq N, \end{aligned}$$
(7.7c)

where

$$\begin{aligned} &\overline{G}_{p \pm 1,q} = g(x_{p \pm 1},y_{q,}z_{p \pm 1,q},v_{p \pm 1,q},w_{p \pm 1,q}, \bar{z}_{x_{p \pm 1,q}},\bar{v}_{x_{p \pm 1,q}}, \\ &\hphantom{\overline{G}_{p \pm 1,q} = }\bar{w}_{x_{p \pm 1,q}}, \bar{z}_{y_{p \pm 1,q}},\bar{v}_{y_{p \pm 1,q}},\bar{w}_{y_{p \pm 1,q}}), \end{aligned}$$
(7.8a)
$$\begin{aligned} &\overline{G}_{p,q \pm 1} = g(x_{p},y_{q \pm 1,}z_{p,q \pm 1},v_{p,q \pm 1},w_{p,q \pm 1}, \bar{z}_{x_{p,q \pm 1}},\bar{v}_{x_{p,q \pm 1}}, \\ &\hphantom{\overline{G}_{p,q \pm 1} =}\bar{w}_{x_{p,q \pm 1}}, \bar{z}_{y_{p,q \pm 1}},\bar{v}_{y_{p,q \pm 1}},\bar{w}_{y_{p,q \pm 1}}), \end{aligned}$$
(7.8b)
$$\begin{aligned} &\bar{\bar{G}}_{p,q} = f(x_{p},y_{q,}z_{p,q},v_{p,q},w_{p,q}, \bar{\bar{z}}_{x_{p,q}},\bar{\bar{v}}_{x_{p,q}},\bar{ \bar{w}}_{x_{p,q}},\bar{\bar{z}}_{y_{p,q}},\bar{\bar{v}}_{y_{p,q}}, \bar{\bar{w}}_{y_{p,q}}), \end{aligned}$$
(7.8c)
$$\begin{aligned} &\hat{\hat{G}}_{p,q} = g(x_{p},y_{q,}z_{p,q},v_{p,q},w_{p,q} \hat{\hat{z}}_{x_{p,q}},\hat{\hat{v}}_{x_{p,q}},\hat{ \hat{w}}_{x_{p,q}},\hat{\hat{z}}_{y_{p,q}},\hat{\hat{v}}_{y_{p,q}}, \hat{\hat{w}}_{y_{p,q}}). \end{aligned}$$
(7.8d)

The approximations associated with (7.8a)–(7.8d) are already defined in Sect. 5.

With the help of boundary values, writing all methods at every interior grid point, one obtains sparse systems of linear algebraic equations for the solution of the multiharmonic equations (7.1), (7.5). A direct solution of these linear systems is impractical because of the large size of the coefficient matrix and enormous storage requirements even for moderate values of the grid size. Classical block iterative methods such as the Gauss–Seidel and successive over relaxation methods are attractive for their low storage requirements as long as convergence is guaranteed.

8 Numerical illustrations

We implement the proposed method on ten benchmark problems of 2D elliptic equations. The exact solution of each problem is given. The right hand side function and Dirichlet boundary conditions are determined from the exact solution in each problem. The system of linear difference equations are solved using the Gauss–Seidel iteration method and that of nonlinear difference equations by the Newton–Raphson iteration method [36, 37]. The iterations are terminated once the absolute error tolerance ≤10−12 is reached and the initial guess made in each case is \(\mathbf{z} = \mathbf{0}\). All computations were performed using MATLAB codes.

Problem 8.1

(Convection–diffusion equation)

$$ z_{xx} + z_{yy} = \beta z_{x},\quad 0 < x,y < 1. $$
(8.1)

The exact solution is given by \(z(x,y) =e^{\beta x/2}\frac{\sin (\pi y)}{\sinh (\sigma )}[2e^{ - \beta /2}\sinh (\sigma x) + \sinh (\sigma (1 - x))]\), \(\sigma = \sqrt{\pi^{2} + \frac{\beta^{2}}{4}}\).

The maximum absolute errors (MAEs) in z are listed in Table 1. Figures 2(a) and (b) give the plots of the exact and numerical solutions for \(h=1/64\) and \(\beta=50\).

Figure 2
figure 2

(a) Exact solution of Problem 8.1 for \(h = 1/64\) and \(\beta =50\). (b) Numerical solution of Problem 8.1 for \(h = 1/64\) and \(\beta= 50\)

Table 1 Problem 8.1: The maximum absolute errors

Problem 8.2

(Poisson’s equation in rθ plane)

$$ z_{rr} + \frac{\alpha}{r}z_{r} + \frac{1}{r^{2}}z_{\theta \theta} = G(r,\theta ),\quad 0 < r,\theta < 1. $$
(8.2)

The exact solution is \(z(r,\theta ) = r^{2}\cos (\pi \theta )\). The MAEs in z are listed in Table 2 for \(\alpha=1\) and 2. Figures 3(a) and (b) give the plots of the exact and numerical solutions for \(h=1/32\) and \(\alpha = 2\).

Figure 3
figure 3

(a) Exact solution of Problem 8.2 for \(h = 1/32\) and \(\alpha = 2 \). (b) Numerical solution of Problem 8.2 for \(h = 1/32\) and \(\alpha = 2\)

Table 2 Problem 8.2: The maximum absolute errors

Problem 8.3

(Poisson’s equation in rw plane)

$$ z_{rr} + \frac{\alpha}{r}z_{r} + z_{ww} = G(r,w),\quad 0 < r,w < 1. $$
(8.3)

The exact solution is \(z(r,w) = \cosh r\cosh w\). The MAEs in z are listed in Table 3 for \(\alpha = 1\) and 2. Figures 4(a) and (b) give the plots of the exact and numerical solutions for \(h=1/64\) and \(\alpha = 2\).

Figure 4
figure 4

(a) Exact solution of Problem 8.3 for \(h= 1/64\) and \(\alpha=2\). (b) Numerical solution of Problem 8.3 for \(h = 1/64\) and \(\alpha= 2\)

Table 3 Problem 8.3: The maximum absolute errors

Problem 8.4

(Burgers’ equation)

$$ \varepsilon (z_{xx} + z_{yy}) = z(z_{x} + z_{y}) + g(x,y),\quad 0 < x,y < 1. $$
(8.4)

The exact solution is \(z(x,y) = e^{x}\sin ( \frac{\pi y}{2} )\). The MAEs in z are listed in Table 4 for \(\varepsilon = 0.1,0.01\) and 0.001. Figures 5(a) and (b) give the plots of the exact and numerical solutions for \(h=1/32\) and \(\varepsilon = 0.01\).

Figure 5
figure 5

(a) Exact solution of Problem 8.4 for \(h = 1/32\) and \(\varepsilon=0.01\). (b) Numerical solution of Problem 8.4 for \(h = 1/32\) and \(\varepsilon=0.01\)

Table 4 Problem 8.4: The maximum absolute errors

Problem 8.5

(2D steady-state Navier–Stokes model equations in rectangular coordinates)

$$\begin{aligned} &\frac{1}{\boldsymbol{R}_{\boldsymbol{e}}}(z_{xx} + z_{yy}) = zz_{x} + vz_{y} + f(x,y),\quad 0< x, y< 1, \end{aligned}$$
(8.5a)
$$\begin{aligned} &\frac{1}{\boldsymbol{R}_{\boldsymbol{e}}}(v_{xx} + v_{yy}) = zv_{x} + vv_{y} + g(x,y),\quad 0< x, y< 1. \end{aligned}$$
(8.5b)

The exact solutions are \(z(x,y) = \sin (\pi x)\sin (\pi y)\), \(v(x,y) = \cos (\pi x)\cos (\pi y)\).

The MAEs in z, v are tabulated in Table 5 for various values of the Reynolds number \(R _{e}\). Figures 6(a) and (b) give the plots of the exact and numerical solutions of z and Figs. 6(c) and (d) give the plots of the exact and numerical solutions of v for \(h=1/64\), \(R _{e}=10\).

Figure 6
figure 6figure 6

(a) Exact solution of Problem 8.5 for \(z(x,y)\) for \(h = 1/64\) and \(\operatorname{Re} = 10\). (b) Numerical solution of Problem 8.5 for \(z(x,y)\) for \(h = 1/64\) and \(\operatorname{Re} = 10\). (c) Exact solution of Problem 8.5 for \(v(x,y)\) for \(h = 1/64\) and \(\operatorname{Re} = 10\). (d) Numerical solution of Problem 8.5 for \(v(x,y)\) for \(h = 1/64\) and \(\operatorname{Re} = 10\)

Table 5 Problem 8.5: The maximum absolute errors

Problem 8.6

(2D steady-state Navier–Stokes model equations in cylindrical polar coordinates in rw plane)

$$\begin{aligned} &\frac{1}{\boldsymbol{R}_{\boldsymbol{e}}} \biggl( z_{rr} + \frac{1}{r}z_{r} + z_{ww} - \frac{1}{r^{2}}z \biggr) = zz_{r} + vz_{w} + H(r,w), \quad 0< r, w< 1, \end{aligned}$$
(8.6a)
$$\begin{aligned} &\frac{1}{\boldsymbol{R}_{\boldsymbol{e}}} \biggl( v_{rr} + \frac{1}{r}v_{r} + v_{zz} \biggr) = zv_{r} + vv_{w} + I(r,w), \quad 0< r, w< 1. \end{aligned}$$
(8.6b)

The exact solutions are given by \(z(r,w) = r^{3}\sinh w\), \(v(r,w) = - 4r^{2}\cosh w\). The MAEs in z, v are listed in Table 6 for various values \(R _{e}\). Figures 7(a) and (b) give the plots of the exact and numerical solutions of z, and Figs. 7(c) and (d) give the plot of the exact and numerical solutions of v for \(h=1/32\) and \(R _{e}=10\).

Figure 7
figure 7figure 7

(a) Exact solution of Problem 8.6 for \(z(x,y)\) for \(h = 1/32\) and \(\operatorname{Re} = 10\). (b) Numerical solution of Problem 8.6 for \(z(x,y)\) for \(h = 1/32\) and \(\operatorname{Re} = 10\). (c) Exact solution of Problem 8.6 for \(v(x,y)\) for \(h = 1/32\) and \(\operatorname{Re} = 10\). (d) Numerical solution of Problem 8.6 for \(v(x,y)\) for \(h = 1/32\) and \(\operatorname{Re} = 10\)

Table 6 Problem 8.6: The maximum absolute errors

Problem 8.7

(Nonlinear elliptic equation)

$$ \bigl(1 + x^{2}\bigr)z_{xx} + \bigl(1 + y^{2} \bigr)z_{yy} = \alpha z(z_{x} + z_{y}) + f(x,y), \quad 0 < x,y < 1. $$
(8.7)

The exact solution is \(z(x,y) = e^{x}\cos (\pi y)\). The MAEs in z are tabulated in Table 7 for various values of α. Figures 8(a) and (b) give the plot of the exact and numerical solutions of z for the values of \(h=1/64\) and \(\alpha=1\).

Figure 8
figure 8

(a) Exact solution of Problem 8.7 for \(h = 1/64\) and \(\alpha =1\). (b) Numerical solution of Problem 8.7 for \(h = 1/64\) and \(\alpha= 1\)

Table 7 Problem 8.7: The maximum absolute errors

Problem 8.8

(Quasilinear elliptic equation)

$$ z_{xx} + \bigl(1 + z^{2}\bigr)z_{yy} = \alpha z(z_{x} + z_{y}) + f(x,y),\quad 0 < x,y < 1. $$
(8.8)

The exact solution is \(z(x,y) = e^{x}\cos (\pi y)\). The MAEs in z are tabulated in Table 8 for various values of α. Figures 9(a) and (b) give the plots of the exact and numerical solutions of z for \(h=1/64\) and \(\alpha = 1\).

Figure 9
figure 9

(a) Exact solution of Problem 8.8 for \(h = 1/32\) and \(\operatorname{Re} = 10\). (b) Numerical solution of Problem 8.8 for \(h = 1/64\) and \(\alpha = 1\)

Table 8 Problem 8.8: The maximum absolute errors

Problem 8.9

(Nonlinear biharmonic equation)

$$ \nabla^{4}z = \alpha z\bigl(z_{x} + z_{y} + \nabla^{2}z_{x} + \nabla^{2}z_{y}\bigr) + f(x,y),\quad 0 < x,y < 1. $$
(8.9)

The exact solution is \(z(x,y) = \sin (\pi x)\cos (\pi y)\). The MAEs in z are tabulated in Table 9 for various values of α. Figures 10(a) and (b) give the plots of the exact and numerical solutions of z for \(h=1/64\) and \(\alpha = 1\).

Figure 10
figure 10

(a) Exact solution of Problem 8.9 for \(h = 1/32\) and \(\alpha = 10\). (b) Numerical solution of Problem 8.9 for \(h = 1/32\) and \(\alpha = 10\)

Table 9 Problem 8.9: The maximum absolute errors

Problem 8.10

(Nonlinear triharmonic Equation)

$$ \nabla^{6}z = \alpha z\bigl(z_{x} + z_{y} + \nabla^{2}z_{x} + \nabla^{2}z_{y} + \nabla^{4}z_{x} + \nabla^{4}z_{y}\bigr) + f(x,y),\quad 0 < x,y < 1. $$
(8.10)

The exact solution is \(z(x,y) = \cos (\pi x)\sin (\pi y)\). The MAEs in z are tabulated in Table 10 for various values of α. Figures 11(a) and (b) give the plots of the exact and numerical solutions of z for \(h=1/64\) and \(\alpha = 1\).

Figure 11
figure 11

(a) Exact solution of Problem 8.10 for \(h = 1/64\) and \(\alpha = 10\). (b) Numerical solution of Problem 8.10 for \(h = 1/64\) and \(\alpha = 10\)

Table 10 Problem 8.10: The maximum absolute errors

9 Conclusions

We have derived a new 9-point compact fourth order numerical method in exponential form for the numerical solution of the system of 2D quasilinear elliptic partial differential equations. The advantage of the exponential method is that it gives higher accuracy results than existing methods in the literature. The theoretical convergence of elliptic equations with constant coefficients has been established. Many benchmark problems like the diffusion–convection equation, Poisson’s equation in polar coordinates, the Navier–Stokes equations of motion both in rectangular and polar coordinates, nonlinear bi- and triharmonic equations, and quasilinear elliptic equations have been solved and compared with the results of existing methods. The theoretical rate of convergence was corroborated by the experimentally observed rate of convergence using the formula \(\rho = \log (e_{h_{1}} / e_{h_{2}}) / \log (h_{1} / h_{2})\), where \(e_{h_{1}}\) and \(e_{h_{2}}\) are the maximum absolute errors for two mesh sizes \(h_{1}\) and \(h_{2}\) respectively. Computational results exhibit that the proposed methods avoid spurious oscillation and give stable results for high Reynolds number even in the vicinity of the singularity. From the results we conclude that the proposed method is competitive in solving the two-dimensional problems and it can be extended to solving three-dimensional problems.

References

  1. Jain, M.K., Jain, R.K., Mohanty, R.K.: Fourth order difference methods for the system of 2D nonlinear elliptic partial differential equations. Numer. Methods Partial Differ. Equ. 7, 227–244 (1991)

    Article  Google Scholar 

  2. Jain, M.K., Jain, R.K., Mohanty, R.K.: A fourth order difference method for elliptic equations with nonlinear first derivative terms. Numer. Methods Partial Differ. Equ. 5, 87–95 (1989)

    Article  MathSciNet  Google Scholar 

  3. Li, M., Tang, T., Fornberg, B.: A compact fourth-order finite difference scheme for the steady state incompressible Navier–Stokes equations. Int. J. Numer. Methods Fluids 20, 1137–1151 (1995)

    Article  Google Scholar 

  4. Erturk, E., Gökcöl, C.: Fourth-order compact formulation of Navier–Stokes equations and driven cavity flow at high Reynolds numbers. Int. J. Numer. Methods Fluids 50, 421–436 (2006)

    Article  Google Scholar 

  5. Liu, J., Wang, C.: A fourth order numerical method for the primitive equations formulated in mean vorticity. Commun. Comput. Phys. 4, 26–55 (2008)

    MathSciNet  MATH  Google Scholar 

  6. Ito, K., Qiao, Z.: A high order compact MAC finite difference scheme for the Stokes equations: augmented variable approach. J. Comput. Phys. 227, 8177–8190 (2008)

    Article  MathSciNet  Google Scholar 

  7. Spotz, W.F., Carey, G.F.: High order compact scheme for the steady stream function vorticity equations. Int. J. Numer. Methods Eng. 38, 3497–3512 (1995)

    Article  MathSciNet  Google Scholar 

  8. Carey, G.F.: Computational Grids: Generation, Adaption and Solution Strategies. Taylor & Francis, Washington (1997)

    MATH  Google Scholar 

  9. Yavneh, I.: Analysis of a fourth-order compact scheme for convection diffusion. J. Comput. Phys. 133, 361–364 (1997)

    Article  MathSciNet  Google Scholar 

  10. Zhang, J.: On convergence of iterative methods for a fourth-order discretization scheme. Appl. Math. Lett. 10, 49–55 (1997)

    Article  MathSciNet  Google Scholar 

  11. Birkhoff, G., Lynch, R.E.: The Numerical Solution of Elliptic Problems. SIAM, Philadelphia (1984)

    Book  Google Scholar 

  12. Mohanty, R.K.: Fourth order finite difference methods for the system of 2D nonlinear elliptic equations with variable coefficients. Int. J. Comput. Math. 46, 195–206 (1992)

    Article  Google Scholar 

  13. Böhmer, K.: On finite element methods for fully nonlinear elliptic equations of second order. SIAM J. Numer. Anal. 46(3), 1212–1249 (2008)

    Article  MathSciNet  Google Scholar 

  14. Feng, X., Neilan, M.: Vanishing moment method and moment solution for fully nonlinear second order partial differential equations. J. Sci. Comput. 38, 78–98 (2009)

    Article  MathSciNet  Google Scholar 

  15. Böhmer, K.: Numerical Methods for Nonlinear Elliptic Differential Equations: A Synopsis. Oxford University Press, Oxford (2010)

    Book  Google Scholar 

  16. Jain, M.K., Jain, R.K., Krishna, M.: Fourth order difference method for quasi-linear Poisson equation in cylindrical symmetry. Commun. Numer. Methods Eng. 10, 291–296 (1994)

    Article  Google Scholar 

  17. Ananthakrishnaiah, U., Saldanha, G.: A fourth order finite difference scheme for two-dimensional non-linear elliptic partial differential equations. Numer. Methods Partial Differ. Equ. 11, 33–40 (1995)

    Article  Google Scholar 

  18. Mohanty, R.K.: Order \(h ^{4}\) difference methods for a class of singular two-space dimensional elliptic boundary value problems. J. Comput. Appl. Math. 81, 229–247 (1997)

    Article  MathSciNet  Google Scholar 

  19. Mohanty, R.K., Dey, S.: A new finite difference discretization of order four for (\(\partial u/ \partial n\)) for two-dimensional quasi-linear elliptic boundary value problems. Int. J. Comput. Math. 76, 505–576 (2001)

    Article  MathSciNet  Google Scholar 

  20. Mohanty, R.K., Singh, S.: A new fourth order discretization for singularly perturbed two dimensional non-linear elliptic boundary value problems. Appl. Math. Comput. 175, 1400–1414 (2006)

    MathSciNet  MATH  Google Scholar 

  21. Mohanty, R.K., Setia, N.: A new compact high order off-step discretization for the system of 2D quasi-linear elliptic partial differential equations. Adv. Differ. Equ. 2013, 223 (2013)

    Article  MathSciNet  Google Scholar 

  22. Zhang, J.: On convergence and performance of iterative methods with fourth order compact schemes. Numer. Methods Partial Differ. Equ. 14, 263–280 (1998)

    Article  MathSciNet  Google Scholar 

  23. Saldanha, G.: Technical note: a fourth order finite difference scheme for a system of 2D nonlinear elliptic partial differential equations. Numer. Methods Partial Differ. Equ. 17, 43–53 (2001)

    Article  MathSciNet  Google Scholar 

  24. Arabshahi, S.M.M., Dehghan, M.: Preconditioned techniques for solving large sparse linear systems arising from the discretization of the elliptic partial differential equations. Appl. Math. Comput. 188, 1371–1388 (2007)

    MathSciNet  MATH  Google Scholar 

  25. Mohanty, R.K.: A new high accuracy finite difference discretization for the solution of 2D non-linear biharmonic equations using coupled approach. Numer. Methods Partial Differ. Equ. 26, 931–944 (2010)

    MATH  Google Scholar 

  26. Mohanty, R.K.: Single cell compact finite difference discretizations of order two and four for multi-dimensional triharmonic problems. Numer. Methods Partial Differ. Equ. 26, 1420–1426 (2010)

    MATH  Google Scholar 

  27. Mohanty, R.K., Jain, M.K., Mishra, B.N.: A compact discretization of \(O(h ^{4})\) for two-dimensional non-linear triharmonic equations. Phys. Scr. 84, 025002 (2011)

    Article  Google Scholar 

  28. Mohanty, R.K., Jain, M.K., Mishra, B.N.: A novel numerical method of \(O(h ^{4})\) for three-dimensional non-linear triharmonic equations. Commun. Comput. Phys. 12, 1417–1433 (2012)

    Article  MathSciNet  Google Scholar 

  29. Khattar, D., Singh, S., Mohanty, R.K.: A new coupled approach high accuracy numerical method for the solution of 3D non-linear biharmonic equations. Appl. Math. Comput. 215, 3036–3044 (2009)

    MathSciNet  MATH  Google Scholar 

  30. Singh, S., Khattar, D., Mohanty, R.K.: A new coupled approach high accuracy numerical method for the solution of 2D non-linear biharmonic equations. Neural Parallel Sci. Comput. 17, 239–256 (2009)

    MathSciNet  MATH  Google Scholar 

  31. Singh, S., Singh, S., Mohanty, R.K.: A new high accuracy off-step discretization for the solution of 2D non-linear triharmonic equations. East Asian J. Appl. Math. 03, 228–246 (2013)

    Article  Google Scholar 

  32. Gautschi, W.: Numerical integration of ordinary differential equations based on trigonometric polynomials. Numer. Math. 3, 381–397 (1961)

    Article  MathSciNet  Google Scholar 

  33. Lyche, T.: Chebyshevian multistep methods for ordinary differential equations. Numer. Math. 19, 65–75 (1972)

    Article  MathSciNet  Google Scholar 

  34. Varga, R.S.: Matrix Iterative Analysis. Springer, New York (2000)

    Book  Google Scholar 

  35. Hageman, L.A., Young, D.M.: Applied Iterative Methods. Dover, New York (2004)

    MATH  Google Scholar 

  36. Kelly, C.T.: Iterative Methods for Linear and Non-Linear Equations. SIAM, Philadelphia (1995)

    Book  Google Scholar 

  37. Saad, Y.: Iterative Methods for Sparse Linear Systems. SIAM, Philadelphia (2003)

    Book  Google Scholar 

Download references

Acknowledgements

The authors thank the reviewers for their valuable suggestions, which substantially improved the quality of the paper.

Funding

This work is supported by Maitreyi College, University of Delhi.

Author information

Authors and Affiliations

Authors

Contributions

All authors drafted the manuscript, and they read and approved the final version.

Corresponding author

Correspondence to R. K. Mohanty.

Ethics declarations

Competing interests

The authors declare that they have no competing interests.

Additional information

Publisher’s Note

Springer Nature remains neutral with regard to jurisdictional claims in published maps and institutional affiliations.

Rights and permissions

Open Access This article is distributed under the terms of the Creative Commons Attribution 4.0 International License (http://creativecommons.org/licenses/by/4.0/), which permits unrestricted use, distribution, and reproduction in any medium, provided you give appropriate credit to the original author(s) and the source, provide a link to the Creative Commons license, and indicate if changes were made.

Reprints and permissions

About this article

Check for updates. Verify currency and authenticity via CrossMark

Cite this article

Mohanty, R.K., Manchanda, G. & Khan, A. Operator compact exponential approximation for the solution of the system of 2D second order quasilinear elliptic partial differential equations. Adv Differ Equ 2019, 47 (2019). https://doi.org/10.1186/s13662-019-1968-9

Download citation

  • Received:

  • Accepted:

  • Published:

  • DOI: https://doi.org/10.1186/s13662-019-1968-9

MSC

Keywords