The goal of this post is to create a program to solve a differential equation in
2D. I will present the programming in C++ but the concepts here are easily
translated to other languages. I myself have used Julia a couple times. The
equation that we are going to solve is the Poisson equation with Dirichlet
boundary conditions in a `$1\times 1$`

square. The equation and the conditions
are given by

`\begin{align} \nabla f(x, y) &= u(x, y),\\ f(0, y) &= g(0, y),\\ f(1, y) &= g(1, y),\\ f(x, 0) &= g(x, 0),\\ f(x, 1) &= g(x, 1). \end{align}`

### Internal points

Before coding anything, we need to expand the `$\nabla f(x,y)$`

term using
finite differences on it. To expand the term is quite easy and results in
`\begin{align} \nabla f(x, y) &= \frac{\partial^2 f}{\partial x^2} + \frac{\partial^2 f}{\partial y^2}. \end{align}`

Finite differences is a method that allows us to remove the continuous `$x$`

and
`$y$`

variables to introduce discrete variables `$i$`

and `$j$`

which will be
used to iterate in a loop. This is an approximation and the accuracy increases
with the number of discrete points one wants to use, but makes the
implementation more complicated as well as costs more computing power. To
discretize a second order partial differential term with second order accuracy,
we have the relation
`\begin{align} \frac{\partial^2 f}{\partial x^2}(x, y) &\approx \frac{f(x_{i-1},y_i)-2f(x_{i}, y_i)+f(x_{i+1},y)}{\Delta x^2}, \end{align}`

where `$\Delta x = \mathrm{span}(x)/N$`

and `$N$`

is the number of points that a
given dimension has. Remember that the continuous variable `$x$`

has infinite
points, but when we use `$x_i$`

we use only a finite set of points of size
`$N$`

. The number of points `$N$`

will be same for both dimensions, so we have
`$\Delta y = \Delta x$`

.

Regarding `$u(x,y)$`

we just evaluate it at `$(i\Delta x, j\Delta x)$`

. This
way, we have our discrete equation

`\begin{align} c_{ij} = \sum_{m=-1,1}f(x_{i-m}, y_i) + \sum_{n=-1,1}f(x_i, y_{j-n}),\\ \frac{1}{\Delta x^2}\left(c_{ij} - 4f(x_i, y_i)\right) = u(i\Delta x, j\Delta x).\label{eq:discretizedPoisson} \end{align}`

Notice that Equation `\eqref{eq:discretizedPoisson}`

is our discretized Poisson
equation. As the focus of this post is to present code using an explicit method,
we isolate the term `$f_{ij}$`

like in
`\begin{align} f_{ij} &= \frac{c_{ij}}{4} -\frac{u_{ij}\Delta x^2}{4}. \label{eq:internalPoints} \end{align}`

We can use Equation `\eqref{eq:internalPoints}`

for any point internal to our
`$1\times 1$`

square. Nevertheless, for the boundary points, i.e where
`$x\in \{0,1\}$`

or `$y\in\{0,1\}$`

, we can not use it as we will not be able to
compute `$i-1$`

, `$i+1$`

, `$j-1$`

or `$j+1$`

. I mean, these points would be
outside our mesh.

### Boundary points

Regarding the boundary points, it is super easy with Dirichlet conditions. We
will just do
`\begin{align}f_{ij} = u_{ij}. \label{eq:boundaryPoints}\end{align}`

### Code

With Equations `\eqref{eq:internalPoints}`

and `\eqref{eq:boundaryPoints}`

, we
can finally code a simple solution to solve our Poisson equation. Notice that I
am testing this code with `g++ 6.2.0`

. I will present the main functions as well
as some of the utility functions I did.

Regarding utilities, I used a
timer class,
so that we can compute how long this method takes, an `abs`

function, so that I
know the absolute differences between iterations, and a `getMesh`

function to
allocate memory and give me a pointer. These three functions are presented in

