Using Neumann boundary conditions

Question

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?

Answer 1

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.

Answer 2

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]