In this question, the equation $uu'=\nu u''$, $x\in (-1,1)$, $u(-1)=1+\delta$, $u(1)=-1$, was considered. The question asked about how to solve this differential equation problem in Mathematica. There were two answers. The first one suggested using the option Method -> FiniteElement in NDSolve for version 12.0. The second solution proposed the trapezoidal finite difference method, which was programed from scratch.
Let us analyze the approximations at $x=0.5$, for $\nu=0.05$ and $\delta=0.1$:
test = ParametricNDSolve[{u''[x] - (1/nu)*u[x]*u'[x] == 0,
u[-1] == 1 + delta, u[1] == -1}, u, {x, -1, 1}, {nu, delta},
Method -> FiniteElement];
u[0.05, 0.1][0.5] /. test
(* 1.09905 *)
fdmODE[nu_, delta_, n_] :=
Module[{h, mesh, f, u, v, eqns, sv, froot, sol},
h = 2/n;
mesh = -1 + h*Range[0, n];
f[{u_, v_}] = {v, (1/nu)*u*v};
eqns =
Flatten[Join[{u[0] == 1 + delta, u[n] == -1},
Table[Thread[{u[i], v[i]} == {u[i - 1], v[i - 1]} +
0.5*h*(f[{u[i - 1], v[i - 1]}] + f[{u[i], v[i]}])], {i, 1,
n}]]];
sv = Flatten[Table[{{u[i], 0}, {v[i], 0}}, {i, 0, n}], 1];(*
starting value, initial guess *)
froot = FindRoot[eqns, sv];
sol = Table[u[i], {i, 0, n}] /. froot;
Return@Thread[{mesh, sol}];
];
fdmODEinterp[nu_, delta_, n_, x_] :=
Interpolation[fdmODE[nu, delta, n], InterpolationOrder -> 1][x]
fdmODEinterp[0.05, 0.1, 10000, 0.5]
(* 1.09923 *)
fdmODEinterp[0.05, 0.1, 11000, 0.5]
(* 1.09923 *)
fdmODEinterp[0.05, 0.1, 12000, 0.5]
(* 1.09923 *)
It seems that Method -> FiniteElement is not sufficiently accurate. Is there any better option to tackle this problem in Mathematica?
