# Inhomogeneous Neumann boundary conditions for diffusion equation

I am new to Mathematica and I have a problem specifying Neumann boundary conditions in diffusion equation. The best result I managed to get is this.

n0 = 4*10^15; d = 0.005; diff = 17.647;  d = 0.005;
Manipulate[
Module[{soln =
NDSolve[{ D[n[x, t], x, x] == 1/diff D[n[x, t], t] +
NeumannValue[1.25*10^18, x == -d] -
NeumannValue[-5.2*10^17, x == d],
DirichletCondition[n[x, t] == n0, t == 0]},
n, {x, -d, d}, {t, 0, 10^-5}]},
Plot[n[x, t0] /. soln, {x, -d, d},
PlotRange -> {{-d, d}, {0, 5 10^15}}]], {t0, 0, 10^-5,
Appearance -> "Labeled"}]


It is actually not so bad, but from the plot you can clearly see that Neumann conditions are not satisfied, as the divergence of n[x,t] at -d and d becomes time-dependent, though it should be constant. I have an analytical solution of this equation obtained by Fourier method and it works fine, but I need a numerical solution as well. My problem seems to originate from lack of understanding how NeumannValue works, but the way it is explained in manuals makes me even more confused. Dirichlet condition is far more straightforward, but with NeumannValue I don't really understand why and in what form do I need to sum it with my equation. I guess in my case I need to multiply NeumannValue on something (presumably, time-dependent and having a dimention of inverse length), but I can't figure it out on what exactly. I tried to search this problem on Stack Exchange, but it didn't help me understand what to do to get a result. I am sorry for my non-native English and probably a silly question.

Eited #1

As Alex Trounev noticed, the boundary and initial conditions are inconsistent, in my case as at t==0 the direction of current changes istanteneously, leaving aside how it can be physically (though experiments show that it is a good approximation for this problem). It is actually not important for analytical solution, but seems to have an impact on numerical computation. I tried to correct it artificially (no idea if it is a correct way) by introducing a piecewise-linear initial condition with a small parameter delta so that the derivatives at boundaries match the boundary conditions, and at the same time in the main region concentration is n0. It did not change a lot.

  n0 = 4*10^15; d = 0.005; dif = 17.647; d = 0.005; delta = 10^-4;
Manipulate[
Module[{soln = NDSolve[{D[n[x, t], x, x] == 1/dif D[n[x, t], t] +
NeumannValue[1.25*10^18, x == -d] -
NeumannValue[-5.2*10^17, x == d],
DirichletCondition[
n[x, t] ==
1.25*10^18 (x + d - delta)*UnitStep[delta - x - d] + n0 -
5.2*10^17 (-d + x + delta)*UnitStep[-d + x + delta],
t == 0]}, n, {x, -d, d}, {t, 0, 10^-5}]},
Plot[n[x, t0] /. soln, {x, -d, d},
PlotRange -> {{-d, d}, {0, 5 10^15}}]], {t0, 0, 10^-5,
Appearance -> "Labeled"}, {delta, 2 10^-4, 0,
Appearance -> "Labeled"}]


Next I tried another variant using derivatives directly instead of NeumannValue, but using the same corrected initial condition. The result is much closer to what I expect, as the boundary conditions are now visually satisfied, but there are still several problems: first, the initial condition is somehow smoothened, so at t0==0 I have a smooth function instead of piecewise-linear, an I get a warning that "boundary and initial conditions are inconsistent" anyway. Second, the solution depends on delta quite strongly, which is certainly not what I need to have. That's how it looks like.

