Non-linear PDE coefficients problem

Question

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}]

I get instead the warning

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

Answer 1

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] := 
  MapThread[
   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[NDSolve`FiniteDifferenceDerivative[#, 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["NDSolve`FEM`"]
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 = NDSolve`VariableData[{"DependentVariables", "Space"} -> {{u}, 
       Table[Unique[X], {dim}]}];
   sd = NDSolve`SolutionData[{"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]

Answer 2

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]} . 
   Inactive[Grad][\[Gamma][u, v], {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!

Answer 3

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.