The problem I would like to solve involves diffusion in the following region

reg = ImplicitRegion[-5 <= x <= 5 && -5 <= y <= 5 && x^2 + y^2 >= 1^2, {x, y}];
RegionPlot[reg, AspectRatio -> Automatic]

Region over which to solve PDE

Additionally, there is an interaction potential $v$ acting in this region:

v = Function[{x, y}, -E^(-x^2 - y^2)];
Plot3D[v[x, y], {x, -5, 5}, {y, -5, 5}, RegionFunction -> Function[{x, y, z}, x^2 + y^2 > 1^2], PlotRange -> {-1, 1}, PlotPoints -> 50]

plot of interaction potential

Physically, I would like to model 2D diffusion confined to this region, in the presence of the biasing potential given by $v$. So what I need to do is solve the PDE $\vec{\nabla} \cdot \left[\vec{\nabla}u(x,y) + u(x,y) \vec{\nabla}v(x,y)\right] = 0$ subject to the boundary conditions

$u(x,-5) = 1$

$u(x,5) =0$

to impose an overall flux, and with zero-flux boundary conditions at the boundary where $x^2 + y^2 = 1$, (as well as on the sides where $x = \pm 5$). So I proceed as follows in Mathematica (10.0.0):

op = Div[Grad[u[x, y], {x, y}] + u[x, y] Grad[v[x, y], {x, y}],{x, y}];
bc = {DirichletCondition[u[x, y] == 1, y == -5],DirichletCondition[u[x, y] == 0, y == 5]};

Note that I did not specify boundary conditions for the circular boundary, nor for the sides $x = \pm 5$, because according to the NDSolve/NeumannValue documentation, these should default to zero-flux boundary conditions.

sol = NDSolve[{op == 0, bc}, u, {x, y} \[Element] reg, Method -> {"FiniteElement","MeshOptions" -> {MaxCellMeasure -> 0.005}}]

The solution looks like this:

Plot3D[Evaluate[u[x, y] /. sol], {x, y} \[Element] reg, PlotRange -> All, AxesLabel -> {"x", "y"}]

Solution of the PDE

Now, this is not quite what I expected, for the following reason:if the PDE is solved correctly, with zero-flux boundary conditions, then the flux into the box (at y = -5) should equal the flux out of the box (at y = 5). You can sort of tell that this is not the case from the plot above, but plotting the normal derivative at the boundaries $y = \pm 5$ makes it clear:

Plot[{Evaluate[(Derivative[0, 1][u][x, 5] + u[x, 5] Derivative[0, 1][v][x, 5]) /.sol],
Evaluate[(Derivative[0, 1][u][x, -5] + u[x, -5] Derivative[0, 1][v][x, -5]) /. sol]}, {x, -5, 5}, PlotRange -> {-0.2, 0}]

normal derivatives

Integrating these will give the flux in and the flux out. Clearly the flux in is not equal to the flux out. Since I have checked to make sure that the PDE $\vec{\nabla} \cdot \left[\vec{\nabla}u(x,y) + u(x,y) \vec{\nabla}v(x,y)\right] = 0$ is obeyed in the interior of the region, it seems that I am not correctly setting the Neumann boundary condition at the circular boundary $x^2 + y^2 = 1$. The boundary condition that I want is

$\left[\vec{\nabla}u(x,y) + u(x,y) \vec{\nabla}v(x,y)\right] \cdot \vec{n} = 0$

but I believe that the boundary condition that I'm getting is simply

$\left[\vec{\nabla}u(x,y)\right] \cdot \vec{n} = 0$

