A finite difference method for a numerical solution of elliptic boundary value problems

In this article, we propose a finite difference method for solving system of elliptic boundary value problems in a square region Ω = [a,b]×[a,b],

∂2u∂x2+∂2u∂y2=g(x,y),(x,y)∈Ω1,f(x,y),(x,y)∈Ω2,$$\begin{array}{} \displaystyle \frac{\partial^2u}{\partial x^2}+\frac{\partial^2u}{\partial y^2}=\begin{cases} g(x,y), & { \quad (x,y)\in\Omega_1 }, \\ \displaystyle f(x,y), & {\quad (x,y)\in\Omega_2 }, \end{cases} \end{array}$$(1)

where regions Ω1={[a,b+a2]×[a,b]}∪{[b+a2,b]×[a,b+a2]}andΩ2=[a+b2,b]×[a+b2,b]$\begin{array}{} {\Omega_1= \{[a,\frac{b+a}{2}]\times[a,b]\} \cup \{[\frac{b+a}{2},b]\times[a,\frac{b+a}{2}]\}} \,{\rm and}\, \Omega_2=[\frac{a+b}{2},b]\times[\frac{a+b}{2},b] \end{array}$ are respectively L-shape and square as shown in figure 1, such that Ω = Ω₁∪Ω₂, subject to the boundary conditions

Fig. 1

u(x,y)=g1(x,y),on∂Ω$$\begin{array}{} \displaystyle u(x,y)= g_1(x,y), \quad\qquad on\qquad \partial\Omega \end{array}$$(2)

where ∂Ω is boundary of Ω.

Many physical problems in engineering and physics can be modeled mathematically. These modeled problems in general governed by partial differential equations. We study the Poisson/Laplace equation a specific partial differential equation and these equations describe the behavior of electric, gravitational, and fluid potentials in the study of electrodynamics and electrostatics, gravitation and fluid dynamics respectively. It is difficult to get analytical solution of the most of the partial differential equations that arise in mathematical models of physical phenomena. So we have to use numerical methods to approximate the solution of such problems. Finite-difference, finite element and finite volume method are three important methods to numerically solve partial differential equations. A powerful and oldest method for solving Poisson**** or Laplace*** equation subject to conditions on boundary is the finite difference method, which makes use of finite-difference approximations. The Finite difference method is very simple and effective. The emphasis in this article is on the development of numerical method not to prove theoretical concepts of convergence and existence. Thus existence and uniqueness of the solution to problems (1) is assumed. We further assume that problems (1) is well posed. We will not consider the specific assumption to ensure existence and uniqueness of the solution to problems (1). Though it is important to prove theorems on uniqueness, existence and convergence and that can be developed in the future.

There are different methods depends on the geometry of the domain to solve Poisson or Laplace equation. The finite element and finite difference methods are mostly preferred method respectively for irregular and regular domain. There also exist several studies in the literature in which the standard finite difference method applied to solve elliptic PDEs in irregular domain too [1], [2], [3], [4]. We are considering Poisson and Laplace equation in a regular domain. For example, we know that the electric potential is related to the charge density by Poisson’s equation and in a charge-free region of space i.e as charge density becomes negligible, this becomes Laplace equation. Consider a problem in a domain where these two equations arise simultaneously and we wish to determine electric potential of this problem. This type of the problem may be considered as an obstacle boundary value problem and numerically solved by following the ideas in [5, 6].

In this article, we consider a finite difference method for numerical solution of system of elliptic equations (1) in a square domain. This is the usual second-order accurate scheme. The derivation of the method depends on difference approximation of the derivative on discrete mesh points in region of interest. So much research have been reported on numerical solution of elliptic partial differential equations, many of them are excellent work. But to best of our knowledge, neither similar problems (1) nor numerical method for the solution of problems (1) has been discussed in literature so for.

The present work is organized as follows. In section 2: we present novel finite difference approximation for system of elliptic equations. A novel finite difference method is presented so that the resulting difference equation need satisfies the boundary conditions exactly. A derivation of the present method will discuss in section 3. A local truncation error in the method and convergence will discuss in section 4 and 5. Finally, the application of the developed method will be present together with illustrative numerical results has been produced to show the efficiency of the method in section 6. A discussion and conclusion on the performance of the method will present in section 7.

