Inhomogeneous Neumann boundary conditions for diffusion equation

Question

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"}]

Answer 1

Updated Background

This is a more detailed response to a comment from the OP about the Neumann Value specification. The following is how I like to think of it.

I was initially confused on how to setup Neumann values. After I read the Element Mesh Generation Tutorial, things got simpler for me because I could make the Mathematica workflow similar to other solver workflows by having the ability to define boundaries and regions.

To use the FEM in Mathematica, it is instructive to cast your equation in coefficient form like so

$$m\frac{{{\partial ^2}}}{{\partial {t^2}}}u + d\frac{\partial }{{\partial t}}u + \nabla \cdot\left( { - c\nabla u - \alpha u + \gamma } \right) + \beta \cdot\nabla u + au - f = 0$$

When in this form, I think of the "0" on the Right Hand Side (RHS) as indicating that the system is closed (i.e., there are no flows in or out of the system). If you want to add flows, $q_i$, to the system, then you add them to the RHS in the form of Neumann values.

$$m\frac{{{\partial ^2}}}{{\partial {t^2}}}u + d\frac{\partial }{{\partial t}}u + \nabla \cdot\left( { - c\nabla u - \alpha u + \gamma } \right) + \beta \cdot\nabla u + au - f = \sum_{i}^{n}q_i$$

By convention, Mathematica defines a positive flow as into the system as I try to illustrate below.

This is where building your own mesh and naming the boundaries comes in handy. For example, if I had a pipe, I would want to designate element IDs to an inlet and an outlet surface. Intuitively, I know that the RHS should be inflow - outflow. So, I would expect the RHS to look like

(NeumannValue[flowin, ElementMarker == bounds[inlet]] - NeumannValue[flowout, ElementMarker == bounds[outlet]])

If you want to model a system that is more complex than a line or a rectangle, then you most likely will have to build a mesh and assign named boundaries and regions. In a valid mesh, Mathematica will know the surface normals and what direction is into the system. You should know if the flows are positive or negative by their names.

Original Post

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["NDSolve`FEM`"]
(* 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"}]

Derivatives at boundary at early time (recall that scaled $t_{max}=1.7647$)

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}}]

Answer 2

The analytical solution is a little tedious but is straightforward:

pde = D[n[x, t], x, x] == D[n[x, t], t]/diff

Separation of variables by multiplication.

n[x_, t_] = X[x] T[t]

