# NDSolve error: CoefficientArray:... is not a polynomial

I get strange errors (CoefficientArrays::poly and NDSolveValue::femper) with NDSolveValue: it seems it is not parsing the equations in the corrrect way, but I cannot understand what's wrong.

Here is the code:

v0=-10^-2;
den[x_, y_, z_] := 10;
p[x_, y_, z_] := 1;
h[x_, y_, z_] := 5;
u0[x_, y_, z_] :=
Exp[ϕ[x, y,
z]] Sqrt@(1 + (v0^2 x^2)/
Sqrt[x^2 + y^2 + z^2]^2 Exp[-2 α1[x, y, z]] + (
v0^2 y^2)/Sqrt[x^2 + y^2 + z^2]^2 Exp[-2 β[x, y, z]] + (
v0^2 z^2)/Sqrt[x^2 + y^2 + z^2]^2 Exp[-2 γ[x, y, z]]);
u1[x_, y_, z_] := (v0 x)/Sqrt[x^2 + y^2 + z^2];
u2[x_, y_, z_] := (v0 y)/Sqrt[x^2 + y^2 + z^2];
u3[x_, y_, z_] := (v0 z)/Sqrt[x^2 + y^2 + z^2];

eqn1 = (-Exp[2*γ[x, y, z]]*D[α1[x, y, z], z] -
Exp[2*γ[x, y, z]]*D[β[x, y, z], z]^2 -
Exp[(2*γ[x, y, z])]*
D[α1[x, y, z],
z]*(D[β[x, y, z], z] - D[γ[x, y, z], z]) +
Exp[2*γ[x, y, z]]*D[β[x, y, z], z]*
D[γ[x, y, z], z] +
Exp[2*γ[x, y, z]]*D[α1[x, y, z], {z, 2}] +
Exp[2*γ[x, y, z]]*D[β[x, y, z], {z, 2}] -
Exp[2*β[x, y, z]]*D[α1[x, y, z], y]^2 +
Exp[2*β[x, y, z]]*D[α1[x, y, z], y]*
D[β[x, y, z], y] -
Exp[2*β[x, y, z]]*D[α1[x, y, z], y]*
D[γ[x, y, z], y] +
Exp[2*β[x, y, z]]*D[β[x, y, z], y]*
D[γ[x, y, z], y] -
Exp[2*β[x, y, z]]*D[γ[x, y, z], y]^2 +
Exp[2*β[x, y, z]]*D[α1[x, y, z], {y, 2}] +
Exp[2*β[x, y, z]]*D[γ[x, y, z], {y, 2}] +
Exp[2*α1[x, y, z]]*D[α1[x, y, z], x]*
D[β[x, y, z], x] -
Exp[2*α1[x, y, z]]*D[β[x, y, z], x]^2 +
Exp[2*α1[x, y, z]]*D[α1[x, y, z], x]*
D[γ[x, y, z], x] -
Exp[2*α1[x, y, z]]*D[β[x, y, z], x]*
D[γ[x, y, z], x] -
Exp[2*α1[x, y, z]]*D[γ[x, y, z], x]^2 +
Exp[2*α1[x, y, z]]*D[β[x, y, z], {x, 2}] +
Exp[2*α1[x, y, z]]*D[γ[x, y, z], {x, 2}])/
Exp[2*ϕ[x, y, z]] == (den[x, y, z] h[x, y, z] -
p[x, y, z]) u0[x, y, z]^2 - p[x, y, z] Exp[2 ϕ[x, y, z]];

