5

I need to perform a DensityPlot, but the figure is partially displayed. The codes and figure are shown below.

    d = 10^-8; 
    z[x_, y_] := Exp[-((1/x + 1/y - 1/400)/d)^2];
    DensityPlot[z[x, y], {x, 780, 820}, {y, 780, 820}, PlotPoints -> 100,  PlotRange -> All, ColorFunction -> (Blend[{Hue[2/3], Hue[0]}, #] &)]

enter image description here

The portions on the up-left corner and down-right corner are not displayed. I think the reason is that the system's precision is not high enough. The correct figure should be like this：

enter image description here

Is there any method to overcome this problem? Thank you!

user14634
  • 771
  • 6
  • 17

2 Answers2

6

Note: While this answer shows a good but arcane use of an option, my other answer shows a better use of the arcane internal structure of graphics to transform this problem from one which DensityPlot was not designed to handle into one which it was. IMHO. In other cases, if unable to solve for the transformation, then this approach may be the best option.

From Specific initial sample points for 3D plots and transforming density with Log:

d = 10^-8;
z[x_, y_] := Exp[-((1/x + 1/y - 1/400)/d)^2];
DensityPlot[Log[2^-1000 + z[x, y]],
 {x, 780, 820}, {y, 780, 820},
 PlotPoints -> {
   45,
   Table[{x, y /. First@Solve[(1/x + 1/y - 1/400) == 0]},
    {x, Subdivide[720, 820, 2 + 3 44][[2 ;; ;; 3]]}]},
 MaxRecursion -> 8, PlotRange -> All, 
 ColorFunction -> (Blend[{Hue[2/3], Hue[0]}, #] &), 
 WorkingPrecision -> MachinePrecision]

enter image description here

The trouble is that d is so small, that the fall-off from the ridge to underflow is probably less than a pixel in width. Even with sample points on the ridge, some of the ridge is missed (if you play with, say, PlotPoints, the sample points, transforming with Log, or WorkingPrecision). I couldn't hit on settings that gave a complete ridge line without the Log transformation. One can hide the log transformation with ScalingFunctions:

d = 10^-8;
z[x_, y_] := Exp[-((1/x + 1/y - 1/400)/d)^2];
DensityPlot[$MinMachineNumber + z[x, y],
 {x, 780, 820}, {y, 780, 820},
 PlotPoints -> {
   45,
   Table[{x, y /. First@Solve[(1/x + 1/y - 1/400) == 0]},
    {x, Subdivide[720, 820, 2 + 3 44][[2 ;; ;; 3]]}]},
 MaxRecursion -> 6, PlotRange -> All, 
 ColorFunction -> (Blend[{Hue[2/3], Hue[0]}, #] &), 
 WorkingPrecision -> MachinePrecision,
 ScalingFunctions -> "Log"]
Michael E2
  • 235,386
  • 17
  • 334
  • 747
  • Thanks a lot for this answer. It is very helpful. But this method doesnot work when d = 10^-13; Is it possible to solve this problem? @ Michael E2 – user14634 Jun 16 '23 at 13:57
  • I met the problem of "SystemException["MemoryAllocationFailure"]" when I use d = 10^-13 – user14634 Jun 16 '23 at 14:00
  • @user14634 It's a poorly scaled problem to begin with. My first idea was to plot the curve (1/x + 1/y - 1/400) == 0, which is in effect what you're doing even in the original problem. My second was to rescale the problem, if possible, or to make a log plot of z. But then the above (after several trials) worked. Now d = 10^-13 makes it harder. For instance, z[790.`16, 790.`16] cannot be computed. So all I can suggest is a manual log plot. DensityPlot[Evaluate@ComplexExpand@Log@z[x, y], {x, 780, 820}, {y, 780, 820}] and explain in the text of your paper if necessary. – Michael E2 Jun 16 '23 at 15:39
  • It's rare to find a problem modeling real world phenomena that cannot be done in standard double-precision floating point numbers. I really can't recall ever running into one that can't be done in Mathematica's arbitrary precision numbers, which can be as small as $10^{-10^{15}}$. (Of course, there are no limits on purely mathematical problem. But maybe you can break it down somehow. Scaling d by k is the same as raising z to the power k^2. So an explanation for d = 1 easily translates to an explanation for d = 10^-8 and d = 10^-13.) – Michael E2 Jun 16 '23 at 15:45
  • In my case, I need to calculate DensityPlot[z[x, y], {x, 799, 801}, {y, 799, 801}] for d = 10^-13. This is a real-world problem I encountered in my study of physics. – user14634 Jun 17 '23 at 06:45
  • I have tried DensityPlot[ Evaluate@ComplexExpand@Log@z[x, y], {x, 799, 801}, {y, 799, 801}, PlotRange -> All, ColorFunction -> (Blend[{Hue[2/3], Hue[0]}, #] &), PlotLegends -> Automatic], but the figure is the same for d=10^-13 and d=10^-23. Is it strange?@Michael E2 – user14634 Jun 17 '23 at 06:50
  • With bounds like 799, 801, I figured it was from physical science. You might want to take a minute and consider the formula for $z$. It has the form $z=\exp(u/d^2)=\exp(u)^{1/d^2}$, where $u=u(x,y)$. So $\log z=\log\exp(u/d^2)=u/d^2$. So the plot of $\log z$ has the shape for any $d$, only scaled by $1/d^2$. Add PlotLegends -> Automatic to DensityPlot[] and it will show the z scale to the side of the graphics. You might ask someone in your field about how to represent $z$. – Michael E2 Jun 17 '23 at 13:31
2

I thought this should be a different answer for two reasons: (1) This is a different approach, which shouldn't be buried in another answer and should be evaluated as an alternative on its own merits. (2) The other answer is over a week old and has been evaluated by the community, and it didn't feel right adding something unrelated to it. [Maybe a third reason: I'm rather happy with the result. I hope the OP is, too. :) ]

We can transform the coordinates so the the region of interest is parallel to a coordinate axis. At first I thought this would be hopeless because the exponential function decreases so rapidly. But it works well with plotting $\log z$. Adding a fixed color function, fixed by the smallest $d$, which gives the greatest dynamic range, allows comparison of plots for different values of $d$. Finally, $\log z$ varies over 37 orders of magnitude for $d=10^{23}$, a log-scale color function and bar legend help show the variation.

The coordinate transformation

We will transform from $(x,y)$ to $(u,w)$, create the plot, and transform back to $(x,y)$. The following shows what the transformation looks like.

Solve[{u == 1/x + 1/y - 1/400, w == 1/x - 1/y}, {x, y}]
(*
{{x -> 800/(1 + 400 u + 400 w), y -> 800/(1 + 400 u - 400 w)}}
*)

Solve[{u == 1/x + 1/y - 1/400, w == 1/x - 1/y}, {u, w}]
% /. {{x -> 780, y -> 780}, {x -> 780, y -> 820}, {x -> 820, 
   y -> 780}, {x -> 820, y -> 820}}
(*
{{u -> -(1/400) + 1/x + 1/y, w -> (-x + y)/(x y)}}


{{{u -> 1/15600, w -> 0}},           (* bounds on {u, w} )
 {{u -> 1/639600,  w -> 1/15990}},
 {{u -> 1/639600, w -> -(1/15990)}},
 {{u -> -(1/16400), w -> 0}}}
)


ParametricPlot[{800/(1 + 400 u + 400 w), 800/(1 + 400 u - 400 w)},
 {u, -(1/15600), 1/15600}, {w, -(1/15990), 1/15990}, 
 FrameLabel -> {x, y}]

Creating the plots

As mentioned above, the color scheme is proportion to $\log\log z$. Equal increments of $\log |\log z|$ since $0<z\le1$ correspond to equal changes in the argument to the color function. Actually, equal increments in $\log(1+|\log z|)$ are used so that a color is defined for $z=1$. We use this to create a bar legend to be used uniformly for all values of $d$. To ensure that the bar legend comprises the range of $z$, construct it using the smallest value of $d$ to be considered. Some adjustment of the scale (ticks) would be needed for different values of $d$.

(* color function, legend, and transformation *)
d = 10^-23;
z[x_, y_] := Exp[-((1/x + 1/y - 1/400)/d)^2];
cfs[z_] := 1 - Log10[1 - z]/38;
leg = BarLegend[{Blend[{Hue[2/3], Hue[0]}, Sow[1 - Log10[1 - #]/38]] &,
    Table[N@ComplexExpand@Log@z[x, y],
      {x, 780, 820}, {y, 780, 820}] // MinMax},
   Table[-10.^(5 k), {k, 7}], LegendLabel -> HoldForm[Log[z]]];
uw2yx = Function @@ {{u, w}, {x, y} /. 
     First@Solve[{u == 1/x + 1/y - 1/400, w == 1/x - 1/y}, {x, y}]};

(* the plots *)
GraphicsRow[
 Table[
  Show[
   DensityPlot[
     Evaluate[
      ComplexExpand@Log@z[x, y] /. {x -> 800/(1 + 400 u + 400 w), 
        y -> 800/(1 + 400 u - 400 w)}],
     {u, -(1/15600) (1 + 4/100), 
      1/15600 (1 + 4/100)}, {w, -(1/15990) (1 + 4/100), 
      1/15990 (1 + 4/100)},
     MaxRecursion -> 6,
     PlotRange -> All,
     ColorFunction -> (Blend[{Hue[2/3], Hue[0]}, cfs[#]] &),
     ColorFunctionScaling -> False,
     WorkingPrecision -> 24,
     PlotLegends -> leg, 
     PlotLabel -> 
      Row[{HoldForm[d], " = ", 
        HoldForm[10]^
         Log10[d]}]] /. {pp_List?(MatrixQ[#, NumericQ] &) /; 
       MatchQ[Dimensions[pp], {_, 2}] :> 
      Transpose[uw2yx @@ Transpose@pp]},
   PlotRange -> {{780, 820}, {780, 820}},
   ImageSize -> 300],
  {d, {10^-8, 10^-13, 10^-23}}],
 ImageSize -> 1300]

Mathematica graphics

Note that the line of interest corresponds to $u=0$, the center of the plot domain for u. The default PlotPoints is odd and so the center is automatically included in the sampling.

Michael E2
  • 235,386
  • 17
  • 334
  • 747