Fourier solution

With the Vesuvius example we have used a numerical method to solve a problem that has an ``analytic solution''. Actually, it has an algebraic solution in Fourier space. For practical purposes this is wonderful because it means we can have a solution on a very large mesh. The cost is only about linear (actually $N\log N$) in the number of grid points while iterative solvers are much more costly. Let us now look at the ``analytic'' solution. Our least squares regression (73) takes the form  

$$\bold 0 \quad\approx\quad \bold{\nabla} \bold \phi \ -\ \bold d$$ (83)

The shortcut to the least squares solution is to apply the adjoint operator. The adjoint to the gradient $[\bold{\nabla}]$is the vector divergence $[\bold{\nabla} \cdot]$.The divergence of the gradient is the negative of the Laplacian. So, we get

$$\bold 0 \quad=\quad \nabla^2 \bold \phi \ +\ \bold{\nabla} \cdot \bold d$$ (84)

Let us visualize this equation in physical space and in Fourier space. In physical space $\bold{\nabla}^2$ is a convolutional operator. In Fourier space it is -(k_x²+k_y²). We can solve simply by dividing by -(k_x²+k_y²) whereas inverting the matrix $\bold{\nabla}^2$(which happens implicitly with conjugate gradients) would be a very big production. Thus, the analytic solution is

$$\phi(x,y) \quad=\quad -\ \bold{ FT}^{-1} { { \bold{ FT} \ \ \nabla\cdot \bold d }\over{ k_x^2+k_y^2 } }$$ (85)

where $\bold{ FT}$ denotes 2-dimensional Fourier transform over x and y.

Instead of representing k_x²+k_y² in the most obvious way, let us represent it in a manner consistant with the finite-difference way we expressed the data $\bold d$. Recall that $-i\omega\Delta t \approx -i\hat\omega\Delta t =1-Z= 1-\exp(-i\omega\Delta t)$which is a Fourier domain way of saying that difference equations tend to differential equations at low frequencies. Likewise a symmetric second time derivative has a finite-difference representation proportional to (-2+Z+1/Z) and in a two-dimensional space, a finite-difference representation of the Laplacian operator is proportional to (-4+X+1/X+Y+1/Y) where $X=\exp(ik_x\Delta x)$ and $Y=\exp(ik_y\Delta y)$.

Fourier solutions have their own peculiarities (periodic boundary conditions) which are not always appropriate in practice, but having these solutions available is often a nice place to start from when solving a problem that cannot be solved in Fourier space. For example, suppose we feel some data values are bad and we would like to throw out the regression equations involving the bad data points. We could define a weighting matrix starting from an identity matrix and replacing some of the ones by zeros. This defines $\bold W$.Now our regression (83) becomes

$$\bold 0 \quad\approx\quad \bold W \ (\bold{\nabla} \bold \ph... ...) \eq (\bold W \ \bold{\nabla}) \bold \phi \ -\ \bold W \bold d$$ (86)

This is a problem we know how to solve, a regression with an operator $\bold W\nabla$ and data $\bold W \bold d$.The weighted problem is not solveable in the Fourier domain because the operator $(\bold W\nabla)'\bold W\nabla$ has no simple expression in the Fourier domain. Thus we would use the analytic solution to the unweighted problem as a starting guess for the iterative solution to the real problem.

With the Vesuvius data how might we construct the weight $\bold W$? We have available the signal strength (which we have not used). We could let the weight be proportional to signal strength. We also have available the curl, which should vanish. Its non-vanishing is an indicator of questionable data which could be weighted down relative to other data.

The laboratory exercise is new this year so it may contain some unexpected difficulties. We're not sure it leads to clear solutions either. Anyway, you are given the Vesuvius data and all the programs in the book. Additionally, you are given a Fourier solver that produces the analytic solution. Please inspect both the Fourier solver and the solution it gets. Go to the web to see what pictures you can find of Vesuvius. Notice the radial drainage patterns on the amplitude of the original complex numbers. It is a little disturbing that we don't see these drainage patterns on the phase data (or maybe you can see them a little?). Any thoughts you have on that issue are certainly welcome. Any other thoughts you have on this lab are certainly welcome. This data is fun so we'd like to get this lab better focused for next year.