How can I impose the correct boundary condition? Have I misunderstood the documentation, which seems to say that my desired boundary condition is the default?

  • $\begingroup$ I think I have an idea what the issue might be, could you try this op = Div[Grad[u[x, y], {x, y}] + Grad[v[x, y], {x, y}] u[x, y], {x, y}] - Div[Grad[v[x, y], {x, y}], {x, y}] u[x, y] and see if this solves your problem. (I know the equation looks strange),if this is the case, I have an idea what is going on. $\endgroup$
    – user21
    Aug 29, 2014 at 9:36
  • $\begingroup$ Or alternatively, this equivalent: op = 2*E^(-x^2 - y^2)*y*Derivative[0, 1][u][x, y] + Derivative[0, 2][u][x, y] + 2*E^(-x^2 - y^2)*x*Derivative[1, 0][u][x, y] + Derivative[2, 0][u][x, y] $\endgroup$
    – user21
    Aug 29, 2014 at 9:41
  • $\begingroup$ Hi user21. When I do that, I get equal fluxes into and out of the box, as desired. But then I'm also solving a different PDE, since your op differs from my op by a term - Div[Grad[v[x, y], {x, y}], {x, y}] u[x, y]. $\endgroup$
    – Greg P
    Aug 29, 2014 at 14:57
  • $\begingroup$ I suspect the natural Neumann b.c. is just what you have observed ( zero gradient, which in this case is not zero flux ). That is to say the result is correct and the issue lies in the docs if they claim natural==zero flux. (My guess they only make that claim in the context of examples w/o the potential term ). $\endgroup$
    – george2079
    Aug 29, 2014 at 18:54
  • $\begingroup$ @george2079, the natural Neumann bc depends on the equation. In this case the equation evaluated to something where the natural bc is $\nabla u \cdot n = 0$. If you do find that claim natural==zero flux in the docs, I'd appreciate if you could let me know. $\endgroup$
    – user21
    Aug 29, 2014 at 20:28

2 Answers 2


Here is what I think the issue is: Let's look at what NDSolve parses.

{state} = 
  NDSolve`ProcessEquations[{op == 0, bc}, u, {x, y} ∈ reg, 
   Method -> {"FiniteElement", 
     "MeshOptions" -> {MaxCellMeasure -> 0.005}}];
femData = state["FiniteElementData"];

{{{{0}}, {{{{0}, {0}}}}}, {{{{{-1, 
      0}, {0, -1}}}}, {{{{0}, {0}}}}, {{{{2 E^(-x^2 - y^2) x, 
      2 E^(-x^2 - y^2) y}}}}, {{4 E^(-x^2 - y^2) - 
     4 E^(-x^2 - y^2) x^2 - 4 E^(-x^2 - y^2) y^2}}}, {{{0}}}, {{{0}}}}

So, there is the diffusion term ({{-1,0},{0,-1}}) the convection term ({2 E^(-x^2 - y^2) x, 2 E^(-x^2 - y^2) y}) and a reaction term. This means NDSolve is modeling something like:

∇⃗ ⋅[∇⃗ u(x,y)] + v(x,y)∇⃗ u(x,y) + u(x,y)∇⃗ v(x,y) = 0

If we manually set the PDE coefficients like this, to be a diffusion coefficient and a conservative convection coefficient (with an added minus sign, as noted in the comments)

methodData = femData["FEMMethodData"];
bcData = femData["BoundaryConditionData"];
sd = state["SolutionData"][[1]];
vd = methodData["VariableData"];

pdeData = 
 InitializePDECoefficients[vd, sd, 
  "DiffusionCoefficients" -> {{-IdentityMatrix[2]}}, 
  "ConservativeConvectionCoefficients" -> {{-{2 E^(-x^2 - y^2) x, 
      2 E^(-x^2 - y^2) y}}}];

This should be more like

∇⃗ ⋅[∇⃗ u(x,y) + u(x,y)∇⃗ v(x,y)]= 0

And discretize and solve the PDE with:

dpde = DiscretizePDE[pdeData, methodData, sd];
dbcs = DiscretizeBoundaryConditions[bcData, methodData, sd];
{l, s, d, m} = dpde["SystemMatrices"];
DeployBoundaryConditions[{l, s, d, m}, dbcs]
lsol = LinearSolve[s, l];
mesh = NDSolve`SolutionDataComponent[sd, "Space"]["ElementMesh"];
sol = {u -> ElementMeshInterpolation[{mesh}, lsol]};

We see that the difference between the in and out flux is relatively small.

Plot[{Evaluate[((Derivative[0, 1][u][x, 5] + 
        u[x, 5] Derivative[0, 1][v][x, 5]) - (Derivative[0, 1][u][
         x, -5] + u[x, -5] Derivative[0, 1][v][x, -5])) /. 
    sol]}, {x, -5, 5}]

enter image description here


In 10.0.2 a better way to input the PDE and get the expected result is made available:

