Problem 1 (MRI in the differential form) 40 pts

In this problem we will build a finite difference solver for an MRI problem reduced to the Helmholtz equation $$ \Delta u({\bf r}) + k^2({\bf r}) u({\bf r}) = v({\bf r}). $$

1. Free space solver

As a first attempt, let us construct a finite differences method simulation tool for the general electromagnetic scattering problem in free-space. In this case $k({\bf r}) \equiv k_0$ Discretize the equation on a box $[0,1]^2$ box using second order finite difference scheme and zero Dirichlet boundary conditions. Let $v$ be an impulse, located at $x = 0.5$ m and $y = 0.5$ m, discretized to a vector of zeros with a single element of $1/h^2$ at the appropriate row.

• (5 pts) Solve the scattering problem using the $\verb|scipy.sparse.linalg.solve|$ direct solver for $f = 21.3$ MHz and $f = 298.3$ MHz, corresponding to the older $0.5$ T permanent magnet MRI and the modern $7$ T superconductor magnet MRI respectively. Show your results as 2D images. (do not forget to use sparse arithmetics in python, e.g. $\verb|scipy.sparse.kron|$ or $\verb|scipy.sparse.spdiags|$ functions

• (10 pts) Using the discrete sine transform, implement two fast solvers for the Dirichlet problem: a second-order accurate solver, and a spectral solver. Using your two fast solvers and your direct $\verb|scipy.sparse.linalg.solve|$ solver, simulate scattering for $f = 298.3$ MHz. Using $\verb|pandas|$ dataframe produce a table with times to solve a system using these 3 solvers. Using exact solution of the equation produce a table with relative errors of all these solvers. Notes: Both tables should contain results for $n=128$, $n=256$ and $n=512$. To obtain the right results, you must use nested grids, i.e. a grid at some scale should contain all the points of a grid at some coarser scale.

2. EM Scattering of a head

• (5 pts) Upload MRI data files of a head with $257^2$ and $1025^2$ points. Solve the scattering problem using a second-order accurate scheme with the $\verb|scipy.sparse.linalg.solve|$ solver for $f = 21.3$ MHz and $f = 298.3$ MHz, corresponding to the $0.5$ T and the $7$ T MRI systems respectively. Show your results as 2D images cropped to the head.

• (2 pts) Can fast DST solver be used in case of nonconstant $k$? Explain the answer.

• (8 pts) The backslash operator becomes very computationally intensive for big $n$, and iterative methods are generally used instead. Here is a simple iterative scheme that makes use of our DST solver is the following equation: $$ \Delta u^{(m+1)} + k_0^2 u^{(m+1)} = v - (k^2 - k_0^2) u^{(m)}. $$ Implement this iterative method with the starting vector set to a vector of zeros, and with the inversion performed using your second-order accurate DST solver. Plot error with resepct to iteration number for $f = 21.3$ MHz and $f = 298.3$ MHz for different $n$. Does the rate of convergence depend on the frequency of excitation?

• (10 pts) Represent the latter iterative process in the form $$ u^{(m+1)} = u^{(m)} - \tau B^{-1} (v - A u^{(m)}), $$ where $A = \Delta + k^2 \mathcal{I}$. Specify relaxation parametr $\tau$ and preconditioner $B$. Give a condition on $k({\bf r})$, given a discretization spacing $h$, that guarantees the iterations will converge. Does this condition hold for both $f = 21.3$ MHz and $f = 298.3$ MHz and both grid sizes?

Problem 2 (Multigrid) (40 pts)

Consider a Poisson equation $$ \begin{align*} -\Delta u &= f, \\ u|_\Gamma &= 0, \end{align*} $$ where $\Omega = [0,1]^d$, $d = 2, 3$ and $\Gamma=\partial\Omega$.

Although the multigrid method has optimal complexity, there are some precomputations to be done. Therefore, it may have a bigger constant than, for instance, the spectral method or even sparse LU. In this problem you will be asked to find out in which case (2D, 3D, grid sizes) multigrid method is more appropriate.

• (20 pts) Implement V- and W-cycles of multigrid method for $d=2$ and $d=3$. Assume that the equations are discretized with finite differences with the second order scheme. Note: all operations should be done with linear grid size complexity. Take a smoother of your choice.

• (10 pts) Set $N=128$. Choose $f$ such that you know analytical solution. Plot errors as a function of the cycle number for the implemented V-,W-cycles of GMG and for AMG (from PyAMG), $d=2$ and $d=3$. Which approach requires fewer number of iterations? Which approach is faster? Note: when doing plottings do not forget to inculde title and axes names

• (10 pts) Set $f\equiv1$. Investigate when (for what cases) GMG is more appropriate than the sparse LU, the spectral method or the DST solver. To do so compare timings for different $N$, $d$ and accuracies. Make conclusions.