Two Full Multigrid Algorithms in Cartesian and Polar

Coordinate Systems to Solve an Elliptic Partial Differential

Equation with a Mixed Derivative Term in a Quarter of
a unit Circle
Osama El-Giara∗

Abstract
In this paper we introduce two full multigrid algorithms in a Cartesian and a Polar coor-
dinate system that have been modified to treat the elliptic partial differential equations with
a mixed derivative term, in the two dimensional space in a quarter of a unit circle, special
modification for the points in the neighborhood of the curved boundary has been done. An
advantage of the two algorithms is that the extension from two to three space dimensions is
straightforward. numerical examples have been given.
Index terms— Multigrid method, Numerical Analysis, Elliptic PDE, Cartesian and Polar Co-
ordinate Systems.
1 Full Multigrid Algorithms in Cartesian Coordinates
1.1 Introduction
Consider the general linear elliptic partial differential equation:
2
LU = AUxx + 2BUxy + CUyy + F = 0 ∈Ω⊂R (1.1a)
2
with Dirichlet boundary conditions U = F ∈ ∂Ω ⊂ R (1.1b)
where ∂Ω is a quarter of a unit circle. p A difference approximation Lh of L in equation 1.1a
is said to be of order p if Lh(U) = O(h ), where h is the mesh size. In our case p = 2, so we
shall use the 6-point difference equation to the mixed derivative term Uxy to avoid the absent
points in our grids because of the curved boundary see [3] and [4]

1.2 Finite Difference Discretizaton

For each grid we replace each term in equation 1.1a by a second order finite difference approx-
imation for two kinds of points as follows:
1. Normal points: These points are interior points are not laying near the curved boundary:
for a grid with mesh size is h, we have the following discretizaton:
Uxx = h−2 ( ui+1,j + ui−1,j − 2ui,j )
Uyy = h−2 ( ui,j+1 + ui,j−1 − 2ui,j ) (1.2)
Uxy = h−2 ( ui,j − ui−1,j − ui,j−1 + ui−1,j−1 )
substitute these values in equation 1.1a we get the following system of linear equations for normal
points:
 - 2(a  c - b) c 0 0 
 ( c  2b ) - 2(a  c - b) c
...
0   u 1   b1 
 ...  u2   b2 
 0 (c - 2b) - 2(a  c - b) 0    
 c        (1.3)
 ... ... ...
- 2(a  c - b)
0      
 0 ... (c - 2b) c    
 (c - 2b)   u N 1   b N 1 
 0 ... 0 - 2(a  c - b) 

∗aDepartment of Basic Science, Modern Academy for Engineering and Technology in Maadi, Cairo, Egypt.
osama996@hotmail.com

1
 2. Curved boundary points: These points are points that laying near the curved boundary:
see [4] and Figure 1 When the boundary is curved and intersects the rectangular grid at
points that are not grid points, so we shall

Figure 1: Curved mesh

give the finite difference approximation to the derivatives at the points in the neighborhood
of the curved boundary thus:
 2 2 2 
U  h 2  u j 1, j  u i 1, j  u
xx
  1    (1   )  i , j 
 2 2 2  (1.4)
U  h 2  u i , j 1  u i , j 1  u i, j 
yy
  (1   ) 1     
−2
Uxy = h (ui,j − ui−1,j − ui,j−1 + ui−1,j−1)

substituting 1.4 into 1.1a we get a system of linear equations for the points in the
neighborhood of the curved boundary given by:

 2 ai , j   a   2 ci, j 
  u i  1 , j  2  i , j  b i , j  u i  1 , j    u i , j 1
  1     1     1    
 ci, j
 2  bi, j
  a c
 u i , j  1  2  i , j  i , j  b i , j

 u i , j  2 b i , j u i  1 , j  1  h 2 f i , j
(1.5)
1      

where 1  i,j  N – 1 , h  1
, N is the number of interior points in each direction
N 1
 h,  h are the distances from the point in the neighborhood of the boundary and the
curve boundary in direction of the x-coordinate and the y-coordinate respectively.

1.3 Multigrid Method

