NeumannBC1pt.nb
Neumann Boundary conditions imposed with 1 boundary point
John Neuberger calls this the Edge Case

This is a GOOD notebook, for handling one neighboring boundary point.  (in the first quadrant).

This is for John Neuberger’s problem of finite difference approximation to the Laplacian on arbitrary regions in the plane with zero Neumann boundary conditions.  

We have a square grid with spacing h = 1.  The center vertex is surrounded by N, S, and W, vertices (North, South, East) at points (0,0), (0,1), (0, -1) and (-1, 0)  The values there are u0, uN, uS, and uW.
The right point (a, 0) is on the boundary, and there is an outward normal vector n = (1, b).  A vector tangent to the boundary is t = (-b, 1).

Note:  For grid spacing h > 0, the vertices are at (0, 0), (0, h), etc.
The right point is at (ah, 0) and the normal vector and tangent vector are the same.  The Laplacian is the “h = 1” Laplacian divided by h^2.

Version Notes:  The folder snow/mathematica/Local Laplacian at Boundary/2016-12-17 has a version of this where I use more points in the interior (namely the point at (-h, -h) giving uSW).  In this file
I stick to the 4 points from the usual 5 that are in the region.

Compute zeros of Bessel Prime.  They used to be built-in commands in the good-old-days (pre-version 6).

besselJprime[n_, z_] := (BesselJ[n-1, z] - BesselJ[n+1, z])/2;
zp01 = z/.FindRoot[besselJprime[0,z] == 0, {z, 3.8}]
zp11 = z/.FindRoot[besselJprime[1,z] == 0, {z, 2}]

(*
pJ0a =Plot3D[BesselJ[0, zp01 Sqrt[x^2 + y^2]], {x, -1,1}, {y, -1,1}, Exclusions -> {x^2 + y^2 == 1}];
pJ0b = ContourPlot[BesselJ[0, zp01 Sqrt[x^2 + y^2]], {x, -1,1}, {y, -1,1}, Exclusions -> {x^2 + y^2 == 1}];
GraphicsRow[{pJ0a, pJ0b}, ImageSize -> 600];
pJ1a = Plot3D[BesselJ[1, zp11 Sqrt[x^2 + y^2]]Cos[ArcTan[x,y]], {x, -1,1}, {y, -1,1}, Exclusions -> {x^2 + y^2 == 1}];
pJ1b = ContourPlot[BesselJ[1, zp11 Sqrt[x^2 + y^2]]Cos[ArcTan[x,y]], {x, -1,1}, {y, -1,1}, Exclusions -> {x^2 + y^2 == 1}];
GraphicsRow[{pJ0a, pJ0b,pJ1a, pJ1b}, ImageSize -> 1000]
*)

Saved figures.

This next cell has what we need to define, with the exception of the h value.

