# LinearSolve problem occurs as solving PDEs with NDSolve

I am new to Mathematica, please forgive me asking naive questions. I tried to solve PDEs numerically using NDSolve, but failed to go through due to errors. Two of three PDEs are time-dependent and another one is time-independent (or implicitly time-dependent). Here is the code:

Needs["NDSolveFEM"]
dr = 0.5; rx = 2.5/2; xs = -(dr + rx); ys = 0.0;
ec = RegionUnion[Disk[{-rx, 0}, dr], Rectangle[{-rx, -dr}, {rx, dr}],
Disk[{rx, 0}, dr]];
Region[ec]
bmesh = ToBoundaryMesh[ec, "MaxBoundaryCellMeasure" -> .05];
bmesh["Wireframe"]
f = Function[{vertices, area},
area > 0.004 (0.9 - 0.5 Norm[Mean[vertices]])];
(mesh = ToElementMesh[bmesh, MeshRefinementFunction -> f])["Wireframe"]
u0 = 11.2; v0 = 1; d0 = 0.97; V2 = -0.2;
Eqs = {D[u[t, x, y], t] ==
d0 (Laplacian[u[t, x, y], {x, y}] + w[t, x, y]),
D[v[t, x, y], t] == Laplacian[v[t, x, y], {x, y}] - w[t, x, y],
Laplacian[
w[t, x, y], {x, y}] == -4 \[Pi] (v[t, x, y] - u[t, x, y])};
ics = {u[0, x, y] == u0, v[0, x, y] == v0, w[0, x, y] == -0.0};
bcs = DirichletCondition[w[t, x, y] == V2, {x, y} \[Element] bmesh];
{usol, vsol, wsol} =
NDSolveValue[{Eqs, ics, bcs}, {u, v, w}, {t, 0,
10^6}, {x, y} \[Element] mesh];


The result shows errors

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.
Set::shape: Lists {usol,vsol,wsol} and NDSolveValue[<<1>>] are not the same shape.
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.


Totally have no ideas where errors come from. Thanks for any instructions and suggestions.

• For the second error message, I remove the {usol,vsol,wsol}= and the ; from your last line of code so I can see what your NDSolveValue returns. That looks like the NDSolveValue failed and just returned the original line. So that error message will likely go away if you can fix your "Zero pivot" error. If I use the little search box at the top of the Stackexchange page and enter "zero pivot" with those quotes then I get lots of results and answers related to that. Can you scan through those and see if any look like they might be helpful?
– Bill
Dec 25, 2023 at 3:52
• I took the liberty of taking a modified version of this example and added it to the parpiv message ref page. I hope that's OK? Jan 1 at 15:16

There are several problems discussed here and here, and we can resolve these problems by adding small regularization term to the last equation as follows

Needs["NDSolveFEM"]
dr = 0.5; rx = 2.5/2; xs = -(dr + rx); ys = 0.0;
ec = RegionUnion[Disk[{-rx, 0}, dr], Rectangle[{-rx, -dr}, {rx, dr}],
Disk[{rx, 0}, dr]];
Region[ec]
bmesh = ToBoundaryMesh[ec, "MaxBoundaryCellMeasure" -> .05];
bmesh["Wireframe"]
f = Function[{vertices, area},
area > 0.004 (0.9 - 0.5 Norm[Mean[vertices]])];
(mesh = ToElementMesh[bmesh, MeshRefinementFunction -> f])["Wireframe"]
u0 = 11.2; v0 = 1; d0 = 0.97; V2 = -0.2;
Eqs = Inactivate[{D[u[t, x, y], t] -
d0 (Laplacian[u[t, x, y], {x, y}] + w[t, x, y]),
D[v[t, x, y], t] - Laplacian[v[t, x, y], {x, y}] + w[t, x, y],
0.001 D[w[t, x, y], t] - Laplacian[w[t, x, y], {x, y}] -
4 \[Pi] (v[t, x, y] - u[t, x, y])}, Laplacian | D];
ics = {u[0, x, y] == u0, v[0, x, y] == v0, w[0, x, y] == 0.0};
bcs = DirichletCondition[w[t, x, y] == V2, True];
{usol, vsol, wsol} =
NDSolveValue[{Activate[Eqs] == {NeumannValue[0, True],
NeumannValue[0, True], 0}, ics, bcs}, {u, v,
w}, {x, y} \[Element] mesh, {t, 0, 10}];


Note, in this case we don't need to integrate up to t=10^6, even t=5 could be fine. Visualization

