# How to improve this plot?

In general, the quality of Mathematica graphics is beyond praise. However, the output of Plot3D command is somewhat unexpected in the following code. Let's solve the Dirichlet problem for the Laplace equation.

NDSolve[{-Laplacian[u[x, y], {x, y}] == 0, DirichletCondition[u[x, y] ==
Boole[y >= 0], True]}, u, {x,y} ∈ Disk[]];Plot3D[Evaluate[u[x, y] /. %],
{x, y} ∈ Disk[],PlotPoints -> 100, PerformanceGoal -> "Quality"]

One sees the superfluous two peaks about (-1,0) and (1,0): it is well known that $u[x,y] \ge 0$ and $u[x,y] \le 1$. Also the exact solution through the Poisson integral formula https://en.wikipedia.org/wiki/Poisson_kernel confirms it.

What causes these peaks? How to get rid of the ones?

• The problem is that NDSolve returns an InterpolatingFunction, which will not be the exact solution to your PDE. Oct 23, 2016 at 11:10
• @ kevin : Thank you. The question is still open and the problem seems to be serious. Oct 23, 2016 at 11:35
• With a strong exposition and a relatable protagonist :-) Oct 23, 2016 at 16:11

I suspect the problem arises from too coarse a mesh when the Disk region is discretized. A better result is obtained if the region is created explicitly with a finer mesh at the edge.

region = DiscretizeRegion[Disk[],
MeshRefinementFunction -> Function[{vertices, area},
area > 0.005 (1 - Norm[Mean[vertices]]^2)]]

sol = NDSolveValue[{-Laplacian[u[x, y], {x, y}] == 0,
DirichletCondition[u[x, y] == Boole[y >= 0], True]},
u, {x, y} ∈ region]

Plot3D[sol[x, y], {x, y} ∈ region, PlotPoints -> 100, PerformanceGoal -> "Quality"]

• @ Simon Woods : Thank you for your constructive answer which solves the problem. Is your region the unit disk? If it is not, your answer is a fake. Oct 23, 2016 at 12:14
• @user64494 It's obviously the unit disk. Look up Disk[], which you used in your answer. Oct 23, 2016 at 12:16
• @user64494, it's an approximation to the unit disk of course, because it's a mesh of triangles. The same applies when you use Disk[] as the region, but in that case the approximation is worse. Oct 23, 2016 at 12:20
• @ Michael E2 : Looking in reference.wolfram.com/language/ref/DiscretizeRegion.html , I am not sure about that. Oct 23, 2016 at 12:20
• @user64494 NDSolve discretizes Disk[] using the same underlying algorithm, but with automatically chosen parameters. What Simon did was alter a parameter to give a finer mesh near the boundary. The discretization in Simon's code is more accurate than the automatic discretization used in yours. Neither, strictly speaking, is exactly the unit disk. Oct 23, 2016 at 14:24

## Interpolation error

The overshoots are the unavoidable result of interpolation. NDSolve computes the values of u of the DirichletCondition to machine-precision accuracy and approximates the values of u at other points in the mesh via the finite element method. Values at intermediate points are interpolated by polynomials that equal the values of u at the mesh points and approximate the derivatives of u at those points. Consequently, there will be some error. Perhaps surprisingly, the maximum error is about 1/8 in both the OP's and Simon Woods' solutions.

When u[x, y] is smooth, the error should be small. In the OP's case, there is a discontinuity, and approximating a discontinuity with a polynomial has limitations, reminiscent of Gibbs phenomenon for trigonometric polynomials (Fourier series). Mathematica has built-in methods for dealing with interpolation of discontinuities in 1D, but as far as I know, not for 2D or higher.

Here is a 1D illustration of the phenomenon, that is similar to what happens in any numerical solution of the OP's PDE produced by NDSolve[]. We interpolate the UnitStep[] function. Using Plot[] will result in interpolating the function between interpolation nodes. However ListLinePlot[] uses just the nodes. We can see the overshoot in the Plot[] output.