Consider the square domain Ω = [a,b]×[a,b] for the solution of problems (1). Let h =(b-a)÷ N be the uniform mesh size in the x and y directions of the Cartesian coordinate system parallel to coordinate axes. Generate mesh poinIn this article, we consider a finite differencets (x_i, y_j), x_i = a + ih,i = 0,1,2,⋅⋅⋅⋅ N and y_j = a + jh, j = 0,1,2,⋅⋅⋅⋅ N. Let denote the mesh point (x_i, y_j) by (i,j). Let us denote the exact, approximate values of the solution of the problems (1) and source function f(x, y) at mesh point (i,j) by U_i,j, u_i,j and f_i,j respectively. Also we define mesh point (xi+12,yj+12)$\begin{array}{} (x_{i+\frac{1}{2}},y_{j+\frac{1}{2}}) \end{array}$ as (xi+h2,yj+h2)$\begin{array}{} (x_{i}+\frac{h}{2},y_{j}+\frac{h}{2}) \end{array}$ and denote (i+12,j+12).$\begin{array}{} (i+\frac{1}{2},j+\frac{1}{2}). \end{array}$ Similarly we can define other notations that we have used in this article. A novel finite difference method following the ideas in[7] for the the solution of problems (1) is defined as follows,

16u0,0+36u0,12+20u0,32+36u12,0+20u32,0−<!-- - -->144u12,12+16u32,32=9h2(5F12,12−<!-- - -->F12,0−<!-- - -->F0,12),i=j=1,$$\begin{array}{} \displaystyle 16u_{0,0}+36u_{0,\frac{1}{2}}+20u_{0,\frac{3}{2}}+36u_{\frac{1}{2},0}+20u_{\frac{3}{2},0}-144u_{\frac{1}{2},\frac{1}{2}}+16u_{\frac{3}{2},\frac{3}{2}}={9h^2}(5F_{\frac{1}{2},\frac{1}{2}}-F_{\frac{1}{2},0}-F_{0,\frac{1}{2}}),\quad { i=j=1 }, \end{array}$$

5ui−32,12+16ui−12,0+ui−32,32+6ui−12,32−34ui−12,12+5ui+12,12+ui+12,32=2h2(4Fi−12,12−Fi−12,0),2≤i<N,j=1,$$\begin{array}{} \displaystyle 5u_{i-\frac{3}{2},\frac{1}{2}}+16u_{i-\frac{1}{2},0}+u_{i-\frac{3}{2},\frac{3}{2}}+6u_{i-\frac{1}{2},\frac{3}{2}}-34u_{i-\frac{1}{2},\frac{1}{2}}+5u_{i+\frac{1}{2},\frac{1}{2}}+u_{i+\frac{1}{2},\frac{3}{2}}\\ \displaystyle \qquad\qquad\qquad\qquad\qquad\qquad\qquad\qquad\qquad\qquad\qquad\qquad\qquad\qquad={2h^2}(4F_{i-\frac{1}{2},\frac{1}{2}}-F_{i-\frac{1}{2},0}), \quad{ 2\leq i\lt N , j=1}, \end{array}$$

20uN−32,0+16uN−32,32+36uN−12,0−144uN−12,12+16uN,0+36uN,12+20uN,32=9h2(5FN−12,12−FN,12−FN−12,0),i=N,j=1,$$\begin{array}{} \displaystyle 20u_{N-\frac{3}{2},0}+16u_{N-\frac{3}{2},\frac{3}{2}}+36u_{N-\frac{1}{2},0}-144u_{N-\frac{1}{2},\frac{1}{2}}+16u_{N,0}+36u_{N,\frac{1}{2}}+20u_{N,\frac{3}{2}}\\ \displaystyle \qquad \qquad\qquad\qquad\qquad\qquad\qquad\qquad\qquad\qquad\qquad\qquad\quad ={9h^2}(5F_{N-\frac{1}{2},\frac{1}{2}}-F_{N,\frac{1}{2}}-F_{N-{1}{2},0}), \qquad { i= N, j=1 }, \end{array}$$

16u0,j−12+5u12,j−32−34u12,j−12+5u12,j+12+u32,j−32+6u32,j−12+u32,j+12=2h2(4F12,j−12F0,j−12),i=1,2≤j<N,$$\begin{array}{} \displaystyle 16u_{0,j-\frac{1}{2}}+5u_{\frac{1}{2},j-\frac{3}{2}}-34u_{\frac{1}{2},j-\frac{1}{2}}+5u_{\frac{1}{2},j+\frac{1}{2}}+u_{\frac{3}{2},j-\frac{3}{2}}+6u_{\frac{3}{2},j-\frac{1}{2}}+u_{\frac{3}{2},j+\frac{1}{2}}\\ \displaystyle \qquad\qquad\qquad\qquad\qquad\qquad\qquad\qquad\qquad\qquad\qquad\qquad\qquad\quad ={2h^2}(4F_{\frac{1}{2},j-\frac{1}{2}}F_{0,j-\frac{1}{2}}),\quad i=1,2\leq j \lt N, \end{array}$$