reg = 
  ImplicitRegion[-5 <= x <= 5 && -5 <= y <= 5 && x^2 + y^2 >= 1^2, {x,
v = Function[{x, y}, -E^(-x^2 - y^2)];
bc = {DirichletCondition[u[x, y] == 1, y == -5], 
   DirichletCondition[u[x, y] == 0, y == 5]};
α = Grad[v[x, y], {x, y}];
op = Inactive[Div][
   Inactive[Plus][Inactive[Grad][u[x, y], {x, y}], 
    Inactive[Times][α, u[x, y]]], {x, y}];
solI = NDSolve[{op == 0, bc}, u, {x, y} ∈ reg];
    (Derivative[0, 1][u][x, y] + 
       u[x, y] Derivative[0, 1][v][x, y]) /. 
     solI] /@ {5, -5}], {x, -5, 5}, PlotRange -> {-0.2, 0}]
  • $\begingroup$ Sorry, this is my first time going into the depths of NDSolve. Can you specify where you got vd,sd,methodData, and bcData? $\endgroup$
    – Greg P
    Aug 29, 2014 at 17:20
  • $\begingroup$ @GregP, sorry my bad. Forgot some of the code. I am not 100% sure this should be considered a bug in NDSolve, It has to do with 'premature' evaluation of the PDE. I'll have to think a bit about this. Let me know what you think. $\endgroup$
    – user21
    Aug 29, 2014 at 18:12
  • $\begingroup$ @GregP, this might be of interest to you. $\endgroup$
    – user21
    Aug 29, 2014 at 18:15
  • $\begingroup$ Many thanks, @user21. It seems to work quite well, although I did change the sign of ConservativeConvectionCoefficients compared to your code. That is, I had "ConservativeConvectionCoefficients" ->{{{-2 E^(-x^2 - y^2) x,-2 E^(-x^2 - y^2) y}}}. This is because the force is minus the gradient of the potential. Alternatively, I suppose I could have set the diffusion tensor to the identity matrix rather than minus the identity matrix. $\endgroup$
    – Greg P
    Aug 29, 2014 at 19:43
  • 1
    $\begingroup$ @GregP, yes, the sgin change makes sense. I did not think too much about this. Would you mind if I use this example or some variation of it in the documentation. I think it makes for a nice example. I think in essence one would like to specify something like Inactive[Div[Grad[u[x,y],{x,y}]+Grad[v[x,y],{x,y}]*u[x,y],{x,y}] but that does not work currently. $\endgroup$
    – user21
    Aug 29, 2014 at 19:45

I would say that you misread the documentation for NeumannValue. I found the documentation fairly difficult to fathom. This was mainly my fault because what they do is substantially different than what I expected. I was trying to cram the round peg that is NeumannValue into the square hole that is my brain.

If you can write the PDE as

Div J = 0

(the form your PDE is in) and you use NeumannValue with its first argument set to zero you are then getting n.J=0. That is if you write your pde as

 Div J = NeumannValue[q - r u,pred]

where u is the dependent variable you are solving for and where pred is true for coordinates that are in the boundaries and false for all other coordinate values then NeumannValue is imposing n.J = q-ru on the boundaries for which pred is true. Thus if you set q =r=0 so that the first argument to NeumannValue is zero then you would in fact be getting the boundary conditions you wanted. If your PDE can't be written as the divergence of a flux (plus time derivatives possibly) or if you want some other boundary conditions then you can't use NeumannValue as far as I know.

  • $\begingroup$ I'm confused. As far as I can tell, your description of the behavior of NDSolve for the PDE $Div J = 0$ and with NeumannValue is exactly what I wanted/hoped for. Given the default NeumannValue[0,pred] behavior when no B.C. is specified, this (desired behavior) would imply no net flux of $J$ out of the box. However, the last plot in my post shows that there is in fact net flux (of $J$). As far as I understand, what you describe above is how I wished NDSolve would work, but didn't seem to. $\endgroup$
    – Greg P
    Nov 14, 2014 at 2:24
  • $\begingroup$ The parsing of the PDE was an issue. If you want a more complicated Neumann value to be satisfied then you'd need to describe your PDE with Inactive $\endgroup$
    – user21
    Dec 16, 2014 at 11:11
  • $\begingroup$ If can suggest a little more precisely what you found difficult to fathom, then it might be possible to improve that. $\endgroup$
    – user21
    Dec 16, 2014 at 13:36

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service and acknowledge you have read our privacy policy.

Not the answer you're looking for? Browse other questions tagged or ask your own question.