[Back to main] [Printable version] [Leave a comment]

Solving a 2D Poisson equation with Neumann boundary conditions through discrete Fourier cosine transform [Leave a comment]

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.

Comments on page '/code/misc-notes/neumann-cosine'

Name:
Email (opt):
Message (no HTML support):
Comment type: Public Private (for page author only)
Spam check: + 2 = 8
Email is optional, not visible to other users, and only used for possible replies.
Please note that the comments are not guaranteed to stay on-line. Trivial questions may be removed after answering them. Corrective and additional notes are sometimes integrated into the main article. Thanks and critisicm are welcome, but are pruned away from time to time. Offensive, abusive, distasteful or aggressive messages are of course removed without a second thougt. Messages sometimes disappear for technical reasons and sometimes just because the administrator felt like removing them.