{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]
"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) in 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:
Simply re-write the equation with
Inactive[Grad]andInactive[Div]is often inefficient.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}]]];
{afunc, bfunc, phifunc} =
ListInterpolation[#, grid] & /@ Partition[solrule[[All, -1]], points];
ListLinePlot[{afunc, bfunc, phifunc}, PlotRange -> All]