Clear[c00, c10, c01, c20, c02, c11, u0, uE, uN, uW, uS, uSW, uNW, x, y,a, b];
fNum = 0; (*0 or 1 for psi01 or psi11 *)
eNum = 0;  (* 0 or 1: eNum = -0 for c11 = 0, eNum = 1 for rotation of axes of the quadratic form so eigenvectors are normal (n) and tangent (t) vectors *)
aVal = 0.9; (* How far the boundary point is from inside grid point *)
yb = 0.5;
hFigs = 0.01;
(*All parameters are above here *)
Which[fNum == 0, bess[x_,  y_] := BesselJ[0, zp01 Sqrt[x^2 + y^2]]; eval = zp01^2; fString = "psi_01",
fNum ==1, bess[x_,  y_] := BesselJ[1, zp11 Sqrt[x^2 + y^2]]Cos[ArcTan[x,y]]; eval = zp11^2; fString = "psi_11"];
f[x_, y_] = c00  + c10 x + c01 y + c20 x^2 + c02 y^2 + c11 x y;
Which[
eNum == 0, eString = " Df(a, 0).(1,b) = 0, c11 = 0";
eqns = {f[0,0] == u0, f[-1,0] == uW, f[0,1] == uN, f[0, -1] == uS, c10 + 2 a c20 +b c01 + a b c11== 0, c11 == 0},

eNum == 1, eString = "Df(a,0).n = 0, D^2(n, t) = 0";
eqns = {f[0,0] == u0, f[-1,0] == uW, f[0,1] == uN, f[0, -1] == uS, c10 + 2 a c20 +b c01 + a b c11== 0, 2 b(c02 - c20) + (1-b^2) c11 == 0}
];
soln = Part[Solve[eqns, {c00, c10, c01, c20, c02, c11}], 1] ;(* There is only one solution, so take the first and only part *)
laplacianWithParams = Collect[Simplify[2(c02 + c20)/.soln], {u0, uN, uS, uW, uSW, uNW}];
Print["Using equations with ", eString, ", the formula for the Laplacian is \n", laplacianWithParams];
a  = aVal;  
xb = Sqrt[1-yb^2];  (* (xb, yb) is a boundary point *)
b = yb/xb;
soln = Part[Solve[eqns, {c00, c10, c01, c20, c02, c11}], 1] ;(* There is only one solution, so take the first and only part *)
laplacian = Simplify[2(c02 + c20)/.soln];
Print["With a = ", a, " and yb = ", yb, ", the boundary point is ", {xb, yb}, " and the normal vector is ",{1,b}];
h = hFigs;
Print["Figures of eigenfunction ", fString, ", which has eigenvalue ", eval, ", for h = ", h];
points = {{xb, yb}, {xb - a h, yb}, {xb - a h - h, yb}, {xb- a h, yb + h}, {xb - a h, yb - h}(*, {xb - a h - h, yb-h}*)};
line = Table[{Sqrt[1 - y^2], y}, {y, yb - h, yb + h, h/5}];
background = Graphics[{PointSize[Large],Point[points], Line[line]}, AspectRatio -> Automatic];
uVals  = {u0 -> bess[xb-a h, yb], uW -> bess[xb - a h - h, yb], uN -> bess[xb - a h, yb + h], uS -> bess[xb - a h, yb - h], uSW -> bess[xb -a h - h, yb - h], uNW -> bess[xb -a h - h, yb + h]};
fit[x_, y_]= f[(x-xb +a h)/h, (y - yb)/h]/.soln/.uVals;
linFit[x_, y_] := (u0 + (c10/h)(x-(xb -a h)) + (c01/h)(y - yb))/.soln/.uVals;
cBessel =ContourPlot[bess[x,y], {x, xb- a h - h, xb - a h + h}, {y, yb-h, yb+h}];
cQuad1 = ContourPlot[bess[x, y]-linFit[x,y], {x, xb- a h - h, xb - a h + h}, {y, yb-h, yb+h}];
cQuad2 = ContourPlot[fit[x, y]-linFit[x,y], {x, xb- a h - h, xb - a h + h}, {y, yb-h, yb+h}];
cDiff =ContourPlot[fit[x,y]-bess[x,y], {x, xb- a h - h, xb - a h + h}, {y, yb-h, yb+h}(*, Contours -> Table[10.^-8 i, {i, -10, 10}]*)];
p1 = Show[cBessel, background, PlotLabel -> "Bessel Function near boundary"];
p2 = Show[cQuad1, background, PlotLabel -> "Bessel Function - linear fit"];
p3 = Show[cQuad2, background, PlotLabel -> "Quadratic fit - linear fit"];
p4 = Show[cDiff, background, PlotLabel ->"Quadratic fit - Bessel function"];
GraphicsRow[{p1, p2, p3, p4}, ImageSize -> 1000]
hExponentMax = 15;
ubVal = bess[xb, yb];
data = {};
Do[h = 1/2^hExponent;
uVals  = {u0 -> bess[xb-a h, yb], uW -> bess[xb - a h - h, yb], uN -> bess[xb - a h, yb + h], uS -> bess[xb - a h, yb - h], uSW -> bess[xb -a h - h, yb - h], uNW -> bess[xb -a h - h, yb + h]};
fit[x_, y_]= f[(x-xb +a h)/h, (y - yb)/h]/.soln/.uVals;
evalEst = -(laplacian/(h^2 u0))/.uVals;
evalError = evalEst - eval;
ubEst = fit[xb, yb];
ubError = ubEst - ubVal;
AppendTo[data, {1/h, evalEst, evalError, evalError/h^1, evalError/h^2, ubError, ubError/h, ubError/h^2, ubError/h^3}];
, {hExponent, 2, hExponentMax}];
TableForm[data, TableHeadings->{None, {"1/h", "evalEst", "evalError", "evalError/h^1", "evalError/h^2", "ubError", "ubError/h^1", "ubError/h^2", "ubError/h^3"}}]

Below are saved results.

Case 1: fNum = 0, eNum = 0
psi_01 and c11 = 0
Local eval estimate is O(1), that is, it converges but to wrong limit
u on boundary estimate is O(h^2)

Case 2: fNum = 0, eNum = 1
psi_01 and rotated axes
Local eval estimate is O(h)
u on boundary estimate is O(h^3)

Case 3: fNum = 1, eNum = 0
psi_11 and c11 = 0
Local eval estimate is O(1)
u on boundary estimate is O(h^2)

Case 4: fNum = 1, eNum = 1
psi_11 and rotated axes
Local eval estimate is O(1), that is estimate converges, but to wrong number
u on boundary estimate is O(h^2)

Here is a summary of the 4 conclusions:

Case 1: fNum = 0, eNum = 0
psi_01 and c11 = 0
Local eval estimate is O(1), that is, it converges but to wrong limit
u on boundary estimate is O(h^2)

Case 2: fNum = 0, eNum = 1
psi_01 and rotated axes
Local eval estimate is O(h)
u on boundary estimate is O(h^3)

Case 3: fNum = 1, eNum = 0
psi_11 and c11 = 0
Local eval estimate is O(1)
u on boundary estimate is O(h^2)

Case 4: fNum = 1, eNum = 1
psi_11 and rotated axes
Local eval estimate is O(1), that is estimate converges, but to wrong number
u on boundary estimate is O(h^2)