# Neuman boundary conditions for system of PDE in 2D

The system of equations for three variables: n(x,y,t), b(x,y,t) and a(x,y,t)

The rest of letters are coefficients. System is defined over a square region with the lower left corner at the origin.

Boundary conditions:

Bold n is the normal vector to the boundary (Neuman condition).

And initial conditions:

My code for solution is the following

(*region*)
\[CapitalOmega] = Rectangle[{0, 0}, {1, 1}];

(*constants*)
\[Mu]n = 0.001;
\[Mu]b = 0.001;
\[Chi] = 0.1;
bHat = 1.5;
\[Lambda]0 = 100;
\[Lambda]1 = 100;
\[Lambda]2 = 10;
\[Lambda]4 = 100;
\[Lambda]5 = 10;
nHat = 1;
\[Lambda]7 = 10;
\[Delta] = 0.01;
\[Alpha] = 2.5;

aNeuman =
NeumannValue[\[Lambda]7*bHat*a[t, x, y], {x, y} \[Element]RegionBoundary[\[CapitalOmega]]];

(*equations*)
eq = {-D[n[t, x, y], t] + \[Mu]n*Laplacian[n[t, x, y], {x, y}] - \[Chi] *
Div[n[t, x, y]*Grad[a[t, x, y], {x, y}], {x, y}] + \[Lambda]1*
a[t, x, y]*b[t, x, y] - \[Lambda]2*n[t, x, y] - \[Lambda]0*(n[t, x, y])^2 == 0,

-D[a[t, x, y], t] + Laplacian[a[t, x, y], {x, y}] +
\[Lambda]4 /2*(1 + Tanh[(1 - b[t, x, y])/\[Delta]]) - (\[Lambda]4 + b[t, x, y])*a[t, x, y] == aNeuman,

-D[b[t, x, y], t] + \[Mu]b*Div[n[t, x, y]*Grad[b[t, x, y], {x, y}], {x, y}] -
Norm[\[Mu]n*Grad[n[t, x, y], {x, y}] - \[Chi] *n[t, x, y]*Grad[a[t, x, y], {x, y}]] == 0};

(*initial conditions*)
incs = {n[0, x, y] == If[x == 0 || y == 0 || x == 1 || y == 1, nHat, 0],
a[0, x, y] == 0,
b[0, x, y] == If[x == 0 || y == 0 || x == 1 || y == 1, bHat, 0]};

(*boundary conditions*)
bcs = {DirichletCondition[n[t, x, y] == nHat*E^(-\[Alpha]*t), True],
DirichletCondition[b[t, x, y] == bHat, True]};

(*solution*)
{nfun, afun, bfun} = NDSolveValue[Join[eq, bcs, incs], {n, a, b}, {t, 0, 2}, {x, y} \[Element] \[CapitalOmega]];


This code gives the error (apparently because of the wrong way of adding Neumann condition)

LinearSolve::parpiv: Zero pivot was detected during the numerical factorization
or there was a problem in the iterative refinement process.
It is possible that the matrix is ill-conditioned or singular.


I actually have two questions.

1. How to properly implement the boundary condition for a(x,y,t)?

2. When I replace this Neuman condition with the Dirichlet one (which is wrong, but for testing puprposes), the solver starts to work, but then I get another error:

NDSolveValue::ndsz: At t == 1.8768713852597982, step size is effectively zero; singularity or stiff system suspected.

This error goes away when I set all coeffitienst of the equations equal to 1. But could you, please, help me with solving the system with original coeffitients?

P.S. My Mathematica version is 12.1

P.P.S. Sorry for that Greek letters. I tried make them consistent with the original equations. If they really distract you, please, tell, I will rename.

• Initial condition looks very artificial. Can you explain physical background of this problem? Commented Aug 14, 2021 at 1:52
• @AlexTrounev this is the modeling of wound healing. b represents vessel density, n represents density of capillary tip density, and a is the concentration of some chemical substance. The initial conditions was motivated by the fact, that on the edge of the wound b and n should have some values, but there are no vessels or tips inside, so b=0 and n=0. And a=0 is also a resonable assumption. Do you think that the case might be in intial conditions discontinuity? Commented Aug 14, 2021 at 8:44

We can solve this system with using FEM. But we need to make some modification concerning Norm. This code has some message for the last point t=2

 (*region*)\[CapitalOmega] = Rectangle[{0, 0}, {1, 1}];