ui−32,j−12+ui+12,j−12+ui−12,j−32+ui−12,j+12−4ui−12,j−12=h2Fi−12,j−12,2≤i<N,2≤j<N,$$\begin{array}{} \displaystyle u_{i-\frac{3}{2},j-\frac{1}{2}}+u_{i+\frac{1}{2},j-\frac{1}{2}}+u_{i-\frac{1}{2},j-\frac{3}{2}}+u_{i-\frac{1}{2},j+\frac{1}{2}}-4u_{i -\frac{1}{2},j-\frac{1}{2}} ={h^2}F_{i-\frac{1}{2},j-\frac{1}{2}}, \quad {2\leq i\lt N , 2\leq j\lt N}, \end{array}$$

uN−32,j−32+6uN−32,j−12+uN−32,j+12+5uN−12,j−32−34uN−12,j−12+5uN−12,j+12+16uN,j−12=2h2(4FN−12,j−12−FN,12),i=N,2≤j<N,$$\begin{array}{} \displaystyle u_{N-\frac{3}{2},j-\frac{3}{2}}+6u_{N-\frac{3}{2},j-\frac{1}{2}}+u_{N-\frac{3}{2},j+\frac{1}{2}}+5u_{N-\frac{1}{2},j-\frac{3}{2}}-34u_{N-\frac{1}{2},j-\frac{1}{2}}+5u_{N-\frac{1}{2},j+\frac{1}{2}}+16u_{N,j-\frac{1}{2}}\\ \displaystyle \qquad \qquad\qquad\qquad\qquad\qquad\qquad\qquad\qquad\qquad\qquad\qquad\qquad ={2h^2}(4F_{N-\frac{1}{2},j-\frac{1}{2}}-F_{N,\frac{1}{2}}), \quad { i= N, 2\leq j \lt N }, \end{array}$$

20u0,N−32+36u0,N−12+16u0,N−144u12,N−12+36u12,N+16u32,N−32+20u32,N=9h2(5F12,N−12−F0,N−12−F12,N),i=1,j=N,$$\begin{array}{} \displaystyle 20u_{0,N-\frac{3}{2}}+36u_{0,N-\frac{1}{2}}+16u_{0,N}-144u_{\frac{1}{2},N-\frac{1}{2}}+36u_{\frac{1}{2},N}+16u_{\frac{3}{2},N-\frac{3}{2}}+20u_{\frac{3}{2},N}\\ \displaystyle \qquad\qquad\qquad\qquad\qquad\qquad\qquad\qquad\qquad\qquad\qquad\qquad\quad ={9h^2}(5F_{\frac{1}{2},N-\frac{1}{2}}-F_{0,N-\frac{1}{2}}-F_{\frac{1}{2},N}), \quad {i=1,j=N }, \end{array}$$

ui−32,N−32+5ui−32,N−12+16ui−12,N+6ui−12,N−32−34ui−12,N−12+ui+12,N−32+5ui+12,N−12=2h2(4F1−12,N−12−Fi−12,N),2≤i<N,j=N,$$\begin{array}{} \displaystyle u_{i-\frac{3}{2},N-\frac{3}{2}}+5u_{i-\frac{3}{2},N-\frac{1}{2}}+16u_{i-\frac{1}{2},N}+6u_{i-\frac{1}{2},N-\frac{3}{2}}-34u_{i-\frac{1}{2},N-\frac{1}{2}}+u_{i+\frac{1}{2},N-\frac{3}{2}}+5u_{i+\frac{1}{2},N-\frac{1}{2}}\\ \qquad \qquad \qquad \qquad \qquad \qquad \qquad \qquad \qquad \qquad \qquad \qquad \qquad \displaystyle =2h^2(4F_{1=\frac12,N-\frac12}-F_{i-1\frac12,N}),\quad 2\leq i \lt N,j=N, \end{array}$$

16uN−32,N−32+20uN−32,N+36uN−12,N−144uN−12,N−12+20uN,N−32+36uN,N−12+16uN,N=9h2(5FN−12,N−12−FN,N−12−FN−12,N),i=j=N,$$\begin{array}{} \displaystyle 16u_{N-\frac{3}{2},N-\frac{3}{2}}+20u_{N-\frac{3}{2},N}+36u_{N-\frac{1}{2},N}-144u_{N-\frac{1}{2},N-\frac{1}{2}}+20u_{N,N-\frac{3}{2}}+36u_{N,N-\frac{1}{2}}+16u_{N,N}\\ \displaystyle \qquad\qquad\qquad\qquad\qquad\qquad\qquad\qquad\qquad\qquad\qquad\quad={9h^2}(5F_{N-\frac{1}{2},N-\frac{1}{2}}-F_{N,N-\frac{1}{2}}-F_{N-\frac{1}{2},N}), \quad { i= j= N }, \end{array}$$(3)

where F_i,j = (∂2u∂x2+∂2u∂y2)i,j,$\begin{array}{} (\frac{\partial^2u}{\partial x^2}+\frac{\partial^2u}{\partial y^2})_{i,j}, \end{array}$ 1≤ i,j ≤ N. Finally we have,

