Multigrid

We consider the Poisson equation in 1-D:

 -u'' = f in (0,1)          and             u(0)=u(1)=0.              (1)

The tri-diagonal matrix using sparse matrix can be generated by

e = ones(n,1)/h;
A = spdiags([-e 2*e -e], -1:1, n, n);

1. Implement the Richardson, Jacobi, and Gauss-Seidal iteration.

   Hint: you can implement the matrix form of these
   methods. Check these related functions in MATLAB:
        diag,  tril,  triu
   

2. Plot the error (in the energy norm) vs iteration steps. You will get a figure similar to this.

   Hint: You can use direct solver to compute the exact solution for a
   given right hand side (rhs) or plug in an exact solution to get a
   rhs. The error in the energy norm can be computed by matrix
   mutiplication.  
   

3. Show the smoothing effect of the above iterative methods by plotting error vectors in the first few steps. For Richardson method, you will get a figure similar to this. Do you see the smoothing affect of Jacobi method?

   Hint: To include all frequency into the error vector, you can set the
   rhs as a constant vector and choose a random initial guess.
   

4. Implement V-cycle multigrid using Richardson or G-S as smoothers. Chose the stopping criterion as 1e-6 and list the iteration steps for different sizes in a table or plot the error (in the energy norm) vs iteration steps.

   Hint: To stop the iteration, you can compute the relative residual in l2 norm: 
   norm(b-A*u)/norm(b)
   or the relative error in the energy norm:
   sqrt((uexact-u)'A*(uexact-u))/sqrt(uexact'*A*uexact) provide the exact solution
   is avaiable.
   

5. Implement PCG method and use diagonal preconditioner and V-cycle multigrid preconditioner. Compare iteration steps for two preconditioners.

   Hint: Code PCG with diagonal preconditioner first and make your PCG works
   well. Then try to use V-cycle MG as the preconditioner.
   

6. (optional) Implement nonlinear Multigrid method (FAS) to solve the nonlinear equation (using finite difference method)

   -u''+ u^3 = 1 in (0,1), u'(0) = u'(1) = 0.
   

   Compare the convergence and computaitonal cost with Newton's method.