4_erdem.dvi


EUROPEAN JOURNAL OF PURE AND APPLIED MATHEMATICS

Vol. 6, No. 1, 2013, 30-43

ISSN 1307-5543 – www.ejpam.com

Comparative Analysis of the Modified SOR and BGC Methods

Applied to the Poleness Conservative Finite Difference Scheme

Arzu Erdem

Department of Mathematics, Kocaeli University, Umuttepe Campus, Izmit-Kocaeli, 41380, Turkey

Abstract. The poleness conservative finite difference scheme based on the weak solution of Poisson

equation in polar coordinates is studied. Due to the singularity at r = 0 in the considered polar domain

Ωrϕ , a special technique of deriving the finite difference scheme in the neighbourhood of the pole point

r = 0 is described. The constructed scheme has the order of approximation O
�
(h2

r
+ h2

ϕ)/r
�

. In the

second part of the paper the structure of the corresponding non-symmetric sparse block matrix is

analyzed. A special algorithm based on SOR-method is presented for the numerical solution of the

corresponding system of linear algebraic equations. The theoretical result are illustrated by numerical

examples for continuous as well as discontinuous source function.

2010 Mathematics Subject Classifications: 65M06, 65F50, 65F10

Key Words and Phrases: Finite difference method, elliptic problem, polar coordinates, sparse matrix

1. Introduction

In this paper we consider the following Dirichlet problem for the Poisson equation in polar

coordinates (r,ϕ):





Au := −1

r

∂

∂ r

�
r ∂ u

∂ r

�
− 1

r2
∂ 2u

∂ ϕ2 = F(r,ϕ), (r,ϕ) ∈ ΩR,

u(R,ϕ) = 0, (r,ϕ) ∈ Γϕ,

u(r, 0) = u(r,β), r ∈ (0,R),

. (1)

where ΩR := {(r,ϕ) ∈ R2 : r ∈ [0,R), ϕ ∈ [0,β)}, Γϕ := {(R,ϕ) : ϕ ∈ [0,β)}, and

β ∈ (0,2π].

This problem is a mathematical model of various physical and engineering problems aris-

ing in steady state flow of an incompressible viscous fluid in a duct of circular cross-section

[15, 13], in the determination of a potential in electrostatics [3] and in the elasticity theory

[8]. The two circumstances may lead to singularities: the geometrical singularity related to

Email address: erdem.arzu�gmail.
om

http://www.ejpam.com 30 c© 2013 EJPAM All rights reserved.



A. Erdem / Eur. J. Pure Appl. Math, 6 (2013), 30-43 31

the corner β ∈ (0,2π) of the polar domain and the pole point r = 0. The first type of singular-

ity means that for some values of the parameter β ∈ (0,π) the second (and higher) derivative

of the solution u(r,ϕ) with respect to r ∈ [0,R) is singular at r = 0. As show examples below

the regularity of the weak solution u ∈ H2(Ω) ∩ �H1(Ω) of problem (1) depends on the value

of the angle β ∈ (0,2π). This behaviour is usual for second order elliptic problems and is

related to the C2-regularity property of the boundary ∂ΩRβ of the considered domain [1, 11].

The singularity at the re-entrant corner makes the numerical solution of these problems chal-

lenging. The second type of singularity is a reason of many difficulties in constructing the

standard finite difference (FD) schemes, when the pole r = 0 is treated as a computational

boundary. To avoid these difficulties various FD and pseudo-spectral (PS) methods have been

suggested in literature [see 4, 14]. These schemes include the necessity of special boundary

closures, which leads to undesirable clustering of grid points in PS schemes [see 4, 6]. The

treatment of the singularities related to the situations r → 0 and sinϕ → 0 have been given

in [9, 10].

This paper is devoted to fill in the lack of result for conservative finite difference scheme

for problem (1) and nonsymmetric sparse block matrices related to finite difference equations

in polar coordinates. We present a conservative finite difference scheme for this problem

and prove its convergence. Our approach is based on the Lax-Wondroff theorem [7], which

guarantees convergence of a conservative FD schemes in the class of weak solutions, as the

polar mesh is refined. Note that a similar technique was used in [12] for problem (1), where

