October 25, 2018 version

This is for John Neuberger’s problem of finite difference approximation to the Laplacian with zero Neumann boundary conditions.  This is for the case were there are two boundary points,
the “East” point at (aE h, 0) and the “North” point at (0, aN h).  We are estimating the Laplacian at (0, 0), called the Local Laplacian, and h is the spacing of the square grid.

u0 = u(0,0)
uS = u(0, -h)
uW = u(-h, 0)

The region is always the unit disk in this notebook.  Given aE, aN and h in (0, 1),  there is a unique point (xc, yc) in the unit disk such that (xc + aE h, yc) and (xc, yc + aN h) are on the unit circle.  The function findx0y0 below finds this center point with the notation
{xc, yc} = findx0y0[aE h, aN h]
========================================================

Note:  When writing the equations, we have a square grid with spacing h = 1.  The center vertex at (0, 0) is surrounded by S and W vertices (South and West).  The values there are u0, uS, and uW.
The East point (aE h, 0) and North point (0, aN h) are on the boundary, and there is an approximate outward normal vector n = (aN, aW) and CCW tangent vector t = (-aW, aN).  

The estimated point on the boundary is the closest point on the line joining  h(aE, 0) and h(0, aN), namely (x*, y*) = h(aN, aW)*(aN aE)/(aN^2 + aW^2).  (This needs proof!)

The Laplacian is the “h = 1” Laplacian divided by h^2.

Compute zeros of Bessel Prime.  Built-in commands were defined 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

Run this cell first to define the function.  usage is {x0, y0} = findx0y0[h aE, h aN];

findx0y0[aE_, aN_]:= Module[{x0, y0, aEE, aNN, soln},  (* Returns grid Poin (x0, y0) with 2 nbrs oustide unit disk *)
soln = Part[Simplify@Solve[{(x0 + aEE)^2 + y0^2 == 1, x0^2 + (y0 + aNN)^2 == 1}, {x0, y0}], 2];  (* The first solution gives a point in Quadrant III outside the circle, the second solution gives a point in Quadrant I inside the circle *)
{x0, y0}/.soln/.{aEE -> aE, aNN -> aN} ]

This next cell runs everything.  If it does not work it may be that the cell above, that computes the zeros of Jprime, has not been run.

Note: Currently I have eNum = 5 defined, but it only has 5 equations, not 6.  there are some commented out attempts to get an equation with arbitrary c11.  Perhaps I can try a pseudo-inverse.  Q: Do I really have to do the Solve twice?  Once before h is set and once after it is set?