ui,j=14(ui−12,j−12+ui+12,j−12+ui−12,j+12+ui+12,j+12),i,j=1,2,...,N−1.$$\begin{array}{} \displaystyle u_{i,j}=\frac{1}{4}(u_{i-\frac{1}{2},j-\frac{1}{2}}+u_{i+\frac{1}{2},j-\frac{1}{2}}+u_{i-\frac{1}{2},j+\frac{1}{2}}+u_{i+\frac{1}{2},j+\frac{1}{2}}),\quad { i, j=1,2,...,N-1}. \end{array}$$(4)

In this section we discuss the development of the method (3). So, let consider the development of the first equation i.e. i = j = 1 in (3). Consider

a1ui−1,j−1+a2ui−1,j−12+a3ui−1,j+12+a4ui−12,j−1+a5ui+12,j−1+a6ui−12,j−12+a7ui+12,j+12=h2(b1(uxx+uyy)i−12,j−12+b2(uxx+uyy)i−12,j−1+b3(uxx+uyy)i−1,j−12)$$\begin{array}{} \displaystyle a_1u_{i-1,j-1}+a_2u_{i-1,j-\frac{1}{2}}+a_3u_{i-1,j+\frac{1}{2}}+a_4u_{i-\frac{1}{2},j-1}+a_5u_{i+\frac{1}{2},j-1}+a_6u_{i-\frac{1}{2},j-\frac{1}{2}}+a_7u_{i+\frac{1}{2},j+\frac{1}{2}}\\ \displaystyle \qquad\qquad\qquad\qquad\qquad\qquad\qquad=h^2(b_1(u_{xx}+u_{yy})_{i-\frac{1}{2},j-\frac{1}{2}}+b_2(u_{xx}+u_{yy})_{i-\frac{1}{2},j-1}+b_3(u_{xx}+u_{yy})_{i-1,j-\frac{1}{2}}) \end{array}$$(5)

where ai′sandbi′s$\begin{array}{} a_i's \,{\rm and}\, b_i's \end{array}$ are real constant to be determined. To determine these constants in (5), let write the terms in (5) in Taylor series about a point i−12,j−12.$\begin{array}{} {i-\frac{1}{2},j-\frac{1}{2}}. \end{array}$ Comparing the coefficients of h^p,p = 0,1,2,3 in so obtained Taylor series expansion for (5) and application of method of undetermined coefficients, we have a system of linear equations in unknown ai′sandbi′s$\begin{array}{} a_i's \,{\rm and}\, b_i's \end{array}$. Thus solving the system of equations, we will obtain following,

(a1,a2,a3,a4,a5,a6,a7,b1,b2,b3)=(16,36,20,36,20,−144,16,45,−9,−9)$$\begin{array}{} \displaystyle (a_1, a_2, a_3, a_4, a_5, a_6, a_7, b_1, b_2, b_3)=(16, 36, 20, 36, 20, -144, 16, 45, -9, -9) \end{array}$$(6)

Let substitute the value of constants from (6) in (5). Thus we have,

16ui−1,j−1+36ui−1,j−12+20ui−1,j+12+36ui−12,j−1+20ui+12,j−1−144ui−12,j−12+16ui+12,j+12=9h2(5Fi−12,j−12−Fi−12,j−1−Fi−1,j−12)+Ti−12,j−12,$$\begin{array}{} \displaystyle 16u_{i-1,j-1}+36u_{i-1,j-\frac{1}{2}}+20u_{i-1,j+\frac{1}{2}}+36u_{i-\frac{1}{2},j-1}+20u_{i+\frac{1}{2},j-1}-144u_{i-\frac{1}{2},j-\frac{1}{2}}+16u_{i+\frac{1}{2},j+\frac{1}{2}}\\ \displaystyle \qquad\qquad\qquad\qquad\qquad\qquad\qquad\qquad\qquad\qquad\qquad \quad={9h^2}(5F_{i-\frac{1}{2},j-\frac{1}{2}}-F_{i-\frac{1}{2},j-1}-F_{i-1,j-\frac{1}{2}})+T_{i-\frac{1}{2},j-\frac{1}{2}}, \end{array}$$(7)

where Fi−12,j−12=(uxx+uyy)i−12,j−12andTi−12,j−12$\begin{array}{} F_{i-\frac{1}{2},j-\frac{1}{2}}= (u_{xx}+u_{yy})_{i-\frac{1}{2},j-\frac{1}{2}}{\rm and}\, T_{i-\frac{1}{2},j-\frac{1}{2}} \end{array}$ terms truncated in Taylor series expansion of (5) or truncating error term. Neglecting the truncation error term Ti−12,j−12$\begin{array}{} T_{i-\frac{1}{2},j-\frac{1}{2}} \end{array}$ in (7), we will get the first equation i.e. i = j = 1 in (3). Similarly we can derive other equations for different values of i and j in (3).