n0 = 4*10^15; d = 0.005; dif = 17.647; d = 0.005; delta = 10^-4;
Manipulate[
Module[{soln =
NDSolve[{dif Derivative[2, 0][n][x, t] ==
Derivative[0, 1][n][x, t],
Derivative[1, 0][n][-d, t] == 1.25*10^18,
Derivative[1, 0][n][d, t] == -5.2*10^17,
n[x, 0] ==
1.25*10^18 (x + d - delta)*UnitStep[delta - x - d] + n0 -
5.2*10^17 (-d + x + delta)*UnitStep[-d + x + delta]},
n, {x, -d, d}, {t, 0, 10^-5}]},
Plot[n[x, t0] /. soln, {x, -d, d},
PlotRange -> {{-d, d}, {0, 5 10^15}}]], {t0, 0, 10^-5,
Appearance -> "Labeled"}, {delta, 2 10^-4, 0,
Appearance -> "Labeled"}]

New contributor
lygeon is a new contributor to this site. Take care in asking for clarification, commenting, and answering. Check out our Code of Conduct.
• What is set at the boundaries - heat flow or derivative? – Alex Trounev 2 days ago
• It is not clear why the initial data is given a constant, and the non-zero derivative is given at the boundaries. The initial and boundary data must be consistent. – Alex Trounev 2 days ago
• Physically it is an equation of ambipolar diffusion of electrons and holes in a diode base with boundaries at -d and d. To obtain an initial carrier distribution we first apply direct current to fill the base with charge carriers by injection (I simplified the equation so that carrier lifetime is infinite). At time moment t=0 we instantaneously change the direction of a current, and the boundary values dn/dx become proportional to reverse current density j. That is why my boundary and initial values are inconsistent, and I believe it shouldn't be a problem (at least for analytical solution). – lygeon 2 days ago
• Analytical solution is given in an article "Reverse recovery processes in silicon power rectifiers" by Benda, Spenke, 1967. – lygeon 2 days ago
• It is impossible physically instantaneously change the direction of a current. There must be a transition process. In the numerical model, it is also desirable to smoothly include the boundary conditions. Describe exactly how the boundary derivatives are given. – Alex Trounev 2 days ago

If you need to use inconsistent boundary conditions, then you could try supplying your own computational mesh where you provide refinement in the neighborhood of the inconsistencies to mitigate their effects. Details on how to build your own mesh can be found in the Element Mesh Generation Tutorial.

You have a lot of parameters that range over many orders of magnitude. I always find it useful to make the equations dimensionless to give the model scales that range closer to unity. The subscript d indicates the variable has dimensions. Unsubscripted implies dimensionless.

$$\frac{1 }{{{D}}}\frac{{\partial n_d}}{{\partial t_d}}- \frac{\partial^2 n_d}{\partial x_{d}^{2}}=q_l+q_r\left\| {\frac{{\varDelta^2}}{n_0}} \right.$$ $$\frac{\varDelta^2 }{{{n_0 D}}}\frac{{\partial n_d}}{{\partial t_d}}- {\frac{{\varDelta^2}}{n_0}}\frac{\partial^2 n_d}{\partial x_{d}^{2}}={\frac{{\varDelta^2}}{n_0}}\left( q_l+q_r\right )$$ The dimensionless form of your equation below. $$\boxed{\frac{{\partial n}}{{\partial t}} - \frac{{{\partial ^2}n}}{{\partial {x^2}}} = \frac{{{\Delta ^2}}}{{{n_0}}}\left( {{q_l} + {q_r}} \right)}$$ Where,

$$\Delta= 2 d$$

$${\tau} = \frac{{\Delta ^2}}{D}$$

$$t = \frac{{{t_d}}}{{{\tau }}};{t_d} = {\tau }t$$

$$x = \frac{{{x_d}}}{\Delta};{x_d} = \Delta x$$

$$n = \frac{{{n_d}}}{{{n_0}}};{n_d} = {n_0}n$$

A 1d transient diffusion problem can sometimes be thought of as convection-diffusion problem in two dimensions. In our case, the "y-axis" will be the time axis. The following code will create a mesh with refinement in the inconsistent regions and solve the dimensionless PDE.

