# The PDE coefficient does not evaluate to a numeric scalar....at the coordinate...{1/2,1/2}; it evaluated to {-0.1} instead. What's wrong in my code?

I have a nonlinear system of two PDEs both of first order to be solved on the rectangle region (0,1)x(0,1). The unknowns must satisfy Dirichlet conditions. The structure of the problem is simple but I always get errors. May be I should need some sophisticated method for NDSolve, but I don't know which one could work better than others. The kind of error is that in the subject, that is "The PDE coefficient does not evaluate to a numeric scalar....at the coordinate...{1/2,1/2}; it evaluated to {-0.1} instead." I know that this kind of error is generated by possible singular points in the coefficients, but I dont see the reason in my case: the in the denominator in the R.H.S. is always strictly positive. Any help would be greatly appreciated. Thanks

fixed1 = { a -> 6.71, b -> -3.94, c -> -1.46, d -> -0.77, e -> -0.3,
h -> 5, k -> 1, m -> 0.05};
fixed2 = { a -> 0.42, b -> -0.25, c -> -0.09, d -> -0.048, e -> -0.02,
h -> 5, k -> 1, m -> 0.05};

F[Y_, Z_, x_] =
Y (a + b x + Y (c + (d + e Y) Y)) (h + k  Max[Z, 0])^(-1) /. fixed1;

G[Y_, Z_, x_] =
Y (a + b x + Y (c + (d + e Y) Y)) (h + k  Max[Z, 0])^(-1)} -
m Max[Z, 0] /. fixed2;

Needs["NDSolveFEM"]

eqn = {D[Y[x, t], x] == F[Y[x, t], Z[x, t], x],
D[Z[x, t], t] == G[Y[x, t], Z[x, t], x]}

sol = NDSolveValue[{eqn,
DirichletCondition[Y[x, t] == 1.0, x == 0 && 0 <= t <= 1],
DirichletCondition[Z[x, t] == 0.0, 0 <= x <= 1. && t == 0]}, {Y,
Z}, {x, 0, 1}, {t, 0, 1},
Method -> {"PDEDiscretization" -> "FiniteElement"}]


Here is the updated code with some modifications about the parameters and symbols . The rest of the code is the same you have just corrected except the use of ParametricNDSolve and Manipulate. My aim is to see the effect of free parameters on the solutions.

    fixed1 = {a -> 6.71, b -> -3.94, c -> -1.46, d -> -0.77, e -> -0.3,
h -> 5};
fixed2 = {a -> 0.42, b -> -0.25, c -> -0.09, d -> -0.048, e -> -0.02,
h -> 5};

params = {m, k}

appro = With[{s = 5. 10^1}, Tanh[s #]/2 + 1/2 &];

F[Y_, Z_, x_] :=
Y (a + b x + Y (c + (d + e Y) Y)) (h + k Z UnitStep[Z])^(-1) /.
fixed1 /. UnitStep -> appro;

G[Y_, Z_, x_] :=
Y (a + b x + Y (c + (d + e Y) Y)) (h + k Z UnitStep[Z])^(-1) -
m Z UnitStep[Z] /. fixed2 /. UnitStep -> appro;

Needs["NDSolveFEM"]
reg = Rectangle[{0, 0}, {1, 1}]; mesh =
ToElementMesh[reg, MaxCellMeasure -> .001]
eqn = {D[Y[x, t], x] == F[Y[x, t], Z[x, t], x],
D[Z[x, t], t] == G[Y[x, t], Z[x, t], x]};

sol = ParametricNDSolveValue[{eqn,
DirichletCondition[Y[x, t] == 1.0, x == 0 && 0 <= t <= 1],
DirichletCondition[Z[x, t] == 0.0, 0 <= x <= 1. && t == 0]}, {Y,
Z}, Element[{x, t}, mesh], params,
Method -> {"PDEDiscretization" -> "FiniteElement"}]

Manipulate[
Plot3D[sol[m, k][[1]][x, t], Element[{x, t}, reg], Mesh -> False,
MeshStyle -> White, AxesLabel -> {x, t, Y}], {m, 0, 1}, {k, 0, 1}]

Manipulate[
Plot3D[sol[m, k][[2]][x, t], Element[{x, t}, reg], Mesh -> False,
MeshStyle -> White, AxesLabel -> {x, t, Y}], {m, 0, 1}, {k, 0, 1}]


• The version of the code you've provided doesn't run on my machine, probably because of a stray extra ] in the second-to-last line. When I remove it, I get the error "There are more dependent variables, {Y[x,y],Z[x,y]}, than equations, so the system is underdetermined" instead. Jun 24, 2021 at 14:10
• Thanks. I removed the extra ] but I don't see why is undetermined: there are two equations and two unknowns. Assigning the two boundary conditions would suffice to solve the problem. I think that the error is related to the discretization method. I don't really know which one could work in this case. Jun 24, 2021 at 15:46
• As an aside, the curly braces in your definitions of F and G mean that F and G return lists with one element. I assume this was not your intention. Jun 24, 2021 at 15:53
• You are right but Mathematica has no problem with this. If I place directly the expressions of F and G into NDSolve, the output does not change and the error is the same, i. e. "the PDE coefficient ..." Jun 24, 2021 at 16:10