If we write system of equations given by (3) at each mesh point, we will obtain N × N linear system of equations if source functions f(x, y) and g(x, y) are linear otherwise nonlinear system of equations in ui−12,j−12,$\begin{array}{} u_{i-\frac{1}{2},j-\frac{1}{2}}, \end{array}$i,j = 1,2,...,N. To obtain a satisfactory solution of this system of equations with reasonable degree of accuracy, we apply an iterative method to solve the system of equation (3).

Let Ti−12,j−12$\begin{array}{} T_{i-\frac{1}{2},j-\frac{1}{2}} \end{array}$ be truncation error at mesh point (xi−12,yj−12),$\begin{array}{} (x_{i-\frac{1}{2}},y_{j-\frac{1}{2}}), \end{array}$ 1 ≤ i,j ≤ N as defined in [8]. Consider Ti−12,j−12$\begin{array}{} T_{i-\frac{1}{2},j-\frac{1}{2}} \end{array}$, and i = j = 1 truncating error term in first equation of (3). By Taylor series expansion of u′s about mesh point (xi−12,yj−12)$\begin{array}{} (x_{i-\frac{1}{2}},y_{j-\frac{1}{2}}) \end{array}$ and from (7), we have

Ti−12,j−12=16ui−1,j−1+36ui−1,j−12+20ui−1,j+12+36ui−12,j−1−144ui−12,j−12+20ui+12,j−1+16ui+12,j+12−9h2(5Fi−12,j−12−Fi−12,j−1−Fi−1,j−12)=3h4320(15∂4u∂x4+4∂4u∂x3∂y+48∂4u∂x2∂y2+4∂4u∂x∂y3+15∂4u∂y4)i,j$$\begin{array}{} \displaystyle T_{i-\frac{1}{2},j-\frac{1}{2}}=16u_{i-1,j-1}+36u_{i-1,j-\frac{1}{2}}+20u_{i-1,j+\frac{1}{2}}+36u_{i-\frac{1}{2},j-1}-144u_{i-\frac{1}{2},j-\frac{1}{2}}+20u_{i+\frac{1}{2},j-1}+16u_{i+\frac{1}{2},j+\frac{1}{2}}\\ \displaystyle \qquad \qquad\qquad\qquad\qquad\qquad\qquad\quad-{9h^2}(5F_{i-\frac{1}{2},j-\frac{1}{2}}-F_{i-\frac{1}{2},j-1}-F_{i-1,j-\frac{1}{2}})\\ \displaystyle \qquad\qquad\qquad =\frac{3h^4}{320}(15\frac{\partial^4u}{\partial x^4}+4\frac{\partial^4u}{\partial x^3\partial y}+48\frac{\partial^4u}{\partial x^2\partial y^2}+4\frac{\partial^4u}{\partial x\partial y^3}+15\frac{\partial^4u}{\partial y^4})_{i,j} \end{array}$$

Thus we have obtained an expression of local truncation error term in case i = j = 1 i.e. first equation of (3). Similarly we can compute local error terms for different values and cases of i and j i.e. in other equations in (3) as follows,