the classical solution of the boundary value problem (1) is considered. Since we are interested

in bounded weak solutions u ∈ H1(ΩR), we require that the solution of problem (1) satisfies

the boundedness at r = 0 condition

lim
r→0

r
∂ u

∂ r
= 0. (2)

Hence one needs to approximate not only the elliptic equation (1), but also condition (2).

This condition will be used for obtaining the conservative finite difference scheme in the

neighbourhood of the pole point r = 0.

Further, the FD approximations of problem (1)-(2) lead to large sparse system of linear

equations, with the special nonsymmetric matrix, due to the periodicity condition

u(r, 0) = u(r,β). These type of linear systems require time-consuming algorithms for their

effective numerical solution [2, 5]. We use the special structure of the obtained matrix [5]

and construct a fast iteration algorithm, based on SOR-method.

The paper is organized as follows. In section 2 the weak solution of problem (1)-(2) is

defined. The piecewise uniform polar mesh and the conservative FD scheme for the problem

is constructed in Section 3. In Section 4 the structure of the corresponding sparse matrix and

iteration algorithm is discussed. Numerical solution of problem (1)-(2) for different types of

source function F(r,ϕ) and results are presented in Section 5.



A. Erdem / Eur. J. Pure Appl. Math, 6 (2013), 30-43 32

2. The Weak Solution and its Regularity Depending on the Parameter β ∈ (0, 2π)

Let v ∈ �H1(ΩRβ) be an arbitrary function, where

�H1(ΩRβ) := {v ∈ H1(ΩRβ) : u(R,ϕ) = 0,ϕ ∈ (0,β); u(r, 0) = u(r,β), r ∈ (0,R)}
and H1(ΩRβ ) is the Sobolev space of functions v = v(r,ϕ) with the norm

‖u‖1 :=





∫ ∫

ΩRβ

�
u2 + |∇u|2
�

rd rdϕ





1/2

.

Let us multiply the both sides of equation (1) to r v(r,ϕ), r ∈ (0,β) and integrate on ΩRβ :

−
∫ ∫

ΩRβ

�
∂

∂ r

�
r
∂ u

∂ r

�
+

1

r

∂ 2u

∂ ϕ2

�
vdrdϕ =

∫ ∫

ΩRβ

F(r,ϕ)v rdrdϕ.

Applying here by part integration, using the boundary and periodicity conditions (2) we get

∫ ∫

ΩRβ

∇u(r,ϕ) · ∇v(r,ϕ)rd rdϕ =

∫ ∫

ΩRβ

F(r,ϕ)v rdrdϕ, ∀v ∈ �H1(ΩRβ). (3)

Here ∇u is the gradient vector in polar coordinates,

∇u=
∂ u

∂ r
e1+

1

r

∂ u

∂ ϕ
e2,

�
e1

e2

�
=

�
Cosϕ, Sinϕ

Sinϕ, Cosϕ

��
i1
i2

�
,

and i1, i2 are unit coordinate vector in Cartesian coordinates. The solution u ∈ �H1(ΩRβ) of the

integral identity (3) is defined as a weak solution of the boundary value problem (1)-(2).

The integral identity (3) shows that if the function u ∈ C2(ΩRβ)∩ C1(ΩRβ ) is the solution

of the boundary value problem (1)-(2), then for all v ∈ �H1(ΩRβ ) this identity holds. Hence the

weak solution of problem (1)-(2) can be defined as function u ∈ �H1(ΩRβ) satisfying this inte-

gral identity for all v ∈ �H1(ΩRβ). According to the general theory for linear elliptic boundary

value problems the regular weak solution of problem (1)-(2) belongs to H2(ΩRβ) ∩ �H1(ΩRβ),

if the boundary ∂ΩRβ is of class C2 [4, 11]. The following example shows that depending on

the values β ∈ (π, 2π) of the parameter β the weak solution of the boundary value problem

(1)-(2) may not belong to the class H2(ΩRβ ).

Example 1. The function

u(r,ϕ) = rπ/βSin(πϕ/β), (r,ϕ) ∈ ΩRβ (4)

satisfies the Laplace equation in polar coordinates and the Dirichlet condition