(*constants*)
\[Mu]n = 0.001;
\[Mu]b = 0.001;
\[Chi] = 0.1;
bHat = 1.5;
\[Lambda]0 = 100;
\[Lambda]1 = 100;
\[Lambda]2 = 10;
\[Lambda]4 = 100;
\[Lambda]5 = 10;
nHat = 1;
\[Lambda]7 = 10;
\[Delta] = 0.01;
\[Alpha] = 2.5;
eps = 10^-5;
aNeuman =
NeumannValue[-\[Lambda]7*bHat*a[t, x, y], {x, y} \[Element]
RegionBoundary[\[CapitalOmega]]];

(*equations*)
eq = {-D[n[t, x, y], t] + \[Mu]n*
Laplacian[n[t, x, y], {x, y}] - \[Chi]*
Div[n[t, x, y]*Grad[a[t, x, y], {x, y}], {x, y}] + \[Lambda]1*
a[t, x, y]*b[t, x, y] - \[Lambda]2*
n[t, x, y] - \[Lambda]0*(n[t, x, y])^2 ==
0, -D[a[t, x, y], t] +
Laplacian[
a[t, x, y], {x, y}] + \[Lambda]4/
2*(1 + Tanh[(1 - b[t, x, y])/\[Delta]]) - (\[Lambda]4 +
b[t, x, y])*a[t, x, y] ==
aNeuman, -D[b[t, x, y], t] + \[Mu]b*
Div[n[t, x, y]*Grad[b[t, x, y], {x, y}], {x, y}] -
Sqrt[(\[Mu]n*Grad[n[t, x, y], {x, y}] - \[Chi]*n[t, x, y]*
Grad[a[t, x, y], {x, y}]) . (\[Mu]n*
Grad[n[t, x, y], {x, y}] - \[Chi]*n[t, x, y]*
Grad[a[t, x, y], {x, y}]) + eps^2] == 0};

(*initial conditions*)
incs = {n[0, x, y] ==
If[x == 0 || y == 0 || x == 1 || y == 1, nHat, 0],
a[0, x, y] == 0,
b[0, x, y] == If[x == 0 || y == 0 || x == 1 || y == 1, bHat, 0]};

(*boundary conditions*)
bcs = {DirichletCondition[n[t, x, y] == nHat*E^(-\[Alpha]*t), True],
DirichletCondition[b[t, x, y] == bHat, True]};

(*solution*)
{nfun, afun, bfun} =
NDSolveValue[
Join[eq, bcs, incs], {n, a, b}, {t, 0,
2}, {x, y} \[Element] \[CapitalOmega],
Method -> {"FiniteElement",
InterpolationOrder -> {n -> 2, a -> 2, b -> 2},
"MeshOptions" -> {"MaxCellMeasure" -> 0.001}},
PrecisionGoal -> 5, AccuracyGoal -> 5] // Quiet;


Visualization

Table[{Row[{"t =", t}],
Plot3D[nfun[t, x, y], {x, y} \[Element] \[CapitalOmega],
PlotRange -> All, AxesLabel -> Automatic, PlotLabel -> "n",
Mesh -> None, PlotPoints -> 50],
Plot3D[afun[t, x, y], {x, y} \[Element] \[CapitalOmega],
PlotRange -> All, AxesLabel -> Automatic, PlotLabel -> "a",
Mesh -> None, PlotPoints -> 50],
Plot3D[bfun[t, x, y], {x, y} \[Element] \[CapitalOmega],
PlotRange -> All, AxesLabel -> Automatic, PlotLabel -> "b",
Mesh -> None, PlotPoints -> 50]}, {t, .1, 1.9, .3}]


• Thank you for answer! Sorry, but there is still this error when I tried your code: NDSolveValue::ndsz: At t == 0.020275884575089772, step size is effectively zero; singularity or stiff system suspected. And what did you mean by the message for the last point? Commented Aug 14, 2021 at 14:20
• This code has been tested with v.12.2 and 12.3. But version 12.1 was mostly unstable (due to covide 19?) and I removed it from my computers. Commented Aug 15, 2021 at 4:17
• Ok, I will install v.12.2 or 12.3. Thank you very much again! Commented Aug 15, 2021 at 10:48
• @АндрейКокорев You are welcome! Commented Aug 15, 2021 at 14:48