We can use UnitStep[Z] instead of Max[Z,0] and some approximation of UnitStep compatible with FEM

fixed1 = {a -> 6.71, b -> -3.94, c -> -1.46, d -> -0.77, e -> -0.3,
h -> 5, k -> 1, m -> 0.05};
fixed2 = {a -> 0.42, b -> -0.25, c -> -0.09, d -> -0.048, e -> -0.02,
h -> 5, k -> 1, m -> 0.05};
appro = With[{k = 5. 10^1}, Tanh[k #]/2 + 1/2 &];
F[Y_, Z_, x_] :=
Y (a + b x + Y (c + (d + e Y) Y)) (h + k Z UnitStep[Z])^(-1) /.
fixed1 /. UnitStep -> appro;

G[Y_, Z_, x_] :=
Y (a + b x + Y (c + (d + e Y) Y)) (h + k Z UnitStep[Z])^(-1) -
m Z UnitStep[Z] /. fixed2 /. UnitStep -> appro;

Needs["NDSolveFEM"]
reg = Rectangle[{0, 0}, {1, 1}]; mesh =
ToElementMesh[reg, MaxCellMeasure -> .001]
eqn = {D[Y[x, t], x] == F[Y[x, t], Z[x, t], x],
D[Z[x, t], t] == G[Y[x, t], Z[x, t], x]};

sol = NDSolveValue[{eqn,
DirichletCondition[Y[x, t] == 1.0, x == 0 && 0 <= t <= 1],
DirichletCondition[Z[x, t] == 0.0, 0 <= x <= 1. && t == 0]}, {Y,
Z}, Element[{x, t}, mesh],
Method -> {"PDEDiscretization" -> "FiniteElement"}]

{Plot3D[sol[x, t][[1]], Element[{x, t}, reg], ColorFunction -> Hue,
MeshStyle -> White, AxesLabel -> Automatic],
Plot3D[sol[x, t][[2]], Element[{x, t}, reg], ColorFunction -> Hue,
MeshStyle -> White, AxesLabel -> Automatic]}


• Thanks very much. I think this is a good step forward for my problem. At least there are no evaluation errors. Jun 24, 2021 at 17:18
• Hi Alex. I bother you for a slight modification of the code. Assume that I wish to leave one or more fixed parameters "params" free to be modified by means of the Manipulate command. Of course I need to use ParametricNDSolve[....,params,...]. and the insert the ParametricFunction into Manipulate. I noticed that the errorInitializePDECoefficients::femcnsd appears again. Is there any way to overcome also these issue? Thanks Jun 24, 2021 at 19:03
• @Fabio Could you provide code you have problem with? Jun 24, 2021 at 19:52
• In the invocation of NDSolve, you ask Mathematica to integrate over {y,0,1}. This should be {t,0,1}.
• Using the curly braces in your definitions for F and G means that F and G return a list with one element when evaluated. This means that instead of (for example) returning the numeric scalar -0.1 when called by NDSolve, they return the single-element list {-0.1}. To fix this, replace ^{-1} with ^(-1) in the definitions.
Fixing these two issues leads to a new and exciting error, which is related to the failure of $$F$$ and $$G$$ to have well-defined derivatives $$\partial F/\partial Z$$ and $$\partial G/\partial Z$$ when $$Z = 0$$. (Specifically, the error messages involve the derivatives of Max[Z,0] with respect to Z.) I am not sure how to address this effectively. However, if you replace Max[Z,0] with Z, the code returns a solution (while warning that the equations are "convection-dominated and may be unstable.")