# Neumann Boundary with ADI-method (FDM); how do I implement this?

I am modelling the Heat equation in 2D in Python. I am using finite difference methods, more specifically the Alternating Direction Implicit method. The model works quite well with Dirichlet boundary conditions, but I can't figure out how to implement Neumann boundary conditions into the tridiagonal matrix.

The equation I'm modelling is: $$\frac{\partial u}{\partial t} = \alpha \left ( \frac{\partial u}{\partial x} + \frac{\partial u}{\partial y}\right )$$

The discretized result in ADI form looks like this:

$$-\gamma_x u^{l+\frac{1}{2}}_{i-1,j} + 2\left ( 1+\gamma_x \right )u^{l+\frac{1}{2}}_{i,j} - \gamma_x u^{l+\frac{1}{2}}_{i+1,j} = 2 u^{l}_{i,j} + \gamma_y \left (u^{l}_{i,j-1} -2u^{l}_{i,j} +u^{l}_{i,j+1} \right )$$

and

$$-\gamma_y u^{l+1}_{i,j-1} + 2\left ( 1+\gamma_y \right )u^{l+1}_{i,j} - \gamma_y u^{l+1}_{i,j+1} = 2 u^{l+\frac{1}{2}}_{i,j} + \gamma_x \left (u^{l+\frac{1}{2}}_{i-1,j} -2u^{l+\frac{1}{2}}_{i,j} +u^{l+\frac{1}{2}}_{i+1,j} \right )$$

And $\gamma$ is $\frac{\alpha \Delta t}{2 \Delta x^2}$. This gives tridiagonal matrices, which look like Ax = b. In the x-direction, this looks like:

$$A = \begin{bmatrix} 1+\gamma &-\gamma & & & & &\\ -\gamma& 2+2\gamma& -\gamma & & & &\\ &-\gamma & 2+2\gamma & -\gamma & & &\\ & &\ddots & \ddots &\ddots & &\\ & & & &-\gamma &2+2\gamma &-\gamma\\ & & & & &-\gamma &1+\gamma \end{bmatrix}$$

$$x = \begin{bmatrix} u^{l+\frac{1}{2}}_{0,j} \\ u^{l+\frac{1}{2}}_{1,j}\\ u^{l+\frac{1}{2}}_{2,j}\\ \vdots \\ u^{l+\frac{1}{2}}_{nx-1,j}\\ u^{l+\frac{1}{2}}_{nx,j} \end{bmatrix}$$

And finally,

$$b = \begin{bmatrix} g_{0,j}^{l} \\ 2u^l_{1,j} + \gamma_y \left ( u^l_{1,j-1} - 2u^l_{1,j} +u^l_{1,j+1} \right )\\ 2u^l_{2,j} + \gamma_y \left ( u^l_{2,j-1} - 2u^l_{2,j} +u^l_{2,j+1} \right )\\ \vdots \\ 2u^l_{nx-1,j} + \gamma_y \left ( u^l_{nx-1,j-1} - 2u^l_{nx-1,j} +u^l_{nx-1,j+1} \right )\\ g_{nx,j}^{l} \end{bmatrix}$$

These are the matrices in the X-sweep direction; the y-direction matrices and equations are the same except for different sub and superscripts.

$g_{0,j}^{l}$ and $g_{nx,j}^{l}$ are (in the current model) constants for the Dirichlet boundary conditions. The first row of the A-matrix has been modified to work for the Dirichlet boundaries, but I have no idea how to make the Neumann boundaries work. I tried the ghost point method, but for some reason this didn't work; it added a negative flux into the model.

I'm wondering if my discretization and matrices are correct. It might be really easy to implement Neumann boundaries, but for some reason I cant manage to. Thanks for reading, and sorry for the wall of text.

• I think that your first row should be $(1, 0, 0, \cdots, 0)$, and the last one $(0, \cdots, 0, 0, 1)$, for your boundary conditions to be satisfied. – nicoguaro May 29 '18 at 15:20
• This doesn't seem to work, it gives a negative flux inwards (the boundary cools the surrounding) while i'd like to have zero flux on some boundaries. – Jeroen Reurink May 30 '18 at 7:53
• I meant "for Dirichlet boundary conditions [...]". – nicoguaro May 30 '18 at 9:23