# Non-linear PDE coefficients problem

I'm trying to solve the following PDE

NDSolveValue[{Derivative[1, 0][\[Gamma]][u, v]*(6*Cot[\[Gamma][u, v]]*Derivative[0, 1][\[Gamma]][u, v]^2 - 2*Derivative[0, 2][\[Gamma]][u, v] + 3*Cos[\[Gamma][u, v]]*Derivative[1, 1][\[Gamma]][u, v]) ==
Derivative[0, 1][\[Gamma]][u, v]*(2*(2 + Cos[2*\[Gamma][u, v]])*Csc[\[Gamma][u, v]]*Derivative[1, 0][\[Gamma]][u, v]^2 + 3*Derivative[1, 1][\[Gamma]][u, v] - 2*Cos[\[Gamma][u, v]]*Derivative[2, 0][\[Gamma]][u, v]),\[Gamma][0, v] == 1 + Sqrt[v], \[Gamma][u, 0] == 1 + u*Sin[u]}, \[Gamma], {u, 0, 1/2}, {v, 0, 1/2}]


however, I get the warning

  The maximum derivative order of the nonlinear PDE coefficients for the Finite Element Method is larger than 1. It may help to rewrite  the PDE in inactive form


which does not make sense, since the coefficients are all of derivative order 1, not larger than it. Is Mathematica not able to get the solutions of these types of quasilinear PDEs ?

When written in inactive form as

       NDSolve[{Inactive[Plus][Inactive[Times][Inactive[Derivative][1, 0][\[Gamma]][u, v], Inactive[Plus][Inactive[Times][6, Inactive[Cot][Inactive[\[Gamma]][u, v]],
Inactive[Power][Inactive[Derivative][0, 1][\[Gamma]][u, v], 2]], Inactive[Times][-1, Inactive[Times][2, Inactive[Derivative][0, 2][\[Gamma]][u, v]]],
Inactive[Times][3, Inactive[Cos][Inactive[\[Gamma]][u, v]], Inactive[Derivative][1, 1][\[Gamma]][u, v]]]], Inactive[Times][-1, Inactive[Times][Inactive[Derivative][0, 1][\[Gamma]][u, v],
Inactive[Plus][Inactive[Times][2, Inactive[Plus][2, Inactive[Cos][Inactive[Times][2, Inactive[\[Gamma]][u, v]]]], Inactive[Csc][Inactive[\[Gamma]][u, v]],
Inactive[Power][Inactive[Derivative][1, 0][\[Gamma]][u, v], 2]], Inactive[Times][3, Inactive[Derivative][1, 1][\[Gamma]][u, v]],
Inactive[Times][-1, Inactive[Times][2, Inactive[Cos][Inactive[\[Gamma]][u, v]], Inactive[Derivative][2, 0][\[Gamma]][u, v]]]]]]] == 0, DirichletCondition[{\[Gamma][0, v] == 1 + Sqrt[v]}, u == 0], DirichletCondition[{\[Gamma][u, 0] == 1 + u*Sin[u]}, v == 0]}, \[Gamma], {u, 0, 1/2}, {v, 0, 1/2}]


 There are more dependent variables than equations, so the system is underdetermined

• Derivative[0, 2] and Derivative[2,0 ] are higher-than-first-order derivatives. Commented Apr 1 at 20:50
• Do you want to solve this problem or just speculate about the warnings? Commented Apr 1 at 21:26
• @AlexTrounev I would solve it if the warnings were clearly related to the problem and gave clear indication of how to fix it. Commented Apr 2 at 14:04
• @azerbajdzan Yes, but it's a bout the coefficients, not the highest derivatives themselves. The equation is semilinear. Commented Apr 2 at 14:04
• @DanielCastro Upfront it is not clear what type of equation we have here. If it is wave equation then we need two boundary conditions at u=0 or v=0. With this updates we can make several steps using method of lines. Commented Apr 2 at 17:45

First, let consider a minimal working example of numerical solution at some boundary conditions to use NDSolve as a solver