ifn = Interpolation[Transpose@{#, UnitStep[# - 1/4]} &@Range[-1., 1., 1/2^4]];
GraphicsRow@{Plot[ifn[x], {x, -1, 1}], ListLinePlot@ifn}

Moreover, the overshoot does not go away or diminish in magnitude as the mesh is refined. All that happens is the width of the overshoot gets smaller:

Plot[Evaluate@ Table[Interpolation[
Transpose@{#, UnitStep[# - 1/4]} &@Range[-1., 1., 1/2^n]][x], {n, 3, 6}],
{x, -1, 1}]

At some point, the bump can be so narrow that Plot[] misses it (unless you increase PlotPoints accordingly):

GraphicsRow@{
Plot[Evaluate@ Interpolation[
Transpose@{#, UnitStep[# - 1/4]} &@Range[-1., 1., 1./2^14]][x],
{x, -1, 1}],
Plot[Evaluate@ Interpolation[
Transpose@{#, UnitStep[# - 1/4]} &@Range[-1., 1., 1./2^14]][x],
{x, -1, 1}, PlotPoints -> 512]
}

There are other ways to control the overshoot in the 1D case, but I think those options do not exist for higher dimension, finite-element solutions.

At this point the question might be divided in two: Is the concern with the values of u computed by NDSolve[], or with the interpolation error?

If the latter, the interpolation error can be constrained to a small region by a finer mesh around the discontinuities, as Simon Woods shows. The error is still there (about 1/8 in magnitude), but Plot3D[], which produces a graphical approximation to the numerical approximation, misses it. Thus the plot looks good. This confines the error to a smaller area and can be considered a significant improvement.

If the former, there is a way to visualize the solution values, similar to ListLinePlot[] in the 1D case. The FEM function NDSolveFEMElementMeshPlot3D plots a 2D ElementMesh interpolating function. In the code for the graphics just below, solOP and solSW represent the InterpolatingFunction solutions from the OP and Simon Woods respectively. (Full code is given further down.) The function addnormals[u] adds the VertexNormals to the surface defined by z == u[x, y]. We can see that there is no overshoot in the plots of the computed values of u[x, y] of either solution.

Needs["NDSolveFEM"];

With[{plot = Show[ElementMeshPlot3D[solOP], Axes -> True]},

With[{plot = Show[ElementMeshPlot3D[solSW], Axes -> True]},

(* add  VertexNormals  to  plot  assumed to be the plot of the function  u  *)
addnormals[u_][plot_Graphics3D /; ! FreeQ[plot, GraphicsComplex]] :=
plot /. GraphicsComplex[p_, rest___] :> GraphicsComplex[
p, rest,
VertexNormals ->
With[{dx = Derivative[1, 0][u], dy = Derivative[0, 1][u], xy = p[[All, {1, 2}]]},
Transpose@{-dx @@@ xy, -dy @@@ xy, ConstantArray[1., Length@p]}]
]

## Further analysis of the interpolation error

We made some claims about the error above, for which we can give some evidence now.

First the code for the solutions: Simon Woods' refinement function may be passed to the underlying mesher via the Method option of NDSolve as shown below. It produces a slightly different mesh; but the mesh has similar properties to Simon's region and has virtually the same numerical properties shown below for solSW. I think this method should be preferred, since the FEM algorithms are designed for solving PDEs.

solOP = u /. First@NDSolve[{-Laplacian[u[x, y], {x, y}] == 0,
DirichletCondition[u[x, y] == Boole[y >= 0], True]},
u, {x, y} ∈ Disk[]];
solSW = NDSolveValue[{-Laplacian[u[x, y], {x, y}] == 0,
DirichletCondition[u[x, y] == Boole[y >= 0], True]},
u, {x, y} ∈ region,
Method -> {"FiniteElement",
"MeshOptions" -> {MeshRefinementFunction ->
Function[{vertices, area}, area > 0.005 (1 - Norm[Mean[vertices]]^2)]}}];

We can see that the range of computed values of u are within the theoretical range:

MinMax[solOP["ValuesOnGrid"]]
MinMax[solSW["ValuesOnGrid"]]
(*  {0., 1.}  *)
(*  {0., 1.}  *)

Next let's check the boundary for interpolation error. We get the mesh nodes and evaluate the solution at the midpoints between them. This leads to an estimate of the error. Note that its magnitude is about the same as the UnitStep function above, around 1/8 in magnitude in both the OP's and Simon's solutions.

coordsOP = SortBy[ArcTan @@ # &]@               (* sort into path order *)
MeshCoordinates@ RegionBoundary@ MeshRegion[solOP["ElementMesh"]];
midpointsOP = Mean /@ Partition[coordsOP, 2, 1];
valuesOP = solOP @@@ midpointsOP;
MinMax[valuesOP] - {0, 1}
(*  {-0.124918, 0.124918}  *)

coordsSW = SortBy[ArcTan @@ # &]@               (* sort into path order *)
MeshCoordinates@ RegionBoundary@ MeshRegion[solSW["ElementMesh"]];
midpointsSW = Mean /@ Partition[coordsSW, 2, 1, 1];
valuesSW = solSW @@@ midpointsSW;
MinMax[valuesSW] - {0, 1}
(*  {-0.124871, 0.124888}  *)
• @ Michael E2: Many Mathematica users derive benefit and pleasure from your excellent answer. I'd like to notice that is known as the Gibbs phenomenon ( see en.wikipedia.org/wiki/Gibbs_phenomenon). Oct 23, 2016 at 19:16
• @user64494 Thanks, for the compliment and for pointing out the typo. Oct 23, 2016 at 19:20
• +1 This is a really good answer. I must admit I hadn't realised that the magnitude of the interpolation error was constant as the mesh is refined. Oct 23, 2016 at 21:25
• Great answer! And (slightly belated) congratulations on reaching 100k! That represents a tremendous amount of your time spend on behalf of the community. Thank you! Oct 24, 2016 at 7:07
• @Mr.Wizard, Thanks, x2! Oct 24, 2016 at 22:47