Äîêóìåíò âçÿò èç êýøà ïîèñêîâîé ìàøèíû. Àäðåñ îðèãèíàëüíîãî äîêóìåíòà : http://num-meth.srcc.msu.ru/english/zhurnal/tom_2001/ps/art1_1.ps Äàòà èçìåíåíèÿ: Mon Dec 16 17:36:04 2002 Äàòà èíäåêñèðîâàíèÿ: Mon Oct 1 22:59:25 2012 Êîäèðîâêà:

Numerical Methods and Programming, 2001, Vol. 2 1
UDC 519.632.4
A ROBUST MILTIGRID TECHNIQUE FOR SOLVING SYSTEMS OF PARTIAL
DIFFERENTIAL EQUATIONS
S. I. Martynenko 1
A robust multigrid technique is used for solving systems of partial differential equations. The tech­
nique unites the processes of adaption of equations to numerical methods, their discretization by
the control volume method, and application of multigrid iterations. The technology is efficient for
solving linear and nonlinear systems on unstaggered and staggered grids. This paper shows that our
multigrid decoupled solver is competitive with the coupled one. New modifications that reduce the
problems into those acceptable for the multigrid technique are developed.
1. Introduction. Recently, a robust multigrid technique (R M T) based on adaption of problems to
numerical methods, control volume discretization and multigrid iterations has been proposed [1]. This paper is
concerned with the development of R M T for handling systems of PDEs. At first, we study four solvers for linear
systems to choose the most suitable of them. Further, a robust solver for CFD­like systems on staggered grids
will be considered in detail.
The paper assumes a basic knowledge of R M T (multigrid structure, vertex­and­cell­centered coarsening,
one­to­one mapping of indices fg, problem­independent transfer operators, sawtooth cycles, equivalent number
of iterations “
š, \Sigma­modification, separate discretizations, approximation boundary conditions on coarse grids,
etc.) [1].
2. Basic multigrid algorithms for systems. Let us consider the application of R M T for solving the
model system of linear elliptic equations
8
? ? ? ? ? ? ? !
? ? ? ? ? ? ? :
@ 2 U
@x 2
+ @ 2 U
@y 2
+ fl V (x; y) + fl Z(x; y) = F (x; y)
@ 2 V
@x 2 + @ 2 V
@y 2 + fl U (x; y) + fl Z(x; y) = F (x; y)
@ 2 Z
@x 2 + @ 2 Z
@y 2 + fl U (x; y) + fl V (x; y) = F (x; y)
(1)
U
fi
fi\Omega = 0 ; V
fi
fi\Omega = 0 ; Z
fi
fi\Omega = 0
in\Omega = (0; 1) \Theta (0; 1). This system is considered only for simplicity of theoretical analysis. The \Sigma­modification
of the solution
U (x; y) = C u (x; y) + “
U (x; y)
V (x; y) = C v (x; y) + “
V (x; y)
Z(x; y) = C z (x; y) + “
Z(x; y)
makes it possible to rewrite the model system in the \Sigma­modified form
8
? ? ? ? ? ? ? !
? ? ? ? ? ? ? :
@ 2 Cu
@x 2 + @ 2 Cu
@y 2 + fl C v (x; y) + fl C z (x; y) = r u (x; y)
@ 2 C v
@x 2
+ @ 2 C v
@y 2
+ fl Cu (x; y) + fl C z (x; y) = r v (x; y)
@ 2 C z
@x 2 + @ 2 C z
@y 2 + fl Cu (x; y) + fl C v (x; y) = r z (x; y)
C u
fi
fi\Omega = \Gamma “
U
fi
fi\Omega ; C v
fi
fi\Omega = \Gamma “
V
fi
fi\Omega ; C z
fi
fi\Omega = \Gamma “
Z
fi
fi\Omega
1 Central Institute of Aviation Motors, Aviamotornaya 2, 111250 Moscow, Russian Federation,
e­mail: martyn s@mail.ru
c
fl Research Computing Center, Moscow State University, 119899 Moscow, Russian Federation

