Solving a 2D Poisson equation with Neumann boundary conditions through discrete Fourier cosine transform

by JARNO ELONEN (elonen@iki.fi), 21.12.2004

The book NUMERICAL RECIPIES IN C, 2ND EDITION (by PRESS, TEUKOLSKY, VETTERLING & FLANNERY) presents a recipe for solving a discretization of 2D Poisson equation $ \frac{\partial^2 u}{\partial x^2} + \frac{\partial^2 u}{\partial y^2} = \rho ( x, y )$ numerically by Fourier transform ("rapid solver"). While it shows the explicit solution for the problem with several other boundary conditions, Neumann condition $ \nabla u = 0$ is handled quite briefly. The following demonstrates in detail how to derive an equation for $ \hat{u}_{m n}$ using the definition of the inverse Fourier cosine transform. The solution turns out, perhaps not very surprisingly, to be exacly the same as for the sine transform.

The finite difference equation is:

$\displaystyle u_{j + 1, l} + u_{j - 1, l} + u_{j, l + 1} + u_{j, l - 1} - 4 u_{j, l}$ $\displaystyle =$ $\displaystyle \rho_{j, l}$  

...from which we need to solve $ u_{j, l}$. First, some abbreviations:
$\displaystyle S_b^a$ $\displaystyle =$ $\displaystyle \sin ( \frac{a \pi b}{A} )$  
$\displaystyle C_b^a$ $\displaystyle =$ $\displaystyle \cos ( \frac{a \pi b}{A} )$  

Now, the inverse 2D discrete cosine transform is:
$\displaystyle f_{j, l}$ $\displaystyle =$ $\displaystyle \frac{4}{JL} \sum^{J - 1}_{m = 0} \sum^{L - 1}_{n = 0} \hat{f}_{m, n} C_m^j C^l_n$  

Substitute this to both sides of the above finite difference equation and remove the summation to obtain:
$\displaystyle \hat{u}_{m, n} ( C_m^{j + 1} C^l_n + C_m^{j - 1} C^l_n + C_m^j C^{l + 1}_n + C_m^j C^{l - 1}_n - 4 )$ $\displaystyle =$ $\displaystyle \Delta^2 \hat{\rho}_{m, n} C_m^j C^l_n$  

Using the following lemmas...
$\displaystyle \cos ( a + x )$ $\displaystyle =$ $\displaystyle \cos ( a ) \cos ( x ) - \sin ( a ) \sin ( x )$  
$\displaystyle \cos \left( \frac{( a + x ) \pi b}{A} \right)$ $\displaystyle =$ $\displaystyle \cos \left( \frac{a \pi b}{A} + \frac{x \pi b}{A} \right)$  
  $\displaystyle \Rightarrow$    
$\displaystyle C_b^{a + x}$ $\displaystyle =$ $\displaystyle C_b^a C_b^x - S^a_b S^x_b$  
$\displaystyle C_b^{- a}$ $\displaystyle =$ $\displaystyle C_b^a$  
$\displaystyle S_b^{- a}$ $\displaystyle =$ $\displaystyle - S_b^a$  

...we simplify and solve for $ \hat{u}_{m, n}$:
$\displaystyle C_m^{j + 1} C^l_n + C_m^{j - 1} C^l_n$ $\displaystyle =$ $\displaystyle ( C_m^j C_m^1 - S_m^j S_m^1 ) C^l_n + ( C_m^j C_m^{- 1} - S_m^j S_m^{- 1} ) C^l_n$  
  $\displaystyle =$ $\displaystyle C^l_n ( ( C_m^j C_m^1 - S_m^j S_m^1 ) + ( C_m^j C_m^1 + S_m^j S_m^1 ) )$  
  $\displaystyle =$ $\displaystyle 2 C^l_n C_m^j C_m^{1 ( j )}$  
  $\displaystyle \operatorname{and}$    
$\displaystyle C_m^j C^{l + 1}_n + C_m^j C^{l - 1}_n$ $\displaystyle =$ $\displaystyle C_m^j ( C_n^l C_n^1 - S_n^l S_n^1 ) + C_m^j ( C_n^l C_n^{- 1} - S_n^l S_n^{- 1} )$  
  $\displaystyle =$ $\displaystyle C_m^j ( ( C_n^l C_n^1 - S_n^l S_n^1 ) + ( C_n^l C_n^1 + S_n^l S_n^1 ) )$  
  $\displaystyle =$ $\displaystyle 2 C_m^j C_n^l C_n^{1 ( l )}$  
  $\displaystyle \Rightarrow$    
$\displaystyle 2 \hat{u}_{m, n} ( C^l_n C_m^j C_m^{1 ( j )} + C_m^j C_n^l C_n^{1 ( l )} - 2 )$ $\displaystyle =$ $\displaystyle \Delta^2 \hat{\rho}_{m, n} C_m^j C^l_n \Vert / ( C_m^j C^l_n ) \Rightarrow$  
$\displaystyle 2 \hat{u}_{m, n} ( C_m^{1 ( j )} + C_n^{1 ( l )} - 2 )$ $\displaystyle =$ $\displaystyle \Delta^2 \hat{\rho}_{m, n} \Rightarrow$  
$\displaystyle \hat{u}_{m, n}$ $\displaystyle =$ $\displaystyle \frac{\Delta^2 \hat{\rho}_{m, n}}{2 ( C_m^{1 ( j )} + C_n^{1 ( l )} - 2 )}$  

Hence, to solve the Poisson equation, first compute the 2D cosine transform $ \hat{\rho}_{m, n}$, then calculate...
$\displaystyle \hat{u}_{m, n} = \frac{\hat{\rho}_{m, n}}{2 ( \cos ( \frac{\pi m}{J} ) + \cos ( \frac{\pi n}{L} ) - 2 )} $
...and finally do an inverse cosine transform to obtain $ u_{j, l}$. In practice, when the denominator becomes 0, substitute it with some small $ \epsilon$ to avoid division by zero.