u(R,ϕ) = Rπ/βSin(πϕ/β), ϕ ∈ (0,β).



A. Erdem / Eur. J. Pure Appl. Math, 6 (2013), 30-43 33

This function belongs to the Sobolev space H2(ΩRβ) ∩ H1(ΩRβ ), for all β ∈ (0,π). However for

the values β ∈ (π, 2π) this solution doesn’t belong to H2(ΩRβ ), since rα /∈ H2(ΩRβ) for α < 1.

The reason, as shown in [11], is that for β ∈ (π, 2π) the boundary ∂ΩRβ doesn’t belong to the

class C2, as requires the Agmon-Nirenberg regularity theorem [1].

Note that the above loss of regularity of the boundary ∂ΩRβ is a result of the introduced

angle α = πϕ/β . When the parameter β changes in (0,2π) the angle α always remains in

(0,π). This also implies that for the function u(r,ϕ), given by (4), the periodicity condition

∂ u(r, 0)

∂ ϕ
=
∂ u(r,β)

∂ ϕ
, r ∈ (0,R) (5)

doesn’t hold. The next example shows that by introducing the parameter α = nπϕ/β , n = 2,

the fulfilment of this condition can be achieved.

Example 2. Consider now the function

u(r,ϕ) = r2π/βSin(2πϕ/β), (r,ϕ) ∈ ΩRβ , (6)

which evidently belongs to H2(ΩRβ ), ∀β ∈ (0,2π). This function satisfies the Laplace equation

and the boundary conditions (1). Observe that the angle α= 2πϕ/β always remains in (0,2π),

when the parameter β changes in (0,2π). This moment removes the lack of smoothness of the

boundary and as a result the solution u(r,ϕ) ∈ H2(ΩRβ )∩ H1(ΩRβ), ∀β ∈ (0,π), given by (6),

also satisfies the periodicity condition (5).

These two solutions show the main distinguished features of the boundary value problem

(1)-(2), and they will be used for testing of the presented finite difference scheme.

3. The Conservative FD Scheme on a Piecewise Uniform Polar Mesh

We assume here β = 2π and introduce the following uniform meshes with respect to

variables r and ϕ

w r := {rn = (n− 0.5)hr : n= 1,2, . . . , N + 1, hr = (2R)/(2N + 1)},
wϕ := {ϕm = (m− 1)hϕ : m= 1,2, . . . , M + 1, hϕ = 2π/M},

with mesh steps hr ,hϕ > 0. Then we obtain the piecewise uniform polar mesh w rϕ := wr×wϕ
(Fig. 1):

w rϕ := {(rn,ϕm) ∈ ΩRϕ : rn ∈ w r , ϕm ∈ wϕ}, dim wrϕ = (N + 1)× (M + 1),

where w rϕ := wrϕ ∪ cγϕ ∪ γ0, and wrϕ := {(rn,ϕm) ∈ ΩRϕ : n = 2, N , m = 2, M}. The

boundary mesh points are defined as follows

γϕ := {(R,ϕm) ∈ Γϕ : m= 1, M}, γ0 := {(rn, 0) ∈ Γ0 : n= 1, N}.



A. Erdem / Eur. J. Pure Appl. Math, 6 (2013), 30-43 34

1 

23
4

5

6
7

8 

10

11

12

13

14

15

16

9

18

19

20

21

22

23

24

17β
0
 

γ
0
 

γφ 

rotated mesh 

Figure 1: Geometry of the poleness polar mesh and its rotated form

Due to the periodicity condition u(r, 0) = u(r, 2π) we will not include the values

u(r, 0), r ∈ [0,R] to the list of unknowns in the discrete problem. The introduced mesh

w rϕ doesn’t include the pole point r = 0, that is in the presented discrete model the domain

ΩRϕ is approximated by the circular disc
o
ΩRϕ. Thus our discrete model does not deal with

the singularity at r = 0, that is usual for the differential problem. Instead we will derive an

approximation of the boundedness condition (2) at the central circle with radius

r = r1, r1 = 0.5hr .

Denote by enm = {(r,ϕ) ∈ ΩRβ : rn ≤ r ≤ rn+1, ϕm ≤ ϕ ≤ ϕm+1} the polar finite element

