UC Berkeley Math 228B, Spring 2024, Problem Set 8

Per-Olof Persson

Optional, not graded

This content is protected and may not be shared, uploaded, or distributed.

  1. In this problem, you will write a multigrid solver for the linear system of equations generated by fempoi and pmesh from previous problem sets. Note that in 2-D it is very hard to be faster than the built-in backslash function, at least without using a compiled language. However, for large problems in 3-D, any multigrid solver should be superior to Gaussian elimination, so here we are more concerned about getting the right convergence behavior rather than a fast solver. To begin with, add two output arguments to the fempoi function to get access to the matrices \(A,b\) in the linear system:

       function [u, A, b] = fempoi(p, t, e)
    1. Write a MATLAB function

         function [u, res] = gauss_seidel(A, b, u, niter)

      that makes niter iterations using the Gauss-Seidel method for \(Au=b\): \[\begin{aligned} u_{m+1} = u_m + (D-L)^{-1} (b - Au_m), \end{aligned}\] starting from the input u and returning the last iterate u. \(D-L\) is the lower triangular part of \(A\), including the diagonal (see the tril command). The output res is a vector of length niter+1 with the infinity norms of the residuals \(b-Au_m\) at each iteration (including the initial and the final iterates). Try the function using the commands:

      pv = [0,0; 2,0; 1.5,1; .5,1; 0,0];
      [p, t, e] = pmesh(pv, 0.5, 3);
      [u0, A, b] = fempoi(p, t, e);
      [u, res] = gauss_seidel(A, b, 0*b, 1000);
      semilogy(0:1000, res)
    2. Write a MATLAB function

         function data = mginit(pv, hmax, nref)

      that computes all the required arrays for a multigrid solution of the Poisson problem using the mesh parameters pv,hmax,nref. Start from the pmesh function, and make appropriate modifications and additions.

      • data(i).p, data(i).t, data(i).e contain the mesh arrays \(p,t,e\) after \(i-1\) refinements, for \(i-1=0,\ldots,n_\mathrm{ref}\).

      • data(i).T contains the interpolation matrix \(T^{(i)}\) from grid \(i\) to grid \(i+1\), for \(i=1\,\ldots,n_\mathrm{ref}\). Use linear interpolation for all the new midpoints (that is, averaging of the neighboring nodes). The second output argument of unique might be useful.

      • data(i).R contains the restriction matrix \(R^{(i)}\) from grid \(i+1\) to grid \(i\). Use the transpose of \(T^{(i)}\), but with the rows scaled to have sums of 1.

      • data(nref+1).A, data(nref+1).b contain \(A,b\) for the finest grid (the actual linear system)

      • data(i).A contains the projected matrices \(A^{(i)}=R^{(i)}A^{(i+1)}T^{(i)}\) for \(i=1,\ldots,n_\mathrm{ref}\).

    3. Write a MATLAB function

         function [u, res] = mgsolve(data, vdown, vup, tol)

      that solves the problem precomputed in data, using multigrid V-cycles with vdown/vup pre/post-smoothing iterations using Gauss-Seidel, until the infinite norm of the residual is less than tol. The outputs are the solution u and the residuals res after each V-cycle (including the residual for the initial solution \(u=0\)).

      Test the function using the commands

      pv = [0,0; 2,0; 1.5,1; .5,1; 0,0];
      for iref = 1:5
          data = mginit(pv, 0.5, iref);
          [u, res] = mgsolve(data, 2, 2, 1e-10);
          semilogy(res), hold on
      hold off

      If everything works correctly, you should see a very fast convergence compared to pure Gauss-Seidel. More importantly, the number of iterations should not increase much when the grid is refined. This leads to the optimal \(\mathcal{O}(n)\) computational cost of the algorithm.