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 ) $$


$$ -\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.

  • $\begingroup$ 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. $\endgroup$ – nicoguaro May 29 '18 at 15:20
  • $\begingroup$ 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. $\endgroup$ – Jeroen Reurink May 30 '18 at 7:53
  • $\begingroup$ I meant "for Dirichlet boundary conditions [...]". $\endgroup$ – nicoguaro May 30 '18 at 9:23

Your Answer

By clicking "Post Your Answer", you acknowledge that you have read our updated terms of service, privacy policy and cookie policy, and that your continued use of the website is subject to these policies.

Browse other questions tagged or ask your own question.