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

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.

• "nonlinear PDE solver" -> "nonlinear FEM solver". Apr 20, 2019 at 5:20
• Are you asking about DSolve or NDSolve ? If both, I think be better to ask separate question on each since these are really different functions and have different purposes. If not, then it is not clear which one are you asking about. Apr 20, 2019 at 6:11
• This question is far too vague as asked. Apr 20, 2019 at 11:32
• If subject header is changed to "...numeric differential equation solver" then it would seem to be a perfectly reasonable question. Apr 20, 2019 at 14:26
• I still think the question too broad. It could perhaps be left open as a CW question, which has been done before with posts that have a multiplicity of questions to be answered. The opening line already shows it is not a Q&A post, and the breadth of the example questions further confirm this. Apr 21, 2019 at 13:04

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= {{0.066, 68.7316}}*)

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

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

(*Out= {{0.8543, 34.6607}}*)

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

(*Out= {{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= 4.4880*)

(*BenchmarkCase=4.519*)

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

(*Out= 4.5274*)

(*BenchmarkCase=4.509*)

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

(*Out= 4.5252*)

(*BenchmarkCase=4.519*)

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

(*Out= {{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[x_, y_] := 0;
VY[x_, y_] := 0;
\[CapitalRho][x_, y_] := 0;
TK[x_, y_] := 0; TX[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= {{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= {{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= {{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= 4.4880*)
In:= (*BenchmarkCase=4.519*)

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

(*Out= 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= 4.5252*)

(*BenchmarkCase=4.519*)

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

(*Out= {{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.

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.

Re question 3: There are NBodySimulation and NBodySimulationData.

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

• I'd be curious if someone could tell about the motivation behind these functions. Seems like a domain-specific wrapper around NDSolve might be better as a package than built-in. Apr 21, 2019 at 13:21
• @ChrisK, the motivation was Stephen wanted such a function. Apr 21, 2019 at 19:56
• @user21 Would such specialized functions also follow as LotkaVolterraSimulation, MalthusGrowth, NSpeciesCoexistence, GravitationalSlingshotSimulation, ExponentialReplicationSimulation, etc.? Apr 22, 2019 at 7:26
• @IstvánZachar, I do not understand your question. Could you please try to reformulate? Apr 22, 2019 at 7:36
• I think I'd prefer a more modular design that would separate the equation-setting stage (like NBodyEquations) from the solving stage (which could be NDSolve). Between this and SystemModels there are now three distinct clusters of functions for solving ODEs. Ah well, I found a couple of twitch videos on the topic to watch: twitch.tv/videos/229615178 and twitch.tv/videos/329322158 Apr 22, 2019 at 16:11