# Problem with PDE boundary and initial conditions

I am trying to numerically solve the following PDE

pde =
With[{d = .25, \[Rho] = .5, M = 10, T = 2,
L = 10, \[Alpha] = 0.2},
{
D[n[x, t], {t, 1}] ==
Piecewise[
Table[{(1 - x^2/(3 d (M T - i T)))*
n[x, t]* g / Sqrt[3 d (M T - i T)] ,
i T < t < i T + \[Alpha] T && -\[Sqrt](3 d*(M T - i T)) <
x < \[Sqrt](3 d*(M T - i T)) }, {i, 0, M - 1}]
]
+
Piecewise[
Table[{d * D[n[x, t], {x, 2}] - \[Rho] * n[x, t],
i T + \[Alpha] T < t < (i + 1) T}, {i, 0, M - 1}]
],
n[x, 0] == InitialCondition,
(D[n[x, t], x] /. x -> -L) == (D[n[x, t], x] /. x -> L) == 0
}
];


However, when I try to get the solution through

sln = With[{M = 10, L=10},
ParametricNDSolveValue[
pde /. InitialCondition -> 0.5  (*uniform*),
n, {x, -L, L}, {t, 0, M T}, {g},
AccuracyGoal -> 20, PrecisionGoal -> 10]
];


and then write, for example, sln[.5][0,5], MMA v13.2 gives an error NDSolve::ibcinc stating the initial and boundary conditions are inconsistent, but it also produces a result. Could someone help me figure out what's happening?

I want to keep g as a parameter to then be able to manipulate the plots by varying it, say

Manipulate[Plot[{sln[g][0, t]}, {t, 0, 10}, PlotRange -> All],{g, 0.1, 10}]

However, in that case MMA throws out the error and then stops the evaluation.

PS: The reason I have added PrecisionGoal is that without it, at longer times the solution seems to have a lot of numerical error. For example, after removing PrecisionGoal and AccuracyGoal, MMA solves the pde with no error, but when I look at Plot[{sln[.1][x, 14]}, {x, -10, 10}, PlotRange -> All], the plot does not make sense to me.

• PS2: I realized with only PrecisionGoal there is no error, but the final results/plots don't make sense. However I got things to work with AccuracyGoal -> 15, PrecisionGoal -> 10. This seems like an accident, since with a different initial condition n[x,0] = Sin[ Pi x /(2L)]^2 the error persists even with these accuracy and precision goals. Commented May 24, 2023 at 20:01
• @AlexTrounev the Neumann boundary condition is on the derivative of the function, not its value Commented May 25, 2023 at 13:01
• What do you try to solve? It looks like delay equation. Commented May 25, 2023 at 15:18
• @AlexTrounev It's a reaction-diffusion system, where at specified time intervals there is particle creation (the first piecewise), and in other intervals there is diffusion and particle annihilation Commented May 26, 2023 at 3:59
• Thank you, please, see my answer. Commented May 26, 2023 at 4:12

## 1 Answer

This problem can be solved with using Mathematica FEM as follows

Needs["NDSolveFEM"]

{d = .25, \[Rho] = .5, M = 10, T = 2, L = 10, \[Alpha] = 0.2}; mesh =
ToElementMesh[ImplicitRegion[-L <= x <= L, {x}],
MaxCellMeasure -> 1/100];
pde = {D[n[x, t], {t, 1}] == (f1[t, x, g] + f3[t]) n[x, t] +
f2[t] D[n[x, t], x, x], n[x, 0] == InitialCondition};

f1[t_, x_, g_] :=
g Piecewise[
Table[{(1 - x^2/(3 d (M T - i T)))/Sqrt[3 d (M T - i T)],
i T < t < i T + \[Alpha] T && -\[Sqrt](3 d*(M T - i T)) <
x < \[Sqrt](3 d*(M T - i T))}, {i, 0, M - 1}]];
f2[t_] :=
Piecewise[
Table[{d, i T + \[Alpha] T < t < (i + 1) T}, {i, 0, M - 1}]];
f3[t_] :=
Piecewise[
Table[{-\[Rho], i T + \[Alpha] T < t < (i + 1) T}, {i, 0, M - 1}]];
sol = NDSolveValue[pde /. {InitialCondition -> .5, g -> .5},
n, {t, 0, M T}, {x} \[Element] mesh];


Visualization

{Plot3D[sol[x, t], {t, 0, T}, {x, -L, L}, ColorFunction -> Hue,
AxesLabel -> Automatic, MeshStyle -> White],
Plot3D[sol[x, t], {t, 0, M T}, {x, -L, L}, ColorFunction -> Hue,
AxesLabel -> Automatic, MeshStyle -> White]}


Update 1. In a case of initial condition in a form of n[x,0] = Sin[ Pi x /(2L)]^2 we also have numerical solution without errors. It looks like this one

• Do we not need to specify boundary conditions? Commented May 26, 2023 at 12:53
• @SaMaSo In a case of Mathematica FEM when using (D[n[x, t], x] /. x -> -L) == (D[n[x, t], x] /. x -> L) == 0 there is automatic NeumannValue[0, True]  applied. So you don't need to take into account these boundary conditions. Commented May 26, 2023 at 16:31
• Ok thank you. Just to make sure I understand correctly, with FEM I would only need to specify the BC if they are Dirichlet or mixed, or otherwise Neumann is automatically enforced. Is that right? Commented May 30, 2023 at 13:14
• Yes you are right, NeumannValue[0, True] is automatically enforced. Commented May 30, 2023 at 15:15