(* Load FEM Package *)
Needs["NDSolveFEM"]
(* Parameters *)
(* tmax = 1.7647 is maximum scaled time t/tau *)
(* Scaled x coordinate is x/delta *)
n0 = 4*10^15; delta = 0.01; d =
delta/2; diff = 17.647; tau = delta^2/diff; tmax = 10^-5/tau;
(* Scaled Fluxes *)
ql = - (delta^2/n0) 1.25*10^18;
qr = -(delta^2/n0) 5.2*10^17;
(* Define Domain vertical is time axis *)
ht = tmax;
len = 1;
top = ht;
bot = 0;
left = -len/2;
right = len/2;

(* create associations for easy BC assignment later *)
bounds = <|l -> 1, r -> 2, init -> 3|>;
regs = <|domain -> 10|>;

(* Create a boundary mesh *)
bmesh = ToBoundaryMesh[
"Coordinates" -> {{left, bot}(*1*), {right, bot}(*2*), {right,
top}(*3*), {left, top}(*4*)},
"BoundaryElements" -> {LineElement[{{1, 2}(*bottom edge 1*), {4,
1}(*left edge 2*), {2, 3}(*right edge 3*), {3,
4}(*top edge 4*)}, {bounds[init], bounds[l], bounds[r], 4}]}];

(* Remove ";" to view mesh *)
bmesh["Wireframe"["MeshElementMarkerStyle" -> Blue,
"MeshElementStyle" -> {Green, Red, Orange, Black},
ImageSize -> Large]];

(* Create 2D mesh with refinement in the r/l bottom corners *)
buf = 0.015;
mesh = ToElementMesh[bmesh,
"RegionMarker" -> {{{(left + right)/2, (bot + top)/2},
regs[domain], 0.001}};
MeshRefinementFunction ->
Function[{vertices, area},
Block[{x, y}, {x, y} = Mean[vertices];
If[(y < bot + buf && x < left + buf) || (y < bot + buf &&
x > right - buf), area > 0.000000625/4, area > 0.00025]]]];

(* Remove ";" to view mesh *)
mesh["Wireframe"["MeshElementStyle" -> {FaceForm[Green]},
ImageSize -> Large]];

(* Set up BCs and Equations *)
nvleft = NeumannValue[ql, ElementMarker == bounds[l]];
nvright = NeumannValue[qr, ElementMarker == bounds[r]];
dcinit = DirichletCondition[n[x, t] == 1,
ElementMarker == bounds[init]];
eqn = D[n[x, t], t] - D[n[x, t], x, x] == nvleft + nvright;

(* Solve the PDE *)
ifun = NDSolveValue[{eqn, dcinit}, {n}, {x, t} \[Element] mesh];

(* Plot the boundary fluxes at early time *)
Plot[{Evaluate[D[ifun[[1]][x, t], x] /. {t -> tp, x -> 1/2}],
Evaluate[D[ifun[[1]][x, t], x] /. {t -> tp, x -> -1/2}]}, {tp,
0.001, tmax/200}, PlotRange -> All]

(* A bunch of plots *)
Plot[Evaluate[
Table[ifun[[1]][x, t0], {t0, 0, tmax, tmax/50}]], {x, -1/2, 1/2},
PlotRange -> {{-1/2, 1/2}, {0.91, 1.001}}, ImageSize -> Large]

(* Flux Differences *)
Abs[(qr - D[ifun[[1]][x, t], x] /. {t -> 0.005, x -> 1/2})/
qr] (* 0.0022975495339616677 *)
Abs[(qr - D[ifun[[1]][x, t], x] /. {t -> 0.1, x -> 1/2})/
qr] (* 0.0003545675754438614 *)
Abs[(ql + D[ifun[[1]][x, t], x] /. {t -> 0.005, x -> -1/2})/
ql] (* 0.005748069660621691 *)
Abs[(ql + D[ifun[[1]][x, t], x] /. {t -> 0.1, x -> -1/2})/
ql] (* 1.5152106243032648*^-6 *)