```
class Timer
{
public:
Timer() : beg_(clock_::now()) {}
void reset() { beg_ = clock_::now(); }
double elapsed() const {
return std::chrono::duration_cast<second_>
(clock_::now() - beg_).count(); }
private:
typedef std::chrono::high_resolution_clock clock_;
typedef std::chrono::duration<double, std::ratio<1> > second_;
std::chrono::time_point<clock_> beg_;
};
double abs(double x) {
if(x < 0) return -x;
return x;
}
// Remember we are considering number of points in x == y
double ** getMesh(int n) {
// we use calloc to initialize everything with zeros
double ** mesh = (double**) calloc (n, sizeof(double*));
for (int i = 0; i < n; i++) {
mesh[i] = (double*) calloc (n, sizeof(double));
}
return mesh;
}
```

The boundary conditions need to be computed only once and it depends on the
function `$g(x,y)$`

that we pass as an argument to another function:

```
/**
g is a function of (x,y) that is applied on the boundaries
*/
void computeBoundaries(double ** mesh, int n, function<double(double, double)> g) {
int i, j;
double dx = 1 / (n - 1);
i = 0;
for (j = 0; j < n; j++) mesh[i][j] = g(i * dx, j * dx);
i = n - 1;
for (j = 0; j < n; j++) mesh[i][j] = g(i * dx, j * dx);
j = 0;
for (i = 0; i < n; i++) mesh[i][j] = g(i * dx, j * dx);
j = n - 1;
for (i = 0; i < n; i++) mesh[i][j] = g(i * dx, j * dx);
return;
}
```

For the internal points, we do something similar by accepting a function as argument

```
/**
f is a function that expects mesh and i, j and computes the value of
that point. It also takes the absolute change reference and adds to it.
The return of this function is the absolute change
*/
double computeInternalPoints(double ** mesh, int n, function<double(double**, int, int, double*)> f) {
double absoluteChange = 0.0;
// only internal points, that's why there is the loop goes from 0 to n - 2
for (int i = 1; i < n - 1; i++) {
for (int j = 1; j < n - 1; j++) {
mesh[i][j] = f(mesh, i, j, &absoluteChange);
}
}
return absoluteChange;
}
```

Now, having these functions, the main idea is simple: run
`computeInternalFunctions`

as long as `absoluteChange`

is smaller than a
threshold. In math representation, this can be represented by
`\begin{align} \sum_{i}\sum_{j} |\mathrm{mesh}_{ij} - \mathrm{mesh}^*_{ij}| < \epsilon, \end{align}`

where `$\epsilon$`

can be pretty small but increases the computation time for
the function to converge, i.e. the time for the program to get a solution. The
variable `$\mathrm{mesh}^*$`

is the variable `$\mathrm{mesh}$`

after all its
internal points were computed again by the function passed to
`computeInternalPoints`

.

The final code in C++ as well as code to plot it (in python) is available here. If you look at the code, I hardcoded the Dirichlet conditions there, but you can change it to your liking. To compile and plot it, I do the following:

```
$ g++ -std=c++14 poissonDirichlet.cpp -o solver_dirichlet
# run with a mesh of size 100 and output file is at out.bin
$ ./solver_dirichlet 100 out.bin
Finished computing! Latest absolute change was: 9.99971e-05
Total number of iterations on mesh: 9497
Total time spent: 4.71717 seconds
Total time per iteration: 0.000496701 seconds
File size: 80004 bytes
File successfully saved at: out.bin
# use python to plot the result on the screen
$ python plot.py --filename out.bin
```

The final result is

### Notes

This is not the most efficient implementation to solve but it is a good starting
point for whoever is starting. If you were curious about how time goes on as we
increase `$N$`

, you can just check at this table:

Mesh size (N) | Time in seconds | Number of Iterations | Time per Iteration (s) |
---|---|---|---|

20 | 0.0107726 | 350 | 3.07788e-05 |

50 | 0.2805 | 2326 | 0.000120593 |

100 | 4.72382 | 9497 | 0.000497401 |

200 | 80.427 | 38377 | 0.00209571 |

Thanks for reading! If you have any questions or suggestions, let me know! You can comment, send me an e-mail, whatever suits you!