Ti−12,j−12=3h4320(15∂4u∂x4+4∂4u∂x3∂y+48∂4u∂x2∂y2+4∂4u∂x∂y3+15∂4u∂y4)i−12,j−12,i=j=1,h48(4∂4u∂x4+6∂4u∂x2∂y2+5∂4u∂y4)i−12,j−12,2≤i<N,j=1,3h4320(15∂4u∂x4−4∂4u∂x3∂y+48∂4u∂x2∂y2−4∂4u∂x∂y3+15∂4u∂y4)i−12,j−12,i=N,j=1,h48(5∂4u∂x4+6∂4u∂x2∂y2+4∂4u∂y4)i−12,j−12,i=1,2≤j<N,h412(∂4u∂x4+∂4u∂y4)i−12,j−12,2≤i<N,2≤j<N,h48(5∂4u∂x4+6∂4u∂x2∂y2+4∂4u∂y4)i−12,j−12,i=N,2≤j<N,3h4320(15∂4u∂x4−4∂4u∂x3∂y+48∂4u∂x2∂y2−4∂4u∂x∂y3+15∂4u∂y4)i−12,j−12,i=1,j=N,h48(4∂4u∂x4+6∂4u∂x2∂y2+5∂4u∂y4)i−12,j−12,2≤i<N,j=N,3h4320(15∂4u∂x4+4∂4u∂x3∂y+48∂4u∂x2∂y2+4∂4u∂x∂y3+15∂4u∂y4)i,j,i=j=N.$$\begin{array}{} \displaystyle T_{i-\frac{1}{2},j-\frac{1}{2}}=\begin{cases} \frac{3h^4}{320}(15\frac{\partial^4u}{\partial x^4}+4\frac{\partial^4u}{\partial x^3\partial y}+48\frac{\partial^4u}{\partial x^2\partial y^2}+4\frac{\partial^4u}{\partial x\partial y^3}+15\frac{\partial^4u}{\partial y^4})_{i-\frac{1}{2},j-\frac{1}{2}}\,, & { i=j=1 }, \\ \displaystyle \frac{h^4}{8}(4\frac{\partial^4u}{\partial x^4}+6\frac{\partial^4u}{\partial x^2\partial y^2}+5\frac{\partial^4u}{\partial y^4})_{i-\frac{1}{2},j-\frac{1}{2}}\,, & { 2\leq i\lt N , j=1},\\ \frac{3h^4}{320}(15\frac{\partial^4u}{\partial x^4}-4\frac{\partial^4u}{\partial x^3\partial y}+48\frac{\partial^4u}{\partial x^2\partial y^2}-4\frac{\partial^4u}{\partial x\partial y^3}+15\frac{\partial^4u}{\partial y^4})_{i-\frac{1}{2},j-\frac{1}{2}}\,, & { i= N, j=1 },\\ \displaystyle \frac{h^4}{8}(5\frac{\partial^4u}{\partial x^4}+6\frac{\partial^4u}{\partial x^2\partial y^2}+4\frac{\partial^4u}{\partial y^4})_{i-\frac{1}{2},j-\frac{1}{2}}\,, & { i=1,2\leq j\lt N },\\ \displaystyle \frac{h^4}{12}(\frac{\partial^4u}{\partial x^4}+\frac{\partial^4u}{\partial y^4})_{i-\frac{1}{2},j-\frac{1}{2}}\,, & { 2\leq i\lt N , 2\leq j\lt N},\\ \displaystyle \frac{h^4}{8}(5\frac{\partial^4u}{\partial x^4}+6\frac{\partial^4u}{\partial x^2\partial y^2}+4\frac{\partial^4u}{\partial y^4})_{i-\frac{1}{2},j-\frac{1}{2}}\,, & { i= N, 2\leq j \lt N },\\ \displaystyle \frac{3h^4}{320}(15\frac{\partial^4u}{\partial x^4}-4\frac{\partial^4u}{\partial x^3\partial y}+48\frac{\partial^4u}{\partial x^2\partial y^2}-4\frac{\partial^4u}{\partial x\partial y^3}+15\frac{\partial^4u}{\partial y^4})_{i-\frac{1}{2},j-\frac{1}{2}}\,, & { i=1,j=N },\\ \displaystyle \frac{h^4}{8}(4\frac{\partial^4u}{\partial x^4}+6\frac{\partial^4u}{\partial x^2\partial y^2}+5\frac{\partial^4u}{\partial y^4})_{i-\frac{1}{2},j-\frac{1}{2}}\,, & { 2\leq i\lt N , j= N},\\ \displaystyle \frac{3h^4}{320}(15\frac{\partial^4u}{\partial x^4}+4\frac{\partial^4u}{\partial x^3\partial y}+48\frac{\partial^4u}{\partial x^2\partial y^2}+4\frac{\partial^4u}{\partial x\partial y^3}+15\frac{\partial^4u}{\partial y^4})_{i,j}\,,& { i= j= N }. \end{cases} \end{array}$$(8)

So we see that the order of local error terms Ti−12,j−12$\begin{array}{} T_{i-\frac{1}{2},j-\frac{1}{2}} \end{array}$ in (8) are of O(h⁴). Thus the order of propose method (3) is O(h²).

We next discuss the convergence of the proposed method and under suitable conditions [9]. We prove that the order of the convergence is O(h²). For each i,j = 1(1)N, let us define

