Bug introduced in 10.3

I recently installed Mathematica v10.3 and checked a project I was working on few months ago under Mathematica v10.1.

In this project I solve a PDE with NDSolve and Finite Element Method over a 2D anulus region and I and ContourPlot that solution.

Under v10.1 (and v10.2) the solution appear like this:

10.1 solution

Under v10.3 the solution appear like this:


The source code is exactly the same.

At the beginning I was convinced the problem was with NDSolve (and this question reflected that, I later edited). After more investigation I discovered the solution in the two releases is nearly the same.

You can reproduce the problem running:


I searched on the documentation, news section, but I didn't find anything that could explain this or restore the previous behavior.

Any chance to understand the reason of the changed behavior and eventually restore the old?


The problem is not related to ContourPlot but to the InterpolatingFunction returned by NDSolve (thanks to JasonB). Probably the InterpolatingFunction is built under the cover by a call to ElementMeshInterpolation. Indeed the following code in Mathematica 10.3 returns result completely different from the results in Mathematica 10.2.


In 10.3:

Mathematica 10.3

In 10.2 (after som InterpolatingFunction::dmval message):

Mathematica 10.2

As you can see, in 10.2, the largest differences in "interpolation errors" are negligible, while in 10.3 they are even greather than the actual expected value.

The problem appear more specifically related with the "scale" of the mesh. Indeed if you change the sizes:


all apear fine:

Mathematica graphics

  • 5
    $\begingroup$ Unless you can post some kind of example it is hard to figure it out, sorry $\endgroup$
    – Jason B.
    Feb 18, 2016 at 10:51
  • 1
    $\begingroup$ @JasonB Of course, just asking if someone is aware of some change in the FEM framework between 10.2 and 10.3... I'll try to prepare a reduced example if possible... $\endgroup$
    – unlikely
    Feb 18, 2016 at 11:11
  • $\begingroup$ What loading/deformation? Plane stress 2 point, 4 point forces deformation/ stress plots? Solves an isotropic 2D stress function pde? $\endgroup$
    – Narasimham
    Feb 18, 2016 at 14:16
  • $\begingroup$ @Narasimham Yes, the plots shoud represents the modulus and horizontal component of the displacement vector; there is a load on the volume, horizontal outward on the x axis and vertical inward on the y axis, that gradually canges in-between; there is an essential homogeneus boundary condition on the inner border, and natural homogeneus on the outer. The material is composed of many isothropic anulus with differents Lamé coefficients. $\endgroup$
    – unlikely
    Feb 18, 2016 at 16:09
  • $\begingroup$ The stress/deformation results are essentially the same. Extra features introduced elsewhere perhaps reduce the display/appearance .. that you can ignore until WRI would correct later on. $\endgroup$
    – Narasimham
    Feb 18, 2016 at 16:21

2 Answers 2


Too long to make as a comment, so I'll post this non-answer as an answer.

I just want to point out that this definitely is not a ContourPlot problem, and it is absolutely a problem with the interpolating function.

So let's evaluate the following code in versions 10.2 and 10.31

{sol1, sol2} = testSolution;
grid = testSolution[[1]]["Grid"];
vals1 = testSolution[[1]]["ValuesOnGrid"];
vals2 = testSolution[[2]]["ValuesOnGrid"];
norm = Norm /@ Transpose[{vals1, vals2}];
normgrid = Transpose[Append[Transpose[grid], norm]];

So I have extracted the grid points from the interpolating functions, as well as the values on those grid points. I have checked, and they are the same across both versions to within an accuracy of 10^-17. So we should be dealing with the same interpolating functions. But we aren't.

Here is the nail in the coffin,

ListPlot[sol1 @@@ grid - vals1]
ListPlot[sol2 @@@ grid - vals2]

Evaluating this code in version 10.2, we get

enter image description here

which looks great. In 10.3.1 we get

enter image description here

which is just plain horrible. So what we have here is some backsliding in the interpolating function resulting from an NDSolve call using a nonrectangular region.