g = NDSolveValue[{Derivative[1, 0][\[Gamma]][u,
v]*(6*Cot[\[Gamma][u, v]]*Derivative[0, 1][\[Gamma]][u, v]^2 -
2*Derivative[0, 2][\[Gamma]][u, v] +
3*Cos[\[Gamma][u, v]]*Derivative[1, 1][\[Gamma]][u, v]) ==
Derivative[0, 1][\[Gamma]][u,
v]*(2*(2 + Cos[2*\[Gamma][u, v]])*Csc[\[Gamma][u, v]]*
Derivative[1, 0][\[Gamma]][u, v]^2 +
3*Derivative[1, 1][\[Gamma]][u, v] -
2*Cos[\[Gamma][u, v]]*
Derivative[2, 0][\[Gamma]][u, v]), \[Gamma][0, v] ==
1 + v, \[Gamma][u, 0] == 1 + u,
Derivative[0, 1][\[Gamma]][u, 0] == 1,
Derivative[1, 0][\[Gamma]][1/2, v] == 0}, \[Gamma], {u, 0,
1/2}, {v, 0, 1/2},
Method -> {"MethodOfLines",
"SpatialDiscretization" -> {"TensorProductGrid",
"MinPoints" -> 33, "MaxPoints" -> 41, "DifferenceOrder" -> 4}}]


There are 2 messages but final solution looks good

Plot3D[g[u, v], {u, 0, 1/2}, {v, 0, 1/2}, Mesh -> None,
ColorFunction -> "Rainbow"]


Update 1. This is FDM solution, not stably but we can see how it should be if we force Mathematica FEM to solve this problem. Note, we solve Goursat problem

