# Solution to the Poisson Equation in C++ using Successive Over-relaxation

Note: The first section of this blog post is mathematical and useful to understand, however, if you are just looking for an implementation you can skip through to the first code block.

Due to the interest in my previous blog post: Solution to the Laplace Equation in C++ using successive over-relaxation, I have decided to write a follow-up post on how to extend the method to allow for solutions to the Poisson Equation.

The Poisson Equation is the general case of the Laplace Equation and takes the following form:

$$\nabla^{2} \phi = f(\vec{r})$$

Working in $\mathbb{R}^{2}$ as before, we can rewrite this:

$$\frac{\partial^{2} \phi}{\partial x^{2}} + \frac{\partial^{2} \phi}{\partial y^{2}} = f(x,y)$$

Using the same central-difference method as before: $f^{\prime \prime}(x) \approx \frac{f(x+h) – 2f(x) + f(x-h)}{h^{2}}$, we get the following relation:

$$\frac{\phi(x+\delta x,y) – 2\phi(x,y) + \phi(x-\delta x,y)}{(\delta x)^{2}}+\frac{\phi(x,y+\delta y) – 2\phi(x,y) + \phi(x,y-\delta y)}{(\delta y)^{2}}=f(x,y)$$

If we sample at every $\delta x$ in the $x$-direction and $\delta y$ in the $y$-direction as before we have:

$$\frac{\phi_{i+1,j} – 2\phi_{i,j} + \phi_{i-1,j}}{(\delta x)^{2}} + \frac{\phi_{i,j+1} – 2\phi_{i,j} + \phi_{i,j-1}}{(\delta y)^{2}} = f_{i,j}$$

Rearranging this we get:

$$\phi_{i,j} = \frac{1}{2\left((\delta x)^{-2} + (\delta y)^{-2}\right)}\left[\frac{\phi_{i+1,j}+\phi_{i-1,j}}{(\delta x)^{2}}+ \frac{\phi_{i,j+1}+\phi_{i,j-1}}{(\delta y)^{2}} – f_{i,j}\right]$$

Setting $\delta x = \delta y = \delta$ and using this in the iterative formula from the previous blog post, we get:

$$\phi_{i,j}^{m+1}=\phi_{i,j}^{m} + \frac{\alpha}{4}\left(\phi_{i+1,j}^{m+1} + \phi_{i-1,j}^{m} + \phi_{i,j+1}^{m+1} + \phi_{i,j-1}^{m} – 4 \phi_{i,j}^{m} – \delta^{2} f_{i,j}\right)$$

We begin by declaring the template function prototype:

template <typename T>
std::size_t PoissonSolver( std::vector<std::vector<T>>& vecGrid,
const std::map<std::pair<std::size_t, std::size_t>, T>& mapPoisson,
T tDelta = T( 1.0 ), std::size_t sMaxNumIterations = 30,
T tEps = std::numeric_limits<T>::epsilon() );


The modification to the body of the LaplaceSolver is only four lines of code in the while loop, giving:

while( sCurrIteration <= sMaxNumIterations && !bPrecisionAchieved )
{
bPrecisionAchieved = true;
for( std::size_t sY = sHeight - 2U; sY >= 1U; --sY )
{
for( std::size_t sX = sWidth - 2U; sX >= 1U; --sX )
{
auto it = mapMask.find( std::make_pair( sX, sY ) );
if( it != mapMask.end() && it->second == false )
continue;

T tResidual = vecGrid[sX][sY + 1U] + vecGrid[sX][sY - 1U] + vecGrid[sX - 1U][sY] + vecGrid[sX + 1U][sY]
- T( 4.0 ) * vecGrid[sX][sY];

auto it2 = mapPoisson.find( std::make_pair( sX, sY ) );
if( it2 != mapPoisson.end() )
{
tResidual -= tDelta * tDelta * it2->second;
}

vecGrid[sX][sY] += tOverRelaxationParameter * tResidual;

bPrecisionAchieved = bPrecisionAchieved && std::fabs( tResidual ) < tEps;
}
}

++sCurrIteration;
}


This checks to see if the current co-ordinate has a non-zero function-value and then subtracts $\delta f_{i,j}$ from the residual. We can check that this code works by considering a point charge $Q$ in $\mathbb{R}^{2}$ in a grounded box, which obeys the Poisson Equation $\nabla^{2} V = -\frac{\rho}{\varepsilon_{0}}$. We know the analytic solution is given by:

$$V(x,y) = \begin{cases}-\frac{Q}{4\pi \varepsilon_{0} \sqrt{(x-0.5)^{2} + (y-0.5)^{2}}} & (x,y) \in [0,1] \times [0,1] \\ 0 & \text{otherwise}\end{cases}$$

Using the PoissonSolver code above, we have the following comparison ($25 \times 25$ grid on the left and analytic solution on the right):

This sort of numeric solution allows us to easily visualize more complex situations involving Poisson’s Equation such as those shown at the top of the page.

Designed and Produced by Thomas Russell © 2014-2017