run := Module[{},
Print["fNum = ", fNum, ", eNum = ", eNum, ", aN = ", aNval, ", aE = ", aEval, ", and hFigs = ", hFigs];
Clear[c00, c10, c01, c20, c02, c11, u0, uE, uN, uW, uS, uSW, uNW, x, y,aN, aE];
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 ==1, eString = "Df(aE,0).n = 0, Df(0, aN).n = 0, c11 = 0";
eqns = {f[0,0] == u0, f[-1,0] == uW, f[0, -1] == uS, aN c10 +2 aE aN c20 +aE c01 + aE^2 c11== 0, aE c01 + 2 aE aN c02 +aN c10 + aN^2 c11== 0, c11 ==  0},

eNum ==2, eString = "Df(aE, 0).xhat = 0, Df(0, aN).yhat = 0, c11 = 0";
eqns = {f[0,0] == u0, f[-1,0] == uW, f[0, -1] == uS, c10 + 2 c20 aE == 0, c01 + 2 c02 aN == 0,c11 == 0},

eNum ==3, eString = "Df(aE, 0).n = 0, Df(0, aN).n = 0, f''(n,t) = 0";
eqns = {f[0,0] == u0, f[-1,0] == uW, f[0, -1] == uS, aN c10 +2 aE aN c20 +aE c01 + aE^2 c11== 0, aE c01 + 2 aE aN c02 +aN c10 + aN^2 c11== 0, 2 aE aN (c02-c20) + c11(aN^2-aE^2) == 0},

eNum ==4, eString = "Df(x*, y*).n = 0, D2f(n,t) = 0, c11 = 0"; 
eqns = {f[0,0] == u0, f[-1,0] == uW, f[0, -1] == uS, aN c10 +aE c01 + (aN aE)/(aN^2 + aE^2)(2 c11 aN aE + 2(c20 aN^2 + c02 aE^2)) == 0, 2 aE aN (c02-c20) + c11(aN^2-aE^2) == 0, c11 ==0},

eNum ==5, eString = "Df(x*, y*).n = 0, D2f(n,t) = 0, and fitting uSW"; 
eqns = {f[0,0] == u0, f[-1,0] == uW, f[0, -1] == uS, aN c10 +aE c01 + (aN aE)/(aN^2 + aE^2)(2 c11 aN aE + 2(c20 aN^2 + c02 aE^2)) == 0, 2 aE aN (c02-c20) + c11(aN^2-aE^2) == 0, f[-1, -1] == uSW},

eNum ==6, eString = "Df(aE, 0).n = 0, D2f(n,t) = 0, and fitting uSW"; 
eqns = {f[0,0] == u0, f[-1,0] == uW, f[0, -1] == uS, aN c10 +2 aE aN c20 +aE c01 + aE^2 c11== 0, 2 aE aN (c02-c20) + c11(aN^2-aE^2) == 0, f[-1, -1] == uSW},

eNum ==7, eString = "Df(aE, 0).n = 0, D2f(n,t) = 0, and fitting uNW"; 
eqns = {f[0,0] == u0, f[-1,0] == uW, f[0, -1] == uS, aN c10 +2 aE aN c20 +aE c01 + aE^2 c11== 0, 2 aE aN (c02-c20) + c11(aN^2-aE^2) == 0, f[-1, 1] == uNW},

eNum == 8, eString = "D2f(n, t) = 0, Df(aE, 0).(1, bE) = 0, Df(0, aN).(bN, 1) = 0";
eqns = {f[0,0] == u0, f[-1,0] == uW, f[0, -1] == uS, 2 aE aN (c02-c20) + c11(aN^2-aE^2) == 0, (c10 + 2 c20 aE) + bE(c01 + c11 aE) == 0, bN(c10+c11 aN ) + (c01+ 2 c02 aN) == 0},

True, Print["value of eNum is not defined"]
];
soln = Part[Solve[eqns, {c00, c10, c01, c20, c02, c11}], 1]; (* There is only one solution, so take the first and only part *)
(*soln = Part[Solve[eqns, {c00, c10, c01, c20, c02}], 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, uSW2, uW2, uS2}];
Print["Equation number eNum = ", eNum, " uses ", eString];Print["Laplacian h^2 = \n", laplacianWithParams];

soln = Part[Solve[eqns, {c00, c10, c01, c20, c02, c11}], 1] ;
(*soln = Part[Solve[eqns, {c00, c10, c01, c20, c02}], 1] ;*)(* There is only one solution, so take the first and only part *)
laplacian = Simplify[2(c02 + c20)/.soln];
aE = aEval;
aN = aNval;
(*Print["With aE = ", aEval, ", aN = ", aNval", the formula for the Laplacian is ", Simplify[2(c02 + c20)/.soln]];
*)
h = hFigs;
{xc, yc} = findx0y0[h aE, h aN];  (* The grid point in disk with 2 neighbors outside disk *)
{xb, yb} = {xc, yc}/Norm[{xc, yc}]; (* the point on boundary closest to {xc, yc} *)
points = {{xc, yc}, {xc - h, yc}, {xc, yc - h},{xc + h aE, yc}, {xc, yc + h aN}(*, {xc + h, yc}, {xc , yc+ h} *),
{xb, yb}};
line = Table[{Sqrt[1 - y^2], y}, {y, yc - h, yc + h, h/5}];
background = Graphics[{PointSize[Large],Point[points], Line[line]}, AspectRatio -> Automatic];

