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?

Question

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

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

   

Answer 1

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

Answer 2

You have a couple of syntax errors in your code that are giving you problems, and a deeper mathematical issue.

• 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.")