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

Per-Olof Persson
persson@berkeley.edu

Optional, not graded

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

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
end
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.