with four nodes. We derive an error approximation for each element. Introducing the half-

nodes r±n = rn ± hr/2, ϕ±m = ϕm ± hϕ/2 and integrating equation (1) on the finite element

eenm := [r−n , r+n ]× [ϕ−m,ϕ+m] we obtain the following balance equation:

∫ ϕ+m

ϕ−m

∫ r+n

r−n

∂

∂ r

�
r
∂ u

∂ r

�
drdϕ+

∫ r+n

r−n

∫ ϕ+m

ϕ−m

1

r

∂ 2u

∂ ϕ2
dϕdr = −
∫ r+n

r−n

∫ ϕ+m

ϕ−m

rF(r,ϕ)dϕdr. (7)

Let us transform the first left integral I r
nm.

I r
nm =

∫ ϕ+m

ϕ−m

�
r
∂ u

∂ r

�r=r+n

r=r−n
dϕ ∼= hϕ

��
r
∂ u

∂ r

�
(r+n ,ϕm)−
�

r
∂ u

∂ r

�
(r−n ,ϕm)

�

We use here the central finite difference formula for approximation of derivatives on the right

hand side, by using the mesh points rn, rn+1/2, rn+1, with mesh step hr/2 = hr/2:

�
r
∂ u

∂ r

�
(r+n ,ϕm)

∼= r+n
u(rn+1,ϕm)− u(rn,ϕm)

hr

,

�
r
∂ u

∂ r

�
(r−n ,ϕm)

∼= r−n
u(rn,ϕm)− u(rn−1,ϕm)

hr

.



A. Erdem / Eur. J. Pure Appl. Math, 6 (2013), 30-43 35

Then we have the following variational finite difference approximation of the integral operator

I r
nm:

I r
nm
∼= hϕ
�

r+n ur,nm− r−n ur ,nm

�
.

By the same way we can derive an approximation of the second integral operator I
ϕ
nm on the

left hand side of (7):

Iϕnmu ∼= hr

rn

�
∂ u

∂ ϕ
(rn,ϕ+m)−

∂ u

∂ ϕ
(rn,ϕ−m)
�

∼= hr

rn

�
uϕ(rn,ϕm)− uϕ(rn,ϕm)

�
.

Applying to the right hand side of (7) the numerical integration (rectangle) formula, finally

we have

−hϕ[r
+
n yr,nm− r−n yr ,nm]−

hr

rn

[yϕ,nm − yϕ,nm] = hrhϕ rnF(rnϕm)

Dividing by hrhϕ rn > 0 we obtain the following finite difference equation

A nmu := −
�

1

r
(r yr)

�

r,nm

−
�

1

r2
yϕϕ

�

nm

= F(rn,ϕm), (rn,ϕm) ∈ωrϕ, n 6= 1. (8)

The finite dimensional operators

A r
nm y :=

�
1

r
(r yr)

�

r,nm

, A ϕ
nm y

�
1

r2
yϕϕ

�

nm

are the finite difference approximations of the differential operators

Aru :=
1

r

∂

∂ r

�
k(r)

∂ u

∂ r

�
, Aϕu :=

1

r2

∂ 2u

∂ ϕ2
, (r,ϕ) ∈ ΩRβ ,

correspondingly. Note that the same approximations can also be obtained from the direct

finite difference approximation of the Poisson equation (1).

The finite difference equation corresponding to the layer r = r1 = hr/2 can be derived by

using the same balance equation (7), substituting r−n = ǫ, r+n = hr :

∫ ϕ+m

ϕ−m

∫ hr

ǫ

∂

∂ r

�
r
∂ u

∂ r

�
drdϕ+

∫ ϕ+m

ϕ−m

∫ hr

ǫ

1

r

∂ 2u

∂ ϕ2
drdϕ = −
∫ ϕ+m

ϕ−m

∫ hr

ǫ

rF(r,ϕ)drdϕ.

Going to the limit ǫ→ 0 here and using condition (2) we obtain

∫ ϕ+m

ϕ−m

∂ u(hr ,ϕ)

∂ r
dϕ+

∫ hr