Now here is something interesting, if you take those grid points and the values on the grid, and you make an interpolation function using Interpolation, then you get basically the same plot in both versions.

 Interpolation[Transpose[{grid, norm}], InterpolationOrder -> 1][x, 
  y], {x, y} ∈ testMesh]

enter image description here

  • $\begingroup$ Paging @user21 - here is some more interpolation madness! $\endgroup$
    – Jason B.
    Feb 19, 2016 at 12:40
  • $\begingroup$ Thanks for your investigation. I was near to confirm the same results. So practically speaking, with ListPlot[sol1 @@@ grid - vals1] test you are proving that the InterpolatingFunction doesn't even interpolate on the grid points? $\endgroup$
    – unlikely
    Feb 19, 2016 at 13:43
  • $\begingroup$ @unlikely - yes!! The most basic thing you can expect from an interpolating function. $\endgroup$
    – Jason B.
    Feb 19, 2016 at 13:47
  • 1
    $\begingroup$ Thanks Jason; I think paging only works if I had a contribution on this particular site already; in any case I did not get a notification. But I certainly will have a look what goes wrong. $\endgroup$
    – user21
    Feb 21, 2016 at 3:01

Update: This is fixed in Version 11.0:

<< NDSolve`FEM`
{Ri, Ro} = {3000000., 6371000.};
rgn = RegionDifference[Disk[{0., 0.}, Ro], Disk[{0., 0.}, Ri]];
mesh = ToElementMesh[rgn];
grid = mesh["Coordinates"];
vog = (Exp[-Norm[#]^2/Ri^2]) & /@ grid;

if = ElementMeshInterpolation[{mesh}, vog];
TakeLargest[(if @@@ grid - vog), 20]

{4.996*10^-16, 4.85723*10^-16, 4.44089*10^-16, 3.88578*10^-16, 
 3.33067*10^-16, 3.33067*10^-16, 2.498*10^-16, 2.498*10^-16, 
 2.22045*10^-16, 2.22045*10^-16, 2.22045*10^-16, 2.22045*10^-16, 
 2.22045*10^-16, 1.94289*10^-16, 1.66533*10^-16, 1.66533*10^-16, 
 1.66533*10^-16, 1.66533*10^-16, 1.66533*10^-16, 1.66533*10^-16}

Old message:

This is a bug; you can workaround that by setting

SetSystemOptions[ "FiniteElementOptions" -> {"CacheInterpolationElements" -> False}]


It seems the default tolerance is not high enough for mesh over a large space compared to the default bounding box:

 "FiniteElementOptions" -> {"InterpolationToleranceFactor" -> 10^-8}]

This would allow you to still make use of cached elements which is faster for structured queries (which *Plot is not)

Sorry about that.

  • $\begingroup$ Sorry? Wait, is this your fault? Grab the pitchforks. :-) $\endgroup$
    – Jason B.
    Feb 21, 2016 at 18:19
  • $\begingroup$ @JasonB, sometimes things just get out of control.... $\endgroup$
    – user21
    Feb 22, 2016 at 20:28
  • $\begingroup$ Both workaround works to fix the aspect of my plots. But no one, apparently, fix the numerical problem exposed in the @JasonB answer (InterpolatingFunction that doesn't interpolate), at least in my specific case: ifn=Import["http://bit.ly/1ThEi3G"];ifn["ValuesOnGrid"]-ifn@@@ifn["ElementMesh"]["Coordinates"]//Abs//Max is not 0... It's 10^-2 w.r.t. the order of magnitude of ifn values on grid... $\endgroup$
    – unlikely
    May 17, 2016 at 9:57
  • $\begingroup$ @unlikely, what do you get 0?, what OS, and $Version? $\endgroup$
    – user21
    May 17, 2016 at 10:11
  • 1
    $\begingroup$ @user21 - I'm using version 10.4 Linux, and here are the results. You have to use the "InterpolationToleranceFactor" option from this answer to get the correct result. $\endgroup$
    – Jason B.
    May 17, 2016 at 11:27

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service and acknowledge you have read our privacy policy.

Not the answer you're looking for? Browse other questions tagged or ask your own question.