# Using Neumann boundary conditions

I had a semi-related physics problem I needed to solve analytically (which I have already done), but I am now curious how I would go about numerically solving the entire system in Mathematica. Consider a conductive rectangular plate - for simplicity, say it goes from $-1\le x\le 1$ and $-1\le y\le 1$, with wires connected at $(-1,0)$ and $(1,0)$. Some total current $I$ flows through the plate through these connections. I want to be able to numerically solve for the potential anywhere in the plate.

From symmetry, we can assume that $V(0,y)=0$, and the current flowing through the wires establish two Neumann boundary conditions along these points. How would I go about solving this system in Mathematica?

My current attempt at solving this is

soln = NDSolveValue[{Laplacian[u[x, y], {x, y}] ==
NeumannValue[-1, y == -1 && x == 0] +
NeumannValue[1, y == 1 && x == 0], u[0, y] == 0},
u, {x, y} \[Element]
Rectangle[{-1, -1}, {1, 1}]];
Plot3D[soln[x, y], {x, y} \[Element] Rectangle[{-1, -1}, {1, 1}],
AxesLabel -> {"x", "y", "V"}]


However, when I run this it gives me two identical error messages corresponding to each use of NeumannValue:

NDSolveValue::bcnop: No places were found on the boundary where Coordinate was True, so NDSolveFEMBoundaryCondition[{Neumann,{1,1},{CompiledFunction[{10,10.4,5568},{_Real,_Real},{{3,0,0},{3,0,1},{3,2,0}},{{{{-1.}},{3,2,0}}},{0,0,2,0,1},{{1}},Function[{x,y},{{-1.}},Listable],Evaluate]}},Coordinate,CompiledFunction[<<1>>],NeumannValue[-1,y==-1&&x==0]] will effectively be ignored. >>

Playing around with this, I've noticed Mathematica doesn't like point-like Neumann conditions (if I set the conditions such that the conditions are (incorrectly) NeumannValue[-1, y == -1] then it at least solves the equations, but then it also doesn't give that u[0,y]==0 over all y...). Does anyone have any suggestions on what might be causing the problem or a better way to do this in Mathematica?

• DirichletConditions apply at points, NeumannValues at edges(2D)/faces(3D); also only having NeumannValues for a stationary PDE will result in non-unique solutions. – user21 Apr 22 '16 at 4:00

## 2 Answers

A point has measure 0 in any numeric approximation possible, so it's not only mathematica who doesn't like point-like conditions.

soln = NDSolveValue[{Laplacian[u[x, y], {x, y}] ==
NeumannValue[-1, y == -1. && Abs[x] < 0.1] +
NeumannValue[1, y == 1 && Abs[x] < 0.1], u[0, y] == 0},
u, {x, y} \[Element] Rectangle[{-1., -1.}, {1., 1.}]]
Plot3D[soln[x, y], {x, y} \[Element] Rectangle[{-1, -1}, {1, 1}],
AxesLabel -> {"x", "y", "V"}]


gives you something, but smaller values for the boundary condition don't. I think it has to do with the mesh FEM is using (when you use region fem is called automatically, I think). So you need to either define the mesh yourself of create a larger boundary condition.

• Ah thank you, this lead me to using Method -> {"FiniteElement", "MeshOptions" -> {MaxCellMeasure -> dA}} in NDSolveValue (or some small value such that the "width" of the connection was less than the order of $\sqrt{dA}$. – Ben Bartlett Apr 22 '16 at 6:43

An alternative is to use a numerical approximation to DiracDelta:

dirac[x_, a_] := Exp[-(x/a)^2]/a/Sqrt[Pi]
soln = NDSolveValue[{Laplacian[u[x, y], {x, y}] ==
NeumannValue[-dirac[x, 0.01], y == -1.] +
NeumannValue[dirac[x, 0.01], y == 1], u[0, y] == 0},
u, {x, y} \[Element] Rectangle[{-1., -1.}, {1., 1.}],
Method -> {"FiniteElement", "MeshOptions" -> {MaxCellMeasure -> dA}}]
Plot3D[soln[x, y], {x, y} \[Element] Rectangle[{-1, -1}, {1, 1}],
AxesLabel -> {"x", "y", "V"}]
StreamPlot[{-D[soln[x, y], x], -D[soln[x, y], y]}, {x, y} \[Element]
Rectangle[{-1, -1}, {1, 1}], AspectRatio -> Automatic,
Evaluated -> True]