n
A sequence of 5 grids Ωkh have been used with mesh sizes kh, k = 2 , n = 0, 1, · · · 4,
1
where h = 64 . The elliptic equation have been solved using full multigrid V (γ1, γ2)-cycle
method with the coarse grid correction (CGC), see Figure 2 and has the following components:
1. Coarse to fine interpolation with operator I
h
2h
: within the V-cycle, we will use correction
interpolation where corrections are to be interpolated from coarser grids and added to finer
grids and the full multigrid interpolation where current values of the approximate solution
on the coarser grid are interpolated as the first approximation to the solution (on the first
visit to the grid) on the finer grids, see [2], we choose for both cases of interpolations the
bilinear interpolation. Now we have two kinds of points see figure 3 :

2
 Figure 2: FMG V(2, 1)

(a) Normal points:

h 2h
u 2 i , 2 j  u i , j for common points in the two grids, (point 13)
u
h
2 i 1 , 2 j

1
2
u 2h
i, j
 u
2h
i 1 , j
 for odd points in the x direction, (point 12)
u
h
2 i , 2 j 1
  u  u  for odd points in the y direction, (point 8)
1
2
2h
i, j
2h
i , j 1

 u u  for remaining points in, (point 7)

h 1 2h 2h 2h 2h
(1.6)
u 2 i 1 , 2 j 1
4
u u
i, j i 1 , j i , j 1 i 1 , j 1

(b) Curved boundary points:

Since we have used red- black gauss-Seidel iteration, then we need interpolation only
for odd points in the x and y directions respectively defined by

u
h
2 i 1 , 2 j

1
 h
 u
2h
i, j
 hu
2h
i 1 , j
for odd points in the x direction, (point 14)

u
h

1 
 u
2h
 hu
2h  for odd points in the  direction, (point 19) (1.7)
2 i , 2 j 1
 h i, j i , j 1


Figure 3: Numbering two grids

2h
2. Fine to coarse restriction operator: The restriction operator denoted by I h takes the
residual equation computed on the fine-grid and transfer it to the coarse-grid according

3
 to the rule: I h2 h R h  R 2 h where the components of Rij2 h of R 2 h are given by the following
two schemes- see [2] and [3]:
(a) Half-weight: If B = 0 in equation 1.1a, we use the red-black Gauss-Seidel relaxation
it is better to use the half-weight operator since the residual at the black (odd) points
must vanish and the stencil of the half-weight operator takes the form:

2h
0 1 0
1 ( 1.8 )
  1 1 
2h
I h
8
4
 0 1 0  h

(b) Full-weight: If B = 0 in equation 1.1a, we use the red-black Gauss-Seidel relaxation,

so we may use either the half-weight operator or the full-weight operator, or it is
better to use the half-weight operator for normal points and full-weight operator for
curved boundary points. The stencil of the full-weight operator takes the form:

2h
1 2 1
1 
2  ( 1.9 )
2h
I h

16 
2 4
 1 2 1  h

1.4 Test problems, Results and Performance

In this section we report the results and performance of some test problems All our implemen-
tations in Fortran were executed on a pentium 4 PC using FORTRAN90 workstation compiler.
The accuracy of the method is measured by both L2 norm of both defect and error and the
max norm of both defect and the error. The general formula for the test problem is given by:

2
AUxx + 2BUxy + CUyy + F = 0 ∈Ω⊂R

Consider the following test problems:

1. Test Problem 1: Let A = C = 1, B = 0 using half-weight restriction:

Uexact = e2xy , F (x,y) = 4(x2 + y2)e2xy (1.10a)

Uexact = sin(3x + y), F (x,y) = −10 sin(3x + y) (1.10b)

2. Test Problem 2: Let A = C = B = 1, using half-weight restriction:

1 2xy
Uexact = e2xy , F (x,y) = 4(x2 + y2 + xy + )e (1.11a)
2
Uexact = sin(3x + y), F (x, y) = −13 sin(3x + y)
(1.11b)
3. Test Problem 3: Let A = C = B = 1, using full-weight restriction:

1 2xy
Uexact = e2xy , F (x,y) = 4(x2 + y2 + xy + )e (1.12a)
Uexact = sin(3x + y), F (x, y) = −13 sin(3x + y) 2 (1.12b)

4
 Test Problem 1: Equation 1.10a and equation 1.10b

