# Problems with Neumann (zero flux) boundary conditions

I have the following code:

(*Parameters*)eps = 1.4434; m = 0.3; c11 = 0.1732; maxCellMeasure = \
0.1;
(*PDEs*)
pde11 :=
D[pp[t, x], t] ==
0.05*Laplacian[pp[t, x], {x}] +
pp[t, x]*(1 - c11*pp[t, x] - z[t, x]/(1 + pp[t, x]^2));
pde21 := D[z[t, x], t] ==
0.05*Laplacian[z[t, x], {x}] +
z[t, x]*(eps*pp[t, x]/(1 + pp[t, x]^2) - m);
(*Initial conditions*)
lo = 48;
hi = 52;
domlen = 100;
ic11[x_] := Which[x > lo && x < hi, 6, True, 0];
ic21[x_] := Which[x < hi && x > lo, 0.5, True, 1/c11];
(*Numerical approximation using NDSolve with zero-flux boundary \
conditions*)
{solp, solz} =
Monitor[NDSolveValue[{pde11, pde21, z[0, x] == ic11[x],
pp[0, x] == ic21[x], (D[pp[t, x], x] /. x -> 0) ==
0, (D[z[t, x], x] /. x -> 0) ==
0, (D[pp[t, x], x] /. x -> domlen) ==
0, (D[z[t, x], x] /. x -> domlen) == 0}, {pp, z}, {t, 29,
30}, {x, 0, domlen},
Method -> {"FiniteElement",
MeshOptions -> MaxCellMeasure -> maxCellMeasure},
EvaluationMonitor :> (monitor = Row[{"t = ", CForm[t]}])], monitor]


To be honest, I don't really understand the usage of NeumannValue. So, I implemented zero flux (zero Neumann) boundary conditions myself. However, I get the following error, which I don't understand:

NDSolveValues: The dependent variable in pp^(0,1)[t,0]==0 in the boundary condition DirichletCondition[pp^(0,1)[t,0]==0,x==0.] needs to be linear.


Why is Mathematica writing something about DirichletConditions here?

Thank you for the help.

• 1. "FiniteElement" method cannot parse b.c.s like (D[z[t, x], x] /. x -> 0) == 0, at least now. In other words, if you want to use non-zero Neumann condition in "FiniteElement", you have to use NeumannValue. 2. However, the default setting for "FiniteElement" method is zero NeumannValue, in other words, you can just omit all of the zero Neumann condition in your code when using "FiniteElement". – xzczd Sep 6 '19 at 5:37

Have a look at the reference page of NeumannValue, there you will find a note in the details section (and examples) that address your issue:

When no boundary condition is specified on a part of the boundary ∂Ω, then the flux term ∇·(-c ∇u-α u+γ)+… over that part is taken to be f=f+0=f+NeumannValue[0,…], so not specifying a boundary condition at all is equivalent to specifying a Neumann 0 condition.

This means that in order to specify a 0 flux you need to: nothing. So this does what you want:

(*Parameters*)eps = 1.4434; m = 0.3; c11 = 0.1732; maxCellMeasure = \
0.1;
(*PDEs*)
pde11 :=
D[pp[t, x], t] ==
0.05*Laplacian[pp[t, x], {x}] +
pp[t, x]*(1 - c11*pp[t, x] - z[t, x]/(1 + pp[t, x]^2));
pde21 := D[z[t, x], t] ==
0.05*Laplacian[z[t, x], {x}] +
z[t, x]*(eps*pp[t, x]/(1 + pp[t, x]^2) - m);
(*Initial conditions*)
lo = 48;
hi = 52;
domlen = 100;
ic11[x_] := Which[x > lo && x < hi, 6, True, 0];
ic21[x_] := Which[x < hi && x > lo, 0.5, True, 1/c11];
(*Numerical approximation using NDSolve with zero-flux boundary \
conditions*)
{solp, solz} =
Monitor[NDSolveValue[{pde11, pde21, z[0, x] == ic11[x],
pp[0, x] == ic21[x]}, {pp, z}, {t, 29,
30}, {x, 0, domlen},
Method -> {"FiniteElement",
MeshOptions -> MaxCellMeasure -> maxCellMeasure},
EvaluationMonitor :> (monitor = Row[{"t = ", CForm[t]}])], monitor]


Also, note that there is a difference between a NeuamannValue and the derivative you specified. You can read up on that in the documentation here. The finite element method parses your input as a DirichletCondition and lets you know that this is not a valid DirichletCondition`.

Let me suggest that you read some of the finite element documentation. As an introduction you can read Solving Partial Differential Equations with Finite Elements.