uVals  = {u0 -> bess[xc, yc], uW -> bess[xc - h, yc], uN -> bess[xc, yc + h], uS -> bess[xc, yc - h], uSW -> bess[xc - h, yc - h], uSW2 -> bess[xc - 2h, yc - 2h], uW2 -> bess[xc - 2h, yc ], uS2 -> bess[xc , yc -2h ], uNW -> bess[xc - h, yc + h ], xE ->xc + h aE, yE -> yc, xN -> xc, yN -> yc + h aN, bN -> xc/(yc+h aN), bE -> yc/(xc + h aE)};
fit[x_, y_]= f[(x-xc)/h, (y - yc)/h]/.soln/.uVals;
linFit[x_, y_] := (u0 + (c10/h)(x-xc) + (c01/h)(y - yc))/.soln/.uVals;
bp = ContourPlot[bess[x, y], {x, xc- h, xc + h}, {y, yc-h, yc+h}, PlotLabel -> "Bessel Function near boundary"];
cp = ContourPlot[bess[x, y]-linFit[x,y], {x, xc- h, xc + h}, {y, yc-h, yc+h}, PlotLabel -> "Bessel Function - linear fit"];
fp = ContourPlot[fit[x, y], {x, xc- h, xc + h}, {y, yc-h, yc+h}, PlotLabel -> "Quadratic fit"];
flp = ContourPlot[fit[x, y]-linFit[x,y], {x, xc- h, xc + h}, {y, yc-h, yc+h}, PlotLabel -> "Quadratic fit - linear fit"];
cd =ContourPlot[fit[x,y]-bess[x,y],{x, xc- h, xc + h}, {y, yc-h, yc+h}, PlotLabel ->"Quadratic fit - Bessel function"];
bp = Show[bp, background];
cp = Show[cp, background];
fp = Show[fp, background];
flp = Show[flp, background];
cd = Show[cd, background];
(*GraphicsRow[{bp,cp,fp, cd}, ImageSize -> 1000];*)
fig = GraphicsRow[{bp, fp,cd}, ImageSize ->750];
Print["Figures are of eigenfunction ", fString, " (fNum = ", fNum,"), which has eigenvalue ", eval, ", for h = ", h];
hExponentMax = 15;
(*ubVal = bess[xb, yb];*)
data = {};
Do[h = 1/2^hExponent;
{xc, yc} = findx0y0[h aE, h aN];
{xb, yb} = {xc, yc}/Norm[{xc, yc}];
ubVal = bess[xb, yb];
uVals  = {u0 -> bess[xc, yc], uW -> bess[xc - h, yc], uN -> bess[xc, yc + h], uS -> bess[xc, yc - h], uSW -> bess[xc - h, yc - h], uSW2 -> bess[xc - 2h, yc - 2h], uW2 -> bess[xc - 2h, yc ], uS2 -> bess[xc , yc -2h ], uNW -> bess[xc - h, yc + h ], xE ->xc + h aE, yE -> yc, xN -> xc, yN -> yc + h aN, bN -> xc/(yc+h aN), bE -> yc/(xc + h aE)};
fit[x_, y_]= f[(x-xc)/h, (y - yc)/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, ubError, ubError/h, ubError/h^2, ubError/h^3}];
, {hExponent, 2, hExponentMax}];
Print@TableForm[data, TableHeadings->{None, {"1/h", "evalEst", "evalError", "evalError/h",  "ubError", "ubError/h", "ubError/h^2", "ubError/h^3"}}];
Show[fig]
]

eNum =2; (* Equation number *)
fNum = 0; (* function number: 0 (psi 01) or 1 (psi11) *)
aEval = .5;aNval = .9;hFigs = 0.125;  (* value of h for the figures *)
(*All parameters are above here *)
run

eNum = 1, 
Equation number eNum = 1 uses Df(aE,0).n = 0, Df(0, aN).n = 0, c11 = 0
Laplacian*h^2 = -((4 (aE+aN) u0)/(aE+aN+2 aE aN))+((4 aE uS)/(aE+aN+2 aE aN))+((4 aN uW)/(aE+aN+2 aE aN))
psi_11: evalError = 0.57, ubError = 0.014 h^2
psi_01: evalError = 7.1, ubError = 0.15 h^2

eNum = 2,    This is the Laplacian built into MATLAB. 
Equation number eNum = 2 uses Df(aE, 0).xhat = 0, Df(0, aN).yhat = 0, c11 = 0
Laplacian*h^2 = (-(2/(1+2 aE))-(2/(1+2 aN))) u0+((2 uS)/(1+2 aN))+((2 uW)/(1+2 aE))
psi_01: evalError = 2.5 h, ubError = -0.20 h^2
psi_11: evalError -> -infinity, ubError = -0.014 h

eNum = 3,
Equation number eNum = 3 uses Df(aE, 0).n = 0, Df(0, aN).n = 0, f’’(n,t) = 0
These three equations are linearly dependent, so this does not give a solution

eNum = 4,
Equation number eNum = 4 uses Df(x*, y*).n = 0, f’’(n,t) = 0, c11 = 0
Laplacian*h^2 = -((4 (aE+aN) u0)/(aE+aN+2 aE aN))+((4 aE uS)/(aE+aN+2 aE aN))+((4 aN uW)/(aE+aN+2 aE aN))
Same as eNum = 1