ϕi−12,j−12=h245Fi−12,j−12+ values on boundary,i=j=1,8Fi−i−12,j−12+ values on boundary,2≤i<N,j=1,45FN−12,j−12+ values on boundary ,i=N,j=1,8Fi−12,j−12+ values on boundary ,i=1,2≤j<N,Fi−12,j−12,2≤i<N,2≤j<N,8Fi−12,j−12+ values on boundary ,i=N,2≤j<N,45Fi−12,j−12+ values on boundary ,i=1,j=N,8F1−12,N−12+ values on boundary ,2≤i<N,j=N,45FN−12,N−12+ values on boundary ,i=j=N.$$\begin{array}{} \displaystyle \phi_{i-\frac{1}{2},j-\frac{1}{2}}={h^2}\begin{cases} 45F_{\frac{i-1}{2},\frac{j-1}{2}}+ \text{ values on boundary}, \quad { i=j=1 }, \\ \displaystyle 8F_{i-\frac{i-1}{2},\frac{j-1}{2}}+ \text{ values on boundary}, \quad { 2\leq i\lt N , j=1},\\ \displaystyle 45F_{N-\frac{1}{2},\frac{j-1}{2}}+ \text{ values on boundary }, \quad { i= N, j=1 },\\ \displaystyle 8F_{\frac{i-1}{2},j-\frac{1}{2}}+ \text{ values on boundary }, \quad { i=1,2\leq j\lt N },\\ \displaystyle F_{i-\frac{1}{2},j-\frac{1}{2}}, \quad\qquad\qquad\qquad\qquad\qquad { 2\leq i\lt N , 2\leq j\lt N},\\ \displaystyle 8F_{i-\frac{1}{2},j-\frac{1}{2}}+ \text{ values on boundary }, \quad { i= N, 2\leq j \lt N },\\ \displaystyle 45F_{\frac{i-1}{2},j-\frac{1}{2}}+ \text{ values on boundary }, \quad { i=1,j=N },\\ \displaystyle 8F_{1-\frac{1}{2},N-\frac{1}{2}}+ \text{ values on boundary }, \quad { 2\leq i \lt N, j= N },\\ \displaystyle 45F_{N-\frac{1}{2},N-\frac{1}{2}}+ \text{ values on boundary }, \quad { i= j= N }. \end{cases} \end{array} $$(9)

and

Ei−12,j−12=ui−12,j−12−Ui−12,j−12,1≤i,j≤N.$$\begin{array}{} \displaystyle E_{i-\frac{1}{2},j-\frac{1}{2}}=u_{i-\frac{1}{2},j-\frac{1}{2}}-U_{i-\frac{1}{2},j-\frac{1}{2}} ,\quad 1\leq i,j\leq N. \end{array}$$

Let define matrix S_{(N–1)² × 1},

S=[S12,12,S32,12,..........,SN−12,12,S12,32,S32,32,.........,SN−12,32,.......,S12,N−12,S32,N−12,.......,SN−12,N−12]T$$\begin{array}{} \displaystyle \mathbf{S}=[S_{\frac{1}{2},\frac{1}{2}},S_{\frac{3}{2},\frac{1}{2}},..........,S_{\frac{N-1}{2},\frac{1}{2}},S_{\frac{1}{2},\frac{3}{2}},S_{\frac{3}{2},\frac{3}{2}},........., S_{\frac{N-1}{2},\frac{3}{2}},.......,S_{\frac{1}{2},\frac{N-1}{2}},S_{\frac{3}{2},\frac{N-1}{2}},.......,S_{\frac{N-1}{2},\frac{N-1}{2}}]^T \end{array} $$

The the difference method (3) represents a system of linear equations in the unknown ui−12,j−12,1≤i,j≤N.$\begin{array}{} u_{i-\frac{1}{2},j-\frac{1}{2}},1\leq i,j\leq N. \end{array}$ Using S = ϕ,u,U and T, we write (3) in matrix form as,

−Mu+ϕ=0$$\begin{array}{} \displaystyle -\mathbf{Mu} + \boldsymbol{\phi} = \boldsymbol{0} \end{array}$$(10)

where M_{(N–1)² ×(N–1)²} is block tridiagonal matrix and defined as:

M=AB0CDCCDC......0BA$$\begin{array}{} \displaystyle \mathbf{M}= \left( \begin{array}{ccccc} A & B& & & 0 \\ C & D & C & & \\ &C & D & C &\\ &..&..&..&\\ 0 & & & B & A \end{array} \right) \nonumber \end{array}$$

in which blocks are

A=14400−534−5−534−5......0144(N−1)×(N−1),B=0−160−1−6−1−1−6−1......0−160(N−1)×(N−1)C=−5−100−100−10......0−1−5(N−1)×(N−1),D=34−60−14−1−14−1......0−634(N−1)×(N−1)$$\begin{array}{} \displaystyle {A}= \left( \begin{array}{ccccc} 144 & 0 & & & 0 \\ -5 & 34 & -5 & & \\ &-5 & 34 & -5 &\\ &..&..&..& \\ 0 & & & & 144 \end{array} \right)_{(N-1)\times(N-1)}~~ ,\quad B =\left( \begin{array}{ccccc} 0&-16& & &0 \\ -1&-6&-1& &\\ &-1 & -6 & -1 & \\ &..&..&..& \\ 0& & &-16&0 \end{array} \right)_{(N-1)\times(N-1)}\\ {C}= \left( \begin{array}{ccccc} -5 & -1 & & & 0 \\ 0 & -1 & 0 & & \\ & 0 & -1 & 0 &\\ &..&..&..& \\ 0 & & & -1 & -5 \end{array} \right)_{(N-1)\times(N-1)} \quad,\quad D =\left( \begin{array}{ccccc} 34 &-6& & &0 \\ -1&4&-1& &\\ &-1 & 4 & -1 & \\ &..&..&..& \\ 0& & &-6& 34\end{array} \right)_{(N-1)\times(N-1)} \end{array} $$

