How do I use the new nonlinear finite element in Mathematica 12 for this equation?

Question

With Mathematica 12 we get new technology for nonlinear finite elements. Out of curiosity, I just wanted to solve the following equation

$$ \frac{d}{dx} \left( c(x) \left[\frac{d}{dx} u(x)\right]^p \right) = r(x) \ , \quad u(0) = u(1) = 0 $$

for some given positive coefficient $c(x) > 0$ and right hand side $r(x)$ in the domain $\Omega = [0,1]$ for some odd integer power $p=1,3,5,\dots$ (I am coming from mechanics and the odd power ensures that the energy is convex such that a minimum exists, ensuring the existence of the solution $u(x)$).

For $p=1$, you have the standard linear DE and you can solve it with finite elements since Mathematica 10.

p = 1;
Omega = Line[{{0}, {1}}];
c[x_] := x^2 + 3;
r[x_] := Sin@x;
eq = D[c[x]*D[u[x], x]^p, x] == r[x];
bc = DirichletCondition[u[x] == 0, True];
usol = NDSolveValue[{eq, bc}, u, Element[{x}, Omega]];
Plot[usol[x], {x, 0, 1}]

But for $p \geq 3$ you get a nonlinear DE. How do you use the new nonlinear finite element method to solve the equation for $p \geq 3$? I tried to use Inactive on D, but I could not figure out where to use it properly. Can you help me out? When you use $p=3$ in the code above, you get the following error

Answer 1

OK, there are a few things going on here. Let me explain them in turn. First, as the message suggests, this should be written in Inactive form (we'll get to the why later). If you click on the three dots in front of the error message

and follow the link to the reference page you will find some information on this error message.

To write the equation in inactive form it's helpful to remember what that form looks like. In the nonlinear case:

$$ \nabla \cdot (-c(t,X,u,\nabla _Xu) \nabla u-\alpha (t,X,u,\nabla _Xu) u $$ $$ + \gamma (t,X,u,\nabla _Xu)) + \beta (t,X,u,\nabla _Xu)\cdot \nabla u+a(t,X,u,\nabla _Xu) u$$ $$ - f(t,X,u,\nabla _Xu)=0.$$

Details are here. And you have to get your equation into that form otherwise you are out of FEM luck. So I rewrote this as:

$$ \frac{d}{dx} \left( \left[ c(x) \left(\frac{d}{dx} u(x)\right)^\left(p-1\right) \right] \frac{d}{dx} u(x) \right) = r(x) \ , \quad u(0) = u(1) = 0 $$

Omega = Line[{{0}, {1}}];
c[x_] := x^2 + 3;
r[x_] := Sin@x;
eq[p_] := 
 Inactive[Div][(c[x]*D[u[x], x]^(p - 1)) Inactive[Grad][
     u[x], {x}], {x}] == r[x]
bc = DirichletCondition[u[x] == 0, True];

Note how I used D[u[x], x]^(p - 1) and Inactive[Grad][u[x],{x}] to separate the original equation.

Now, when you use that:

NDSolveValue[{eq[3], bc}, u, Element[{x}, Omega]];

FindRoot::nosol: Linear equation encountered that has no solution.

FindRoot::sszero: The step size in the search has become less than the tolerance prescribed by the PrecisionGoal option, but the function value is still greater than the tolerance prescribed by the AccuracyGoal option.

So NDSolve (or better FindRoot) can not find a solution. The first thing to try if this happens is to use a less non-linear problem (in this case the linear PDE) as a starting value like so:

usolP1 = NDSolveValue[{eq[1], bc}, u, Element[{x}, Omega]];

usolP3 = NDSolveValue[{eq[3], bc}, u, Element[{x}, Omega], 
   InitialSeeding -> {u[x] == usolP1[x]}];
Plot[usolP3[x], {x, 0, 1}]

usolP5 = NDSolveValue[{eq[5], bc}, u, Element[{x}, Omega], 
   InitialSeeding -> {u[x] == usolP3[x]}];
Plot[usolP5[x], {x, 0, 1}]

The scope of the nonlinear FEM solver is given here:

• The coefficients can be functions of space, time, parameters, dependent variables and first order derivatives of dependent variables.

Concerning the use of Inactive. This comes up when you have derivatives of the dependent variable as a nonlinear diffusion coefficient. This is explained in detail in the last part of the section Formal Partial Differential Equations.

You can find more examples of nonlinear PDEs in the usual places (FEMDocumentation) and in the nonlinear FEM verification tests that are now part of the FEMDocumentation. Direct your help system to FEMDocumentation/tutorial/NonlinearFiniteElementVerificationTests (The web version looks bad and needs to be fixed but is here)

All details of the implementation including code for the nonlinear FEM solver is documented here.

Update:

You could also use ParamatricNDSolveValue to automate the task somewhat:

Set up the equation, bcs and region:

eqn[p_] := 
 Inactive[Div][((x^2 + 3)*D[u[x], x]^(p - 1)) Inactive[Grad][
     u[x], {x}], {x}] - Sin[x]
bc = DirichletCondition[u[x] == 0, True];
\[CapitalOmega] = Line[{{0}, {1}}];

Create an initial seeding:

Clear[seeding]
seeding[x_?NumericQ] := 0

Set up a ParamatricNDSolveValue function with parameter p.

psol = ParametricNDSolveValue[{eqn[p] == 0, bc}, u, 
  Element[{x}, \[CapitalOmega]], {p}, 
  InitialSeeding -> {u[x] == seeding[x]}]

Solve for p=1:

s1 = psol[1]

Clear the seeding and reset it to point to the previous solution s1

Clear[seeding]
seeding[x_?NumericQ] := s1[x]

Solve for p=3:

s2 = psol[3]