The solution of the 2D Poisson equation:
$$
\frac{\partial^2 V}{\partial x^2}+\frac{\partial^2 V}{\partial y^2}=\delta(x-x_s)\delta(y-y_s) \tag{*}
$$
subjected to the boundary condition $V(\infty)=0$ is:
$$
G(x,y)=\frac{1}{2\pi}\ln|r-r_s|
$$
Question: What is the solution of (*) subjected to the following boundary conditions?
$$
V(x\to\pm\infty)=0\\
V(y\to-\infty)=0\\
\frac{\partial V}{\partial y}\left|_{y=0}\right.=0
$$
[Math] Green’s function of 2D Poisson equation in the lower half-plane with Neumann boundary
greens functionpartial differential equationspoisson's equation
Related Solutions
Finding the general solution of the PDE (below) is not the more difficult part of the task :
As usual, the arduous part is to determine, among the infinite number of solutions which one fits with the boundary conditions :
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.
Best Answer
I started working to this answer after seeing the comment of Dylan and the following comment of A Slow Learner, therefeore I was a little misplaced by the answer by Dylan himself (which however focuses more on the image method for constructing Green's functions): however I think it could be useful to see how the answer is obtained, and so I've decided to provide a(n almost) complete development below.
The Green's function $\mathscr{G}_N(x-x_s,y-y_s)$ for Neumann's problem for the Poisson equation (*) in the lower half plane is the solution of the following Dirichlet boundary value problem: $$ \begin{cases} \Delta \mathscr{G}_N(x-x_s,y-y_s)=\delta(x-x_s,y-y_s),\quad (x,y),(x_s,y_s)\in \mathbb{R}\times\mathbb{R}^-\\ \frac{\partial\mathscr{G}_N}{\partial y}(x-x_s,-y_s)=0 \end{cases}.\tag{1}\label{1} $$ Note that $\frac{\partial\mathscr{G}_N}{\partial\mathbf{n}}=\frac{\partial\mathscr{G}}{\partial y}$ on the boundary $\{y=0\}$ of the lower half plane.
According to Vladimirov ((1983) ยง29.1, p. 368)), the general solution of problem \eqref{1} has the following form $$ \mathscr{G}_N(x-x_s,y-y_s)=\frac{1}{2\pi}\ln|(x,y)-(x_s,y_s)|+g(x-x_s,y-y_s) $$ where the first term on the right side is the fundamental solution of the laplacian while the second one is a function harmonic in the whole $\mathbb{R}\times\mathbb{R}^-$. In the problem posed by the OP, $$ \mathscr{G}_N(x-x_s,y-y_s)=\frac{1}{2\pi}\big(\ln|(x,y)-(x_s,y_s)|+\ln|(x,y)-(x_s,-y_s)|\big).\tag{2}\label{2} $$ Since $\frac{\partial\mathscr{G}_N}{\partial\mathbf{n}}=\frac{\partial\mathscr{G}_N}{\partial y}$ for problem \eqref{1}, by nothing that $$ \begin{split} \frac{\partial}{\partial y}\frac{1}{2\pi}\ln|(x,y)-(x_s,y_s)|& =\frac{\partial}{\partial y}\frac{1}{2\pi}\ln\sqrt{(x-x_s)^2+(y-y_s)^2}\\ &=\frac{1}{\pi}\frac{y-y_s}{(x-x_s)^2+(y-y_s)^2}, \end{split} $$ it is easy to verify that the function \eqref{2} solves problem \eqref{1} and , by Green's formula it solves the general Neumann problem for the 2D Poisson's equation on the half plane.
[1] Vladimirov, V. S. (1983)[1970], Equations of mathematical physics, Moscow: Mir Publishers, 2nd ed., pp. 464, MR0764399, Zbl 0207.09101 (the Zbmath review refers to the first English edition).