This is a derivation of the 2D Laplacian finite difference approximation on 2D grid with Neumann boundary conditions for solving the elliptic PDE. The PDE is\begin {align*} -\nabla ^{2}u & =f(x,y)\\ -\left ( \frac {\partial ^{2}u}{\partial x^{2}}+\frac {\partial ^{2}u}{\partial y^{2}}\right ) & =f(x,y) \end {align*}
The derivation is given for the general case of non-uniform grid (meaning \(h_{x}\) is not equal to \(h_{y}\)). At the end it is simplified to uniform grid by setting \(h_{x}=h_{y}\).
Notation used is based on the following grid where it is assumed that the grid size is \(n\times n\)
Starting with centered difference scheme\[ -\left ( \frac {U_{i,j-1}-2U_{i,j}+U_{i,j+1}}{h_{x}^{2}}+\frac {U_{i-1,j}-2U_{i,j}+U_{i+1,j}}{h_{y}^{2}}\right ) =f_{i,j}\] And solving for \(U_{i,j}\) gives\begin {align*} 2\left ( \frac {1}{h_{x}^{2}}+\frac {1}{h_{y}^{2}}\right ) U_{ij}+\frac {-U_{i,j-1}-U_{i,j+1}}{h_{x}^{2}}+\frac {-U_{i-1,j}-U_{i+1,j}}{h_{y}^{2}} & =f_{i,j}\\ 2\left ( \frac {h_{x}^{2}+h_{y}^{2}}{h_{x}^{2}h_{y}^{2}}\right ) U_{ij}+\frac {-h_{y}^{2}\left ( U_{i,j-1}+U_{i,j+1}\right ) -h_{x}^{2}\left ( U_{i-1,j}+U_{i+1,j}\right ) }{h_{x}^{2}h_{y}^{2}} & =f_{i,j}\\ \frac {2\left ( h_{x}^{2}+h_{y}^{2}\right ) }{h_{x}^{2}h_{y}^{2}}U_{ij}-\frac {h_{y}^{2}\left ( U_{i,j-1}+U_{i,j+1}\right ) +h_{x}^{2}\left ( U_{i-1,j}+U_{i+1,j}\right ) }{h_{x}^{2}h_{y}^{2}} & =f_{i,j} \end {align*}
Hence\begin {equation} U_{i,j}=\frac {1}{2\left ( h_{x}^{2}+h_{y}^{2}\right ) }\left [ h_{y}^{2}\left ( U_{i,j-1}+U_{i,j+1}\right ) +h_{x}^{2}\left ( U_{i-1,j}+U_{i+1,j}\right ) \right ] +\frac {h_{x}^{2}h_{y}^{2}}{2\left ( h_{x}^{2}+h_{y}^{2}\right ) }f_{i,j}\tag {1A} \end {equation} Therefore, for node \(i,j\) the equation \(-\nabla ^{2}u=f\left ( x,y\right ) \) becomes\begin {align*} -\left ( \frac {h_{y}^{2}\left ( U_{i,j-1}+U_{i,j+1}\right ) +h_{x}^{2}\left ( U_{i-1,j}+U_{i+1,j}\right ) -2\left ( h_{x}^{2}+h_{y}^{2}\right ) U_{i,j}}{h_{x}^{2}h_{y}^{2}}\right ) & =f_{i,j}\\ -2\left ( h_{x}^{2}+h_{y}^{2}\right ) U_{i,j}+h_{y}^{2}\left ( U_{i,j-1}+U_{i,j+1}\right ) +h_{x}^{2}\left ( U_{i-1,j}+U_{i+1,j}\right ) & =-h_{x}^{2}h_{y}^{2}f_{i,j} \end {align*}
The above formula reduces to the standard one when \(h_{x}=h_{y}\)\[ U_{i,j}=\frac {U_{i,j-1}+U_{i,j+1}+U_{i-1,j}+U_{i+1,j}+h^{2}f_{i,j}}{4}\] Therefore, the equation \(-\nabla ^{2}u=f\left ( x,y\right ) \) becomes\[ -\left ( -4U_{i,j}+U_{i,j-1}+U_{i,j+1}+U_{i-1,j}+U_{i+1,j}\right ) =h^{2}f_{i,j}\]
Below is the derivation of the discretization for the case when Neumann boundary conditions are used.
Given a 2D grid, if there exists a Neumann boundary condition on an edge, for example, on the left edge, then this implies that \(\frac {\partial u}{\partial x}\) in the normal direction to the edge is some function of \(y\).
Assuming that \(\frac {\partial u}{\partial x}=g\left ( y\right ) \) on the left edge is as shown in the following diagram
The value of the unknown \(u\left ( x,y\right ) \) is not given on the left boundary, only its derivative (in the normal direction to the edge, pointing outwards) is given. Therefore it is not possible to use the standard 5-point Laplacian on the first column \(j=1\), since that requires \(u\) be known on the left edge.
Therefore, when \(j=1\), then \(U_{i,1}\) is not known (if this was a Dirichlet boundary conditions, then \(U_{i,1}\) will be known).
Neumann boundary conditions are handled as follows. Taking the case of the left edge, an imaginary boundary located \(h_{x}\) to the left of the actual left edge is added as shown below
Then \(\frac {\partial u}{\partial x}\) is approximated on the line \(j=1\) by using \(\frac {\partial u}{\partial x}\approx \frac {U_{i,2}-U_{i,0}}{2h_{x}}\). Therefore for \(j=1\) the following results\begin {equation} \frac {U_{i,2}-U_{i,0}}{2h_{x}}=g_{i,1}\tag {2} \end {equation} From eq (1A) above, given below again\[ U_{i,j}=\frac {1}{2\left ( h_{x}^{2}+h_{y}^{2}\right ) }\left [ h_{y}^{2}\left ( U_{i,j-1}+U_{i,j+1}\right ) +h_{x}^{2}\left ( U_{i-1,j}+U_{i+1,j}\right ) \right ] +\frac {h_{x}^{2}h_{y}^{2}}{2\left ( h_{x}^{2}+h_{y}^{2}\right ) }f_{i,j}\] When \(j=1\) the above becomes\begin {equation} U_{i,1}=\frac {1}{2\left ( h_{x}^{2}+h_{y}^{2}\right ) }\left [ h_{y}^{2}\left ( U_{i,0}+U_{i,2}\right ) +h_{x}^{2}\left ( U_{i-1,1}+U_{i+1,1}\right ) \right ] +\frac {h_{x}^{2}h_{y}^{2}}{2\left ( h_{x}^{2}+h_{y}^{2}\right ) }f_{i,j}\tag {3} \end {equation} But using Eq. (2) gives \[ U_{i,0}=U_{i,2}-2h_{x}g_{i,1}\] Substituting the above into Eq. (3) results in\begin {equation} U_{i,1}=\frac {1}{2\left ( h_{x}^{2}+h_{y}^{2}\right ) }\left [ h_{y}^{2}\left ( U_{i,2}-2h_{x}g_{i,1}+U_{i,2}\right ) +h_{x}^{2}\left ( U_{i-1,1}+U_{i+1,1}\right ) \right ] +\frac {h_{x}^{2}h_{y}^{2}}{2\left ( h_{x}^{2}+h_{y}^{2}\right ) }f_{i,j}\tag {3A} \end {equation} An expression for \(U\) on the left edge is obtained. This will be added to the set of the equations to solve for.
Hence, For Neumann boundary conditions, the process starts by looking for the edge which have this condition specified on it, and for each such edge a set of equations as the above are added to the current set of equations for the internal points.
Therefore, for node \(\left ( i,1\right ) \) the equation \(-\nabla ^{2}u=f\) \(\ \)becomes\begin {align*} \frac {2\left ( h_{x}^{2}+h_{y}^{2}\right ) U_{i,1}-h_{y}^{2}\left ( 2U_{i,2}-2h_{x}g_{i,1}\right ) -h_{x}^{2}\left ( U_{i-1,1}+U_{i+1,1}\right ) }{h_{x}^{2}h_{y}^{2}} & =f_{i,j}\\ 2\left ( h_{x}^{2}+h_{y}^{2}\right ) U_{i,1}-h_{y}^{2}\left ( 2U_{i,2}-2h_{x}g_{i,1}\right ) -h_{x}^{2}\left ( U_{i-1,1}+U_{i+1,1}\right ) & =h_{x}^{2}h_{y}^{2}f_{i,j} \end {align*}
The edge points which have Neumann boundary conditions are initialized the same way as the internal units before the iterative process begins.
Eq. (3A) above reduces to the standard one when \(h_{x}=h_{y}\) giving\[ U_{i,1}=\frac {1}{4}\left ( 2U_{i,2}-2h_{x}g_{i,1}+U_{i-1,1}+U_{i+1,1}+h^{2}f_{i,j}\right ) \] Therefore, the equation \(-\nabla ^{2}u=f\) for such a node becomes\[ 4U_{i,1}-2U_{i,2}-U_{i-1,1}-U_{i+1,1}=h^{2}f_{i,j}-2h_{x}g_{i,1}\] The term \(2h_{x}g_{i,1}\) is moved to the RHS since it is a known quantity and will go into the force vector since only unknowns are kept in the LHS.
Now the case for right edge is considered. Assuming that \(\frac {\partial u}{\partial x}=q\left ( y\right ) \) is given.
\(\frac {\partial u}{\partial x}\) is approximated on the line \(j=end\) by using \(\frac {\partial u}{\partial x}\approx \frac {U_{i,end-1}-U_{i,end+1}}{2h_{x}}\), hence for the last column (\(j=n)\) this results in\begin {equation} \frac {U_{i,n-1}-U_{i,n+1}}{2h_{x}}=q_{i,n}\tag {1} \end {equation} From Eq. (1A) above, shown below again\[ U_{i,j}=\frac {1}{2\left ( h_{x}^{2}+h_{y}^{2}\right ) }\left [ h_{y}^{2}\left ( U_{i,j-1}+U_{i,j+1}\right ) +h_{x}^{2}\left ( U_{i-1,j}+U_{i+1,j}\right ) \right ] +\frac {h_{x}^{2}h_{y}^{2}}{2\left ( h_{x}^{2}+h_{y}^{2}\right ) }f_{i,j}\] For \(j=n\) it becomes\begin {equation} U_{i,n}=\frac {1}{2\left ( h_{x}^{2}+h_{y}^{2}\right ) }\left [ h_{y}^{2}\left ( U_{i,n-1}+U_{i,n+1}\right ) +h_{x}^{2}\left ( U_{i-1,n}+U_{i+1,n}\right ) \right ] +\frac {h_{x}^{2}h_{y}^{2}}{2\left ( h_{x}^{2}+h_{y}^{2}\right ) }f_{i,n}\tag {2} \end {equation} But from Eq. (1) \[ U_{i,n+1}=U_{i,n-1}-2h_{x}q_{i,n}\] Substituting the above into Eq. (2) results in\begin {equation} U_{i,n}=\frac {1}{2\left ( h_{x}^{2}+h_{y}^{2}\right ) }\left [ h_{y}^{2}\left ( 2U_{i,n-1}-2h_{x}q_{i,n}\right ) +h_{x}^{2}\left ( U_{i-1,n}+U_{i+1,n}\right ) \right ] +\frac {h_{x}^{2}h_{y}^{2}}{2\left ( h_{x}^{2}+h_{y}^{2}\right ) }f_{i,n}\tag {3} \end {equation} Therefore, for node \(\left ( i,n\right ) \) the equation \(-\nabla ^{2}u=f\) \(\ \)becomes\[ \frac {2\left ( h_{x}^{2}+h_{y}^{2}\right ) U_{i,n}-h_{y}^{2}\left ( 2U_{i,n-1}-2h_{x}q_{i,n}\right ) -h_{x}^{2}\left ( U_{i-1,n}+U_{i+1,n}\right ) }{h_{x}^{2}h_{y}^{2}}=f_{i,n}\] Eq. (3) reduces to the standard one when \(h_{x}=h_{y}\) giving\[ U_{i,n}=\frac {1}{4}\left ( 2U_{i,n-1}-2h_{x}q_{i,n}+U_{i-1,n}+U_{i+1,n}+h^{2}f_{i,n}\right ) \] Therefore, equation \(-\nabla ^{2}u=f\) for such a node becomes\[ 4U_{i,n}-2U_{i,n-1}-U_{i-1,n}-U_{i+1,n}=h^{2}f_{i,n}-2hq_{i,n}\]
Assuming \(\frac {\partial u}{\partial y}=\beta \left ( x\right ) \) is given for the top edge. \(\frac {\partial u}{\partial x}\) is approximated on the line \(i=1\) by using \(\frac {\partial u}{\partial y}\approx \frac {U_{2,j}-U_{0,j}}{2h_{y}}\), resulting in\begin {equation} \frac {U_{2,j}-U_{0,j}}{2h_{y}}=\beta _{1,j}\tag {1} \end {equation} Using Eq. (1A) above, shown below again\[ U_{i,j}=\frac {1}{2\left ( h_{x}^{2}+h_{y}^{2}\right ) }\left [ h_{y}^{2}\left ( U_{i,j-1}+U_{i,j+1}\right ) +h_{x}^{2}\left ( U_{i-1,j}+U_{i+1,j}\right ) \right ] +\frac {h_{x}^{2}h_{y}^{2}}{2\left ( h_{x}^{2}+h_{y}^{2}\right ) }f_{1,j}\] For the \(i=1\) it becomes\begin {equation} U_{1,j}=\frac {1}{2\left ( h_{x}^{2}+h_{y}^{2}\right ) }\left [ h_{y}^{2}\left ( U_{1,j-1}+U_{1,j+1}\right ) +h_{x}^{2}\left ( U_{0,j}+U_{2,j}\right ) \right ] +\frac {h_{x}^{2}h_{y}^{2}}{2\left ( h_{x}^{2}+h_{y}^{2}\right ) }f_{1,j}\tag {2} \end {equation} But from Eq. (1) \[ U_{0,j}=U_{2,j}-2h_{y}\ \beta _{1,j}\] Substituting the above into Eq. (2) results in\begin {equation} U_{1,j}=\frac {1}{2\left ( h_{x}^{2}+h_{y}^{2}\right ) }\left [ h_{y}^{2}\left ( U_{1,j-1}+U_{1,j+1}\right ) +h_{x}^{2}\left ( 2U_{2,j}-2h_{y}\ \beta _{1,j}\right ) \right ] +\frac {h_{x}^{2}h_{y}^{2}}{2\left ( h_{x}^{2}+h_{y}^{2}\right ) }f_{1,j}\tag {3} \end {equation} Therefore, for node \(\left ( 1,j\right ) \) equation \(-\nabla ^{2}u=f\) \(\ \)becomes\begin {align*} \frac {2\left ( h_{x}^{2}+h_{y}^{2}\right ) U_{1,j}-h_{y}^{2}\left ( U_{1,j-1}+U_{1,j+1}\right ) -h_{x}^{2}\left ( 2U_{2,j}-2h_{y}\ \beta _{1,j}\right ) }{h_{x}^{2}h_{y}^{2}} & =f_{1,j}\\ 2\left ( h_{x}^{2}+h_{y}^{2}\right ) U_{1,j}-h_{y}^{2}\left ( U_{1,j-1}+U_{1,j+1}\right ) -h_{x}^{2}\left ( 2U_{2,j}-2h_{y}\ \beta _{1,j}\right ) & =h_{x}^{2}h_{y}^{2}f_{1,j} \end {align*}
Eq. (3) reduces to the standard one when \(k=0\) and when \(h_{x}=h_{y}\) giving\[ U_{1,j}=\frac {U_{1,j-1}+U_{1,j+1}+2U_{2,j}-2h\ \beta _{1,j}+h^{2}f_{1,j}}{4}\] And equation \(-\nabla ^{2}u=f\) reduces to\[ 4U_{1,j}-U_{1,j-1}-U_{1,j+1}-2U_{2,j}=h^{2}f_{1,j}-2h\ \beta _{1,j}\]
Assuming \(\frac {\partial u}{\partial y}=\alpha \left ( x\right ) \) is given on the bottom edge.
\(\frac {\partial u}{\partial x}\) is approximated on the line (\(i=n)\) by using \(\frac {\partial u}{\partial y}\approx \frac {U_{n-1,j}-U_{n+1,j}}{2h_{y}}\) giving\begin {equation} \frac {U_{n-1,j}-U_{n+1,j}}{2h_{y}}=\alpha _{n,j}\tag {1} \end {equation} From Eq. (1A) above, shown below again\[ U_{i,j}=\frac {1}{2\left ( h_{x}^{2}+h_{y}^{2}\right ) }\left [ h_{y}^{2}\left ( U_{i,j-1}+U_{i,j+1}\right ) +h_{x}^{2}\left ( U_{i-1,j}+U_{i+1,j}\right ) \right ] +\frac {h_{x}^{2}h_{y}^{2}}{2\left ( h_{x}^{2}+h_{y}^{2}\right ) }f_{n,j}\] For \(i=n\) it becomes\begin {equation} U_{n,j}=\frac {1}{2\left ( h_{x}^{2}+h_{y}^{2}\right ) }\left [ h_{y}^{2}\left ( U_{n,j-1}+U_{n,j+1}\right ) +h_{x}^{2}\left ( U_{n-1,j}+U_{n+1,j}\right ) \right ] +\frac {h_{x}^{2}h_{y}^{2}}{2\left ( h_{x}^{2}+h_{y}^{2}\right ) }f_{n,j}\tag {2} \end {equation} But from Eq. (1) \[ U_{n+1,j}=U_{n-1,j}-2h_{y}\ \alpha _{n,j}\] Substituting the above into Eq. (2) results in\begin {equation} U_{n,j}=\frac {1}{2\left ( h_{x}^{2}+h_{y}^{2}\right ) }\left [ h_{y}^{2}\left ( U_{n,j-1}+U_{n,j+1}\right ) +h_{x}^{2}\left ( 2U_{n-1,j}-2h_{y}\ \alpha _{n,j}\right ) \right ] +\frac {h_{x}^{2}h_{y}^{2}}{2\left ( h_{x}^{2}+h_{y}^{2}\right ) }f_{n,j}\tag {3} \end {equation} Therefore, for node \(\left ( n,j\right ) \) equation \(-\nabla ^{2}u=f\) \(\ \)becomes\[ \frac {\left ( 2\left ( h_{x}^{2}+h_{y}^{2}\right ) -h_{x}^{2}h_{y}^{2}k^{2}\right ) U_{n,j}-h_{y}^{2}\left ( U_{n,j-1}+U_{n,j+1}\right ) -h_{x}^{2}\left ( 2U_{n-1,j}-2h_{y}\ \alpha _{n,j}\right ) }{h_{x}^{2}h_{y}^{2}}=f_{n,j}\] Eq. (3) reduces to the standard one when \(k=0\) and when \(h_{x}=h_{y}\) giving\[ U_{n,j}=\frac {U_{n,j-1}+U_{n,j+1}+2U_{n-1,j}-2h\ \alpha _{n,j}+h^{2}f_{n,j}}{4}\] Therefore, equation \(-\nabla ^{2}u=f\) reduces to\[ 4U_{n,j}-U_{n,j-1}-U_{n,j+1}-2U_{n-1,j}=h^{2}f_{n,j}-2h\ \alpha _{n,j}\] The following diagram shows a typical stencil for the edges under Neumann boundary conditions
From Eq. (1A) it was found that\begin {equation} U_{i,j}=\frac {1}{2\left ( h_{x}^{2}+h_{y}^{2}\right ) }\left [ h_{y}^{2}\left ( U_{i,j-1}+U_{i,j+1}\right ) +h_{x}^{2}\left ( U_{i-1,j}+U_{i+1,j}\right ) \right ] +\frac {h_{x}^{2}h_{y}^{2}}{2\left ( h_{x}^{2}+h_{y}^{2}\right ) }f_{i,j}\tag {1} \end {equation} \(-\nabla ^{2}u=f\) is the following\[ -\left ( \frac {h_{y}^{2}\left ( U_{i,j-1}+U_{i,j+1}\right ) +h_{x}^{2}\left ( U_{i-1,j}+U_{i+1,j}\right ) -2\left ( h_{x}^{2}+h_{y}^{2}\right ) U_{i,j}}{h_{x}^{2}h_{y}^{2}}\right ) =f_{i,j}\] Hence\begin {align*} U_{i,j}^{k+1} & =\frac {1}{2\left ( h_{x}^{2}+h_{y}^{2}\right ) }\left [ h_{y}^{2}\left ( U_{i,j-1}^{k}+U_{i,j+1}^{k}\right ) +h_{x}^{2}\left ( U_{i-1,j}^{k}+U_{i+1,j}^{k}\right ) \right ] +\frac {h_{x}^{2}h_{y}^{2}}{2\left ( h_{x}^{2}+h_{y}^{2}\right ) }f_{i,j}\\ & =\frac {h_{y}^{2}\left ( U_{i,j-1}^{k}+U_{i,j+1}^{k}\right ) +h_{x}^{2}\left ( U_{i-1,j}^{k}+U_{i+1,j}^{k}\right ) +h_{x}^{2}h_{y}^{2}f_{i,j}}{2\left ( h_{x}^{2}+h_{y}^{2}\right ) } \end {align*}
The residue is\begin {align} r_{ij} & =f_{ij}+\nabla ^{2}u\nonumber \\ & =f_{ij}+\left ( \frac {h_{y}^{2}\left ( U_{i,j-1}+U_{i,j+1}\right ) +h_{x}^{2}\left ( U_{i-1,j}+U_{i+1,j}\right ) -2\left ( h_{x}^{2}+h_{y}^{2}\right ) U_{i,j}}{h_{x}^{2}h_{y}^{2}}\right ) \tag {2} \end {align}
Eq. (2) reduces to the following when \(h_{x}=h_{y}\) Where\[ r_{ij}=f_{ij}-\frac {4U_{ij}+U_{i,j-1}+U_{i,j+1}+U_{i-1,j}+U_{i+1,j}}{h^{2}}\]
When iteration runs from top to bottom, left to right, the following can be used\[ U_{i,j}^{k}=\frac {h_{y}^{2}\left ( U_{i,j-1}^{k}+U_{i,j+1}^{k}\right ) +h_{x}^{2}\left ( U_{i-1,j}^{k}+U_{i+1,j}^{k}\right ) +h_{x}^{2}h_{y}^{2}f_{i,j}}{2\left ( h_{x}^{2}+h_{y}^{2}\right ) }\] Some books write the above as\[ U_{i,j}^{k+1}=\frac {h_{y}^{2}\left ( U_{i,j-1}^{k+1}+U_{i,j+1}^{k}\right ) +h_{x}^{2}\left ( U_{i-1,j}^{k+1}+U_{i+1,j}^{k}\right ) +h_{x}^{2}h_{y}^{2}f_{i,j}}{2\left ( h_{x}^{2}+h_{y}^{2}\right ) }\] It leads to the same result in software, as long as iterations is run left to right, top to down, using the index notation given here.
Another option is to run the iterations bottom to top, right to left.
As long as one is consistent, same result will be obtained. The idea is to use the most recent updated cell in the matrix all the time.
The residue is\begin {align} r_{ij} & =f_{ij}+\nabla ^{2}u\nonumber \\ & =f_{ij}+\left ( \frac {h_{y}^{2}\left ( U_{i,j-1}+U_{i,j+1}\right ) +h_{x}^{2}\left ( U_{i-1,j}+U_{i+1,j}\right ) -2\left ( h_{x}^{2}+h_{y}^{2}\right ) U_{i,j}}{h_{x}^{2}h_{y}^{2}}\right ) \tag {2} \end {align}
It is important in the above, to find the residual before updating the cell. In the code, the following should be used\begin {align*} r_{ij} & =f_{ij}+\left ( \frac {h_{y}^{2}\left ( U_{i,j-1}^{k}+U_{i,j+1}^{k}\right ) +h_{x}^{2}\left ( U_{i-1,j}^{k}+U_{i+1,j}^{k}\right ) -2\left ( h_{x}^{2}+h_{y}^{2}\right ) U_{i,j}^{k}}{h_{x}^{2}h_{y}^{2}}\right ) \\ U_{i,j}^{k} & =\frac {h_{y}^{2}\left ( U_{i,j-1}^{k}+U_{i,j+1}^{k}\right ) +h_{x}^{2}\left ( U_{i-1,j}^{k}+U_{i+1,j}^{k}\right ) +h_{x}^{2}h_{y}^{2}f_{i,j}}{2\left ( h_{x}^{2}+h_{y}^{2}\right ) } \end {align*}