pde/n[x, t]
(*X''[x]/X[x] == T'[t]/(diff T[t])*)

Each side equal to the same constant. The constant could be 0, positive, or negative. We need only the 0 and negative constants.

The 0 constant.

X0 = X[x] /. DSolve[X''[x]/X[x] == 0, X[x], x][[1]] /.{C[1] -> c1, C[2] -> c2};

T0 = T[t] /. DSolve[T'[t]/(diff T[t]) == 0, T[t], t][[1]] /. C[1] -> 1;

The negative constant.

X1 = X[x] /. DSolve[X''[x]/X[x] == -α^2, X[x], x][[1]] /. {C[1] -> c3, C[2] -> c4}

T1 = T[t] /. DSolve[T'[t]/(diff T[t]) == -α^2, T[t], t][[1]] /. C[1] -> 1

Typically for the diffusivity equation, we have to add separation of variables by addition.

n[x_, t_] = X[x] + T[t]

pde
(*X''[x] == T'[t]/diff*)

Setting each side to a single constant will do in this case.

X2 = X[x] /. DSolve[X''[x] == β, X[x], x][[1]] /.{C[1] -> c5, C[2] -> c6};

T2 = T[t] /. DSolve[T'[t]/diff == β, T[t], t][[1]] /. C[1] -> c7

Combine all solutions.

Clear[n]

    nk[x_, t_] = X0 T0 + X1 T1 + X2 + T2
(*c1 + c2 x + E^(α^2 (-diff) t) (c3 Cos[α x] + c4 Sin[α x]) + c5 + c6 x + c7 + β diff t + (β x^2)/2*)

We only need one constant term and only one coefficient of x, so simplify by:

c6 = 0;
c5 = 0;
c7 = 0;

We get a simpler solution by working the problem from x = 0 to 2d so we will do that for now. At the end, we will back shift the x-coordinate by d. The BC at x = 0

(D[nk[x, t], x] /. x -> 0) == m1
(*c2 + α c4 E^(α^2 (-diff) t) == m1*)

Satisfy this by:

c2 = m1;
c4 = 0;

BC at x = 2d

(D[nk[x, t], x] /. x -> 2 d) == m2
(*-(α c3 Sin[2 α d] E^(α^2 (-diff) t)) + 2 β d + m1 == m2*)

From this

β = β /. Solve[m1 + 2 β d == m2, β][[1]];

Also require Sin[2 α d] == 0 to get α

α = (k π)/(2 d);

with

$Assumptions = k ∈ Integers && d > 0

nk[x_, t_] = nk[x, t] // Simplify

Apply the initial condition.

nk[x, 0] == n0 // Simplify
(*c1 + c3 Cos[α x] + (x^2 (m2 - m1))/(4 d) + m1 x == n0*)

We can get there with

c1 = n0

and

c3 Cos[(π k x)/(2 d)] == -((x^2 (m2 - m1))/(4 d) + m1 x)

Use the orthogonality integral to find c3.

c3*Integrate[Cos[(π*k*x)/(2*d)]^2, {x, 0, 2*d}] == -Integrate[((x^2*(m2 - m1))/(4*d) + m1*x)*Cos[(π*k*x)/(2*d)], 
    {x, 0, 2*d}]

c3 = c3 /. Solve[%, c3][[1]] // Simplify
(*(4 d ((-1)^(k + 1) m2 + m1))/(π^2 k^2)*)

Since we have a Cos fourier series we need a k = 0 term.

c30*Integrate[1^2, {x, 0, 2*d}] == -Integrate[((x^2*(m2 - m1))/(4*d) + m1*x)*1, {x, 0, 2*d}]

c30 = c30 /. Solve[%, c30][[1]] // Simplify
(*-(1/3) d (2 m1 + m2)*)

Add c30 to the solution and shift the x-coordinate by -d to get the boundaries from -d to d.

nk[x_, t_] = nk[x, t] + c30 /. x -> x + d

We get an infinite sum as part of the solution by summing k from 1 to infinity. For computing and plotting we will use a finite sum by using mm as the number of terms in the sum.

n[x_, t_, mm_] := (diff*t*(m2 - m1))/(2*d) + ((d + x)^2*(m2 - m1))/(4*d) - (1/3)*d*(2*m1 + m2) + m1*(d + x) + n0 + 
   ((4*d)/Pi^2)*Sum[(((-1)^(k + 1)*m2 + m1)*Cos[(Pi*k*(d + x))/(2*d)])/(E^((Pi^2*diff*k^2*t)/(4*d^2))*k^2), 
     {k, 1, mm}]

Plug in the numbers

n0 = 4*10^15;
d = 0.005;
diff = 17.647;
m1 = 1.25 10^18;
m2 = -5.2 10^17;
tmax = 10^-5;

The exponential term will give us underflow warnings, so turn them off.

Off[General::munfl]

tb = Table[
   Plot[Evaluate[n[x, t, 50]], {x, -d, d}, 
    PlotRange -> {-10 10^15, 5 10^15}, PlotLabel -> "t" N[t]], {t, 0, tmax/2, tmax/100}];
ListAnimate[tb]

And for the final time so everyone can check their numerical answers.

Plot[Evaluate[n[x, tmax, 20]], {x, -d, d}]

It is interesting to note that the Exp portion causes the Sum term of the solution to damp out after a fairly short time, leaving only the polynomial terms contributing to the solution. The series converges very rapidly except at t = 0 when the exponential term does not contribute.

Answer 3

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]}

Add a comparison with the analytical solution that found Bill Watts

n[x_, t_, 
  mm_] := (dif*t*(m2 - m1))/(2*d) + ((d + x)^2*(m2 - m1))/(4*d) - (1/
     3)*d*(2*m1 + m2) + m1*(d + x) + 
  n0 + ((4*d)/Pi^2)*
   Sum[(((-1)^(k + 1)*m2 + m1)*
       Cos[(Pi*k*(d + x))/(2*d)])/(E^((Pi^2*dif*k^2*t)/(4*d^2))*
       k^2), {k, 1, mm}]
n0 = 4*10^15; d = 0.005; dif = 17.647; tm = 10^-5; m1 = 1.25 10^18;
m2 = -5.2 10^17;

soln = ParametricNDSolveValue[{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}, {k}];


Table[{Row[{"mm=", mm}], 
  Table[Plot[
    Evaluate[
     Table[(soln[k][x, t]*n0 - n[x, t, mm])/n0, {t, 0, 
       tm, .2*tm}]], {x, -d, d}, PlotRange -> All, 
    PlotLegends -> Automatic, 
    AxesLabel -> {"x", "\[CapitalDelta]n/n0"}, 
    PlotLabel -> Row[{"k=", k}]], {k, {500, 1000, 10^4}}]}, {mm, {50, 
   100, 200}}]

The best agreement between numerical and analytical solutions is obtained for k = 500, although it seems intuitively that this should be for k = 10,000. It is difficult to explain this discrepancy with an increase in k.

Answer 4

Updated version:

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. The key point, however, is to pay close attention to where to put the diffusion coefficient. If you put it in front of the Laplacian it will have an effect on what the NeumannValue means. To reproduce the analytical solution given in one of the other posts I put it front of the time derivative.

soln = NDSolveValue[{D[n[x, t], x, x] == 
    1/diff D[n[x, t], t] + NeumannValue[m1, x == -d] - 
     NeumannValue[m2, x == d], n[x, 0] == n0}, 
  n, {x} \[Element] Line[{{-d}, {d}}], {t, 0, 10^-5}]

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.

tt = 10^-5;
Plot[{soln[x, tt], n[x, tt, 10]}, {x, -d, d}]

On such a large scale the difference between the two curves will be large, though.