# Numerical Issues Using Results of NDSolveValue

I am solving Laplace's Eqn. (Cartesians) in 2D in a region which is infinite in x, semi-infinite in y, and bounded by a function of x (h[x] in my code) close to y=0. In order to keep things under control, I set a "large" region in which to solve, i.e. |x|<12, y<12, (setting this to Infinity results in error messages) and I am interested only in the region |x|<4, y<2.5. The first NDSolveValue itself seems to work quite well, but I differentiate the result to make a couple of figures. The second results in numerical issues, as you can see.

h[x_] = 1/(1 + x^2);
yMnSlv = 0;
rng = 12;
regSlv = ImplicitRegion[{y > h[x]}, {{x, -rng, rng}, {y, yMnSlv,
rng}}]; psi =
NDSolveValue[{D[u[x, y], x, x] + D[u[x, y], y, y] == 0,
DirichletCondition[u[x, y] == 0, y == h[x]],
DirichletCondition[u[x, y] == y, x == -rng],
DirichletCondition[u[x, y] == y, x == rng],
DirichletCondition[u[x, y] == y, y == rng]},
u, {x, y} \[Element] regSlv, PrecisionGoal -> 20];
{xMn, xMx, yMn, yMx} = {-4, 4, 0, 2.5};
vx = Evaluate[D[psi[x, y], y]];
vy = Evaluate[-D[psi[x, y], x]];
v = {vx, vy};
StreamPlot[v, {x, xMn, xMx}, {y, yMn, yMx}, PlotRangePadding -> None,
psixx = Evaluate[D[psi[x, y], x, x]];
regPl = ImplicitRegion[{y > h[x]}, {{x, xMn, xMx}, {y, yMn, yMx}}];
ContourPlot[psixx == 0, {x, y} \[Element] regPl, PlotRange -> All,
PlotLegends -> Automatic, AspectRatio -> 0.8]


The stream plot is good, the contour plot is poor. This is physics, and the contour plot should be smooth. Moreover, my input function is even, so that its deviations from mirror symmetry (about x=0) are a measure of the inaccuracies due to numerical issues. I have tried adjusting the PrecisionGoal and using AccuracyGoal, but they result only in slowing down the code and no improvement. Can anyone suggest a way around these precision issues?

Thanks.

• Note that the max precision goal achievable using machine precision numbers is a little less than \$MachinePrecision, which is about 15.95. – Michael E2 Jan 22 at 13:24

You need to control the element mesh used to generate the solution. Add the option

Method -> {"FiniteElement",
"MeshOptions" -> {"MaxCellMeasure" -> 0.001}}


The error is reduced:

ContourPlot[psixx == 0, {x, y} \[Element] regPl, PlotRange -> All,
PlotLegends -> Automatic, AspectRatio -> 0.8]


Here is the mesh created with the default options:

psi["ElementMesh"]@"Wireframe"


With "MaxCellMeasure" -> 0.001, it looks solid black (not shown).

• One of the issues is that FEM does not have an error estimator that it can use to refine the mesh. The user can run the integration over two meshes and compare the max difference between the solutions for evidence of convergence. As a rule of thumb, to get a smooth plot over a domain of diameter ~3, I'd expect to have a mesh with a max length around 3/100 = 0.03 or "MaxCellMeasure" -> {"Length" -> 0.03}. That produces a solution that is 200+ MB. While that might be enough for a smooth plot, the solution still suffers from truncation (discretization) error. – Michael E2 Jan 22 at 13:22
• Even if the FEM did have an error estimator that would not help with the derivative recovery (psixx) of order two that is needed here. What might help here, I think, is a MeshOrder->3 but even then this will be hard / need a fine mesh. – user21 Jan 22 at 13:54
• @user21 Oops, I didn't even look at what was being plotted. How embarrassing. Thanks for the comment. :) (I'll update later with your observation later, I hope, when I have time.) – Michael E2 Jan 22 at 13:56
• Thanks @Michael E2 for the answer. I have adopted the approach you suggested and experimented with the values. It is working much better now. – Paul Harrison Jan 22 at 23:14