0

1

r

�
∂ u(r,ϕ+m)

∂ ϕ
− ∂ u(r,ϕ−m)

∂ ϕ

�
dr +

∫ hr

0

r

∫ ϕ+m

ϕ−m

F(r,ϕ)dϕdr = 0.



A. Erdem / Eur. J. Pure Appl. Math, 6 (2013), 30-43 36

Again using the numerical integration formula in the first integral we get

hrhϕ
∂ u(hr ,ϕm)

∂ r
+ 2

�
∂ u

∂ ϕ

�
hr

2
,ϕm+

hϕ

2

�
− ∂ u

∂ ϕ

�
hr

2
,ϕm −

hϕ

2

��
+

h2
rhϕ

2
F

�
hr

2
,ϕm

�
≈ 0

To derive this finite difference equation in canonical form we use the standard definitions

yr,1m :=
1

hr

[y(r1 + hr ,ϕm)− y(r1,ϕm)],

yϕ,1m :=
1

hϕ
[y(r1,ϕm)− y(r1,ϕm−1)],

yϕ,1m :=
1

hϕ
[y(r1,ϕm+1)− y(r1,ϕm)].

Then we have

hrhϕ yr(r1,ϕm) + 2hϕ yϕ,ϕ(r1,ϕm) +
h2

rhϕ

2
F(r1,ϕm) = 0, m = 1,2, . . . , M

Dividing the both sides to h2
rhϕ/2 finally we obtain the finite difference equation correspond-

ing to the layer r1 = hr/2:

− 2

hr

yr(r1,ϕm)−
4

h2
r

yϕ,ϕ(r1,ϕm) = F(r1,ϕm) = 0, m= 1,2, . . . , M . (9)

Equations (8)-(9) represent the finite difference analogue of the Poisson equation (1) in the

constructed polar mesh ωrϕ.

Lemma 1. If u ∈ C4(ΩRβ) then the order of the approximation error of the finite difference

schemes (8)-(9) in C-norm is ψrϕ :ψrϕ = O
�
(h2

r + h2
ϕ)/rn

�
, where

ψrϕ =
h2

r

6rm


∂

3u

∂ r3
+

rn + hr/2

4

∂ 4u(ern,ϕm)

∂ r4
+

rn− hr/2

4

∂ 4u(eern,ϕ)

∂ r4


+

h2
ϕ

12rm

∂ 4u(rn, eϕm)

∂ ϕ4
.

(10)

The explicit form (10) of ψ(hr ,ϕm) and the proof of this result is given in [11].

The lemma shows that as r → r1 the function ψrϕ increases, and at the first layer r = r1

(r1 = hr/2) becomes ψrϕ = O(hr + h2
ϕ/hr).

Note that the boundedness condition (2) is necessary for the above approximation, and

hence for the convergence of the finite difference schemes (8)-(9), although this condition is

not used explicitly on deriving the approximation error. To show this, consider the following:

Example 3. The function

u(r,ϕ) = ln
1

r
Sinϕ, (r,ϕ) ∈ ΩR,ϕ,



A. Erdem / Eur. J. Pure Appl. Math, 6 (2013), 30-43 37

satisfies the Poisson equation (1) in ΩRβ , with for R = 1, β = 2π, and the right hand side

F(r,ϕ) = 1

r2 Sinϕ. Evidently this solution also satisfies the boundary and periodicity conditions

(1), but does not satisfy the boundedness condition (2), since r∂ u/∂ r = Sinϕ 6→ 0, as r → 0.

Calculating the right hand side of (10) for n= 1 we obviously observe that for r = r1

ψ(r1,ϕm) =
h2
ϕ

12r3
1

sinϕm.

This shows that ψ(r1,ϕm) 6→ 0, as r1→ 0, and there is no approximation in the neighbourhood

of the pole point r = 0.

To formulate the discrete problem we need to add to equations (8)-(9), the equations,

obtained from the periodicity condition (1).

4. The Algebraic Problem with Nonsymmetric Sparse Matrix: Iteration

Algorithm

The finite difference equations (8)-(9) compose K := N M number of algebraic equations

with K unknowns