Test problem h M axN (R) L2(R) M axN (E) L2(E)

1
2 0.917481D − 07 0.917481D − 07 0.432282D − 02 0.432228D − 02
1
4 0.185526D − 02 0.796756D − 03 0.220622D − 02 0.130441D − 02
1
8 0.121447D − 02 0.333357D − 03 0.102566D − 02 0.440546D − 03
Equation 1.10a 1
16 0.316902D − 03 0.885846D − 04 0.365635D − 03 0.186466D − 03
1
32 0.847075D − 04 0.208229D − 04 0.147426D − 03 0.737563D − 04
1
64 0.298698D − 04 0.662750D − 05 0.484918D − 04 0.248272D − 04
1
2 0.526033D − 07 0.526033D − 07 0.335108D − 01 0.335108D − 01
1
4 0.535409D − 02 0.238619D − 02 0.137040D − 01 0.889376D − 02
Equation 1.10b 1
8 0.203170D − 02 0.563616D − 03 0.420590D − 02 0.217381D − 02
1
16 0.412457D − 03 0.112770D − 03 0.100574D − 02 0.499224D − 03
1
32 0.123807D − 03 0.221409D − 04 0.247875D − 03 0.117913D − 03
1
64 0.472992D − 04 0.439256D − 05 0.605037D − 04 0.281259D − 04

Table 1: Fortran implementation of test problem 1 using F M V (2, 1) on PC

Test Problem 2: Equation 1.11a and equation 1.11b

Test problem h M axN (R) L2(R) M axN (E) L2(E)

1
2 0.553447D − 07 0.553447D − 07 0.136968D + 00 0.136968D + 00
1
4 0.150842D − 02 0.518837D − 03 0.743628D − 01 0.487242D − 01
1
Equation 1.11a 8 0.115986D − 02 0.359802D − 03 0.387126D − 01 0.211317D − 01
1
16
0.522688D − 03 0.131486D − 03 0.193226D − 01 0.979662D − 02
1 0.128848D − 03 0.228685D − 04
32 0.955943D − 02 0.465064D − 02
1 0.489233D − 04 0.461051D − 05 0.474638D − 02 0.226454D − 02
64
1 0.357284D − 07 0.357284D − 07 0.666329D − 01 0.666329D − 01
2
1 0.172004D − 02 0.595731D − 03 0.454477D − 01 0.240583D − 01
4
Equation 1.11b 1 0.218721D − 02 0.689057D − 03 0.150908D − 01 0.718392D − 02
8
1 0.501558D − 03 0.105073D − 03 0.567628D − 02 0.260250D − 02
16
1 0.174007D − 03 0.172096D − 04 0.242653D − 02 0.114201D − 02
32
1 0.595715D − 04 0.301087D − 05 0.124441D − 02 0.544308D − 03
64

Table 2: Fortran implementation of test problem 2 using F M V (2, 1) on PC

5
 Test Problem 3: Equation 1.12a and equation 1.12b

Test problem h M axN (R) L2(R) M axN (E) L2(E)

1
2 0.553447D − 07 0.553447D − 07 0.136968D + 00 0.136968D + 00
1
4 0.150842D − 02 0.518837D − 03 0.743628D − 01 0.487242D − 01
1
Equation 1.12a 8 0.116307D − 02 0.278797D − 03 0.385639D − 01 0.211171D − 01
1
16 0.452654D − 03 0.931737D − 04 0.198604D − 01 0.100401D − 01
1
32 0.117215D − 03 0.158437D − 04 0.992671D − 02 0.482168D − 02
1
64 0.249014D − 04 0.365150D − 05 0.494005D − 02 0.235517D − 02
1
2 0.357284D − 07 0.357284D − 07 0.666329D − 01 0.666329D − 01
1
4 0.172004D − 02 0.595731D − 03 0.454477D − 01 0.240583D − 01
Equation 1.12b 1
8 0.104549D − 02 0.408385D − 03 0.142263D − 01 0.687403D − 02
1
16 0.355017D − 03 0.855387D − 04 0.577923D − 02 0.265536D − 02
1
32 0.105237D − 03 0.139371D − 04 0.254165D − 02 0.118886D − 02
1
64 0.257567D − 04 0.255912D − 05 0.129553D − 02 0.566869D − 03

