New things and limitations in Version 12 numerical differential equation solver?

Question

This question is intended to be a place to summarize users' exemplary experience in solving differential equation with the NDSolve family in MMA’s latest version 12. Since there are so many differential equations and some are difficult, each time a major update may help resolve early bugs/issues and solve some new types. And certainly, many remain unable to be handled.

I personally have not yet tried much with the new solvers. Let me just list a few possible things of interest to get the ball rolling.

1. At least it is now equipped with a nonlinear FEM solver. But it certainly cannot solve all and contains bugs as well. Any info on that would be helpful to many of us.
2. Is it possible to control DifferenceOrder in NDSolve or NDEigenvalue, say, when solving an ODE?
3. What other new equations become solvable now?

---

Please feel free to edit this question.

Answer 1

I test nonlinear FEM using solutions obtained by other methods. I developed one of these methods for the problem of the nature convection, aerodynamics, and unsteady hydrodynamics, using linear FEM - see Solver for unsteady flow with the use of Mathematica FEM . Let me give an example. After the release of Mathematica version 12, I tested a non-linear FEM on the convection problem and compared it with my method. Surprisingly, the coincidence is one to one. The following code generates the same output as my code implementing the method of the false transient from the paper Vahl Davis, G.de (1983) : Natural convection of air in a square cavity : A bench mark numerical solution. Int. J. Numer. Methods Fluids 3, 249-264.

Code based on nonlinear FEM

Pr0 = .72; Ra = 10^5; R = Ra*Pr0; a = 0;
\[CapitalOmega] = ImplicitRegion[0 <= x <= 1 && 0 <= y <= 1, {x, y}];

{UX, VY, \[CapitalRho], TK, TX} = 
  NDSolveValue[{{Inactive[
          Div][({{-\[Mu], 0}, {0, -\[Mu]}}.Inactive[Grad][
            u[x, y], {x, y}]), {x, y}] + D[p[x, y], x] + 
        u[x, y]*D[u[x, y], x] + v[x, y]*D[u[x, y], y] - 
        R*T[x, y]*Sin[a], 
       Inactive[
          Div][({{-\[Mu], 0}, {0, -\[Mu]}}.Inactive[Grad][
            v[x, y], {x, y}]), {x, y}] + D[p[x, y], y] + 
        u[x, y]*D[v[x, y], x] + v[x, y]*D[v[x, y], y] - 
        R*T[x, y]*Cos[a], D[u[x, y], x] + D[v[x, y], y]} == {0, 0, 
       0} /. \[Mu] -> Pr0, {
     DirichletCondition[{p[x, y] == 0}, y == 1 && x == 1], 
     DirichletCondition[{u[x, y] == 0, v[x, y] == 0}, 
      y == 0 || y == 1]}, {u[x, y]*D[T[x, y], x] + 
       v[x, y]*D[T[x, y], y] - (D[T[x, y], x, x] + 
         D[T[x, y], y, y]) == NeumannValue[0, y == 0 || y == 1], 
     DirichletCondition[{u[x, y] == 0, v[x, y] == 0, T[x, y] == 1}, 
      x == 0], 
     DirichletCondition[{u[x, y] == 0, v[x, y] == 0, T[x, y] == 0}, 
      x == 1]}}, {u, v, p, T, 
\!\(\*SuperscriptBox[\(T\), 
TagBox[
RowBox[{"(", 
RowBox[{"1", ",", "0"}], ")"}],
Derivative],
MultilineFunction->None]\)}, {x, y} \[Element] \[CapitalOmega], 
   Method -> {"FiniteElement", 
     "InterpolationOrder" -> {u -> 2, v -> 2, p -> 1, T -> 2}, 
     "MeshOptions" -> {"MaxCellMeasure" -> 0.0001}}];

Compare with a bench mark numerical solution

MaximalBy[Table[{x, VY[x, .5]}, {x, 0, 1, .0001}], Last]

(*Out[5]= {{0.066, 68.7316}}*)

(*BenchmarkCase={x=0.066,vmax=68.59}*)

MaximalBy[Table[{y, UX[0.5, y]}, {y, 0, 1, .0001}], Last]

(*Out[8]= {{0.8543, 34.6607}}*)

(*BenchmarkCase={y=0.855,umax=34.73}*)
 MaximalBy[Table[{y, -TX[0, y]}, {y, 0, 1, 0.0001}], Last]

(*Out[12]= {{0.0747, 7.73417}}*)

(*BenchmarkCase={y=0.081,Numax=7.717}*)

 NuAv = 
 NIntegrate[UX[x, y]*TK[x, y] - TX[x, y], {x, 0, 1}, {y, 0, 1}, 
  WorkingPrecision -> 5]

(*Out[14]= 4.4880*)

(*BenchmarkCase=4.519*)


Nu0 = NIntegrate[-TX[0, y], {y, 0, 1}, WorkingPrecision -> 5]

(*Out[16]= 4.5274*)

(*BenchmarkCase=4.509*)

Nu05 = 
 NIntegrate[UX[.5, y]*TK[.5, y] - TX[0.5, y], {y, 0, 1}, 
  WorkingPrecision -> 5]

(*Out[18]= 4.5252*)

(*BenchmarkCase=4.519*)

MinimalBy[Table[{y, -TX[0, y]}, {y, 0, 1, 0.0001}], Last]

(*Out[19]= {{1., 0.726786}}*)

(*BenchmarkCase={y=1,Numin=0.729}*)

Code based on linear FEM

