The Multigrid Method is based on two elegant ideas:
1. Nested Iteration
2. Coarse Grid Correction
Both of these ideas are clearly dealt with in the SIAM publication, "A Multigrid Tutorial" by Bill Briggs.
To describe these ideas, let Gh denote a uniform finite difference grid with spacing h on the interval [0,1]. ( if the grid points are labeled x0, x1, x2, ... xN+1, h = 1/(N+1) ) We will assume N+1 is a power of 2.
Gh will be referred to as the fine grid; G2h a coarser grid with grid points separated by 2h; G4h a still coarser grid with grid points separated by 4h; ... ; G1/2 the coarsest nontrivial grid on [0,1].
Let Au = f denote a system of linear equations to be solved by an iterative method, such as the Gauss-Seidel method.
Nested Iteration:
Iterate on Au = f on G1/2 to obtain initial guess for G1/4 Iterate on Au = f on G1/4 to obtain intial guess for G1/8 . . . . . .
Iterate on Au = f on G4h to obtain intial guess for G2h Iterate on Au = f on G2h to obtain initial guess for Gh Iterate on Au = f on Gh
We illustrate the Multigrid Method strategy of Nested Iteration with a simple one-dimesnsional example obtained from Netlib (as an exercise, you should use Netlib to locate the code discussed below):
c Numerical Analysis:
c The Mathematics of Scientific Computing
c D.R. Kincaid & E.W. Cheney
c Brooks/Cole Publ., 1990
c
c Section 9.8
c
c multigrid program used on the equation
c u'' = f, in one variable.
c
c parameter m - the number of mesh
c refinements used in the method
c n - 2**m (should be n+1=2**m)
c k - number of iterations in
c the Gauss-Seidel method
c to be used on each grid.
parameter (m=6, k=3, np1=65)
dimension v(0:np1),w(0:np1)
c the function f is the right-hand side
c of the differential equation.
c the function g is the solution function.
c
real f,g
f(x) = cos(x)
g(x) = -cos(x) + x*(cos(1.0) - 1.0) + * 1.0
c first we define the solution of the
c difference equations on the
c coarsest grid. this grid has three
c points, two of which are on the
c boundary.
c
h = 0.5
n = 1
v(0) = 0.0
v(2) = 0.0
v(1) = -f(0.5)/8.0
print *,' n =',n,' h =',h
print *
c
c now the main iteration (outer loop)
c begins.
c
do 8 i=1,m-1
h = 0.5*h
n = 2*n + 1
print *,' n =',n,' h =',h
c the next part is the interpolation
c phase in which an approximate
c solution on a finer grid is manufactured
c from an approximation on a coarse grid.
do 2 j=0,n-1
w(2*j) = v(j)
2 continue
c
do 3 j=1,n-1
w(2*j-1) = 0.5*(v(j-1) + v(j))
3 continue
c
do 4 j=0,2*n-2
v(j) = w(j)
4 continue
c
c next is the Gauss-Seidel iteration on
c the finer grid (k steps)
c
do 6 p=1,k
do 5 j=1,n
v(j) = 0.5*(v(j-1) + v(j+1) -
* h*h*f(real(j)*h))
5 continue
6 continue
c
do 7 j=0,n+1
w(j) = abs(g(real(j)*h) - v(j))
7 continue
c
print *,' error norm =',vnorm(n+1,w)
8 continue
c
stop
end
function vnorm(n,v)
c
c infinity norm of vector v
c
dimension v(0:n)
real max
c
temp = 0.0
do 2 i=0,n
temp = max(temp,abs(v(i)))
2 continue
vnorm = temp
c
return
end
Multigrid method example
Section 9.8, Kincaid-Cheney
n = 1 h = 0.500000
n = 3 h = 0.250000
error norm = 2.20695E-03
n = 7 h = 0.125000
error norm = 3.75146E-03
n = 15 h = 6.25000E-02
error norm = 4.16399E-03
n = 31 h = 3.12500E-02
error norm = 4.24080E-03
n = 63 h = 1.56250E-02
error norm = 4.27556E-03
Now, let's compare these multigrid results with ordinary Gauss-Seidel.
c gs program used on the equation
c u'' = f, in one variable.
c
c parameter n - the number of interior
c grid points
c k - number of iterations in
c the Gauss-Seidel method to be used
c
c file: gs.f
c
parameter (k=15,n=63, np1=64)
dimension v(0:np1),w(0:np1)
c
c the function f is the right-hand side of the differential equation.
c the function g is the solution function.
c
real f,g
f(x) = cos(x)
g(x) = -cos(x) + x*(cos(1.0) - 1.0) +
* 1.0
c
print *,'GS method example'
print *
c
c first we initialize the solution of
c the difference equations on the grid
c
h = 1.0/float(np1)
do 1 i = 0,np1
v(i) = .5*(g(0)+g(1.0))
1 continue
c
c the Gauss-Seidel iteration on the fine grid (k steps)
c
do 6 p=1,k
do 5 j=1,n
v(j) = 0.5*(v(j-1) + v(j+1) -
* h*h*f(real(j)*h))
5 continue
do 7 j=0,n+1
w(j) = abs(g(real(j)*h) - v(j))
7 continue
print *,'iter',p, ' error
* =',vnorm(n+1,w)
print *
6 continue
end
GS method example
iter 1.00000 error = 0.107414
iter 2.00000 error = 0.107197
iter 3.00000 error = 0.106979
iter 4.00000 error = 0.106761
iter 5.00000 error = 0.106544
iter 6.00000 error = 0.106326
iter 7.00000 error = 0.106109
iter 8.00000 error = 0.105892
iter 9.00000 error = 0.105674
iter 10.00000 error = 0.105457
iter 11.0000 error = 0.105240
iter 12.0000 error = 0.105023
iter 13.0000 error = 1.04806E-01
iter 14.0000 error = 1.04589E-01
iter 15.0000 error = 1.04372E-01
The above examples illustrate that the technique of Nested Iteration can lead to dramatically impvoved results when compared to standard Gauss-Seidel on the fine grid. Not all calculations will be as favorable as this example. What should we do if the final solution produced by Nested Iterataion still has significant errors?
The second idea on which multigrid is based is one we have used before -- we called it iterative refinement. It is the basis for the Coarse Grid Correction strategy.
To review:
Let u denote the exact solution to Au = f
Let v denote an approximate solution to Au = f
Let e = u - v be the error incurred in using v to approx u
Let r be the residual associated with v; r = f - Av
It follows that Ae = r
If we approximately solve for e, then an improved estimate for u is obtained from v + e
With this brief review, we can state the basis for the Coarse Grid Correction Strategy as follows:
Coarse Grid Correction
Iterate on Gh to obtain vh as an approximation to Au = f
Compute the residual r = f - Avh Iterate on G2h to obtain an approximate solution toAe = r
Correct the approximation obtained on Gh by replacing vh with vh + e
A recursive application of Coarse Grid Correction is the basis for the Multigrid V-Cycle:
c
c Numerical Analysis:
c The Mathematics of Scientific Computing
c D.R. Kincaid & E.W. Cheney
c Brooks/Cole Publ., 1990
c
c Section 9.8
c
c complete V-cycle in the multigrid algorithm to solve
c u" = g(x), with u(0)=u(l)=0.
c parameter m gives the number of grids to be employed.
c parameter k gives the number of iterations to be performed
c in each gauss-seidel iteration.
c parameter np1 is simply n+1
c
c
c file: mgrid2.f
c
parameter (m=7, np1=128)
dimension v(m,0:np1), f(m,0:np1), w(0:np1), error(15)
integer p
real g,true
g(x) = cos(x)
true(x) = -g(x) + x*(g(1.0) - 1.0) + 1.0
c
print *
print *,' V-cycle in the multigrid method'
print *,' Section 9.8, Kincaid-Cheney'
print *
c
error(1) = 1.0
do 12 k=2,15
c
c initialize arrays and
c put best available guess into v**m
c
n = 2**m - 1
h = 1.0/real(n+1)
print *,' n =',n,' h =',h
c
do 1 i=1,m
do 1 j=0,np1
f(i,j) = 0.0
v(i,j) = 0.0
1 continue
c
c put data from the bvp into f**m
c
do 2 j=1,n
f(m,j) = g(real(j)*h)
2 continue
c
do 6 i=m,2,-1
c
c gauss-seidel iteration (k steps)
c
do 3 p=1,k
do 3 j=1,n
v(i,j) = 0.5*(v(i,j-1) + v(i,j+1) - (h**2)*f(i,j))
3 continue
c
c residual vector
c
do 4 j=1,n
w(j) = f(i,j) - (v(i,j-1) - 2.0*v(i,j) + v(i,j+1))/h**2
4 continue
c
c apply restriction operator
c
do 5 j=1,(n-1)/2
f(i-1,j) = w(2*j)
5 continue
c
h = 2.0*h
n = (n-1)/2
print *,' n =',n,' h =',h
6 continue
c
c solve coarsest (smallest) system exactly
c
v(1,1) = -g(0.5)/8.0
print *,' Bottom of V-cycle'
c
c end downward phase of V-cycle
c start upward phase of V-cycle
c
do 10 i=2,m
h = 0.5*h
n = 2*n + 1
print *,' n =',n,' h =',h
c
c apply extension operation, i.e.
c interpolation up from coarse grid to fine grid
c
do 7 j=0,(n+1)/2
w(2*j) = v(i-1,j)
7 continue
c
do 8 j=1,(n+1)/2
w(2*j-1) = 0.5*(v(i-1,j-1) + v(i-1,j))
8 continue
c
c add correction to extension operator
c
do 9 j=0,n+1
v(i,j) = v(i,j) + w(j)
9 continue
c
c gauss-seidel iteration (k steps)
c
do 10 p=1,k
do 10 j=1,n
v(i,j) = 0.5*(v(i,j-1) + v(i,j+1) - (h**2)*f(i,j))
10 continue
c
c output phase
c
print *,' k =',k,' m =',m
do 11 j=0,np1
w(j) = true(real(j)*h) - v(m,j)
11 continue
error(k) = vnorm(np1,w)
print *,' maximun error =', error(k),
+ ' ratio =', error(k)/error(k-1)
c
12 continue
stop
end