But U is exact solution of (3), which in matrix form can be written as

−MU+ϕ+T=0$$\begin{array}{} \displaystyle -\mathbf{MU} + \boldsymbol{\phi}+\mathbf{T} = \boldsymbol{0} \end{array}$$(11)

where T is truncation error matrix and each element is of O(h⁴). Let there are no roundoff errors in solution of difference equation (3), so from (10) and (11) we get the error equation as

ME=T$$\begin{array}{} \displaystyle \mathbf{ME}=\mathbf{T} \end{array}$$(12)

The matrix M is row diagonally dominant. Let M be the adjacency matrix of some graph G. We may easily prove that graph G is connected. From this fact it follows that adjacency matrix M is irreducible [10]. By the row sum criterion it follows that M is monotone [11]. Thus positive M^–1 exist. Moreover M is regular, thus system of equations (3) can be solved by Gaussian elimination method by preserving its band structure. Thus from (12), we have

∥E∥∞≤∥M−1∥∞∥T∥∞.$$\begin{array}{} \displaystyle \|\mathbf{E}\|_\infty \leq \|\mathbf{M}^{-1}\|_\infty \|\mathbf{T}\|_\infty . \end{array}$$(13)

With the help of (8) and (13), for sufficiently small h, we have

∥E∥≤O(h2)$$\begin{array}{} \displaystyle \|\mathbf{E}\|\leq O(h^2) \end{array}$$(14)

Thus the proposed difference method (3) converges and the order of the convergence is quadratic.

The validity and effectiveness of the proposed finite difference method on uniform mesh have been tested on linear model problem. In our experiment we took a square as a region of integration and covered it with a uniform mesh of size h. We have used iterative method Gauss Seidel to solve the system of linear equations. In table, we have shown maximum absolute error MAU, computed for different values of N, using following formula,

MAU=max2≤i,j≤N−1u(xi,yj)−ui,j$$\begin{array}{} \displaystyle MAU= \max_{\substack {2\leq{i,j}\leq N-1}}{\left\vert u(x_i,y_j)-u_{i,j}\right\vert } \end{array} $$

All the computations in the experiment were performed on MS Window 2007 professional operating system in the GNU FORTRAN environment version -99 compiler(2.95 of gcc) running on Intel Duo core 2.20 Ghz PC . The stopping condition for iteration was either error of order 10^–10 or number of iterations 10³⋅

Table 1

Maximum absolute errors in u(x, y) in Ω₁ and Ω₂ for problem 1.

Table 2

a0=0.516,a1=1.0−a02.0,$\begin{array}{} a_0=\frac{0.5}{16}, a_1=1.0-\frac{a_0}{2.0}, \end{array}$a₂ = a₀ exp(–1:5)

Table 3

a0=0.5128,a1=1.0−a02.0,$\begin{array}{} a_0=\frac{0.5}{128}, a_1=1.0-\frac{a_0}{2.0}, \end{array}$a₂ = a₀ exp(–1:5)

In this article we have discussed finite difference method for numerical solution of system of elliptic equations. The method is theoretically second order accurate and Dirichlet boundary conditions incorporated into the method. In considered model problem we find that the computational performance of the method and accuracy depends on the constructed exact solution. The results obtained in model problem is not in good agreement to the estimated order of the method. In future we wish to improve proposed method. Though method developed on rectangular domain and used equally spaced grid mesh systems, it has good potential for efficient application to many problems on different geometries and dimensions, work in this specific direction is in progress.

We are interested in computing approximate numerical solution of the elliptic obstacle problem at the mesh points (x_i, y_j) in two dimension space. The development of the method in this article may be applied for problems in higher dimensions and different geometry. Fortunately we able to find some idea about the convergence of finite difference method under an appropriate condition but unable to measure maximum absolute error in terms of the mesh size h in numerical experiment. In numerical experiment we have considered exact solution and its first order partial derivatives are continuous at the boundary points (a+b2,b)and(b,a+b2)$\begin{array}{} (\frac{a+b}{2},b) \,{\rm and}\, (b,\frac{a+b}{2}) \end{array}$ of regions L and square. However in theoretical estimation of error we have considered smooth u(x, y) i.e. fourth order partial derivatives are continuous in considered domain. Thus under an appropriate assumption, we obtain an error in U(x_i, y_j) exact solution and u an approximate solution obtained by proposed finite difference method (3) of the considered problem by using formula |U(x_i, y_j)–u_i,j| < Ch^p, where C is constant i.e. independent of U(x, y), h the mesh size of the discretization and p the order of the method estimated by numerical experiment. In tables we have find that maximum absolute errors in the exact and approximate solutions for the considered model problem for different values of h=1N.$\begin{array}{} h=\frac{1}{N}. \end{array}$ We observe from the tables that as N increases maximum absolute error decreases.