eqn2 = (Exp[2*γ[x, y, z]]*D[β[x, y, z], z]^2 -
Exp[2*γ[x, y, z]]*
D[β[x, y, z],
z]*(D[γ[x, y, z], z] - D[ϕ[x, y, z], z]) -
Exp[2*γ[x, y, z]]*D[γ[x, y, z], z]*
D[ϕ[x, y, z], z] +
Exp[2*γ[x, y, z]]*D[ϕ[x, y, z], z]^2 -
Exp[2*γ[x, y, z]]*D[β[x, y, z], {z, 2}] -
Exp[2*γ[x, y, z]]*D[ϕ[x, y, z], {z, 2}] -
Exp[2*β[x, y, z]]*D[β[x, y, z], y]*
D[γ[x, y, z], y] +
Exp[2*β[x, y, z]]*D[γ[x, y, z], y]^2 -
Exp[2*β[x, y, z]]*D[β[x, y, z], y]*
D[ϕ[x, y, z], y] +
Exp[2*β[x, y, z]]*D[γ[x, y, z], y]*
D[ϕ[x, y, z], y] +
Exp[2*β[x, y, z]]*D[ϕ[x, y, z], y]^2 -
Exp[2*β[x, y, z]]*D[γ[x, y, z], {y, 2}] -
Exp[2*β[x, y, z]]*D[ϕ[x, y, z], {y, 2}] +
Exp[2*α1[x, y, z]]*D[β[x, y, z], x]*
D[γ[x, y, z], x] +
Exp[2*α1[x, y, z]]*D[β[x, y, z], x]*
D[ϕ[x, y, z], x] +
Exp[2*α1[x, y, z]]*D[γ[x, y, z], x]*
D[ϕ[x, y, z], x])/
Exp[2*α1[x, y, z]] == (den[x, y, z] h[x, y, z] -
p[x, y, z]) u1[x, y, z]^2 +
p[x, y, z] Exp[2 α1[x, y, z]];

eqn3 = (Exp[2*γ[x, y, z]]*D[α1[x, y, z], z]^2 -
Exp[2*γ[x, y, z]]*
D[α1[x, y, z],
z]*(D[γ[x, y, z], z] - D[ϕ[x, y, z], z]) -
Exp[2*γ[x, y, z]]*D[γ[x, y, z], z]*
D[ϕ[x, y, z], z] +
Exp[2*γ[x, y, z]]*D[ϕ[x, y, z], z]^2 -
Exp[2*γ[x, y, z]]*D[α1[x, y, z], {z, 2}] -
Exp[2*γ[x, y, z]]*D[ϕ[x, y, z], {z, 2}] +
Exp[2*β[x, y, z]]*D[α1[x, y, z], y]*
D[γ[x, y, z], y] +
Exp[2*β[x, y, z]]*D[α1[x, y, z], y]*
D[ϕ[x, y, z], y] +
Exp[2*β[x, y, z]]*D[γ[x, y, z], y]*
D[ϕ[x, y, z], y] -
Exp[2*α1[x, y, z]]*D[α1[x, y, z], x]*
D[γ[x, y, z], x] +
Exp[2*α1[x, y, z]]*D[γ[x, y, z], x]^2 -
Exp[2*α1[x, y, z]]*D[α1[x, y, z], x]*
D[ϕ[x, y, z], x] +
Exp[2*α1[x, y, z]]*D[γ[x, y, z], x]*
D[ϕ[x, y, z], x] +
Exp[2*α1[x, y, z]]*D[ϕ[x, y, z], x]^2 -
Exp[2*α1[x, y, z]]*D[γ[x, y, z], {x, 2}] -
Exp[2*α1[x, y, z]]*D[ϕ[x, y, z], {x, 2}])/
Exp[2*β[x, y, z]] == (den[x, y, z] h[x, y, z] -
p[x, y, z]) u2[x, y, z]^2 + p[x, y, z] Exp[2 β[x, y, z]];

