The following code plots the difference between two functions in a contour plot.

kmsq = om^2/
   2*( -ma^2/om^2 + n^2 + 1 - Sqrt[(1 ma^2/om^2 + n^2 - 1)^2 + x^2]);

km = Sqrt[kmsq];

kax = Sqrt[om^2 - ma^2];

fp[m_, g_] := 
  Log10[Abs[(kax - km)] /. {om -> 1.0, ma -> m, n -> 2.0, 
     x -> 2.0*194.0*g/1.0*10^(-9.0)}], 500]

ContourPlot[fp[m, g], {m, 10^(-10), 10^1}, {g, 10^(-10), 10^5}, 
 ScalingFunctions -> {"Log", "Log"}, PlotPoints -> 60, 
 PlotLegends -> Automatic]

The result is: enter image description here

For some points of the contour plot the two functions can become very similar, e.g. the difference that I plot can become very small. In some parts of the parameterspace (white areas of the contour plot) the difference becomes zero what should not be the case. I guess this is a floating point problem?

How can I resolve this? I already tries to increase the Precision, but it did not work.

Many thanks in advance!


As the OP seems to suspect, the issue arises from the "rounding plateau" around a relative error of $10^{-16}$ of standard (double) floating-point numbers. Increasing precision is the right move, but you have to do that before computing the function value. The OP's code has SetPrecision wrapped around the expression that computes the function value. The computation will be completed at the precision of the components, some of which are machine-precision numbers, before the precision is raised. You need to set the precision of the components before the computation starts. In particular, the input arguments and the parameter settings all need high precision. (For instance floats in the parameters like 1.0 need to be changed to exact numbers like 1.)

ClearAll[ma, n, om, x, fp];
kmsq = om^2/
    2*(-ma^2/om^2 + n^2 + 1 - Sqrt[(1 ma^2/om^2 + n^2 - 1)^2 + x^2]);
km = Sqrt[kmsq];
kax = Sqrt[om^2 - ma^2];

fp[m0_, g0_] := 
 Block[{m = SetPrecision[m0, 32], g = SetPrecision[g0, 32]},
  Log10[Abs[(kax - km)] /. {om -> 1, ma -> m, n -> 2, 
     x -> 2*194*g/1*10^(-9)}]

ContourPlot[fp[m, g], {m, 10^(-10), 10^1}, {g, 10^(-10), 10^5}, 
 ScalingFunctions -> {"Log", "Log"}, PlotLegends -> Automatic]

enter image description here

  • $\begingroup$ If I just copy your code, I still get the problem. Maybe there is a typo in the code? $\endgroup$ – jan schütte-engel Oct 15 at 12:23
  • $\begingroup$ @janschütte-engel Did you copy the OP's definitions of kmsq etc.? It works on a new kernel for me, when I include the rest of the OP's code. (It seems typical of the site often to show only what's been changed.) $\endgroup$ – Michael E2 Oct 15 at 12:37
  • $\begingroup$ Yes sure. I just substituted the code in the answer to my code. The result looks better, but there a still white holes in the plot. And it still does not go down to 10**-30 $\endgroup$ – jan schütte-engel Oct 15 at 12:42
  • $\begingroup$ I realize, that when you write in Mathematica for example "2.0" and copy it to the code window here you can see just "2" and not "2.0" . Maybe it i due to that issue? In my mathematica code I have written for example "kmsq=om^2.0/2.0*...". But this is not visible in the OP... $\endgroup$ – jan schütte-engel Oct 15 at 12:50
  • $\begingroup$ Resolved now. Thanks a lot for the help! $\endgroup$ – jan schütte-engel Oct 15 at 12:50

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service, privacy policy and cookie policy

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