# Mathematica thinks that an initial condition is a boundary condition

I would like to solve the following reaction-diffusion problem in Mathematica using NDSolve:

DiffConst = 150;
kon = 0.14/60;
koff = 0.18/60;
leftBoundary = 0.;
rightBoundary = 10.;
eqs = {D[L[x, t], t] == DiffConst*D[L[x, t], {x, 2}] + koff*LR[x, t] -
kon*L[x, t]*R[x, t] + NeumannValue[0, x == leftBoundary] +
NeumannValue[0, x == rightBoundary],
D[LR[x, t], t] == -koff*LR[x, t] + kon*L[x, t]*R[x, t] + NeumannValue[0, x == leftBoundary] +
NeumannValue[0, x == rightBoundary],
D[R[x, t], t] == koff*LR[x, t] - kon*L[x, t]*R[x, t] + NeumannValue[0, x == leftBoundary] +
NeumannValue[0, x == rightBoundary],
L[x, 0] == 10,
R[0, 0] == 10,
LR[x, 0] == 0}
NDSolve[eqs, {L, R, LR}, {t, 0, 100}, {x, 0, 10}]


The idea is to apply no-flux boundary conditions (with the NeumannValue functions) and apply the initial conditions in the last three lines (LR is zero all over, L is 10 everywhere, and R is 10 at x=0). I get the "NDSolve::bcedge: Boundary condition R[0,0]==10 is not specified on a single edge of the boundary of the computational domain" error message. R[0,0] is supposed to be an initial condition, but Mathematica believes it is a boundary condition. What's the problem? Thanks for any help.

• Perhaps it should be R[0, t] == 10 ? Sep 16, 2023 at 16:36
• Or rather R[x,0] == 10, since it's an initial condition. Setting R[0,0]=10 clashes with your Neumann conditions for the left boundary at x=0, since they are independent of t. You might also consider using the predicate in the Neumann condition to excise the point {0,0} to avoid any further conflict using an inequality instead of the default True everywhere.
– user87932
Sep 16, 2023 at 17:52

Here is a way to solve it. Note that there is no need for the NeumannValue, as the 0 NeumannValue is the default. I have changed the initial condition on R and gave a method option to use the finite element method:

DiffConst = 150;
kon = 0.14/60;
koff = 0.18/60;
leftBoundary = 0.;
rightBoundary = 10.;
eqs = {D[L[x, t], t] ==
DiffConst*D[L[x, t], {x, 2}] + koff*LR[x, t] -
kon*L[x, t]*R[x, t](*+NeumannValue[0,x==leftBoundary]+
NeumannValue[0,x==rightBoundary]*),
D[LR[x, t], t] == -koff*LR[x, t] +
kon*L[x, t]*R[x, t](*+NeumannValue[0,x==leftBoundary]+
NeumannValue[0,x==rightBoundary]*),
D[R[x, t], t] ==
koff*LR[x, t] - kon*L[x, t]*R[x, t](*+NeumannValue[0,x==
leftBoundary]+NeumannValue[0,x==rightBoundary]*),
L[x, 0] == 10, R[x, 0] == 10, LR[x, 0] == 0};
res = NDSolveValue[eqs, {L, R, LR}, {t, 0, 100}, {x, 0, 10},
Method -> {"PDEDiscretization" -> {"MethodOfLines",
"SpatialDiscretization" -> "FiniteElement"}}]

• Thanks. This is syntactically correct, but not what I wanted to achieve. I wanted R to be 0 everywhere at t=0, except at a certain location. This is the code I replaced R[x,0]==10 with: R[x, 0] == If[x == 0, 10, 0] or R[x, 0] == If[ 0 <= x <= 0.1, 10, 0] Sep 18, 2023 at 15:03