2D Poisson equation with Dirichlet and Neumann boundary conditions

Question

The problem is given as follows

\begin{align} -\Delta u &= f, \text{in} \: \Omega \\ u &= 0, \text{on} \: \delta \Omega_D \\ H(u) &= 0, \text{on} \: \delta \Omega_N, \end{align} where $\Omega$ is the domain to be considered and the Dirichlet and Neumann boundary conditions are applied on $\delta \Omega_D$ and $\delta \Omega_N$, respectively. What does $H(u)$ mean in this case?

This is an exam question, and I have the solution manual available. I'm trying to understand the proposed solution where the 5-point formula is used to obtain a matrix $A$ such that $AU=F$. I've set up the 5-point formula which looks like this $$ -\frac{1}{h^2}(U_{i+1}^{j}+U_{i-1}^{j}+U_{i}^{j+1}+U_{i}^{j-1}-4U_{i}^{j}) = f_i^j,$$ where the indexes $i$ and $j$ are iterated in the $x$ and $y$ direction, respectively, and $U_{i}^{j}$ is the numerical approximation of $u(x_i^j,y_i^j)$.

In the solution manual, the points on the boundary $\delta\Omega_N$ rely on twice the contribution from the points in $\Omega$. E.g. say the point $(X^j_i,Y^j_i) \in \delta\Omega_N$ has a neighbouring point $(X^{j}_i,Y^{j+1}_i) \in \Omega$ which lies right above. Then according to the solution manual, the 5-point formula gives $$-\frac{1}{h^2}(2U_{i}^{j+1}-4U_{i}^{j}) = f_i^j.$$ I suspect my derivation of the 5-point formula above doesn't take into account the correct given boundary conditions. Can you explain how the contribution from the inner points of the domain $\Omega$ is twice as much when evaluating the points on the Neumann boundary $\delta \Omega_N$?

Answer 1

• The meaning of $H(u)$ is a matter of convention. A Neumann boundary condition (BC) writes $\frac{\partial u}{\partial n}(x,y) = k(x,y)$. The derivative with respect to $n$ means "in the direction normal to the boundary" with an agreed convention on the sign of this normal. I'll get back to this later. For now, you can put $k$ on the left-hand side and call that $H(u)$, i.e. $H(u)=\frac{\partial u}{\partial n}(x,y) - k(x,y)$ . In your example I suspect $k\equiv 0$, but let's assume it is not for now.

• The Neumann BC involves a derivative and we need to represent it using finite differences. You should use an order of finite difference that is the same as the one you are using inside of the domain. The 5 points stencil is second order. We thus need to express the Neumann BC using a second order scheme. We have two main options: central differences or forward/backward differences.

• To go further we need to know which boundary we are talking about. It seem your example concern a boundary in a plane $x=cst$ which is the "inlet" of the domain. Such a boundary would typically correspond to the value $j=0$. Note that all the things derived from now on are specific to this boundary. The normal to this boundary, pointing towards the inside of the domain, is thus $n=e_y$. The derivative "in the direction normal to the boundary", will be taken as the derivative along $y$ that is we need to ensure that $\frac{\partial u}{\partial y}(x_i,y_j) = k(x_i,y_j)$, where $j=0$. Now we need to represent that using a second order finite difference scheme: central differences or forward differences (forward since we are at the inlet and we can't really use information backward). Let's take the example of the central differences.

• If we use central differences, we will go backward by one point. I've said we couldn't use information backward, but in this case it introduces an "unknown" that we call a ghost point, and we can solve for its value since we know $k$. The central difference approximation writes: \begin{align} j=0, \quad \frac{\partial u}{\partial y}(x_i,y_j)=k(x_i,y_j) \rightarrow \frac{U_i^{j-1}-U_i^{j+1}}{2 \Delta y} = k(x_i,y_j) \end{align} which can be solved for the ``ghost point'' as \begin{align} j=0,\quad U_i^{j-1} = U_i^{j+1} + 2 k(x_i,y_j)\Delta y \end{align} In your case $\Delta y=h$. Note that if we were at the outlet (say $j=n$), the ghost point would be at $n+1$. At this point, the normal towards the domain is along $-e_y$ and we will have: \begin{align} j=n, \quad U_i^{j+1} =U_i^{j-1} + 2 k(x_i,y_j)\Delta y \end{align} This is a matter of convention between what we put in the variable $k$ and what we choose to be the normal.

• Now we need to ensure that the boundary condition is met for the Poisson equation. We write the Poisson equation at the boundary point itself (that's just the general formula, at $j=0$): $$ j=0,\quad -\frac{1}{h^2}(U_{i+1}^{j}+U_{i-1}^{j}+U_{i}^{j+1}+[U_{i}^{j-1}]-4U_{i}^{j}) = f_i^j $$ we can replace $[U_i^{j-1}]$ in this formula by the value we obtained above: $$j=0,\quad -\frac{1}{h^2}(U_{i+1}^{j}+U_{i-1}^{j}+U_{i}^{j+1}+[ U_i^{j+1} + 2 k(x_i,y_j)h]-4U_{i}^{j}) = f_i^j$$ Gathering the terms: $$j=0,\quad -\frac{1}{h^2}(U_{i+1}^{j}+U_{i-1}^{j}+2U_{i}^{j+1}+ 2 k(x_i,y_j)h-4U_{i}^{j}) = f_i^j$$ You see that the term $U_{i}^{j+1}$ appears twice indeed. Now in your case, I suspect $k=0$ so the equation at the boundary is: $$j=0,\quad -\frac{1}{h^2}(U_{i+1}^{j}+U_{i-1}^{j}+2U_{i}^{j+1}-4U_{i}^{j}) = f_i^j$$

• The formula above is close to the formula you give. Yet, there is no reason for $U_{i+1}^{j}+U_{i-1}^{j}$ to be equal to zero. The only case where this can happen is if these two points are on a boundary where a Dirichlet value of 0 is specified. In this case your domain is just made of three nodes in the $x$ direction, counting the boundary nodes.