y = (y11, y12, . . . , y1M , y21, . . . , yN M )
T , dim y = K .

We can rewrite these equations in the form of the system of linear algebraic equations

A y =F with the following positive band matrix A , dimA = K × K:

A :=




A11 A12 0 0 . . . 0 0 0

A21 A22 A23 0 . . . 0 0 0

0 A32 A33 A34 . . . 0 0 0

. . . . . . . . . . . .
. . . . . . . . . . . .

0 0 0 0 . . . AN−1,N−2 AN−1,N−1 AN−1,N

0 0 0 0 . . . 0 AN ,N−1 AN ,N




.

Here M ×M -dimensional block matrices Ai j are of the following structure:

Aii =




• • 0 0 . . . 0 0 •
• • • 0 . . . 0 0 0

0 • • • . . . 0 0 0

. . . . . . . . . . . .
. . . . . . . . . . . .

0 0 0 0 . . . • • •
• 0 0 0 . . . 0 • •




and

Aii+1 = αi I , Ai−1i = βi I .



A. Erdem / Eur. J. Pure Appl. Math, 6 (2013), 30-43 38

The parameters αi and βi can be defined from the finite difference schemes (8)-(9) as follows:

α1 =
2

h2
r

; αi =
ki

rih
2
r

, i = 2, N ; βi =
ki+1

rih
2
r

, i = 1, N .

Since αi 6= βi, the matrix A is not symmetric one. Moreover, A is a sparse matrix of spe-

cial structure, corresponding to the polar mesh and periodicity condition. Evidently such a

system of linear algebraic equations with non-symmetric sparse matrix needs to be solved by

iteration methods [12, 13, 14]. However, as the computational experiments below show, use

of compact storage for the band matrix A and then an application of any effective iteration

method requires large enough time for the solution of the linear system of algebraic equations

A y =F . The reason is that the bandwidth of the non-symmetric matrix A is

bw = 3M

and there are many null terms in the band. Specifically, the above block matrices Aii, Aii+1,

Ai−1i contain M2 − 3M , M2 −M , M2 −M zero elements, correspondingly. Hence for M = N

the number of zero and non-zero elements of the band are 3N3 − 4N2 − 4N and 5N2 + N ,

respectively. For the mesh with N = 30, this means that the number of non-zero elements of

the band matrix A is about 4.3% of all band elements. Therefore one needs to construct a

special algorithm which can store and operate with only nonzero elements αi and βi, realising

the multiplications U y and Ly, to minimize the time required for the solution of the consid-

ered problem by any iteration method. Here the U and L are the upper and lower triangular

matrices: A = D+U+ L. Note that for some class of systems, arising from the finite-element

discretization, similar algorithm was constructed in [15].

Table 1 illustrates the comparative analysis of the standard SOR method

yk+1 = (D+ωL)−1[(1−ω)D−ωU]yk +ω(D+ωL)−1F , (11)

by using MATLAB code, and SOR method with the constructed here special algorithm. As a

test example the analytical solution given in Example 2, with β = 2π, is used. The iteration

parameter ω ∈ (0,2) in (11) was defined by the formula

ω =
2

1+
p
λmin(2−λmin)

,

where λmin is the minimal eigenvalue of the Laplace operator, and λmin = 2sin2(π/2N), for

the square mesh M = N .

Table 1 shows that direct application of the SOR method is expensive in the sense of

the required CPU time. This time increases as the number K = N M of mesh points increases.

Computational results show that this increase has the character 2n, i.e. n-times increase of the

number of mesh points leads 2n-times increase of the CPU time. This is due to the moderately

ill-conditionedness, according to [12], of the matrix A , as the fourth column of the table

shows. At the same time, as the second column of Table 1 shows, the SOR method with the

special algorithm requires less than 1 second for all considered meshes.



A. Erdem / Eur. J. Pure Appl. Math, 6 (2013), 30-43 39

Table 1: Comparison of the standard and modified SOR algorithms

N ×M SOR with special algorithm Standard SOR Condition number

CPU time(sec.) CPU time(sec.) of the matrix A
20× 20 0 9 1.5× 104

30× 30 0 48 8.0× 104

40× 40 0 126 2.5× 105