eqn4 = (Exp[2*γ[x, y, z]]*D[β[x, y, z], z]*
D[ϕ[x, y, z], z] +
Exp[2*γ[x, y, z]]*
D[α1[x, y, z],
z]*(D[β[x, y, z], z] + D[ϕ[x, y, z], z]) +
Exp[2*β[x, y, z]]*D[α1[x, y, z], y]^2 -
Exp[2*β[x, y, z]]*D[α1[x, y, z], y]*
D[β[x, y, z], y] +
Exp[2*β[x, y, z]]*D[α1[x, y, z], y]*
D[ϕ[x, y, z], y] -
Exp[2*β[x, y, z]]*D[β[x, y, z], y]*
D[ϕ[x, y, z], y] +
Exp[2*β[x, y, z]]*D[ϕ[x, y, z], y]^2 -
Exp[2*β[x, y, z]]*D[α1[x, y, z], {y, 2}] -
Exp[2*β[x, y, z]]*D[ϕ[x, y, z], {y, 2}] -
Exp[2*α1[x, y, z]]*D[α1[x, y, z], x]*
D[β[x, y, z], x] +
Exp[2*α1[x, y, z]]*D[β[x, y, z], x]^2 -
Exp[2*α1[x, y, z]]*D[α1[x, y, z], x]*
D[ϕ[x, y, z], x] +
Exp[2*α1[x, y, z]]*D[β[x, y, z], x]*
D[ϕ[x, y, z], x] +
Exp[2*α1[x, y, z]]*D[ϕ[x, y, z], x]^2 -
Exp[2*α1[x, y, z]]*D[β[x, y, z], {x, 2}] -
Exp[2*α1[x, y, z]]*D[ϕ[x, y, z], {x, 2}])/
Exp[2*γ[x, y, z]] == (den[x, y, z] h[x, y, z] -
p[x, y, z]) u3[x, y, z]^2 +
p[x, y, z] Exp[2 γ[x, y, z]];


Here the NDSolveValue step:

NDSolveValue[{eqn1, eqn2, eqn3, eqn4,
DirichletCondition[ϕ[x, y, z] == 0, True],
DirichletCondition[α1[x, y, z] == 0, True],
DirichletCondition[β[x, y, z] == 0, True],
DirichletCondition[γ[x, y, z] == 0,
True]}, {ϕ, α1, β, γ}, {x, -25,
25}, {y, -25, 25}, {z, -25, 25}]


And here is the error:

I am sorry for the long code, but I couldn't find a MWE.

• If I evaluate your code on a fresh kernel (V12.3 Windows), I get – after a couple of minutes – a different error. Jan 9 at 17:38
• @Domen I use version 11.3 for Linux Jan 9 at 17:46
• Your system looks nonlinear. FEM was enhance to deal with nonlinear systems in V12.0, I think. It'd be nice to add the message names so that others can search for the error names to find answers that might help them. I'm getting the same as @Domen in V13, Mac. The kernel did have a "virtual memory size" of over 400GB, though, and "real memory size" of 9GB. Jan 9 at 19:21
• @MichaelE2 the error names are CoefficientArrays::poly and NDSolveValue::femper. But would it be possible to solve the system with another method? Jan 9 at 19:26
• I think "MethodOfLines" is the only other in NDSolve, but you'd have to reformulate the BCs somewhat, I think. You have to have a variable for the time integration. Others have implemented FDM on this site; you can search for "FDM". I don't do a lot of work with PDEs, so my experience is not extensive. Jan 9 at 19:46

The message you see is because the version of Mathematica that you use does not have the nonlinear finite element solver implemented. When I run your system in Version 13.0 with a coarser then default mesh:

mesh = NDSolveFEMToElementMesh[Cuboid[-25*{1, 1, 1}, 25*{1, 1, 1}],
"MaxCellMeasure" -> 50]

NDSolveValue[{eqn1, eqn2, eqn3, eqn4,
DirichletCondition[\[Phi][x, y, z] == 0, True],
DirichletCondition[\[Alpha]1[x, y, z] == 0, True],
DirichletCondition[\[Beta][x, y, z] == 0, True],
DirichletCondition[\[Gamma][x, y, z] == 0,
True]}, {\[Phi], \[Alpha]1, \[Beta], \[Gamma]}, {x, y,
z} \[Element] mesh]


Which means that the solver could not find a solution. That can have various causes. The PDE is not correct, a solution does not exists, the initial seeding of 0 is not good, the solver can not find a solution.... I'd first double check that the PDE is what you want it to be.