where regions
$\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 Ω = Ω1∪Ω2, subject to the boundary conditions
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.
The Finite Difference Method
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 (xi, yj), xi = a + ih,i = 0,1,2,⋅⋅⋅⋅ N and yj = a + jh, j = 0,1,2,⋅⋅⋅⋅ N. Let denote the mesh point (xi, yj) 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 Ui,j, ui,j and fi,j respectively. Also we define mesh point
$\begin{array}{}
(x_{i+\frac{1}{2}},y_{j+\frac{1}{2}})
\end{array}$
as
$\begin{array}{}
(x_{i}+\frac{h}{2},y_{j}+\frac{h}{2})
\end{array}$
and denote
$\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,
where
$\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
$\begin{array}{}
{i-\frac{1}{2},j-\frac{1}{2}}.
\end{array}$
Comparing the coefficients of hp,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
$\begin{array}{}
a_i's \,{\rm and}\, b_i's
\end{array}$.
Thus solving the system of equations, we will obtain following,
where
$\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
$\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
$\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).
Local Truncation Error
Let
$\begin{array}{}
T_{i-\frac{1}{2},j-\frac{1}{2}}
\end{array}$
be truncation error at mesh point
$\begin{array}{}
(x_{i-\frac{1}{2}},y_{j-\frac{1}{2}}),
\end{array}$
1 ≤ i,j ≤ N as defined in [8]. Consider
$\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
$\begin{array}{}
(x_{i-\frac{1}{2}},y_{j-\frac{1}{2}})
\end{array}$
and from (7), we have
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,
So we see that the order of local error terms
$\begin{array}{}
T_{i-\frac{1}{2},j-\frac{1}{2}}
\end{array}$
in (8) are of O(h4). Thus the order of propose method (3) is O(h2).
Convergence of the method
We next discuss the convergence of the proposed method and under suitable conditions [9]. We prove that the order of the convergence is O(h2). For each i,j = 1(1)N, let us define
$$\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}
$$
and
$$\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}$$
The the difference method (3) represents a system of linear equations in the unknown
$\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,
where M(N–1)2 ×(N–1)2 is block tridiagonal matrix and defined as:
$$\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}$$
where T is truncation error matrix and each element is of O(h4). Let there are no roundoff errors in solution of difference equation (3), so from (10) and (11) we get the error equation as
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
Thus the proposed difference method (3) converges and the order of the convergence is quadratic.
Numerical Experiments
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,
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 103⋅
Problem 1
Consider a linear elliptic boundary value problems in a square region Ω = [0,1]×[0,1] which, when solving consists of
where regions
$\begin{array}{}
{\Omega_1= \{[0,\frac{1}{2}]\times[0,1]\} \cup \{[\frac{1}{2},1]\times[0,\frac{1}{2}]\}}\,{\rm and}\,
\Omega_2=[\frac{1}{2},1]\times[\frac{1}{2},1]
\end{array}$
are respectively L-shape and square such that {Ω = Ω1∪Ω2, subject to the boundary conditions u(x, y)on ∂Ω. The source function f(x, y) is given such that the exact solution is
where a0,a1 and a2 are constants and determined [12] by the considering the continuity of u, ux and uy points on ∂Ω = ∂Ω1 ∩ ∂Ω2. In table 1, we have presented the time Etime in seconds and the computed MAU for a0 = 0.5, a1 = 1.0–
$\begin{array}{}
\frac{a_0}{2.0},
\end{array}
$a2 = a0 exp(–1.5) and different values of N in constructed solutions. We have presented MAU for different values of a0,a1 and a2table 2 and table 3.
Maximum absolute errors in u(x, y) in Ω1 and Ω2 for problem 1.
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.
Remark
We are interested in computing approximate numerical solution of the elliptic obstacle problem at the mesh points (xi, yj) 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
$\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(xi, yj) exact solution and u an approximate solution obtained by proposed finite difference method (3) of the considered problem by using formula |U(xi, yj)–ui,j| < Chp, 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
$\begin{array}{}
h=\frac{1}{N}.
\end{array}$
We observe from the tables that as N increases maximum absolute error decreases.