Multigrid

Multigrid is a method for solving boundary-value problems (ODEs and PDEs) by recursive iteration. Consider solving the BVP . Of course, the solution is y=0. Consider solving this equation numerically, on a grid of N points with an initial guess . Stationary iterative methods like Jacobi and Gauss-Siedel damp the high frequency components quickly. However neither method damps the low frequency waves very much.

Now notice that, on a grid of N/2 points, low frequency components look more oscillatory. Jacobi or G-S on this coarse grid removes additional components. You may be asking How can we transfer the y's to a coarse grid, and what is the correct ODE on that grid? To address these questions, make use of linearity. By discritizing, the ODE can be written

Let be the approximate solution on the grid with gridsize h. The error is e=u-v, which we can't compute. But the residual is r=b-Av, which we can compute. Write

Subtract these to find

Now we can describe the fundamental ideas of a two-grid method schematically as follows:

• Relax Au=b on to obtain approximation
• Compute the residual
• Transfer the residual to the grid , and relax the residual equation Ae=r, to obtain an error {Interpolate the error back onto , and correct

If you think recursively, you can see that, to get a better correction on the coarse grid, you could coarsen again, continuing until, on a grid with one interior point, the equation can be solved exactly.

Details of the method will be demonstrated in class. For now, Here is a sample code available from the MGNet. This code solves a simple 1D problem. (In 1 space dimension, the data structures are easy.)

Copy the code. Note the grid size, governed by `K'. How well is the symmetry of the problem maintained by the code?

For HOMEWORK, change the Gauss-Seidel solver to Jacobi - indeed, to a weighted Jacobi, with weight 2/3.

{

C      _______________________________________________
C     |                                               |
C     |    SOLVE  -X''(T) = 1, X(0) = X(1) = 0        |
C     |          BY THE MULTIGRID METHOD              |
C     |_______________________________________________|
      REAL R(518),X(518),D0,D1,S,T,U
      INTEGER I,IT,J,K,L,LL,M,N,NL
      K = 8
      N = 2**K
      IT = 4
      T = .0001
      U = .7
C     --------------------------
C     |*** INPUT RIGHT SIDE ***|
C     --------------------------
      DO 10 I = 1,N
10         R(I) = 1.
C     --------------------
C     |*** INITIALIZE ***|
C     --------------------
      S = (1./N)**2
      DO 20 I = 1,N
20         R(I) = S*R(I)
      M = N
      L = 1
      LL = N
      NL = N + N + K - 2
      DO 30 I = 1,NL
30         X(I) = 0.0
      D1 = 0.
40    J = 0
C     --------------------------------
C     |*** GAUSS-SEIDEL ITERATION ***|
C     --------------------------------
50    D0 = D1
      D1 = 0.
      J = J + 1
      I = L
60    I = I + 1
           S = .5*(X(I-1)+X(I+1)+R(I))
           D1 = D1 + ABS(S-X(I))
           X(I) = S
      IF ( I .LT. LL ) GOTO 60
      WRITE(6,70) D1
70    FORMAT(' DIF:',F20.10)
      IF ( J .LT. IT ) GOTO 50
      IF ( D1  .LT.  T ) GOTO 100
      IF ( D1/D0 .LT. U ) GOTO 50
      IF ( N .EQ. 2 ) GOTO 50
C     -----------------------------------------
C     |*** COARSER MESH (SLOW CONVERGENCE) ***|
C     -----------------------------------------
      I = LL + 2
80    L = L + 2
      I = I + 1
      IF ( L .GT. LL ) GOTO 90
           X(I) = 0.
           R(I) = 4*(R(L)+X(L-1)-2*X(L)+X(L+1))
      GOTO 80
90    N = N/2
      WRITE(6,*) 'NUMBER OF MESH INTERVALS:',N
      LL = LL + N + 1
      L = L + 1
      GOTO 40
C     ---------------------------------------
C     |*** FINER MESH (FAST CONVERGENCE) ***|
C     ---------------------------------------
100   IF ( N .EQ. M ) GOTO 120
      I = L - 3
      J = LL
110        X(I) = X(I) + X(J)
           X(I+1) = X(I+1) + .5*(X(J)+X(J+1))
           I = I - 2
           J = J - 1
      IF ( J .GT. L ) GOTO 110
      N = N + N
      WRITE(6,*) 'NUMBER OF MESH INTERVALS:',N
      LL = L - 2
      L = I
      X(I+1) = X(I+1) + .5*X(J+1)
      GOTO 40
C     ----------------------
C     |*** PRINT ANSWER ***|
C     ----------------------
120   M = N + 1
      DO 130 I = 1,M
           J = I - 1
           S = J/FLOAT(N)
130        WRITE(6,140) J,S,X(I)
140   FORMAT(I5,F10.5,F15.6)
      END