(* Scaled manipulate  *)
Manipulate[
Plot[Evaluate[ifun[[1]][x, t0]], {x, -1/2, 1/2},
PlotRange -> {{-1/2, 1/2}, {0.91, 1.001}}], {t0, 0, tmax,
Appearance -> "Labeled"}]
(* Unscaled manipulate *)
Manipulate[
Plot[Evaluate[n0 ifun[[1]][x/delta, t0/tau]], {x, -d, d},
PlotRange -> {{-d, d}, {0.91 n0, 1.001 n0}}], {t0, 0, 10^-5,
Appearance -> "Labeled"}]


## 50 Plots from $$0\leq t\leq t_{max}$$

Besides the line at $$t=0$$, the tangent lines at the boundaries appear to be parallel.

## Computational Mesh in Lower Left Corner

We can get a better view of the discretization by zooming into the lower left corner.

Show[mesh[
"Wireframe"["MeshElementStyle" -> {FaceForm[Green]},
ImageSize -> Large]], PlotRange -> {{-1/2, -1/2.2}, {0, 0.04}}]


New contributor
Tim Laska is a new contributor to this site. Take care in asking for clarification, commenting, and answering. Check out our Code of Conduct.

It is necessary to normalize the data to n0, add an option and a smooth transition function.

n0 = 4*10^15; d = 0.005; dif = 17.647; tm = 10^-5;
Manipulate[
Module[{soln =
NDSolve[{dif Derivative[2, 0][n][x, t] ==
Derivative[0, 1][n][x, t],
Derivative[1, 0][n][-d,
t] == (1.25*10^18/n0)*(1 - Exp[-k*t/tm]),
Derivative[1, 0][n][d, t] == -(5.2*10^17/n0)*(1 - Exp[-k*t/tm]),
n[x, 0] == 1}, n, {x, -d, d}, {t, 0, tm},
MaxStepSize -> 0.01]},
Plot[n[x, t0] /. soln, {x, -d, d}, PlotRange -> All]], {t0, 0,
tm/100, Appearance -> "Labeled"}, {k, 1000, 10000,
Appearance -> "Labeled"}]
With[{k = 1000},
soln = NDSolve[{dif Derivative[2, 0][n][x, t] ==
Derivative[0, 1][n][x, t],
Derivative[1, 0][n][-d, t] == 1.25*10^18/n0*(1 - Exp[-k*t/tm]),
Derivative[1, 0][n][d, t] == -.52*10^18/n0*(1 - Exp[-k*t/tm]),
n[x, 0] == 1}, n, {x, -d, d}, {t, 0, tm},MaxStepSize -> 0.01]]
{Plot[Evaluate[
Table[n[x, t] /. soln, {t, 0, tm/1000, .0001*tm}]], {x, -d, d},
PlotRange -> All, PlotLegends -> Automatic],
Plot3D[n[x, t] /. soln, {t, 0, tm/1000}, {x, -d, d},
PlotRange -> All, Mesh -> None, ColorFunction -> Hue]}


• That works just perfect, thank you! Just in case if anyone knows whether it is possible to make it work with NeumannValue and inconsistent boundary and initial conditials as in my first code (which somehow worked not complaining on inconsistency though not exactly correct because I specified NeumannValue in a wrong way), it would be much appreciated. – lygeon yesterday

Here are a few thoughts: I assume you want to solve this as a time dependent problem. For that you should not specify the initial condition as a DirichletCondition. Then, as other have noted the scales have huge differences and this affects the time integration. For that reason I added a StartingStepSize that is small enough to get NDSolve to integrate.

soln = NDSolveValue[{diff D[n[x, t], x, x] ==
D[n[x, t], t] + NeumannValue[1.25*10^18, x == -d] -
NeumannValue[-5.2*10^17, x == d], n[x, 0] == 0},
n, {x} \[Element] Line[{{-d}, {d}}], {t, 0, 10^-5},
StartingStepSize -> 10^-14]


Note also that I used Line[{{-d},{d}} - this forces the Finite Element Method, other wise you'd need to specify that as a method option.

Plot[soln[x, 10^-5], {x, -d, d}]


Also make sure that the signs for the diff and NeumannValues are correct.