5. Numerical Solution of the Elliptic Problem (1)-(2) by the Poleness

Conservative Schemes (8)-(9)

In this section we discuss results of computational experiments related to numerical so-

lution of the Dirichlet problem (1)-(2) by the poleness conservative scheme (8)-(9). In the

first series of the computational experiments, the convergence and accuracy of the numerical

solution of the Dirichlet problem (1)-(2), with the analytical solution

u(r,ϕ) = (r − 1)Sinϕ, (r,ϕ) ∈ ΩRβ , R= 1,

is studied. Note that the corresponding source function

F(r,ϕ) =
1

r2
Sinϕ, (r,ϕ) ∈ ΩRβ , (12)

has singularity at r = 0. The two appropriate iterative methods - SOR method and Bi-

Congugate Gradient (BCG) method - are applied for the iterative solution of the linear system

of algebraic equations, corresponding to the finite difference schemes (8)-(9). In all cases the

MATLAB codes of this methods with the above mentioned special algorithm is used. Results

are presented in the Table 2. For the comparison, the numerical results obtained by the Gauss-

Seidel method, are also presented. Here and below the value of the stopping parameter δ > 0

in

‖uni−1 − uni‖0 ≤ δ,

where ‖ · ‖0 is the L2-norm, is taken δ = 10−5.

Table 2 shows numbers of iterations corresponding to all three iteration methods, with

H1-relative and L∞-absolute errors. Results given in the table show that for relatively coarse

meshes BGC method is more effective than the SOR method. But for the meshes 45 × 45

and higher, SOR method is more effective in the sense of iterations. Moreover, the number of

iterations ni in SOR method increases slowly, as increases the number of mesh points. Thus

ni = 234 and ni = 277, for the meshes 40× 40 and 50× 50, respectively. As show the last

three columns of the table, the H1-relative and L∞-absolute errors are small enough, although

the source function (12) has singularity at the pole point r = 0.

The next series of computational experiments is realized for the smooth continuous source

function

F̃(x , y) = F(r,ϕ) =