XYgrid[dom_List, pts_List] :=
N@Range[Sequence @@ #1, Abs[Subtract @@ #1]/#2] &, {dom,
pts - 1}];
BoundaryIndex[xgridlen_, ygridlen_] :=
Module[{tmp, left, right, bot, top},
tmp = Table[(n - 1) ygridlen + Range[1, ygridlen], {n, 1,
xgridlen}]; {left, right} = tmp[[{1, -1}]]; {bot, top} =
Transpose[{First[#], Last[#]} & /@ tmp]; {top, right[[2 ;; -2]],
bot, left[[2 ;; -2]]}];
FDMat[deriv_, xygrid_, difforder_] :=
Map[NDSolveFiniteDifferenceDerivative[#, xygrid,
"DifferenceOrder" -> difforder]["DifferentiationMatrix"] &, deriv]
{tmax, domain, pts, difforder} = {10, {{0, 1/2}, {0, 1/2}}, {12, 12},
4};
xygrid = XYgrid[domain, pts]; {nx, ny} =
Map[Length, xygrid]; {top, right, bot, left} =
BoundaryIndex[nx, ny]; {dx, dy, dx2, dy2, dxy} =
FDMat[{{1, 0}, {0, 1}, {2, 0}, {0, 2}, {1, 1}}, xygrid,
difforder]; boundaries = Join[top, right, bot, left]; sgrid =
Flatten[Outer[List, Sequence @@ xygrid], 1]; g =
Table[uu[i], {i, nx ny}];

eqs = 6 Cot[g] (dy . g)^2 (dx . g) - 2 (dy2 . g) (dx . g) -
4 Csc[g] (dy . g) (dx . g)^2 -
2 Cos[2 g] Csc[g] (dy . g) (dx . g)^2 - 3 (dy . g) (dxy . g) +
3 Cos[g] (dx . g) (dxy . g) + 2 Cos[g] (dy . g) (dx2 . g);

eqs[[left]] = g[[left]] - (1 + Sqrt[sgrid[[left]][[All, 2]]]);

eqs[[bot]] =
g[[bot]] - (1 + sgrid[[bot]][[All, 1]] Sin[sgrid[[bot]][[All, 1]]]);
gini = (1 + Sqrt[#[[2]]] + #[[1]] Sin[#[[1]]]) & /@
sgrid;
sol1 = FindRoot[Table[eqs[[i]] == 0, {i, Length[eqs]}],
Table[{g[[i]], gini[[i]]}, {i, Length[g]}], MaxIterations -> 10000, Method -> {"Newton", "StepControl" -> "TrustRegion"}];


Visualization

{U1} = Map[Interpolation@Join[sgrid, Transpose@List@#, 2] &,
Partition[g /. sol1, Length[sgrid]]];
Plot3D[U1[u, v], {u, 0, 1/2}, {v, 0, 1/2}, Mesh -> None,
ColorFunction -> "Rainbow", PlotRange -> All,
AxesLabel -> {"u", "v", "\[Gamma]"}, Boxed -> False]


Update 2. This is hand made FEM solution using FEM differentiation matrices proposed by user21 (thanks to him) here

Needs["NDSolveFEM"]
FiniteElementDerivative[order : {__Integer}, mesh_ElementMesh] /;
1 <= Length[order] <= 3 :=
Block[{dim, nr, vd, sd, mdata, ccoef, pos, dcoef, cdata},
dim = Length[order];
nr = ToNumericalRegion[mesh];
vd = NDSolveVariableData[{"DependentVariables", "Space"} -> {{u},
Table[Unique[X], {dim}]}];
sd = NDSolveSolutionData[{"Space"} -> {nr}];
mdata = InitializePDEMethodData[vd, sd];
ccoef = ConstantArray[0, dim];
pos = Flatten[Position[order, 1]];
ccoef[[pos]] = 1;
dcoef = ConstantArray[0, dim];
pos = Flatten[Position[order, 2]];
dcoef[[pos]] = 1;
dcoef = DiagonalMatrix[dcoef];
(*"Pure ConvectionCoefficients" will trigger a warning*)
Quiet[cdata =
InitializePDECoefficients[vd, sd,
"DiffusionCoefficients" -> {{dcoef}},
"ConvectionCoefficients" -> {{ccoef}}], \
{InitializePDECoefficients::femcscd}];
DiscretizePDE[cdata, mdata, sd]];


First we define mesh, matrices and boundary as follows

mesh = ToElementMesh[Rectangle[{0, 0}, {1/2, 1/2}],
MaxCellMeasure -> 1/24^2]

{dx, dy, dx2, dy2, dxy} =
FiniteElementDerivative[#, mesh]["StiffnessMatrix"] & /@ {{1,
0}, {0, 1}, {2, 0}, {0, 2}, {1, 1}};
boundary = mesh["BoundaryElements"][[1, 1]] // Flatten;
coord = mesh["Coordinates"];
bc = coord[[boundary]];
left = Table[
If[coord[[boundary]][[i]][[1]] == 0, boundary[[i]], Nothing], {i,
Length[boundary]}] // DeleteDuplicates;
right = Table[
If[coord[[boundary]][[i]][[1]] == 1/2, boundary[[i]],
Nothing], {i, Length[boundary]}] // DeleteDuplicates;
bot = Table[
If[coord[[boundary]][[i]][[2]] == 0, boundary[[i]], Nothing], {i,
Length[boundary]}] // DeleteDuplicates;
top = Table[
If[coord[[boundary]][[i]][[2]] == 1/2, boundary[[i]],
Nothing], {i, Length[boundary]}] // DeleteDuplicates;


Then we use equations and boundary conditions from FDM shown above

g = Table[uu[i], {i, Length[coord]}];

eqs = 6  Cot[g] (dy . g)^2 (dx . g) - 2  (dy2 . g) (dx . g) -
4  Csc[g] (dy . g) (dx . g)^2 -
2  Cos[2  g] Csc[g] (dy . g) (dx . g)^2 - 3 (dy . g) (dxy . g) +
3  Cos[g] (dx . g) (dxy . g) + 2  Cos[g] (dy . g) (dx2 . g);
eqs[[left]] = g[[left]] - (1 + Sqrt[coord[[left]][[All, 2]]]);

eqs[[bot]] =
g[[bot]] - (1 + coord[[bot]][[All, 1]] Sin[coord[[bot]][[All, 1]]]);
gini = (1 + Sqrt[#[[2]]] + #[[1]] Sin[#[[1]]]) & /@
coord;


Finally we solve equations and plot numerical solution

sol2 = FindRoot[Table[eqs[[i]] == 0, {i, Length[eqs]}],
Table[{g[[i]], gini[[i]]}, {i, Length[g]}],
Method -> {"Newton", "StepControl" -> "TrustRegion"},
MaxIterations -> 10000];

U2 = ElementMeshInterpolation[{mesh}, g /. sol2];
Plot3D[U2[u, v], {u, 0, 1/2}, {v, 0, 1/2}, Mesh -> None,
ColorFunction -> "Rainbow", PlotRange -> All,
AxesLabel -> {"u", "v", "\[Gamma]"}, PlotTheme -> "Scientific"]


It is interesting that FDM solution looks much smoother than FEM solution. This can be explained by the fact that for FDM we use a scheme of 4 orders of accuracy, and for FEM only 2. If we use a mesh with double the number of elements, the solution becomes much smoother.

mesh = ToElementMesh[Rectangle[{0, 0}, {1/2, 1/2}],
MaxCellMeasure -> 1/34^2]


• I see, thank you. So it may be an issue of the boundary conditions ? As per the active-inactive way of the syntax I find it unreasonable that it may work well for some conditions and not for others. Commented Apr 2 at 19:58
• @DanielCastro Mathematica FEM is not universal solver. We can try some hand made code like FDM to see what is the main issue with your boundary conditions. Commented Apr 2 at 20:56
• @DanielCastro Please see Update 1 to my answer with FDM solution. The good news is that boundary conditions you try to use with FEM are quite suitable for solving the problem. Commented Apr 3 at 1:35
• @DanielCastro Please, see Update 2 with FEM solution. Commented Apr 4 at 1:49
• @AlexTrounev Very interesting approch (as is so often the case). In the FEM approach update 2 I'm missing NeumannValue conditions which come from integrating the second order derivatives of the unknown function. It's not clear to me either why it's allowed to factor the FEM-matrices to get the final ode-system. Thanks! Commented Apr 9 at 10:38

When you look at the message:

You will see a blue circle-i at the end. Click on that, it will take you to the message ref page for this issue, and that page explains the cause of this and how to avoid it.

For reference: see here and here.

What you need to do: Write down your equation, then write down the coefficient form of the FEM equation referred to in the message ref page and the finite element best practice tutorial. You need to match your equation to the FEM coefficient form. If you have problems with that you could add the information to your question.

Here is how that procedure would work in principal. We start from an Expand-ed form of your equation:

yourOp =
Derivative[1, 0][\[Gamma]][u,
v]*(6*Cot[\[Gamma][u, v]]*Derivative[0, 1][\[Gamma]][u, v]^2 -
2*Derivative[0, 2][\[Gamma]][u, v] +
3*Cos[\[Gamma][u, v]]*
Derivative[1, 1][\[Gamma]][u, v]) - (Derivative[0,
1][\[Gamma]][u,
v]*(2*(2 + Cos[2*\[Gamma][u, v]])*Csc[\[Gamma][u, v]]*
Derivative[1, 0][\[Gamma]][u, v]^2 +
3*Derivative[1, 1][\[Gamma]][u, v] -
2*Cos[\[Gamma][u, v]]*Derivative[2, 0][\[Gamma]][u, v]));

Expand[yourOp] // InputForm

yourOpExpanded =
6*Cot[\[Gamma][u, v]]*Derivative[0, 1][\[Gamma]][u, v]^2*
Derivative[1, 0][\[Gamma]][u, v] -
2*Derivative[0, 2][\[Gamma]][u, v]*
Derivative[1, 0][\[Gamma]][u, v] -
4*Csc[\[Gamma][u, v]]*Derivative[0, 1][\[Gamma]][u, v]*
Derivative[1, 0][\[Gamma]][u, v]^2 -
2*Cos[2*\[Gamma][u, v]]*Csc[\[Gamma][u, v]]*
Derivative[0, 1][\[Gamma]][u, v]*
Derivative[1, 0][\[Gamma]][u, v]^2 -
3*Derivative[0, 1][\[Gamma]][u, v]*
Derivative[1, 1][\[Gamma]][u, v] +
3*Cos[\[Gamma][u, v]]*Derivative[1, 0][\[Gamma]][u, v]*
Derivative[1, 1][\[Gamma]][u, v] +
2*Cos[\[Gamma][u, v]]*Derivative[0, 1][\[Gamma]][u, v]*
Derivative[2, 0][\[Gamma]][u, v];


Now, we start by looking at one term:

yourOpExpanded =(*6*Cot[\[Gamma][u,v]]*Derivative[0,1][\[Gamma]][u,
v]^2*Derivative[1,0][\[Gamma]][u,
v]*)-2*Derivative[0, 2][\[Gamma]][u, v]*
Derivative[1, 0][\[Gamma]][u,
v](*-4*Csc[\[Gamma][u,v]]*Derivative[0,1][\[Gamma]][u,v]*\
Derivative[1,0][\[Gamma]][u,v]^2-2*Cos[2*\[Gamma][u,v]]*Csc[\[Gamma][\
u,v]]*Derivative[0,1][\[Gamma]][u,v]*Derivative[1,0][\[Gamma]][u,v]^2-\
3*Derivative[0,1][\[Gamma]][u,v]*Derivative[1,1][\[Gamma]][u,v]+3*Cos[\
\[Gamma][u,v]]*Derivative[1,0][\[Gamma]][u,v]*Derivative[1,1][\[Gamma]\
][u,v]+2*Cos[\[Gamma][u,v]]*Derivative[0,1][\[Gamma]][u,v]*Derivative[\
2,0][\[Gamma]][u,v]*);

yourOpExpanded // InputForm
-2*Derivative[0, 2][\[Gamma]][u, v]*
Derivative[1, 0][\[Gamma]][u, v]


Next, we set of an Inactive Div/Grad term:

inactiveOp =
Inactive[
Div][-{{0 , 0}, {0, 2*Derivative[1, 0][\[Gamma]][u, v]}} .
Inactive[Grad][\[Gamma][u, v], {u, v}], {u, v}]


And we compare to your form:

(Activate[inactiveOp] - yourOpExpanded) // InputForm
-2*Derivative[0, 1][\[Gamma]][u, v]*
Derivative[1, 1][\[Gamma]][u, v]


Because we pulled the term -2*Derivative[1, 0][\[Gamma]][u, v] into the inactive Div/Grad term we need to compensate that:

inactiveOp =
Inactive[
Div][-{{0 , 0}, {0, 2*Derivative[1, 0][\[Gamma]][u, v]}} .
Inactive[Grad][\[Gamma][u, v], {u, v}], {u, v}] + {0,
2*Derivative[1, 1][\[Gamma]][u, v]} .

Activate[inactiveOp] - yourOpExpanded
0


You need to do this for all terms. Also, note that when adding a new term sometimes it may be better to re-write a previous addition. Because we do not know where your equation comes from it's hard to say how to approach this. But also note that the Inactive form of the PDE is, very, very general. Good luck!

• Please see edit about the warning when written in inactive form. Commented Apr 1 at 9:44
• @user21 Could you explain how we can force NDSolve to solve Goursat problem for hyperbolic equation discussed here? Commented Apr 6 at 10:53
• @AlexTrounev, I have tried to outline how I this could be tried. Is something not clear? Commented Apr 9 at 7:34
• @user21 Is the Guorsat problem solution algorithm implemented in NDSolve ? Could you explain this algorithm in a case of NDSolve[{Derivative[2, 0][\[Gamma]][u, v] - 2*Derivative[1, 1][\[Gamma]][u, v] - Derivative[0, 2][\[Gamma]][u, v] == 0, [Gamma][0, v] == 1 + Sqrt[v], \[Gamma][u, 0] == 1 + u*Sin[u]}, \[Gamma], {u, 0, 1/2}, {v, 0, 1/2}]? Commented Apr 9 at 9:37
• @AlexTrounev, sorry I do not understand the question. Commented Apr 9 at 14:57

Your NDSolveValue call seems to contain a system of partial differential equations (PDEs) to be solved for the function [Gamma][u, v]. However, there is a minor error in the use of the NDSolveValue function, as well as some syntax issues in the definition of the equations.

Here's the corrected version of your code:

sol = NDSolveValue[{
Derivative[1, 0][\[Gamma]][u, v]*(6*Cot[\[Gamma][u, v]]*Derivative[0, 1][\[Gamma]][u, v]^2 -
2*Derivative[0, 2][\[Gamma]][u, v] + 3*Cos[\[Gamma][u, v]]*Derivative[1, 1][\[Gamma]][u, v]) ==
Derivative[0, 1][\[Gamma]][u, v]*(2*(2 + Cos[2*\[Gamma][u, v]])*Csc[\[Gamma][u, v]]*
Derivative[1, 0][\[Gamma]][u, v]^2 + 3*Derivative[1, 1][\[Gamma]][u, v] -
2*Cos[\[Gamma][u, v]]*Derivative[2, 0][\[Gamma]][u, v]),
\[Gamma][0, v] == 1 + Sqrt[v],
\[Gamma][u, 0] == 1 + u*Sin[u]
},
\[Gamma], {u, 0, 1/2}, {v, 0, 1/2}
];


Here's the explanation:

• NDSolveValue is the function used to solve differential equations numerically.
• The equations are provided as a list inside {}.
• The boundary conditions are specified for [Gamma] at the boundaries u = 0 and v = 0.
• The range of u and v is specified as {u, 0, 1/2} and {v, 0, 1/2}, respectively.
• The solution is stored in the variable sol.

After running this code, sol` will contain an InterpolatingFunction representing the solution for [Gamma][u, v]. You can then use this function for further analysis or plotting.

• Thank you, but it's not helping in anything. Commented Apr 2 at 19:53