Computational Grid

The symmetry of the problem leads us to try a change of coordinate system, from Cartesian to polar: The mapping: $\left\{\begin{array}{l} x=r\cos\theta\\ y=r\sin\theta\\ x^2+y^2=r^2\end{array}\right.$

Hence, in the polar coordinate system we obtain a square lattice $\widetilde {\cal D}$ associated with ${\cal D}$, the cylindrical domain via the mapping above. The square grid $\widetilde {\cal D}$ has four edges, denoted by $\partial \widehat {\cal D}$. The two domains are illustrated in Figure 37.

Figure 37: The Cartesian and Polar domains ${\cal D}$ and $\tilde {\cal D}$, respectively

$$\begin{array}{lll} \mbox {let } & U(r,\theta)=U(r\cos\theta... ...\right.\\ & g(r,\theta)=g(r\cos\theta, \sin\theta) \end{array}$$

$$\nabla^2 = \frac{\partial ^2}{\partial x^2}+\frac{\partial ^2}{\partial y^2}$$$$\mbox {(cartesian)}\Rightarrow \nabla^2=\frac{\partial ^2}{\part... ...eta}+ \frac{1}{r^2}\frac{\partial ^2}{\partial \theta^2}\quad \mbox {(Polar)}$$

$\therefore$ (173) in Polar Coordinates is:

$$\therefore \frac{\partial ^2 u }{\pa... ...2u=g\quad \left\{\begin{array}{l} 0<r<1 0\le\theta\le 2\pi \end{array}\right.$$ (176)

Boundary Conditions: on $D$ we have Dirichlet condition on the edge of the disk. Two other conditions are based on our expectation of the solution. For example, we could demand that the solution be bounded at the center of the disk and that it be periodic in $\theta$.

$\widetilde D$ has 4 edges for which we need 4 conditions, that is, on ${\partial \widetilde D}\quad :\quad r=1, 0\le\theta\le 2\pi$

$$\mbox {the } U(r\cos\theta, r\sin\theta)=\phi(\theta) \mbox {gives}$$

$$\Rightarrow u(1,\theta)=\phi(\theta)\quad 0\le\theta\le 2\pi.$$

$$\mbox {on } 0<r<1\quad \theta =0, \quad \mbox {and } \quad 0<r<1, \theta=2\pi$$

$$u(r,0)=u(r,2\pi)\quad 0<r<1 \quad \mbox {(Periodic)}$$

Finally, we need a condition at $r=0$, for $0\le \theta\le 2\pi$. This whole line corresponds to just a single point in ${D}$, namely $x=0, y=0$: the procedure outlined in problem # 3 HW 9 is a general procedure for determining the condition at the origin. Here we could impose the condition that $u$ be constant along the line $r=0$. This implies that

$$\frac{\partial u}{\partial \theta}(0,\theta)=0\quad 0\le\theta\le 2\pi$$

At this point, we have a well-posed problem. We solve (177) on $\widetilde D$ subject to

$$\begin{array}{ll} u(1,\theta)=\phi & 0\le\theta\le 2\pi\\ u(r,0)... ...rac{\partial u}{\partial \theta} (0,\theta)=0 & 0\le \theta\le 2\pi \end{array}$$

Remark: The coordinate transformation is not only a proper choice of the coordinate system, preserving symmetrices of the solution: it also avoids grid errors due to interpolation which would be required in the discretization. Notice as well, that the coordinate system turned our problem into a ``separable'' PDE (see elementary PDE book).

We can use finite differences to solve the problem and then use the mapping to go from $\widetilde {\cal D}$ to ${\cal D}$ and thus obtain the approximate solution $u(x,y)$. We will use Fourier methods, combined with finite differences instead. Fourier methods are particularly suited to the solution of $L_2$ functions which are periodic. In this case, the periodicity is in $\theta$ hence in the $\theta$ direction we'll use Fourier, and we'll use finite differences for the radial direction The fastest algorithm for a discrete Fourier transform is the FFT (a useful reference for Fourier Methods is ``The DFT'' by the W. Briggs & V.E. Henson. The reason we call this section ``Fast Fourier Solver for the Poisson/Helmholtz Equation'' is because we're using FFT's rather than DFT's. In this sense it is fast.

Since the solution is periodic in $\theta$, we'll apply the Fourier transform:

$$\widehat u_m(r)=\frac{1}{2\pi}\int^{2\pi}_0 u(r,\theta)e^{-im\theta}d\theta\quad m\in {\mathbb{Z}}$$

$\widehat u$ is complex (in the implementation we could use a sine or cosine transform, but we will use the complex transform, for simplicity and generality). The goal is to convert the PDE into an infinite set of ODE's for the Fourier coefficients $$\left\{\hat u_m\right\} \displaystyle ^{\infty}_{m=-\infty}$$.

It's easy to deduce that

$$\widehat{\left(\frac{\partial u}{\partial r}\right)}_u = u'_m(r) \qquad \mbox { and} \widehat{\left(\frac{\partial ^2 u}{\partial r^2}\right)}_m = \widehat u''_{m}(r)\quad m\in {\mathbb{Z}}$$

$$\widehat{\left(\frac{ \partial ^2 u}{\partial \theta^2}\right)}_m = \frac{1}{2\pi}\int^{2\pi}_0\frac{\partial ^2 u}{\partial \theta^2... ...{2\pi}_0\frac{\partial u(r,\theta)}{\partial \theta}e^{-im\theta}d\theta\right]$$

$$= \frac{im}{2\pi}\int^{2\pi}_0\frac{\partial u(r,\theta)}{\partial \theta}e ^{-im\theta}d\theta =$$

$$= \frac{im}{2\pi}\left[e^{-im\theta}u(r,\theta)\Big\vert^{2\pi}_0+im\int^ {2\pi}_0 u(r,\theta)e^{-im\theta}d\theta\right] =-m^2\hat u_m(r)$$

So, we multiply (177) by $e^{-im\theta}$ and integrate from 0 to $2\pi$ and divide by $2\pi$. Using orthoqonality, we obtain the set of ODE's for the Fourier coefficients:

$$\widehat u''_m+\frac{1}{r}\widehat u'_m-\frac{m^2}{r^2}\widehat u_m+\kappa^2\widehat u_m=\widehat g_m \quad \mbox {for } r\in (0,1)\qquad m\in {\mathbb{Z}}$$ (177)

We also transform the boundary conditions:

$$\widehat u_m(1)=\widehat \phi_m\quad m\in {\mathbb{Z}}$$

(Note that solution is already periodic in $2\pi$ by having used Fourier transforms). Now we need to take care of the condition $$\frac{\partial }{\partial \theta}u(0,\theta)=0$$:

$$\mbox {Since } u(r,\theta)=\sum^{\infty}_{m=-\infty}\widehat u_m(r) e^{im\theta}\quad 0\le r\le 1,\quad 0\le\theta\le 2\pi$$

differentiate with respect to $\theta$ and set $r=0$:

$$i\sum^{\infty}_{m=-\infty}m\hat u_m(0)e^{im\theta}\equiv 0,\quad 0\le\theta\le 2\pi$$

Since the $e^{im\theta}$ are linearly independent, for $m\in{\mathbb{Z}}\Rightarrow$ all coefficients must vanish identically:

$$\hat v_m(0)=0\quad m\in {\mathbb{Z}}\backslash \{0\}$$

This leaves us with just a single missing item of information, namely the boundary condition for the zeroth harmonic.

$$\mbox {Since } u(r,\theta)=U(r\cos\theta, r\sin\theta)$$

$$\therefore\frac{\partial u(0,\theta)... ...\theta \frac{\partial U(0,0)}{dx}+\sin\theta \frac{\partial U(0,0)}{\partial y}$$

$$\therefore\hat{u}'(0)=\frac{1}{2\pi}... ...x}+ \frac{1}{2\pi}\int^{2\pi}_0\sin \theta \frac{\partial U(0,0)}{\partial y}=0$$

Since

$$\int^{2\pi}_0\cos\theta d\theta = \int^{2\pi}_0\sin\theta d\theta =0 \quad \therefore$$

$$u'(0)=0$$ (178)

Now, we're ready on the ODE part.

Remark: The key fact is that the fourier transform uncouples the harmonics. We need to solve the ODE system and then use the ``inverse'' Fourier formula

$$u(r,\theta)=\sum^{\infty}_{m=-\infty}\hat u_m(r)e^{im\theta}$$

It turns out that (178) has an analytical solution in terms of Bessel functions. More generally, we would use one of the boundary value techniques presented in this course, i.e. finite differences, shooting method, FEM.

Exercise) Show that, by using center differences, the FD approximation $w(k\Delta r)=u(k\Delta r)+{\cal O} (\Delta r^2)$ is given by the solution of

$$\left[ \begin{array}{l} (1-\frac{1}{... ...}=0 \quad (\mbox {forward difference approx of (\ref{ft7})}) \end{array}\right.$$ (179)

The outcome is a tridiagonal system for every $m\ne 0$ and an almost tridiagonal system for $m=0$. Such systems can be efficiently solved by sparse $LU$ factorization.

Implementation comments: Of course, we need to truncate the infinite set of ODE's by a subset. Say

$$-M+1\le m\le M$$

(the truncation error induced depends on how large $M$ is. See Canutto, Quarteroni, Hussaini, book on Spectral Methods). Provided $\phi$ and $g$ are smooth, good accuracy is attained and the error drops exponentially with $M$!!

(1)

Use an FFT, so pick $M=2^{n-1}\quad n=1, 2, 3, \cdots$ transform $g$ and $\phi$ to get $\widehat g$ and $\widehat \phi$
Pick $d$ and solve (180) for $\widehat w_{m,k}$ for $-M+1\le m\le M$ and $k=1, 2, \cdots d-1$

(2)

Then employ a $d-1$ inverse FFT to produce $w$ on a $d\times(2m)$ square grid.

(3)

find approximation of $u$ on $D$ by using the mapping on $\widetilde{D}$.

Roughly, we need ${\cal O}(M\log_2 M)$ ops for FFT in (175), ${\cal O}(dM)$ for LU solution in FD solution in (175), plus ${\cal O}(dM\log_2M)$ for the reconstruction. This compares very favorably with full finite difference methods.