25
$\begingroup$

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

solution for p=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

error

$\endgroup$
  • $\begingroup$ fyi, You can same output with much simpler input p = 2; Omega = Line[{{0}, {1}}]; eq = D[u[x], x]^p == 1; bc = DirichletCondition[u[x] == 0, True]; usol = NDSolveValue[{eq, bc}, u, Element[{x}, Omega]] It seems it does not like the derivative being raised to higher powers than 1. Nonlinear coefficients are not supported in this version of NDSolve $\endgroup$ – Nasser Apr 18 at 7:49
25
$\begingroup$

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

enter image description here

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

enter image description here

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

enter image description here

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]
$\endgroup$
  • $\begingroup$ Cool! thank you very much, I did not see the dependency of the coefficients on the gradient of the solution. Thank you! $\endgroup$ – Mauricio Fernández Apr 18 at 8:19
  • 2
    $\begingroup$ @MauricioFernández, hi, would you mind if I use this example in the ref/message/NDSolveValue/femnlmdor help page? $\endgroup$ – user21 Apr 26 at 6:23
  • $\begingroup$ I dont mind at all, please use it! Whatever helps to make something more comprehensible is a good thing :D $\endgroup$ – Mauricio Fernández Apr 27 at 7:28
  • $\begingroup$ @MauricioFernández,, done. Thanks. $\endgroup$ – user21 Apr 29 at 7:11
  • $\begingroup$ @xzczd, thanks for the code edit. $\endgroup$ – user21 May 14 at 6:12

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service, privacy policy and cookie policy

Not the answer you're looking for? Browse other questions tagged or ask your own question.