Plot[{usol[t, 0, 0], vsol[t, 0, 0], wsol[t, 0, 0]}, {t, 0, 10},
PlotLegends -> {u, v, w}, Frame -> True, FrameLabel -> Automatic,
PlotRange -> All]


{DensityPlot[usol[10, x, y], {x, y} \[Element] mesh,
ColorFunction -> "Rainbow", PlotLegends -> Automatic,
AspectRatio -> Automatic],
DensityPlot[vsol[10, x, y], {x, y} \[Element] mesh,
ColorFunction -> "Rainbow", PlotLegends -> Automatic,
AspectRatio -> Automatic],
DensityPlot[wsol[10, x, y], {x, y} \[Element] mesh,
ColorFunction -> "Rainbow", PlotLegends -> Automatic,
AspectRatio -> Automatic]}


• Thanks to Alex's help and let me learn how to solve PDEs numerically by regularization. Dec 27, 2023 at 9:06
• (+1), but you do not need the NeumannValue and probably you can also get rid of the Inactivate and Activate a few lines later. Jan 1 at 15:17
• @user21 Thank you. Anyway this answer looks as not finished since real system is nonlinear. Jan 1 at 16:55

Thanks to Alex's prompt answer and instruction about using regularization term. The PDEs in the last question is a simplified set of equations for pilot study because I encountered the same problem in the PDEs to be solved. Then I tried the method of regularization to the following:

Needs["NDSolveFEM"]
dr = 0.5; rx = 2.5/2; xs = -(dr + rx); ys = 0.0;
ec = RegionUnion[Disk[{-rx, 0}, dr], Rectangle[{-rx, -dr}, {rx, dr}],
Disk[{rx, 0}, dr]];
Region[ec]
bmesh = ToBoundaryMesh[ec, "MaxBoundaryCellMeasure" -> .05];
bmesh["Wireframe"]
f = Function[{vertices, area},
area > 0.004 (0.9 - 0.5 Norm[Mean[vertices]])];
(mesh = ToElementMesh[bmesh, MeshRefinementFunction -> f])["Wireframe"]
u0 = 11.2; v0 = 1; d0 = 0.97; V2 = -0.2;
Eqs = Inactivate[{D[u[t, x, y], t] -
d0 (Laplacian[u[t, x, y], {x, y}] +
Div[(u[t, x, y]) Grad[w[t, x, y], {x, y}], {x, y}]),
D[v[t, x, y], t] - Laplacian[v[t, x, y], {x, y}] -
Div[(v[t, x, y]) Grad[w[t, x, y], {x, y}], {x, y}],
0.001 D[w[t, x, y], t] - Laplacian[w[t, x, y], {x, y}] -
4 \[Pi] (v[t, x, y] - u[t, x, y])}, Laplacian | D];
ics = {u[0, x, y] == u0, v[0, x, y] == v0, w[0, x, y] == 0.0};
bcs = DirichletCondition[w[t, x, y] == V2, True];
{usol, vsol, wsol} =
NDSolveValue[{Activate[Eqs] == {NeumannValue[0, True],
NeumannValue[0, True], 0}, ics, bcs}, {u, v,
w}, {x, y} \[Element] mesh, {t, 0, 10}];


You may find out the differences are the divergence of u or v times the gradient of w following the Laplacian term in the 1st and 2nd PDEs. Small regularization term is applied to the 3rd PDE as the prior setting. Soon I encounter the error as:

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


I serached stackexchange to see if it could be solved and modified NDSolve as

  {usol, vsol, wsol} =
NDSolveValue[{Activate[Eqs] == {NeumannValue[0, True],
NeumannValue[0, True], 0}, ics, bcs}, {u, v,
w}, {x, y} \[Element] mesh, {t, 0, 1},
Method -> {"FiniteElement",
"InterpolationOrder" -> {u -> 2, v -> 2, w -> 1},
"MeshOptions" -> {"MaxCellMeasure" -> 0.0005}}];
`

However, it ran over a day and not finished yet. I guess it's not the best way to remove sigularity or stiffness here. I'd like to learn form you or anone if any approach to solve this type of problem? Thanks for any help and instruction.

• You should add this information to the question. That's better than posting it as an answer which it is not. You could also make a follow-up question. Dec 27, 2023 at 9:15
• Sorry that I don't know how to make a follow-up question. I just notice there is an "Edit" buttom after the question. Is that the way to make a follow-up question? Thanks. Dec 27, 2023 at 9:31
• @dopey The system you are considering here is very different from the system you considered in your main question. It could be better to start a new topic, as suggested above by user21. Dec 28, 2023 at 2:25