Convergence of PDE solution using method of lines

I'm afraid that this will turn more into a math question rather than a Mathematica one.

I'm trying to solve the equation

$$\frac {\partial n}{\partial t}=D\frac {\partial^2n}{\partial x^2}$$ $$\partial_x n(0,t)=n(0,t)$$

$$\partial_x n(L,t)=1$$ $$n(x,0)=1+x$$

using the code

sol = NDSolveValue[{D[u[t, x], {t, 1}] - D[u[t, x], {x, 2}] ==
NeumannValue[u[t, x], x == 0] + NeumannValue[1, x == 10],
u[0, x] == 1 + x},
u, {x, 0, 10}, {t, 0, 10}, Method -> {"MethodOfLines",
"SpatialDiscretization" -> {"FiniteElement"}}]


Please note that the initial condition $$n(0,x)$$ is the function that solves the equation in steady state. So I would be expecting that the solution $$n(t,x)$$ would stay the same instead of giving the result that you can obtain by running for example:

plots = Table[
Plot[sol[i, x], {x, 0, 10}, PlotRange -> All],
{i, 0, 10}
];

ListAnimate[plots]


Please also note that at some point the left boundary condition is not even being respected.

There is a simple sign error in the set up I would think; the left NeumannValues needs a negative sign:

sol = NDSolveValue[{D[u[t, x], {t, 1}] - D[u[t, x], {x, 2}] ==
NeumannValue[-u[t, x], x == 0] + NeumannValue[1, x == 10],
u[0, x] == 1 + x}, u, {x, 0, 10}, {t, 0, 10},
Method -> {"MethodOfLines",
"SpatialDiscretization" -> {"FiniteElement"}}];

plots = Table[
Plot[sol[i, x], {x, 0, 10}, PlotRange -> All], {i, 0, 10}];
ListAnimate[plots]
Plot[(D[sol[t, x], x] - sol[t, x]) /. x -> 0, {t, 0, 10},
PlotStyle -> Blue]


This then gives:

For details see the 'details' section of the ref page of NeumannValue and also the section The Relation between NeumannValue and Boundary Derivatives.

• Thanks for pointing out that! So, as far as I understood by those documents, we have to put the minus sign at x=0 because of the way the normal to the boundary is defined?
– AJHC
Commented Oct 31, 2018 at 22:16

It looks like a bug. On the other hand, if we correctly formulate the task and use the automatic method, we get the expected result.

sol = NDSolveValue[{D[u[t, x], {t, 1}] - D[u[t, x], {x, 2}] == 0,
u[0, x] == 1 + x, u[t, 0] ==
\!$$\*SuperscriptBox[\(u$$,
TagBox[
RowBox[{"(",
RowBox[{"0", ",", "1"}], ")"}],
Derivative],
MultilineFunction->None]\)[t, 0],
\!$$\*SuperscriptBox[\(u$$,
TagBox[
RowBox[{"(",
RowBox[{"0", ",", "1"}], ")"}],
Derivative],
MultilineFunction->None]\)[t, 10] == 1}, u, {x, 0, 10}, {t, 0, 10}]

{ContourPlot[sol[t, x], {x, 0, 10}, {t, 0, 10},
PlotLegends -> Automatic, PlotRange -> All],Plot[Table[sol[t, x], {t, 0, 10, 1}], {x, 0, 10},
PlotRange -> {0, 20}],Plot[(D[sol[t, x], x] - sol[t, x]) /. x -> 0, {t, 0, 10},
PlotStyle -> Blue]}