eNum = 5,
Equation number eNum = 5 uses Df(x*, y*).n = 0, D2f(n,t) = 0, and fitting uSW
Laplacian*h^2 = (((-3 (aE^2) aN-(aN^3)-(aE^3) (1+2 aN)-aE (aN^2) (3+2 aN)) u0)/(aE aN (aE+aN+2 aE aN)))+(((3 (aE^2) aN-aE (aN^2)+(aN^3)+2 aE (aN^3)+(aE^3) (1+2 aN)) uS)/(aE aN (aE+aN+2 aE aN)))+((((aE^2) aN+aE (aN^2)-(aN^3)-2 aE (aN^3)-(aE^3) (1+2 aN)) uSW)/(aE aN (aE+aN+2 aE aN)))+(((-(aE^2) aN+3 aE (aN^2)+(aN^3)+2 aE (aN^3)+(aE^3) (1+2 aN)) uW)/(aE aN (aE+aN+2 aE aN)))
psi_01: evalError = -6.3 h, ubError = -0.016 h^3
psi_11: evalError = 0.23, ubError = 0.022h^2

eNum = 6,
Equation number eNum = 6 uses Df(aE, 0).n = 0, D2f(n,t) = 0, and fitting uSW
Same Laplacian as eNum5

eNum = 7
Equation number eNum = 7 uses Df(aE, 0).n = 0, D2f(n,t) = 0, and fitting uNW
Laplacian*h^2 = 
((((aE^3)-(aE^2) aN-2 (aE^3) aN-3 (aN^3)-aE (aN^2) (5+2 aN)) u0)/((1+2 aE) (aN^2) (aE+aN)))+((((aE^3)-(aE^2) aN+2 (aE^3) aN-aE (aN^2)+(aN^3)+2 aE (aN^3)) uNW)/((1+2 aE) (aN^2) (aE+aN)))+(((-(aE^3)+3 (aE^2) aN+2 (aE^3) aN+aE (aN^2)+(aN^3)+2 aE (aN^3)) uS)/((1+2 aE) (aN^2) (aE+aN)))+(((-(aE^3)-(aE^2) aN-2 (aE^3) aN+5 aE (aN^2)+(aN^3)-2 aE (aN^3)) uW)/((1+2 aE) (aN^2) (aE+aN)))
psi_01: evalError = -5.2 h, uBerror = 0.0039 h^3
psi_11: evalError = -0.17, ubError = 0.033 h^2

eNum = 8
Equation number eNum = 8 uses D2f(n, t) = 0, Df(aE, 0).(1, bE) = 0, Df(0, aN).(bN, 1) = 0
Laplacian*h^2 =
-((2 (2 (aE^3)+2 (aE^3) bN-2 (aN^2) (1+aN+aN bE-bE bN)+2 aE (aN^2) (-1+bN+2 bE bN)-2 (aE^2) (-1+bE bN+aN (-1+bE+2 bE bN))) u0)/((aE^3) (2+4 aN)+(aN^2) (-1-2 aN+bE bN)-2 aE (aN^2) (1+2 aN-(1+bE) bN)-(aE^2) (-1+4 (aN^2) (bE-bN)+bE bN+2 aN (-1+bE+bE bN))))-((2 (-(aE^2)-2 (aE^3)+(aN^2)+2 aE (aN^2)+2 (aE^2) aN bE+2 (aN^3) bE+(aE^2) bE bN-(aN^2) bE bN-4 aE (aN^2) bE bN) uS)/((aE^3) (2+4 aN)+(aN^2) (-1-2 aN+bE bN)-2 aE (aN^2) (1+2 aN-(1+bE) bN)-(aE^2) (-1+4 (aN^2) (bE-bN)+bE bN+2 aN (-1+bE+bE bN))))-((2 (-(aE^2)-2 (aE^2) aN+(aN^2)+2 (aN^3)-2 (aE^3) bN-2 aE (aN^2) bN+(aE^2) bE bN+4 (aE^2) aN bE bN-(aN^2) bE bN) uW)/((aE^3) (2+4 aN)+(aN^2) (-1-2 aN+bE bN)-2 aE (aN^2) (1+2 aN-(1+bE) bN)-(aE^2) (-1+4 (aN^2) (bE-bN)+bE bN+2 aN (-1+bE+bE bN))))
psi_01: evalError = -5.6 h, ubError = O(h^3)?  Roundoff error is large
psi_11: evalError diverges, ubError = 0.016 h