# Inactive nonlinear ODE system results in NDSolveValue::femcnmd warning

I am currently trying to solve a set of 1D steady state ODEs using NDSolveValue in the form of Inactive[Grad] and Inactive[Div]. However, I encountered error:

NDSolveValue::femcnmd: The PDE coefficient 0. does not evaluate to a numeric matrix of dimensions {1,1} at the coordinate {0.000025}; it evaluated to 0. instead.

Here is the essential code:

Creating mesh,

ClearAll["Global*"]
k = 1;
m = 5;
d[x_] = 0.5/1000*(1 - Tanh[m*(x + 10)]) +
0.5/1000*(1 - Tanh[-m*(x - 2)]) + 0.00005/k;(*With dip*)
xlist = {0};
xend = xlist[[Length[xlist]]];
coord = {{1}};
i = 1;
While[xend < 4,
i = i + 1;
coord = Join[coord, {{i}}];
xend = xend + d[xend];
If[xend > 4, xend = 4;, xend = xend;];
xlist = Insert[xlist, xend, (Length[xlist] + 1)];
]
(*ReplacePart[xlist,Length[xlist]->4]*)
Needs["NDSolveFEM"]
xmesh = ToElementMesh[Map[{#} &, xlist]];

Defining ODEs and Boundary conditions to solve using NDSolveValue,

testEquations = {-50.5 x^2 a[x] phi[x] (1 - Tanh[20 (-1 + x)]) +
Inactive[
Div][(x^2 a[
x] Inactive[
Grad][(-1 + Log[a[x]] + Tanh[20 (-1 + x)]), {x}]), {x}] ==
NeumannValue[0, x == 0],
50. x^2 a[x] phi[x] (1 - Tanh[20 (-1 + x)]) +
Inactive[
Div][(x^2 b[
x] Inactive[
Grad][(-1 + Log[b[x]] + Tanh[20 (-1 + x)]), {x}]), {x}] ==
NeumannValue[0, x == 0],
0.05 x^2 a[x] phi[x] (1 - Tanh[20 (-1 + x)]) +
0.001 Inactive[
Log[10. (1 - phi[x])], {x}]), {x}] ==
NeumannValue[0, x == 0], a[4] == 1, b[4] == 0, phi[4] == 1.};
testSolution =
NDSolveValue[testEquations , {a, b, phi}, Element[x, xmesh]];

Apart from that, is it possible to perform divergence on a 1D vector in spherical coordinate for my ODEs (to become equivalent to )? Currently, I have to manually converting to this using cartesian form by performing multiplication of r^2 into the vector.

• Can we avoid 1/0 error when b[4]==0 and in equations there is 1/b[x] ? Jul 18, 2022 at 7:49
• @AlexTrounev Do you mean that when we expand the divergence term in second ODE? Jul 18, 2022 at 8:12
• Yes it is, second equation looks like 50 x^2 a[x] phi[x] (1 - Tanh[20 (-1 + x)]) + 2 x b[x] (20 Sech[20 (-1 + x)]^2 + Derivative[1][b][x]/b[x]) + x^2 Derivative[1][b][ x] (20 Sech[20 (-1 + x)]^2 + Derivative[1][b][x]/b[x]) + x^2 b[x] (-800 Sech[20 (-1 + x)]^2 Tanh[20 (-1 + x)] - Derivative[1][b][x]^2/b[x]^2 + (b^\[Prime]\[Prime])[x]/b[x]) == 0 . There are several terms with 1/b[x], 1/b[x]^2 while b[4]==0. Jul 18, 2022 at 12:31
• This system can be solved with collocation method - see my answer. Jul 19, 2022 at 4:49

"I want to solve a differential equation system. I have a quick search, it seems that everybody is talking about a powerful method called FEM, so I decide to use it, too." ← If this is the reason you keep trying solving your nonlinear ODE system with FiniteElement method, please remember, there is no silver bullet. Though it's hot these days, finite element method has its own drawbacks, and you've encountered the most troublesome one (in my view).

The warning can be reproduced with the following sample (you should have put more effort in simplifying your code sample BTW):

eqsample = 1. Inactive[Div][Inactive[Grad][1. phi[x] + 1, {x}], {x}] == 0;
NDSolveValue[{eqsample, phi[0] == 0, phi[4] == 1}, phi, x ∈ xmesh]

NDSolveValue::femcnmd: The PDE coefficient 0. does not evaluate to a numeric matrix of dimensions {1,1} at the coordinate {0.000025}; it evaluated to 0. instead.

Definition of xmesh is the same as yours. I use it to reproduce exactly the same warning as mentioned in the question. If you use {x, 0, 4} instead, the warning will be

NDSolveValue::femcnmd: The PDE coefficient 0. does not evaluate to a numeric matrix of dimensions {1,1} at the coordinate {2.}; it evaluated to 0. instead.

NDSolve shouldn't have any difficulty in handling eqsample, it's just a simple phi''[x] == 0:

eqsample // Activate
(* 1. phi''[x] == 0 *)

Further check proves something is wrong in parsing stage:

NDSolveProcessEquations[{eqsample, phi[0] == 0, phi[4] == 1}, phi, {x, 0, 4}]

NDSolveValue::femcnmd ……

Why does this happen? It's because that, in order to solve differential equation, NDSolve always need to transform the equation to certain standard form. Different standard forms are required by different methods. (For example, standard forms of initial value problem of ODE is discussed in this post. ) FiniteElement method is bounded in the most complicated and restricted standard form (which is called formal form or inactive form in document) of NDSolve. This makes the pre-processing of FiniteElement method rather hard and error-prone. The topic is non-trivial so I'd like not to talk too much about it in this post, please read through Formal Partial Differential Equations section in the document for more info. In short:

1. Simply re-write the equation with Inactive[Grad] and Inactive[Div] is often insufficient.

2. Personally I suggest one always check if the equation is parsed properly with NDSolveFEMGetInactivePDE.

"OK, then how to help FiniteElement parsing the equation system? " I don't know, and I don't bother to. The biggest advantage of FiniteElement (in my view) is its capability for solving PDE(s) on irregular domain, but you're just dealing with an ODE system! Why not turning to the old-good finite difference method (FDM)?

I'll use pdetoae for the task. Definition of testEquations is the same as yours.

{eq, bc} = {testEquations[[;; -4]], {testEquations[[-3 ;;]],
phi'[0] == 0, a'[0] == 0, b'[0] == 0}} /. NeumannValue[__] :> 0 // Activate;

domain = {0, 4};
points = 500;
grid = Array[#1 &, points, domain];
difforder = 4;
(* Definition of pdetoae isn't included in this post,
ptoafunc = pdetoae[{a, b, phi}[x], grid, difforder];
ae = (#1[[2 ;; -2]] &) /@ ptoafunc[eq];
aebc = ptoafunc[bc];
icguess[var_][x_] := 1;
solrule = FindRoot[
Flatten[{ae, aebc}], (Flatten[#1, 1] &)[
Table[{var[x], icguess[var][x]}, {var, {a, b, phi}}, {x, grid}]],
Jacobian -> "FiniteDifference"]; // AbsoluteTiming
(* {0.990207, Null} *)
{afunc, bfunc, phifunc} =
ListInterpolation[#, grid] & /@ Partition[solrule[[All, -1]], points];
ListLinePlot[{afunc, bfunc, phifunc}, PlotRange -> All]

• This is a nice solution (+1). Why it is not correlated with wavelets method? Jul 19, 2022 at 9:10
• @AlexTrounev I'm not knowledgable enough to explain this, but the convergency of the discretized system doesn't seem to be fast, for example, if I decrease the difference order to difforder = 2;, the numeric error is still obvious even for points = 800;. Might be related to the sudden change around x == 1? Jul 19, 2022 at 9:29
• Ok! Finally I got your solution with 512 collocation points - see update to my answer. But my code is terribly slow compare to yours. Jul 19, 2022 at 10:12
• @Alex Adding Jacobian -> "FiniteDifference" to FindRoot speeds up your code. With this option, it takes about 90 seconds on my laptop. Without it, the calculation doesn't finish for several minutes. (We also need to add Exclusions -> None to Plot in this case. ) Jul 19, 2022 at 13:02
• Thank you very much, this option really works. Jul 19, 2022 at 13:32

This problem can be solved with using the Euler wavelets collocation method described here, here, and hear. First we rationalize coefficients and map solution on the unit interval, then we have

Clear["Global*"]

testEquations = {-505/10 L0^2 x^2 a[x] phi[
x] (1 - Tanh[20 (-1 + L0 x)]) +
D[(x^2 a[
x] D[(-1 + Log[a[x]] + Tanh[20 (-1 + L0 x)]), {x}]), {x}] ==
0, 50 L0^2 x^2 a[x] phi[x] (1 - Tanh[20 (-1 + L0 x)]) +
D[x^2 b[x] D[(-1 + Log[b[x]] +
Tanh[20 (-1 + L0 x)]), {x}], {x}] == 0,
5/100  L0^2 x^2 a[x] phi[x] (1 - Tanh[20 (-1 + L0 x)]) +
1/1000 D[(x^2 (1 - phi[x]) D[Log[10 (1 - phi[x])], {x}]), {x}] ==
0};

Let transform system to the collocation method

testEquations /. {a'[x] -> a1[x], a''[x] -> a2[x],
b'[x] -> b1[x], b''[x] -> b2[x], phi'[x] -> phi1[x],
phi''[x] -> phi2[x]}

Out[]= {2 x a[x] (a1[x]/a[x] + 20 L0 Sech[20 (-1 + L0 x)]^2) +
x^2 a1[x] (a1[x]/a[x] + 20 L0 Sech[20 (-1 + L0 x)]^2) -
101/2 L0^2 x^2 a[x] phi[x] (1 - Tanh[20 (-1 + L0 x)]) +
x^2 a[x] (-(a1[x]^2/a[x]^2) + a2[x]/a[x] -
800 L0^2 Sech[20 (-1 + L0 x)]^2 Tanh[20 (-1 + L0 x)]) == 0,
2 x b[x] (b1[x]/b[x] + 20 L0 Sech[20 (-1 + L0 x)]^2) +
x^2 b1[x] (b1[x]/b[x] + 20 L0 Sech[20 (-1 + L0 x)]^2) +
50 L0^2 x^2 a[x] phi[x] (1 - Tanh[20 (-1 + L0 x)]) +
x^2 b[x] (-(b1[x]^2/b[x]^2) + b2[x]/b[x] -
800 L0^2 Sech[20 (-1 + L0 x)]^2 Tanh[20 (-1 + L0 x)]) ==
0, (-2 x phi1[x] - x^2 phi2[x])/1000 +
1/20 L0^2 x^2 a[x] phi[x] (1 - Tanh[20 (-1 + L0 x)]) == 0}

This system we transform to the algebraic system as follows

UE[m_, t_] := EulerE[m, t]
psi[k_, n_, m_, t_] :=
Piecewise[{{2^(k/2) Sqrt[2/Pi] UE[m, 2^k t - 2 n + 1], (n - 1)/
2^(k - 1) <= t < n/2^(k - 1)}, {0, True}}]
PsiE[k_, M_, t_] :=
Flatten[Table[psi[k, n, m, t], {n, 1, 2^(k - 1)}, {m, 0, M - 1}]]
k0 = 7; M0 = 8; With[{k = k0, M = M0},
nn = Length[Flatten[Table[1, {n, 1, 2^(k - 1)}, {m, 0, M - 1}]]]];
dx = 1/(nn);  xl = Table[ l*dx, {l, 0, nn}];
xcol = Table[(xl[[l - 1]] + xl[[l]])/2, {l, 2, nn + 1}]; Psijk =
With[{k = k0, M = M0}, PsiE[k, M, t1]]; Int1 =
With[{k = k0, M = M0}, Integrate[PsiE[k, M, t1], t1]];
Int2 = Integrate[Int1, t1]; Psi[y_] := Psijk /. t1 -> y;
int1[y_] := Int1 /. t1 -> y; int2[y_] := Int2 /. t1 -> y;

var1 = Array[v1, {nn}]; var2 = Array[v2, {nn}]; var3 =
Array[v3, {nn}]; con = Array[c, {6}];
L0 = 4; a2[x_] = var1 . Psi[x]; a1[x_] = var1 . int1[x] + c[1];
a[x_] = var1 . int2[x] + c[1] x + c[2];
b2[x_] = var2 . Psi[x]; b1[x_] = var2 . int1[x] + c[3];
b[x_] = var2 . int2[x] + c[3] x + c[4];
phi2[x_] = var3 . Psi[x]; phi1[x_] = var3 . int1[x] + c[5];
phi[x_] = var3 . int2[x] + c[5] x + c[6]; var =
Join[var1, var2, var3, con];

eq = Flatten[
Table[{2 x a[x] (a1[x]/a[x] + 20 L0 Sech[20 (-1 + L0 x)]^2) +
x^2 a1[x] (a1[x]/a[x] + 20 L0 Sech[20 (-1 + L0 x)]^2) -
101/2 L0^2 x^2 a[x] phi[x] (1 - Tanh[20 (-1 + L0 x)]) +
x^2 a[x] (-(a1[x]^2/a[x]^2) + a2[x]/a[x] -
800 L0^2 Sech[20 (-1 + L0 x)]^2 Tanh[20 (-1 + L0 x)]) == 0,
2 x b[x] (b1[x]/b[x] + 20 L0 Sech[20 (-1 + L0 x)]^2) +
x^2 b1[x] (b1[x]/b[x] + 20 L0 Sech[20 (-1 + L0 x)]^2) +
50 L0^2 x^2 a[x] phi[x] (1 - Tanh[20 (-1 + L0 x)]) +
x^2 b[x] (-(b1[x]^2/b[x]^2) + b2[x]/b[x] -
800 L0^2 Sech[20 (-1 + L0 x)]^2 Tanh[20 (-1 + L0 x)]) ==
0, (-2 x phi1[x] - x^2 phi2[x])/1000 +
1/20 L0^2 x^2 a[x] phi[x] (1 - Tanh[20 (-1 + L0 x)]) == 0}, {x,
xcol}]];

Finally we add boundary conditions and solve with FindRoot

bc = {a[x] == 1, b[x] == 0, phi[x] == 1} /.
x -> 1; bc1 = {a1[x] == 0, b1[x] == 0, phi1[x] == 0} /. x -> 0;

eqs = Join[eq, bc, bc1];

sol = FindRoot[eqs, Table[{var[[i]], 1/10}, {i, Length[var]}],Jacobian -> "FiniteDifference"];

Visualization

Plot[Evaluate[{a[x/L0] /. sol, b[x/L0] /. sol, phi[x/L0] /. sol}], {x,
0, L0}, PlotLegends -> {"a", "b", "phi"}, Exclusions -> None]

We also can solve this problem with linear FEM, using iterative false transient algorithm discussed here, here, and here. We use transformed system of equations with adding iteration parameter dt=1/5 and linearized as follows

Clear["Global*"]

Needs["NDSolveFEM"]
xmesh = ToElementMesh[ImplicitRegion[0 <= x <= 1, {x}],
MaxCellMeasure -> 2 10^-3]
eqn = {-(a[x] - A[i - 1][x])/dt -
808  a[x] P[i - 1][x] (1 - Tanh[20 (-1 + 4 x)]) +
2/x  (80 Sech[20 (-1 + 4 x)]^2 a[x] +
Derivative[1][a][x]) +  (80 Sech[20 (-1 + 4 x)]^2 Derivative[
1][a][x] +
A1[i - 1][x] a'[x]/A[i - 1][x]) +  (-12800 Sech[
20 (-1 + 4 x)]^2 Tanh[20 (-1 + 4 x)] a[x] -
A1[i - 1][x] a'[x]/A[i - 1][x] + (a^\[Prime]\[Prime])[x]) ==
NeumannValue[0, x == 0], -(b[x] - B[i - 1][x])/dt +
800  P[i - 1][x] a[x] (1 - Tanh[20 (-1 + 4 x)]) +
2/x  (80 Sech[20 (-1 + 4 x)]^2 b[x] +
Derivative[1][b][x]) + (80 Sech[20 (-1 + 4 x)]^2 Derivative[1][
b][x] + B1[i - 1][x] b'[x]/B[i - 1][x]) + (-12800 Sech[
20 (-1 + 4 x)]^2 Tanh[20 (-1 + 4 x)] b[x] -
B1[i - 1][x] b'[x]/B[i - 1][x] + (b^\[Prime]\[Prime])[x]) ==
NeumannValue[0, x == 0], (phi[x] - P[i - 1][x])/dt +
4/5  P[i - 1][x] a[
x] (1 - Tanh[20 (-1 + 4 x)]) + (-2/x  Derivative[1][phi][x] - (
phi^\[Prime]\[Prime])[x])/1000 == NeumannValue[0, x == 0]}; bc =
DirichletCondition[{a[x] == 1, b[x] == 0, phi[x] == 1}, x == 1];

Initial guess is very important in this case, but we use the simple one

A[0][x_] := 1; B[0][x_] = 1/100; P[0][x_] = 1; A1[0][x_] := 0;
B1[0][x_] = 0; P1[0][x_] := 0;

L0 = 4; dt =
1/5; nn = 231; Do[{A[i], B[i], P[i], A1[i], B1[i], P1[i]} =
NDSolveValue[{eqn, bc}, {a, b, phi, a', b', phi'},
Element[{x}, xmesh]]; , {i, 1, nn}]

Visualization

Table[Plot[Evaluate[{A[i][x/L0], B[i][x/L0], P[i][x/L0]}], {x, 0, L0},
PlotLegends -> {"a", "b", "phi"}, PlotRange -> All,
PlotLabel -> i], {i, 228, 231}]

Also FEM solution looks similar to above computed with wavelets and computed with pdetoae at xzczd post, the algorithm not stable and solution diverges at nn>231.

• There seem to be simple mistakes in mapping, for example, D[(x^2 a[x] D[… should be D[((L0 x)^2 a[x] D[…. I suggest using DSolveChangeVariables or DChange to avoid possible mistakes. Jul 19, 2022 at 7:40
• Actually there is no mistake in mapping, since we can test my testEquations with bc, bc1 with your code and got absolutely same result as yours. But solution strictly depends on the number of collocation points. May be this is a problem? Jul 19, 2022 at 8:43
• Oh you're right, I forgot about the L0 generated by D… With M0 = 9 our results (almost) matches :) . Jul 19, 2022 at 9:08
• Yes, but in my example there is 98 points and I have tested it with 112 and 128 points as well. May be there are several branches in numerical solution? Jul 19, 2022 at 9:14
• If the system owns multiple solutions, according to my (limited) experience, the initial guess should play an important role, but this problem doesn't seem to be the case… Jul 19, 2022 at 9:34