# How to help NDSolveValue handle the addition of another dependent variable in a set of reaction-diffusion equations

I am trying to use NDSolveValue to solve a system of coupled reaction-diffusion equations in a 2D space (technically 3D cylindrical coordinates, but due to radial symmetry I only have to worry about coordinates $$r$$ and $$z$$). I have a simplified version of what I am trying to do below, but it is still a little complex, so if more detail is needed I can include it.

To summarize the system, I have one species "$$P$$" that can diffuse throughout the entire space (a "rectangular" region defined for $$0\leq r\leq\infty$$ and $$0\leq z\leq H$$ for some finite $$H$$), and then it can react with an immobile species "$$pL$$" that is only present within some strip above the $$z=0$$ plane $$0\leq z\leq s. This forms a complex "$$P\cdot pL$$" which then produces a final species "$$L$$" (which can diffuse, but for testing right now I have omitted its diffusion).

I am currently trying to put in one element at a time to make sure I have put in each part correctly. Everything works fine up until I try to put in species "$$L$$", which then gives me an error about a singular matrix in the differential equation solver.

Working Code

The mesh I am currently using to accentuate the reaction region (which is probably built in an inefficient way) is created by

Needs["NDSolveFEM"]

dr = .01; (*Close to r = 0*)
rFar = 100.; (*To "approximate infinity*)
h = 5.;
s = 1.;

coords = Join[Table[{r, z}, {z, 0, 1.2*s, s/5}, {r, dr, rFar, (rFar - dr)/90}] //Flatten[#, 1] &, Table[{r, z}, {z, 1.2*s + (h - 1.2*s)/5, h, (h - 1.2*s)/5}, {r, dr,rFar, (rFar - dr)/90}] // Flatten[#, 1] &];

elems = Table[If[Mod[a, 91] == 0, Nothing, {a, a + 1, a + 92, a + 91}], {a,Length@coords - 92}];

meshD = ToElementMesh["Coordinates" -> coords, "MeshElements" -> {QuadElement[elems]}];


Everything works fine if I include the production of "$$L$$" (given by terms with the production rate $$k_c$$) without actually keeping track of "$$L$$" in a differential equation. So the differential equations that work are given by:

pSource = (P0*k)/(2.*\[Pi]*\[Sigma]^2)*Exp[-((r - dr)^2/(2.*\[Sigma]^2)) - k*t]*(1 - Exp[-100.*t]);(*Gaussian protease source that turns on quickly and depletes slowly*)

diffeqP = D[P[r, z, t], t] - DP*Laplacian[P[r, z, t], {r, \[Theta], z}, "Cylindrical"] + kon*P[r, z, t]*pL[r, z, t] - (koff + kc)*PpL[r, z, t] == NeumannValue[pSource, z == 0];

diffeqpL = D[pL[r, z, t], t] == -kon*P[r, z, t]*pL[r, z, t] + koff*PpL[r, z, t];

diffeqPpL = D[PpL[r, z, t], t] == kon*P[r, z, t]*pL[r, z, t] - (koff + kc)*PpL[r, z, t];


Along with boundary/initial conditions:

(*Set concentration to 0 far away*)

condrFar = DirichletCondition[P[r, z, t] == 0, r == rFar];

(*Initial Conditions*)

initP = P[r, z, 0.] == 0.;

initpL = pL[r, z, 0.] == If[0 <= z <= s, pL0, 0.];

initPpL = PpL[r, z, 0.] == 0.;


With this code, a solution is obtained using NDSolveValue

params = {P0 -> 1000., k -> 5., \[Sigma] -> 20., DP -> 20., pL0 -> 1.,kon -> 100., koff -> 5., kc -> 5.};

eqns = {diffeqP, diffeqpL, diffeqPpL, condrFar, initP, initpL, initPpL} /. params;

sol = NDSolveValue[eqns, {P, pL, PpL}, {r, z} \[Element] meshD, {t, 0, 10}]


I have been looking at solutions using Plot3D in a manipulate function. i.e.

Manipulate[Plot3D[sol[[1]][r, z, t], {r, dr, rFar}, {z, 0, h}, PlotRange -> {-.00001, .06}], {t, 0, 10, .1}]

Manipulate[Plot3D[sol[[2]][r, z, t], {r, dr, rFar}, {z, 0, h}, PlotRange -> {-.00001, 1.0001}], {t, 0, 10, .1}]

Manipulate[Plot3D[sol[[3]][r, z, t], {r, dr, rFar}, {z, 0, h}, PlotRange -> {-.00001, 1.0001}], {t, 0, 10, .1}]


Broken Code

The issue comes in if I try to keep track of how much of the species "$$L$$" we have over time. This just involves adding in an additional differential equation and initial condition to our set of equations

diffeqL = D[L[r, z, t], t] == kc*PpL[r, z, t];

initL = L[r, z, 0.] == 0.;

eqns = {diffeqP, diffeqpL, diffeqPpL, diffeqL, condrFar, initP, initpL, initPpL, initL} /. params;


Then I run the differential equation solver

sol = NDSolveValue[eqns, {P, pL, PpL, L}, {r, z} \[Element] meshD, {t, 0, 10}]


and I get an error of

LinearSolve::sing: Matrix SparseArray[<<1>>] is singular.

It seems odd to me that adding in an additional species gives an error. The way I am currently handling "$$L$$" is similar to (and I would say simpler than) how I handle "$$pL$$" and "$$P\cdot pL$$", yet those are being handled fine.

Is there an explicit setting I can change in the differential equation solver to get around this "linear solve" error?

NDSolve can not always figure out which boundary condition belong to which equation. The issue you are seeing is described in the documentation. There are two things you can do in this case. One is to rename the dependent variables of your system of PDEs into something that has alphabetical ordering like so:

rename = {P -> aa, pL -> bb, PpL -> cc, L -> dd};
eqns = {diffeqP, diffeqpL, diffeqPpL, diffeqL, condrFar, initP,
initpL, initPpL, initL} /. rename /. params;
sol2 = NDSolveValue[eqns,
rename[[All, 2]], {r, z} \[Element] meshD, {t, 0, 10}]


A second approach is to put the DirichletCondition into the equation it belongs to, much like the NeumannValue:

diffeqP2 =
D[P[r, z, t], t] -
DP*Laplacian[P[r, z, t], {r, \[Theta], z}, "Cylindrical"] +
kon*P[r, z, t]*pL[r, z, t] - (koff + kc)*PpL[r, z, t] ==
NeumannValue[pSource, z == 0] + condrFar;
eqns = {diffeqP2, diffeqpL, diffeqPpL, diffeqL, initP, initpL,
initPpL, initL} /. params;
sol = NDSolveValue[
eqns, {P, pL, PpL, L}, {r, z} \[Element] meshD, {t, 0, 10}]


I understand this is very inconvenient and I will try to find a better (automatic) solution for this. Since this issue is also holding back some other FEM projects I'll have to figure this out at some point in time - but it's quite complicated and tricky.... I'll do my best, sorry for the inconvenience.

• Well that is very interesting. I was reading through that page but hadn't hit that section or even recognized its relevance to me here. Thanks for the information. Hopefully I don't hit any other things like this as I build up to my full model. It seems like I will need to be very careful about how I define my equations, coefficients, etc. Dec 19 '19 at 14:58