Table 3: Fortran implementation of test problem 3 using F M V (2, 1) on PC

2 Full Multigrid Algorithms in Polar Coordinates

2.1 Introduction
A multigrid algorithm has been modified to treat the elliptic boundary value problem for
partial differential equations with a mixed derivative term in two dimensional space in polar
coordinates. A full multigrid algorithm in which W (γ1, γ2) with a line relaxation scheme are
used to deal with the finite difference approximation of the mixed derivative term in polar
coordinates. Full weight restriction of residuals and linear interpolation operators are used.
Numerical examples are given.
Consider the elliptic partial differential equation with a mixed derivative term Urθ in polar
1
coordinates which is equivalent to equation 1.1 with A = C = 1 , B = :
2

1 1
( 1 + sin  cos  ) Urr + ( 1 – sin  cos  ) Ur + 2 ( 1 – sin  cos  ) U 
r r


1
r2
  1
 
sin 2   cos 2  U   cos 2   sin 2  U r   F ( r ,  ),
r
 (2.1)

where   r , ,0  r  1,0  0    , subject to the dirchlit boundary conditions:

 4

U  g   (2.2)

2.2 Finite Difference Discretizaton

For each grid we replace each term in equation 2.1 by a second order finite approximation as
follows - see Figure 4
Urr = h−2 (ui+1,j + ui−1,j − 2ui,j )
Ur = (2h)−1 (ui+1,j − ui−1,j )
Uθθ = (δθ)−2 (ui,j+1 + ui,j−1 − 2ui,j ) (2.3)
−1
Uθ = (2δθ) (ui,j+1 − ui,j−1) and
−1
Urθ = (4hδθ) (ui+1,j+1 + ui−1,j−1 − ui+1,j−1 − ui−1,j+1)

6
 Figure 4 : Descritization of a quarter of a unit circle
substitute these values in equation 2.1 we get the following system of linear equations:
D1 B1 0 ... 0 
C   u1   b1 
 1 D 1 B1 ... 0   u   b 
 0 C  D B1 0   2   2  ( 2.4 )

1 1
       
 ... ... ... ... 0     
 0 ... C  D B1       

1 1
  u N  1   b N  1 
 0 ... 0 C1  D 1 

1 1
where 1  I, j  N – 1 , h ,  = , are the mesh size and
N 1 2 N  1
D1 = 2 1  1 sin 2   1 1 
 1  sin 2   

 2 i  2
 2 
B1= 1  1  cos 2  ( 2.5 )
 1  sin 2   
i  2  2  2 i 2 
C1= 1  1  cos 2 
 1  sin 2   
i  2  2  2 i 2 
where N is the number of interior points in each direction

2.3 Multigrid Method

n
A sequence of 5 grids Ωkh have been used with mesh sizes kh, k = 2 , n = 0, 1, · · · 4,
1
Where h = .The elliptic equation have been solved using full multigrid W-cycle method with the
64
coarse grid correction (CGC), see Figure 5

Figure 5: FMG W-Cycle

and has the following components:

1. Line relaxation: Since the finite difference approximation of the mixed derivative term
contains the four corner points, then a scheme of line-relaxation should be used, (first we
may first relax the solution at even (odd) lines, then relax at the odd (even) lines, this
gives the tridiagonal system of linear equation 2.4.
2. Interpolation operator I h2 h : It takes coarse-grid vectors u2h defined on G2h and produces
fine grid vectors uh defined on Gh using I
h
2h
u 2 h = u h , the components of uh are divided
into 4 groups as follows- see figure 4:
( a ) Common points of the two grids
u2hi , 2 j  ui2, hj ( point 5 ) ( 2.6 )

7
 ( b ) Radial points

u
h
2 i 1, 2 j

1
2
u h
2 i 1 , 2 j
 u
2h
i 1 , j
 (points 4 , 6) (2.7)

( c ) Angle points

u
h
2 i , 2 j 1

1
4
u h
2 i 1, 2 j
h h
 u 2 i 1, 2 j  u 2 i 1, 2 j  2  u 2 i 1, 2 j  2
h
 (points 2 , 8) (2.8)

where the right hand side is obtained from 2.7.

( d ) Intermediate points: They are two kinds:
i. Points are not in the neighborhood of the center of the circle

u
h
2 i 1 , 2 j 1

1
4
u 2h
i, j
2h 2h 2h
 u i  1 , j  u i , j 1  u i 1 , j 1  (points 3, 9) (2.9)

ii. Points are in the neighborhood of the center of the circle

u
h
2 i 1, 2 j 1

1
3
u 2h
i, j
2h
 ui 1, j  ui 1, j 1
2h
 (points 1, 7) (2.10)

3. Restriction operator I h2 h :It takes fine grid vectors Rh defined on Gh and produces coarse
grid vectors R2h defined on G2h using I h2 h R h = R 2 h since we use the line relaxation for
every grid it is better to use the full-weight restriction operator whose stencil takes the
form:
2h
 s 1 

 1 i   2
i  2 
  (2.11)
2h

1 21  1  2  1 
1
 
21 
1 
I h
 1    
2i  i   2


2i 
8  i      
 i    
1 2 1

 i   i  2
2
 h

2.4 Test problems, Results and Performance

In this section we report the results and performance of some test problems All our implementations in
Fortran were executed on a pentium 4 PC using FORTRAN90 workstation compiler. The accuracy of the
method is measured by both L2 norm of both defect and error and the max norm of both defect and the
error. Consider the following test problems: Consider the following test problems:

1. Test Problem 1:
0.5r2 sinθcosθ 2
Uexact = 3e + 2r sin θ cos θ (2.12)

Uexact = e r (sin  cos )

4 2
2. Test Problem 2: (2.13)

2 (2.14)
3. Test Problem 3: Uexact = sin (r sinθ cos θ + 1)

4. Test Problem 4:
6 3 4 2 2
Uexact = r (sin θ cos θ) + 3r (sin θ cos θ) +6r (sin θ cos θ) + 2 (2.15)

5. Test Problem 5:
2 sinθcosθ
Uexact = e r ( ) (2.16)

8
 Test Problem 1: Equation 2.12

Test problem h MaxN (R) L2(R) MaxN (E) L2(E)

1
2 0.487803D − 06 0.487803D − 06 0.823281D − 01 0.823281D − 01
1
4 0.655207D − 04 0.497669D − 04 0.130835D − 01 0.681314D − 02
1
8 0.433554D − 03 0.120637D − 03 0.318050D − 02 0.142206D − 02
Equation 2.12 1
16 0.617784D − 03 0.143102D − 03 0.981808D − 03 0.434354D − 03
1
32 0.390044D − 03 0.622140D − 04 0.471115D − 03 0.168117D − 03
1
64 0.846599D − 03 0.359533D − 04 0.272989D − 03 0.962188D − 04

Table 4: Fortran implementation of test problem 1 using F M W (2, 1) on PC

Test Problem 2: Equation 2.13

Test problem h MaxN (R) L2(R) MaxN (E) L2(E)

1
2 0.340116D − 06 0.340116D − 06 0.355722D − 01 0.3557221D − 01
1
4 0.321983D − 04 0.245546D − 04 0.938702D − 02 0.435724D − 02
1
8 0.346041D − 03 0.985342D − 04 0.240290D − 02 0.103357D − 02
Equation 2.13 1
16 0.226426D − 03 0.430177D − 04 0.701070D − 03 0.274752D − 03
1
32 0.973725D − 04 0.173268D − 04 0.262976D − 03 0.867464D − 04
1
64 0.164582D − 03 0.962157D − 05 0.821352D − 04 0.199558D − 04

Table 5: Fortran implementation of test problem 2 using F M W (2, 1) on PC

Test Problem 3: Equation 2.14

Test problem h MaxN (R) L2(R) MaxN (E) L2(E)

1
2 0.150675D − 06 0.150675D − 06 0.673306D − 02 0.673306D − 02
1
4 0.818226D − 05 0.621353D − 05 0.143027D − 02 0.782289D − 03
1
8 0.118259D − 03 0.332533D − 04 0.363469D − 03 0.194491D − 03
Equation 2.14 1
16 0.124193D − 03 0.245218D − 04 0.150561D − 03 0.490143D − 04
1
32 0.583000D − 04 0.103898D − 04 0.891089D − 04 0.232035D − 04
1
64 0.139612D − 03 0.786312D − 05 0.367165D − 04 0.106076D − 04

Table 6: Fortran implementation of test problem 3 using F M W (2, 1) on PC

Test Problem 4: Equation 2.15

Test problem h MaxN (R) L2(R) MaxN (E) L2(E)

1
2 0.814372D − 07 0.814372D − 07 0.163969D + 00 0163969D + 00
1
4 0.145828D − 03 0.111533D − 03 0.298817D − 01 0.153387D − 01
1
8 0.926281D − 03 0.304502D − 03 0.751638D − 02 0.344485D − 02
Equation 2.15 1
16 0.436931D − 03 0.958402D − 04 0.206471D − 02 0.939828D − 03
1
32 0.184554D − 03 0.331101D − 04 0.545025D − 03 0.260191D − 03
1
64 0.535118D − 03 0.238269D − 04 0.237465D − 03 0.761428D − 04

Table 7: Fortran implementation of test problem 5 using F M W (2, 1) on PC

9
 Test Problem 5: Equation 2.16

Test problem h MaxN (R) L2(R) MaxN (E) L2(E)

1
2 0.413364D − 06 0.413364D − 06 0.282279D − 01 0.282279D − 01
1
4 0.258139D − 04 0.190722D − 04 0.539088D − 02 0.272324D − 02
1
8 0.163754D − 03 0.533917D − 04 0.134575D − 02 0.611430D − 03
Equation 2.16 1
16 0.723312D − 04 0.163671D − 04 0.369430D − 03 0.166279D − 03
1
32 0.587217D − 04 0.633453D − 05 0.935793D − 04 0.460934D − 04
1
64 0.216094D − 03 0.107476D − 04 0.404120D − 04 0.124871D − 04
Table 8: Fortran implementation of test problem 5 using F M W (2, 1) on PC

3 Conclusion
The implementation of an elliptic partial differential equation in the 2-dimensional space in a
quarter of a unit circle have been shown to be an efficient method. The convergence rate of the
method in Cartesian coordinates has been improved as we increase the number of relaxation
steps before coarse grid correction, in the polar coordinates we should use the W-cycle instead
of V-cycle to achieve almost the same or a little bit better convergence rate than that of the
method in Cartesian coordinates. The other advantage of the two algorithms in Cartesian and
in Polar coordinates is that the extension from two to three space dimensions is straightforward.

References
[1] Brandt A., “Multi-level adaptive technique (MLAT) for fast numerical solution to bound-
ary value problem,” Springer, Berlin - 1973.
[2] Briggs L., “A Multigrid Tutorial,” Society for Industrial and applied mathematics - 2000.
[3] Hackbucsh W., “ Multigrid Methods and applications,” Springer Berlin Heidelberg Berlin,-
2003.
[4] G. D. Smith, “ Numerical Solution of Partial Differential Equations: Finite Difference
Methods.”, Oxford Applied Mathematics & Computing Science Series, - January 16, 1986.
[5] McCormick S. F., “Mulitigrid Methods ,” SIAM, Phildelphia-Pennsylvania - 1987
[6] Briggs W. L, “Mulitigrid tutorial ,” SIAM, Phildelphia-Pennsylvania - 1987
[7] John D. Anderson, Jr. “Computational Fluid Dynamics: The Basics with Applications
”, McGraw-Hill Inc, - 1995.
[8] J.E. Dennis and R.B. Schnabel, “Numerical Methods for Unconstrained Optimization
and Nonlinear Equations”, SIAM Books, Philadelphia, - 1996.
[9] J.C. Tannehill, D.A. Anderson, “Computational Fluid Mechanics and Heat Transfer, Sec-
ond Edition ”, Series in Computational and Physical Processes in Mechanics and Thermal
Sciences. , - April 1, 1997.
[10] W. L. Briggs, V. E. Henson, and S. F. McCormick, “A Multigrid Tutorial, Second Edition”
SIAM Books, Philadelphia - 2000
[11] Bobby Philip, “An Introduction to Multigrid Techniques” The Institute for Mathematics
and Its Applications, Minnesota - January 11, 2008.

10