2 Numerical Methods and Programming, 2001, Vol. 2
where
r u (x; y) = F (x; y) \Gamma @ 2 “
U
@x 2 \Gamma @ 2 “
U
@y 2 \Gamma fl “
V (x; y) \Gamma fl “
Z(x; y)
r v (x; y) = F (x; y) \Gamma @ 2 “
V
@x 2 \Gamma @ 2 “
V
@y 2 \Gamma fl “
U (x; y) \Gamma fl “
Z(x; y)
r z (x; y) = F (x; y) \Gamma @ 2 “
Z
@x 2
\Gamma @ 2 “
Z
@y 2
\Gamma fl “
U (x; y) \Gamma fl “
V (x; y)
The functions “
U , “
V , and “
Z will be the approximations to the solutions U , V , and Z, whereas the functions C u ,
C v , and C z will be the corrections used in the subsequent multigrid iterations.
The second step of R M T is the generation of a multigrid structure [1]. The finest computational grids for
the control volume discretization consists of two sets of grid points G v = G v
x [ G v
y and G f = G f
x [ G f
y , where
G v
x = fx v
i : x v
i = \Delta(i \Gamma 1); i = 1; : : : ; H 0 + 1; \Delta = H 0
\Gamma1 g
G f
x = fx f
i : x f
i = 0:5 (x v
i + x v
i+1 ); i = 1; : : : ; H 0 g
G v
y = fy v
j : y v
j = \Delta(j \Gamma 1); j = 1; : : : ; H 0 + 1; \Delta = H 0
\Gamma1 g
G f
y = fy f
j : y f
j = 0:5 (y v
j + y v
j+1 ); j = 1; : : : ; H 0 g
Assume that the functions U , V , and Z are assigned to the grid points (x v ; y v ). Integration of the \Sigma­modified
system over the control
volume\Omega fijg
\Omega fijg =
n
(x;y) j x f
fi\Gamma1g 6 x 6 x f
fig ; y f
fj \Gamma1g 6 y 6 y f
fjg
o
yields the following finite­difference scheme:
L fijg (Cu ) + fl C vfijg + fl C zfijg = Jufijg
L fijg (C v ) + fl C ufijg + fl C zfijg = J vfijg
L fijg (C z ) + fl Cufijg + fl C vfijg = J zfijg
Here L fijg is the standard five­point finite­difference representation of the Laplacian:
L fijg (OE) =
OE fi\Gamma1jg \Gamma 2OE fijg + OE fi+1jg
\Delta 2 3 2L +
OE fij \Gamma1g \Gamma 2OE fijg + OE fij+1g
\Delta 2 3 2L
and the integrals J u , J v and Ju are defined as
J ffl fijg = 1
\Delta 2 3 2L
x f
fig
Z
x f
fi\Gamma1g
y f
fjg
Z
y f
fj\Gamma1g
r ffl (x; y) dy dx; ffl =
0
B B @
u
v
z
1
C C A
The integrals J ffl can be computed on the finest grid (L = 0) as follows:
(R ?
u ) mk = Fmk \Gamma Lmk ( “
U ) \Gamma fl “
Vmk \Gamma fl “
Zmk + o(\Delta 2 )
(R ?
v ) mk = Fmk \Gamma Lmk ( “
V ) \Gamma fl “
Umk \Gamma fl “
Zmk + o(\Delta 2 )
(R ?
z ) mk = Fmk \Gamma Lmk ( “
Z) \Gamma fl “
Umk \Gamma fl “
Vmk + o(\Delta 2 )
The finite­difference model system may be written as
8
? !
? :
ACu + fl C v + fl C z = Ju
AC v + fl Cu + fl C z = J v
AC z + fl Cu + fl C v = J z
(2)
where the matrix A corresponds to the five­point Laplacian (L fijg ). Consider the following algorithms for
solving the discretized model system.

Numerical Methods and Programming, 2001, Vol. 2 3
Test 1: Coupled Algorithm. The first step in the coupled algorithm is to combine all the finite­difference
equations (2) into the system 0
@
A flE flE
flE A flE
flE flE A
1
A
0
@
C u
C v
C z
1
A =
0
@
J u
J v
J z
1
A
where E is the unity matrix. Let the coefficient matrix be split as
0
@
A flE flE
flE A flE
flE flE A
1
A =
0
@
AL +AD 0 0
flE AL +AD 0
flE flE AL +AD
1
A +
0
@
AU flE flE
0 AU flE
0 0 AU
1
A
where AD , AL , and AU represent the diagonal, lower­triangular, and upper­triangular parts of the matrix A,
respectively. Iterations of the Gauss--Seidel method in the coupled algorithm are defined by
0
@
AL +AD 0 0
flE AL +AD 0
flE flE AL + AD
1
A
0
@
C u
C v
C z
1
A
n+1
=
0
@
Ju
J v
J z
1
A \Gamma
0
@
AU flE flE
0 AU flE
0 0 AU
1
A
0
@
C u
C v
C z
1
A
n
In the decoupled algorithms, each unknown is determined with the others fixed at their previous iterative values.
Here we study three decoupled solvers.
Test 2: Classical Decoupled Algorithm. The following program performs iterations in the decoupled manner
and checks the convergence:
C
C Classical Decoupled Algorithm
C
k = 0
Starting guesses: C (0)
u , C (0)
v and C (0)
z
1 k = k + 1
AC (k)
u = J (k)
u \Gamma flC (k\Gamma1)
v \Gamma flC (k\Gamma1)
z ! C (k)
u
AC (k)
v = J (k)
v \Gamma flC (k)
u \Gamma flC (k\Gamma1)
z ! C (k)
v
AC (k)
z = J (k)
z \Gamma flC (k)
u \Gamma flC (k)
v ! C (k)
z
C
C Check convergence; continue (Go To 1) if necessary
C
End
Each equation is solved by R M T until
kR ?
ffl k
kFk ! 10 \Gamma7 ; ffl =
0
B B
@
u
v
z
1
C C
A
where k \Delta k is the l 2 ­norm, R ?
u , R ?
v and R ?
z are the residuals on the finest grid of the first, second and third
equation, respectively, and kFk is the l 2 ­norm of the right­hand side vector.
Test 3: Simplified Multigrid Decoupled Algorithm. This algorithm performs the smoothing iterations on the
multigrid structure as shown in Figure 1.
Test 4: Complicated Multigrid Decoupled Algorithm. This algorithm performs the decoupled iterations on
the coarsest level (L + ) and the smoothing iterations on the finer levels (Figure 2).
Several numerical experiments are performed under the following conditions:
1. the function F is chosen so that the exact solution of the model system be
U e (x; y) = V e (x; y) = Z e (x; y) = f(x)f(y)
where
f(%) = 10
\Gamma
e % + (1 \Gamma e)% \Gamma 1
\Delta

4 Numerical Methods and Programming, 2001, Vol. 2
Fig. 1. Simplified Multigrid Decoupled Algorithm
Fig. 2. Complicated Multigrid Decoupled Algorithm
2. the starting guess is taken to be zero: U (0) = V (0) = Z (0) = 0;
3. the uniform finest 401 \Theta 401 grid is chosen, multigrid structure has five levels: L + = 4;
4. the error of the numerical solution (E) is defined as
E = max(EU ; EV ; EZ ); where
EU = max
ij
j U e (x v
i ; y v
j ) \Gamma U ij j
EV = max
ij
j V e (x v
i ; y v
j ) \Gamma V ij j
EZ = max
ij
j Z e (x v
i ; y v
j ) \Gamma Z ij j
5. the three smoothing iterations of Alternating Line Gauss­Seidel are performed on the finer levels: š = 3 ; 0 6
L ! L + ;
6. the global stopping criterion is defined by max (kR ?
u k ; kR ?
v k ; kR ?
z k) ! 10 \Gamma7 kFk .
In order to compare the rates of convergence of the above algorithms, we define the efficiency of a (de)coupled
solver as
\Xi = “
š \Sigma (0)
“
š \Sigma (fl)
where “ š \Sigma is the equivalent number of smoothing iterations [1]. When the model system is solved as fast as the
three independent equations (fl = 0), the efficiency of the solver tends to unity: \Xi ! 1 as “
š \Sigma (fl) ! “
š \Sigma (0). A
slow convergence of (de)coupled iterations corresponds to \Xi ! 0 as “
š \Sigma (fl) AE “
š \Sigma (0). The convergence rates in
these tests is shown in Figure 3 and Table 1.
Let us estimate the rate of convergence of decoupled iterations. The vector U (1) obtained in the first
iteration satisfies the equation
AU (1) = F \Gamma fl V (0) \Gamma fl Z (0) (3)
Current errors of the numerical solutions can be defined by
ffi (k)
U = U \Gamma U (k) ; ffi (k)
V = V \Gamma V (k) ; ffi (k)
Z = Z \Gamma Z (k)

Numerical Methods and Programming, 2001, Vol. 2 5
Fig. 3. Convergence behavior in Tests 1; 2; 3; 4
Equation (3) can be rewritten as
Affi (1)
U = \Gammafl ffi (0)
V \Gamma fl ffi (0)
Z
The error ffi (1)
U can be estimated as
ffi (1)
U = \Gammafl A \Gamma1 (ffi (0)
V + ffi (0)
Z ) ) kffi (1)
U k 6 fl kA \Gamma1 k kffi (0)
V + ffi (0)
Z k (4)
A similar estimate for the solution U is given by
kAk kUk ? kF \Gamma fl V \Gamma fl Zk (5)
Estimates (4) -- (5) can be united by
kffi (1)
U k
kUk 6 fl kA \Gamma1 k kAk kffi (0)
V + ffi (0)
Z k
kF \Gamma fl V \Gamma fl Zk
(6)
The condition number has the following asymptotic behavior [2]:
kA \Gamma1 k kAk !
`
2
ú\Delta
' 2
as \Delta ! 0
Hence, estimate (6) can be rewritten as
kU \Gamma U (1) k
kUk 6
`
2
ú\Delta
' 2 kffi (0)
V + ffi (0)
Z k
kF \Gamma fl V \Gamma fl Zk fl
Since
kffi (0)
V + ffi (0)
Z k 6 kffi (0)
V k + kffi (0)
Z k 6 (K \Gamma 1)ffi max ; ffi max = max(kffi (0)
V k; kffi (0)
Z k)
estimate (6) can be written down in its final form as
kU \Gamma U (1) k
kUk 6
`
2
ú\Delta
' 2 (K \Gamma 1)fl ffi max
kF \Gamma fl V \Gamma fl Zk
(7)
where K is the number of equations.
The number of the decoupled iterations depends strongly on the starting guess ffi max , the number K of
differential equations in the system, the grid spacing \Delta and the parameter fl. The strong coupling of equations
in the system (fl '') results in the deterioration of convergence.

6 Numerical Methods and Programming, 2001, Vol. 2
Table 1. Convergence rate of the multigrid (de)coupled algorithms in Tests 1; 2; 3; 4
fl 0 0.5 1 3 5 7
Test 1 “
š \Sigma 272 280 288 406 516 949
\Xi 1 0.971 0.944 0.670 0.558 0.287
Test 2 “
š \Sigma 272 660 770 1254 1958 3649
\Xi 1 0.412 0.353 0.217 0.139 0.075
Test 3 “
š \Sigma 272 340 343 482 827 1448
\Xi 1 0.800 0.793 0.564 0.329 0.188
Test 4 “
š \Sigma 272 284 293 417 539 1006
\Xi 1 0.958 0.928 0.652 0.505 0.270
Figure 3 and Table 1 report a slow convergence of the classical decoupled algorithm (Test 2). The first
iteration reduces the current error of the solution U as follows:
kffi (q)
U k 6 ff q kffi (0)
U k ; ff 2 (0;1) ) q 6
1
ln ff \Gamma1
ln kffi (0)
U k
flkA \Gamma1 (ffi (0)
V + ffi (0)
Z )k
Let us denote
q \Lambda =
''
1
ln ff \Gamma1
ln kffi (0)
U k
flkA \Gamma1 (ffi (0)
V + ffi (0)
Z )k
#
+ 1
where the square brackets indicate the integer part. When q 6 q \Lambda , each multigrid iteration of the classical
decoupled algorithm improves the solution, but when q ? q \Lambda , the multigrid iteration deminishes only a residual
of the discretization equation without any improvement of the solution (Figure 3). An appropriate stopping
criterion for the classical decoupled algorithm should depend on the errors ffi max of the right­hand side vector.
Proximity of solutions obtained on the neighboring level leads to a convergence acceleration of decoupled
iterations: ffi max # ) \Xi '' (Test 3).
Further development of the decoupled algorithms is to minimize ffi max (this will accelerate the convergence
of the iterations to the solution (Test 4)). Consider the two­level structure and assume that the coarsest level
solution has been obtained. Accuracy of the solutions can be estimated as
M e (x v
i ;y v
j ) \Gamma M L=0
ij = \Delta 2 S M
ij + o(\Delta 4 ) ;
M e (x v
i ;y v
j ) \Gamma M L=1
ij = 9\Delta 2 S M
ij + o(\Delta 4 ) ;
M =
0
@
U
V
Z
1
A
therefore
M L=0
ij \Gamma M L=1
ij = ffi (0)
M = 8\Delta 2 S M
ij + o(\Delta 4 ) ) kffi (0)
M k 6 8\Delta 2 S M
max
where S M is a bounded function: max ij jS M
ij j 6 Smax . The norm of errors in the right­hand side vector can be
estimated as follows:
kffi (0)
V + ffi (0)
Z k 6 kffi (0)
V k + kffi (0)
Z k 6 8\Delta 2 (K \Gamma 1)Smax
Estimate (7) is expressed as
kU \Gamma U (1) k
kUk 6
32
ú 2
(K \Gamma 1)flS max
kF \Gamma flV \Gamma flZk
Convergence of the algorithm is independent of the mesh size \Delta. The proximity of the solutions obtained on
the neighboring levels and the absence of prolongation errors result in an impressive acceleration of convergence
of the decoupled solver in Test 4.
It should be noted that the coefficient matrix in Tests 1, 2, 3, 4 is not diagonally dominant; therefore,
convergence of the (de)coupled algorithms strongly decreases with increasing fl (Table 1). Consider the second

Numerical Methods and Programming, 2001, Vol. 2 7
Table 2. Convergence rate of the multigrid (de)coupled algorithms in Tests 1 \Lambda , 2 \Lambda , 3 \Lambda , 4 \Lambda
fl 0 0.5 1 3 5 7 10
Test 1 \Lambda “
š \Sigma 272 276 276 288 295 300 312
\Xi 1 0.986 0.986 0.944 0.922 0.907 0.872
Test 2 \Lambda “
š \Sigma 272 638 770 1029 1260 1413 1674
\Xi 1 0.426 0.353 0.264 0.216 0.192 0.162
Test 3 \Lambda “
š \Sigma 272 335 332 378 480 522 560
\Xi 1 0.806 0.819 0.720 0.567 0.521 0.486
Test 4 \Lambda “
š \Sigma 272 276 283 295 308 322 340
\Xi 1 0.986 0.961 0.922 0.883 0.845 0.800
model system with the same exact solution U e (x; y) = V e (x; y) = Z e (x; y) = f(x)f(y):
8
? ? ? ? ? ? ? !
? ? ? ? ? ? ? :
@ 2 U
@x 2 + @ 2 U
@y 2 \Gamma 2flU (x; y) + fl V (x; y) + fl Z(x; y) = F (x; y)
@ 2 V
@x 2
+ @ 2 V
@y 2 \Gamma 2flV (x; y) + fl U (x; y) + fl Z(x; y) = F (x; y)
@ 2 Z
@x 2
+ @ 2 Z
@y 2 \Gamma 2flZ(x; y) + fl U (x; y) + fl V (x; y) = F (x; y)
(8)
The coefficient matrix of the second model system is diagonally dominant for any fl. The computational
procedures in Test 1 \Lambda , 2 \Lambda , 3 \Lambda , 4 \Lambda are similar to those in Tests 1, 2, 3, 4. Table 2 represents the convergence rate
for the second model system.
Our numerical tests and theoretical analysis show that the complicated multigrid decoupled algorithm
(Tests 4, 4 \Lambda ) for solving systems of partial differential equations is competitive in comparison with the coupled
one (Tests 1, 1 \Lambda ). However, the decoupled algorithm does not require a global linearization of nonlinear dis­
cretization equations. In addition, this algorithm is easy for programming. Below, this solver (Tests 4, 4 \Lambda ) will
be named Multigrid Decoupled Algorithm (or, briefly, MDA).
3. A model system of nonlinear equations. Multidimensional fluid flow and heat­transfer phenomena
are governed by a system of coupled parial differential equations that express the transport of mass, momenta,
heat, and other scalar variables. Consider the CFD­like model system of nonlinear equations in the domain
\Omega = (0; 1) \Theta (0; 1)
8
? ? !
? ? :
@(u 2 )
@x + @(vu)
@y = '' @ 2 u
@x 2 + '' @ 2 u
@y 2 + F u (x; y) ;
@(uv)
@x
+ @(v 2 )
@y
= '' @ 2 v
@x 2
+ '' @ 2 v
@y 2
+ F v (x; y) ;
u
fi
fi\Omega = 0 ; v
fi
fi\Omega = 0 (9)
The model system can be rewritten in the \Sigma­modified form as
8
? ? !
? ? :
\Gamma @(c 2
u + 2“uc u )
@x \Gamma @((c v + “
v)c u )
@y
+ '' @ 2 c u
@x 2
+ '' @ 2 c u
@y 2
= r u (x; y)
\Gamma @((c u + “
u)c v )
@x
\Gamma @(c 2
v + 2“vc v )
@y
+ '' @ 2 c v
@x 2
+ '' @ 2 c v
@y 2
= r v (x; y)
where
r u (x; y) = @(“u) 2
@x
+ @(v“u)
@y
\Gamma '' @ 2 “
u
@x 2
\Gamma '' @ 2 “
u
@y 2
\Gamma F u (x; y)
r v (x; y) = @(u“v)
@x + @(“v) 2
@y \Gamma ''
@ 2 “
v
@x 2 \Gamma ''
@ 2 “ v
@y 2 \Gamma F v (x; y)
Staggered grids are often used in computational fluid dynamics, as introduced by Harlow and Welch. The
most obvious way of constructing the control volumes is to place their faces midway between neighboring grid
points [3]. In the staggered grid, the velocity components are calculated for the points that lie on the faces of

8 Numerical Methods and Programming, 2001, Vol. 2
Fig. 4. The staggered mesh pattern and definition of control volumes
the control volumes. The functions u and v are assigned to the grid points (x v
i ; y f
j ) and (x f
i ; y v
j ) as shown in
Figure 4 (i.e., u ij 2 G v
x [ G f
y and v ij 2 G f
x [ G v
y ).
Consider a control volume discretization of the first equation in the model system. Since the discretized
system will be solved in the multigrid decoupled manner, the discretization equations can be obtained by
integrating over the control volumes (Figure 4) as follows:
\Gamma 1
\Delta3 L
0
B
@(cu ) 2
fjg + 2
\Delta3 L
y v
fj+1g Z
y v
fjg
“
u dy \Delta (c u ) fjg
1
C
A
fi
fi
fi fi fi fi fi
x f
fig
x f
fi\Gamma1g
(10)
\Gamma 1
\Delta3 L
0
B
B
@
1
\Delta3 L
x f
fig
Z
x f
fi\Gamma1g
(c v + “
v) dx \Delta (c u ) fig
1
C
C
A
fi
fi
fi
fi
fi
fi
fi
fi
y v
fj+1g
y v
fjg
+ '' L fijg (c u ) = (J u ) fijg
Note that the variable c u is considered for c v , “ v, “
u whose values are fixed from the previous iteration. The
following integrals
2
\Delta3 L
y v
fj+1g Z
y v
fjg
“ u dy ; 1
\Delta3 L
x f
fig
Z
x f
fi\Gamma1g
(c v + “
v) dx ; (J u ) fijg = 1
\Delta 2 3 2L
x f
fig
Z
x f
fi\Gamma1g
y v
fj+1g Z
y v
fjg
r u (x; y) dx dy
should be computed on the finest grid to obtain an accurate formulation of the discrete problem on coarse
grids. The false diffusion can be reduced by a separate discretization of the left­hand side operator in (10) and
right­hand side integral J u [1]. The convection terms in (10) should be approximated by the upwind difference
to obtain a monotone finite­difference scheme. However, the right­hand side integral J u should be evaluated on
the finest grid to obtain a more accurate numerical solution as follows:
(R ?
u ) mk = (“u) 2
e \Gamma (“u) 2
w
\Delta + vn “
un \Gamma v s “ u s
\Delta \Gamma '' Lmk (“u) \Gamma (Fu ) mk + o(\Delta 2 )
Here
“ u e = 0:5 (“u mk + “
um+1k ) + o(\Delta 2 ) ; “
uw = 0:5 (“u mk + “ um\Gamma1k ) + o(\Delta 2 )
“
un = 0:5 (“u mk + “
umk+1 ) + o(\Delta 2 ) ; “
u s = 0:5 (“u mk + “ umk\Gamma1 ) + o(\Delta 2 )
v n = 0:5 (v m\Gamma1k+1 + v mk+1 ) + o(\Delta 2 ) ; v s = 0:5 (v m\Gamma1k + vmk ) + o(\Delta 2 )
Unfortunately, there exist examples for which Newton's method for solving nonlinear discretization equations
fails and the convergence takes place only at a good starting guess. Consider the very simple problem
du 2
dx
= ''
d 2 u
dx 2
+ 2ú sin(4úx) + 4ú 2 '' sin(2úx) ; u(0) = u(1) = 0

Numerical Methods and Programming, 2001, Vol. 2 9
Table 3. Multigrid convergence in Test 5
k '' k q \Xi Eu ae u
q Ev ae v
q
0 2 +0 1 0:67 1:999 \Delta 10 +0 0:008 1:998 \Delta 10 +0 0:008
1 2 \Gamma1 1 0:65 1:026 \Delta 10 +0 0:244 1:024 \Delta 10 +0 0:245
2 2 \Gamma2 1 0:65 5:013 \Delta 10 \Gamma1 0:256 4:994 \Delta 10 \Gamma1 0:258
3 2 \Gamma3 1 0:66 2:494 \Delta 10 \Gamma1 0:245 2:484 \Delta 10 \Gamma1 0:247
4 2 \Gamma4 1 0:66 1:240 \Delta 10 \Gamma1 0:246 1:236 \Delta 10 \Gamma1 0:248
5 2 \Gamma5 1 0:66 6:108 \Delta 10 \Gamma2 0:248 6:084 \Delta 10 \Gamma2 0:250
6 2 \Gamma6 1 0:66 2:958 \Delta 10 \Gamma2 0:247 2:946 \Delta 10 \Gamma2 0:249
7 2 \Gamma7 1 0:66 1:384 \Delta 10 \Gamma2 0:247 1:378 \Delta 10 \Gamma2 0:249
8 2 \Gamma8 1 0:66 5:978 \Delta 10 \Gamma3 0:247 5:950 \Delta 10 \Gamma3 0:249
9 2 \Gamma9 1 0:66 2:067 \Delta 10 \Gamma3 0:247 2:054 \Delta 10 \Gamma3 0:249
10 2 \Gamma10 5 0:81 4:644 \Delta 10 \Gamma8 0:093 4:639 \Delta 10 \Gamma8 0:093
Its exact solution is u e (x) = sin(2úx). Starting with the zero guess u (0) = 0, the first approximation to the
solution u (1) obtained by Newton's method takes the form
u (1) = 1
8ú'' sin(4úx) + sin(2úx)
It follows from the above that Newton's iterations diverge as '' ! 0:
max
[0;1]
ju e (x) \Gamma u (0) j = max
[0;1]
j sin(2úx)j = 1 Ü max
[0;1]
ju e (x) \Gamma u (1) j = 1
8ú'' max
[0;1]
j sin(4úx)j = 1
8ú''
Convergence of Newton's iterations depends critically on '': if '' ! +1 (diffusion dominated problem), then
convergence of Newton's iteration is independent of a starting guess; if '' ! 0 (convection dominated problem),
then convergence depends strongly on a starting guess.
From the R M Tpoint of view, the global convergence of Newton's iterations can be obtained by solving the
model problem with a variable small parameter ''. The computational procedure starts at a sufficiently large
value of '' to provide convergence of Newton's iterations. Then, the value of '' decreases gradually down to the
given value. A solution obtained with a larger value of '' is considered to be a starting guess for a solution with
a smaller value of ''. In what follows, variation of '' will be called the ''­modification.
Test 5: Uniform 411 \Theta 411 finest grid (L + = 4). The functions Fu and F v are chosen so that the exact
solution be u e (x; y) = v e (x; y) = 2'' sin(2úx) sin(2úy), where '' = 2 \Gamma10 . The computational procedure for this
test can be represented by the following program:
C
C Globally convergent modification of Newton's method
C
DO k = 0,10
'' k = 2 \Gammak
IF(k.Eq.0), starting guess: u (0)
ij ('' k ) = v (0)
ij ('' k ) = 0
IF(k.Ne.0), starting guess: u (0)
ij ('' k ) = u ij ('' k\Gamma1 ) ; v (0)
ij ('' k ) = v ij ('' k\Gamma1 )
Solve problem (9) ) u ij ('' k ) and v ij ('' k )
End DO
End
The efficiency of MDA in the sixth test is defined by \Xi = “ š \Sigma =(2“š \Lambda
\Sigma
), where “
š \Lambda
\Sigma
is the equivalent number of
smoothing iterations for solving the single equation
@(u 2 )
@x
+ @(v e u)
@y
= '' @ 2 u
@x 2
+ '' @ 2 u
@y 2
+ Fu (x; y)
Table 3 reports results of our numerical experiment, where ae q is the average reduction factor of the residual [4].
It should be emphasized that R M T is efficient both for unstaggered and staggered grids. The main components
of R M T (multigrid structure, vertex­and­cell­centered coarsening, mapping of indices, problem­independent

10 Numerical Methods and Programming, 2001, Vol. 2
transfer operators, sawtooth cycles, etc.) are independent of definition of the control volumes and assignment
of solutions to the grid points. As opposed to this fact, the efficiency of classical multigrid methods depends on
the choice of control volumes. For example, the extension to a system by the black box multigrid method for
staggered grids is not very satisfactory [4].
1. Problem classification. From the R M Tpoint of view, all problems to be solved can be classified
according to the following definition.
Definition. The problems of order k of complexity require k modifications for the adaption to those
acceptable for R M T.
The Poisson equation, anisotropic equations, equations with discontinuous coefficients [1] and model systems
(1), (8) are examples of problems of first­order complexity. All of them require a single \Sigma­modification for
adaption to R M T. The \Pi­modification seems to be more preferable for some nonlinear equations [1].
Singular perturbation problems have the second­order complexity. In addition to the \Sigma­modification, they
require an additional ''­modification for the adaption to R M T. The Taylor series expansion about '' k yields
u(x; y; '' k ) = u(x; y; '' k\Gamma1 ) + 2'' k
@u
@''
fi fi fi fi
'' k
+ o('' 2
k ) ; '' k = '' k\Gamma1
2
The correction c (0)
k (x; y) in the computational procedure used in Test 5 can be estimated as
max jc (0)
k (x; y)j = max ju(x; y; '' k ) \Gamma u(x; y; '' k\Gamma1 )j = O('' k )
Contrary to the function u(x; y; ''), the correction c (0)
k (x; y) has no singularity of the boundary­layer type.
In general, the Navier--Stokes equations have the third­order complexity. First, these equations should be
modified to obtain the strong coupling between velocity components and pressure. Second, the Navier--Stokes
equation should be adapted for R M T. Third, the ''­modification makes it possible to avoid the divergence of
Newton's iterations.
REFERENCES
1. S.I. Martynenko, ``Robust multigrid technique for solving partial differential equations on structured grids'',
Numerical Methods and Programming, Volume 1, Section 1, pp. 83--100, 2000.
2. L.A. Hageman and D.M. Young, Applied Iterative Methods, New York, 1981.
3. S. Patankar, Numerical Heat Transfer and Fluid Flow, New York, 1980.
4. J. E. Dendy, ``Black box multigrid for systems'', Appl. Math. Comput., 19: 57--74, 1986.
16 October 2000