(
50ex p(− ǫ2

ǫ2−r2 ), 0< r < ǫ

0, ǫ < r < R, R= 1,
(13)



A. Erdem / Eur. J. Pure Appl. Math, 6 (2013), 30-43 40

Table 2: Comparative analysis of the modified SOR and BGC methods applied to the poleness

conservative finite difference scheme (8)-(9)

Methods N ×M Time Number of Rel. error Abs. error Abs. error

(sec.) iterations H1-norm L∞-norm L∞-norm

ni r = r1 r = R/2

30× 30 0 185 8.0× 10−6 3.8× 10−3 1.3× 10−3

SOR 40× 40 1 234 2.7× 10−6 2.5× 10−3 8.1× 10−4

50× 50 2 277 1.4× 10−6 1.9× 10−3 5.9× 10−4

30× 30 0 110 8.0× 10−6 3.4× 10−3 1.2× 10−3

BCG 40× 40 2 224 2.5× 10−6 2.0× 10−3 6.8× 10−4

50× 50 4 385 1.0× 10−6 1.3× 10−3 4.3× 10−4

30× 30 33 473 8.5× 10−6 4.9× 10−3 1.2× 10−3

Gauss- 40× 40 95 762 5.0× 10−6 3.0× 10−3 7.8× 10−4

Seidel 50× 50 219 1105 1.1× 10−5 2.0× 10−3 1.6× 10−3

with ǫ = 1/2, approximating in weak sense the Dirac δ-function (Figure 2a). This function is

taken as a given data for the Dirichlet problem (1)-(2). The right pane, Figure 2b, illustrates

the numerical solution ũh(x , y) = uh(r,ϕ) of problem (1)-(2) by the poleness conservative

schemes (8)-(9), for the mesh size 30× 30.

Finally we consider the weak solution of the Dirichlet problem (1)-(2), when the source

F(r,ϕ) is a discontinuous at r = 1/2 function

F̃(x , y) = F(r,ϕ) =





50ǫ
exp(−ǫ2/4r)p

4πr3
, 0< r < ǫ

50ǫ
exp(−ǫ2/4)p

4π
, ǫ < r < R, R= 1,

given in the left pane of Figure 3. This function with ǫ = 1/2 is taken as a given data

for the Dirichlet problem (1)-(2). The numerical solution ũh(x , y) = uh(r,ϕ) of problem

(1)-(2) obtained for the mesh 30× 30 is plotted in the right pane, Figure 3b. To estimate

an accuracy of the numerical solution, in particular at the discontinuity point r = 1/2, the

numerical solutions u
(1)
h
(r,ϕ) = u

(2)
h
(r,ϕ) corresponding to two different meshes w(1)

rϕ = w(2)rϕ

is compared. The relative error

ǫh =







u

(1)

h
(r,ϕ)− u

(2)

h
(r,ϕ)

0.5
�

u
(1)
h
(r,ϕ) + u

(2)

h
(r,ϕ)
�







∞

,

is about ǫh = 10−3÷ 10−3, including the discontinuity point. This shows high accuracy of the

numerical method in the case of discontinuous source function, also.



A. Erdem / Eur. J. Pure Appl. Math, 6 (2013), 30-43 41

(a) Continuous Source Solution (b) Numerical Solution

Figure 2: Dirichlet problem in polar coordinates

(a) Discontinuous Source Solution (b) Numerical Solution

Figure 3: Dirichlet problem in polar coordinates



REFERENCES 42

6. Conclusion

We studied poleness conservative finite difference scheme for Laplace operator in polar

coordinates. The scheme with the modification of the SOR method allows to construct an

effective numerical method for solving the Dirichlet problem in the polar coordinates, based

on the weak solution approach. Numerical results presented for discontinuous source function

shows high accuracy of the method on acceptable meshes.

Extension of results given here can be made for positive elliptic operators with discontin-

uous coefficients, and for nonlinear monotone operators of Plateau type, as well. This require

some additional techniques that will be done in next studies.

References

[1] S. Agmon, A. Douglis, and L. Nirenberg. Estimates near the boundary for solutions of ellip-

tic partial differential equations satisfying general boundary conditions. Communications

on Pure and Applied Mathematics. 17, 35-92. 1964.

[2] A.M. Bruaset. A Survey of Preconditioned Iterative Methods. Addison-Wesley. 1995.

[3] J.D. Jackson. Classical Electrodynamics. 2nd Ed., Wiley, New York. 1975.

[4] M.D. Griffin, E. Jones, and J.D. Anderson. A computational fluid dynamic technique valid

at the centerline for non-axisymmetric problems in cylindrical coordinates. Journal of Com-

putational Physics. 30, 352-364. 1979.

[5] W. Hackbush. Iterative Solution of Large Sparse Systems of Equations. Springer-Verlag,

Berlin. 1994.

[6] W. Huang and D.M. Sloan. Pole condition for singular problems: The pseudo-spectral ap-

proximation. Journal of Computational Physics. 107, 254-365. 1993.

[7] P.D. Lax and B. Wendroff. Systems of conservation laws. Communications on Pure and

Applied Mathematics. 13, 217-237. 1960.

[8] A.I. Lurie. Three-Dimensional Problems of the Theory of Elasticity. Interscience Publishers,

New York. 1964.

[9] P.E. Merilees. The pseudospectral approximation applied to the shallow water equations on

a sphere. Atmosphere, 13(1), 897-910.1973.

[10] K. Mohseni and T. Colonius. Numerical treatment of polar coordinate singularities. Journal

of Computational Physics. 157, 787-795. 2000.

[11] M. Renardy and J. Rogers. Introduction to Partial Differential Equations. Springer-Verlag,

New York, 1993.



REFERENCES 43

[12] A.A. Samarskii and V.B. Andreev. Difference Methods for Elliptic Problems (in Russian).

Nauka, Moscow. 1976.

[13] F.S. Sharman. Viscous Flow. McGraw-Hill, New York. 1990.

[14] R. Verzicco and P. Orlandi. A finite difference scheme for three dimensional incompressible

flows in cylindrical coordinates. Journal of Computational Physics. 123, 402-415. 1996.

[15] F.M. White. Viscous Fluid Flow. McGraw-Hill, New York. 1991.