k = 100; Pr0 = .72; Ra = 10^5; R = Ra*Pr0; t0 = 1/100; a = 0;
\[CapitalOmega] = ImplicitRegion[0 <= x <= 1 && 0 <= y <= 1, {x, y}];
UX[0][x_, y_] := 0;
VY[0][x_, y_] := 0;
\[CapitalRho][0][x_, y_] := 0;
TK[0][x_, y_] := 0; TX[0][x_, y_] := 0;
Do[
  {UX[i], VY[i], \[CapitalRho][i], TK[i], TX[i]} = 
   NDSolveValue[{{Inactive[
           Div][({{-\[Mu], 0}, {0, -\[Mu]}}.Inactive[Grad][
             u[x, y], {x, y}]), {x, y}] + D[p[x, y], x] + 
         UX[i - 1][x, y]*D[u[x, y], x] + 
         VY[i - 1][x, y]*D[u[x, y], y] + (u[x, y] - UX[i - 1][x, y])/
          t0 - R*T[x, y]*Sin[a], 
        Inactive[
           Div][({{-\[Mu], 0}, {0, -\[Mu]}}.Inactive[Grad][
             v[x, y], {x, y}]), {x, y}] + D[p[x, y], y] + 
         UX[i - 1][x, y]*D[v[x, y], x] + 
         VY[i - 1][x, y]*D[v[x, y], y] - 
         R*T[x, y]*Cos[a] + (v[x, y] - VY[i - 1][x, y])/t0, 
        D[u[x, y], x] + D[v[x, y], y]} == {0, 0, 0} /. \[Mu] -> Pr0, {
      DirichletCondition[{p[x, y] == 0}, y == 1 && x == 1], 
      DirichletCondition[{u[x, y] == 0, v[x, y] == 0}, 
       y == 0 || y == 1]}, {UX[i - 1][x, y]*D[T[x, y], x] + 
        VY[i - 1][x, y]*D[T[x, y], y] + (T[x, y] - TK[i - 1][x, y])/
         t0 - (D[T[x, y], x, x] + D[T[x, y], y, y]) == 
       NeumannValue[0, y == 0 || y == 1], 
      DirichletCondition[{u[x, y] == 0, v[x, y] == 0, T[x, y] == 1}, 
       x == 0], 
      DirichletCondition[{u[x, y] == 0, v[x, y] == 0, T[x, y] == 0}, 
       x == 1]}}, {u, v, p, T, 
\!\(\*SuperscriptBox[\(T\), 
TagBox[
RowBox[{"(", 
RowBox[{"1", ",", "0"}], ")"}],
Derivative],
MultilineFunction->None]\)}, {x, y} \[Element] \[CapitalOmega], 
    Method -> {"FiniteElement", 
      "InterpolationOrder" -> {u -> 2, v -> 2, p -> 1, T -> 2}, 
      "MeshOptions" -> {"MaxCellMeasure" -> 0.0001}}], {i, 1, k}];

Compare with a bench mark numerical solution

MaximalBy[Table[{x, VY[k][x, .5]}, {x, 0, 1, .0001}], Last]

(*Out[9]= {{0.066, 68.7316}}*)

(*BenchmarkCase={x=0.066,vmax=68.59}*)

MaximalBy[Table[{y, UX[k][0.5, y]}, {y, 0, 1, .0001}], Last]

(*Out[12]= {{0.8543, 34.6607}}*)
(*BenchmarkCase={y=0.855,umax=34.73}*)

MaximalBy[Table[{y, -TX[k][0, y]}, {y, 0, 1, 0.0001}], Last]

(*Out[16]= {{0.0747, 7.73417}}*)

(*BenchmarkCase={y=0.081,Numax=7.717}*)

 NuAv = 
 NIntegrate[
  UX[k][x, y]*TK[k][x, y] - TX[k][x, y], {x, 0, 1}, {y, 0, 1}, 
  WorkingPrecision -> 5]

(*Out[18]= 4.4880*)
In[19]:= (*BenchmarkCase=4.519*)


 Nu0 = 
 NIntegrate[-TX[k][0, y], {y, 0, 1}, WorkingPrecision -> 5]

(*Out[20]= 4.5274*)

(*BenchmarkCase=4.509*)

 Nu05 = 
 NIntegrate[UX[k][.5, y]*TK[k][.5, y] - TX[k][0.5, y], {y, 0, 1}, 
  WorkingPrecision -> 5]

(*Out[22]= 4.5252*)

(*BenchmarkCase=4.519*)

MinimalBy[Table[{y, -TX[k][0, y]}, {y, 0, 1, 0.0001}], Last]

(*Out[23]= {{1., 0.726786}}*)

(*BenchmarkCase={y=1,Numin=0.729}*)

Note that both methods nonlinear and linear FEM are consistent with each other and with the third method, published by Davis, G.de (1983) : Natural convection of air in a square cavity : A bench mark numerical solution. Int. J. Numer. Methods Fluids 3, 249-264.

Answer 2

Re question 2 and NDSolve:

How to control "DifferenceOrder" at each step has been shown in this tutorial, which has been available for several versions. Dynamic control of order is usually not simple, since you have to evaluate whether to change the order every step (or every $k^{\rm th}$ step).

One can control the order for certain method families, such as the Runge-Kutta ones (explicit, implicit, and symplectic). The connection between difference order and interpolation order is explored in this answer.

Answer 3

Re question 3: There are NBodySimulation and NBodySimulationData.

Related Q&A: (25562), (31094), (38478), (63571), (81340), (81340), (105342), (121035), (133349), (135857), (164196), (181275)