Inactive nonlinear ODE system results in NDSolveValue::femcnmd warning

Question

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["NDSolve`FEM`"]
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[
        Div][(x^2 (1 - phi[x]) Inactive[Grad][
          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.

Answer 1

"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:

NDSolve`ProcessEquations[{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 NDSolve`FEM`GetInactivePDE.

"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,
   please find it in the link above. *)
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]

Answer 2

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["NDSolve`FEM`"]
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.