How to refine ContourPlot3D result?

Question

From this interesting question How to plot a surface given implicitly by 3 equations in 5 coordinate variables? I have found a subproblem that I would like to clarify:

Often it is advantageous to use ContourPlot3D to solve (visualize) an equation eq[yh,yd,xh]==0

eq = 0 == -1 + (2 (1 + yh))/(xh^2 + (1 + yh)^2) - (
   2 xh (1 + yh) ((-3 + xh)^2 + (yd + yh)^2))/((-3 + 
      xh) (xh^2 + (1 + yh)^2)^2) - 
   Log[200] - (4 yh (-9 xh + 6 xh^2 - xh^3 + 3 yd - xh yd - 
        3 xh^2 yd + xh^3 yd + xh yd^2 + 3 yh - 10 xh yh + 3 xh^2 yh + 
        6 yd yh + xh yd^2 yh + 6 yh^2 - xh yh^2 + 3 yd yh^2 + 
        xh yd yh^2 + 
        3 yh^3) Log[-((800 yh^2 (-9 xh + 6 xh^2 - xh^3 + 3 yd - 
              xh yd - 3 xh^2 yd + xh^3 yd + xh yd^2 + 3 yh - 
              10 xh yh + 3 xh^2 yh + 6 yd yh + xh yd^2 yh + 6 yh^2 - 
              xh yh^2 + 3 yd yh^2 + xh yd yh^2 + 3 yh^3))/((-3 + 
              xh) (yd + yh) (1 + xh^2 + 2 yh + yh^2)^2))])/((-3 + 
        xh) (yd + yh) (1 + xh^2 + 2 yh + yh^2)^2)

ContourPlot3D shows a smooth solution surface

pic = ContourPlot3D[
Evaluate[ eq ] , {yh, 0, 5}, {yd, 0, 5}, {xh, -2.5, 3 - .01}, 
ColorFunctionScaling -> False (*MaxRecursion\[Rule]4*) , 
Mesh -> False , PerformanceGoal -> "Quality" ]
pts = pic[[1, 1]][[1]]; (* save points from pic*)

But checking the solution subsequently

ptsF = Select[pts, ! (eq /. Equal -> Subtract /. {yh -> #[[1]], yd -> #[[2]],xh -> #[[3]]})^2 < .01 &];
Show[{pic , ListPointPlot3D[ptsF, PlotStyle -> Red]},AxesLabel -> {yh, yd, xh}]

shows parts (red points) in the solution surface where eq==0 isn't fullfilled. (To be precise, the squared residual error in the red zone is larger than 0.01. )

I know how to improve the quality of ContourPlot3D using options PlotPoints and MaxRecursion but the basic problem remains.

My question:

How can I force better results from ContourPlot3D?

Thanks!

Do you mean you want to make the "solution" found by ContourPlot3D numerically accurate? — xzczd, Jan 03 '24 at 11:41

xzczd · Answer 1 · 2024-01-03T15:01:17.717

First of all, it's worth pointing out that, setting a larger PlotPoints or MaxRecursion i.e. using more points for visualization improves not only the visualization, but also the numerical accuracy, it's just rather inefficient. And if we don't consider increasing plot points as a choice, then it's not possible to improve the numerical accuracy by adjusting option of ContourPlot3D AFAIK. (I sincerely hope I'm wrong. ) I think it's reasonable: ContourPlot3D is just for visualization, so the points used for creating the surface only need to be visually accurate i.e. it's not necessary to make the "solution" for visualization strictly fullfill the plotted equation.

Still, it's possible to improve the "solution" found by ContourPlot3D. The idea is simple: we solve eq with FindRoot, using points in pic as initial guesses.

The code is as follows. Definition of eq is the same as yours.

(* Modify pic slightly for better illustration. *)
pic = ContourPlot3D[Evaluate[eq], {yh, 0, 5}, {yd, 0, 5}, {xh, -2.5, 3 - .01}, 
  Mesh -> False, PerformanceGoal -> "Quality", ContourStyle -> Opacity[0.5]]
eqfunc[yh_, yd_, xh_] = eq;
(* More general definition for pts, compatible with earlier versions. *)
pts = Cases[pic, GraphicsComplex[a_, __] :> a, ∞][[1]];
ClearAll@calibrate;
calibrate[yh_, yd_, xh_] := {yh, yd, xhvar} /. 
  FindRoot[eqfunc[yh, yd, xhvar], {xhvar, xh}]
ptsbetter = calibrate @@@ pts; // AbsoluteTiming
(* There'll be some warnings, don't worry. )
( {29.324, Null} *)
picbetter = pic /. pts -> ptsbetter

FindRoot spits out warnings, and the picbetter is visually the same as pic. Is ptsbetter really improved? Let's check the error:

error = ptsbetter\[Transpose] // Apply@eqfunc // Last // Abs // Sort;
error[[-10 ;;]]
(* {0.0859047, 0.0859366, 0.106304, 0.108949, 0.119662, 
    0.120376,  0.13304,   0.319378, 0.355189, 0.361554} *)
thre = Function@Evaluate@Simplify`PWToUnitStep@PiecewiseExpand@If[# > 0.01, 1, 0]
picadd = ListPointPlot3D[Pick[ptsbetter, thre@error, 1]];
picbetter~Show~picadd

Blue points are those whose absolute residual error are larger than 0.01. (Notice I haven't added the ^2 i.e. this is amount to setting <.0001 in your code. ) As we can see, Accuracy of ptsbetter is indeed better.

Thanks for your answer, very helpful. Actual I`m trying to improve the surface "in normal direction" but still didn't succeed — Ulrich Neumann, Jan 03 '24 at 15:10

score 3 · Answer 2 · answered Jan 03 '24 at 16:31

Similar to xzczd's helpful answer it's indead possible to improve the accuracy of the solution points found by ContourPlot3D

For this we use the knowledge of VertexNormals vn given by ContourPlot3D and try to improve the solution in this direction.

eqfunc[yh_, yd_, xh_] = eq; 
vn = Cases[pic, 
GraphicsComplex[__, VertexNormals -> q_] :> q, -1][[1]];
improve[ p_, n_] :=p + alfa n /. FindRoot[Apply[eqfunc, p + alfa n], {alfa, 0}]
ptsN = MapThread[ improve[#1 , #2 ] &, {pts, vn}, 1]; // Quiet (* new points *)

Surprisingly no "red points" left

ptsF = Select[ptsN, ! (eq /. Equal -> Subtract /. {yh -> #[[1]], yd ->#[[2]], xh -> #[[3]]})^2 < .01 &];
(* {} *)

All new points ptsN fullfill the equation accurately

eqo = eq /. Equal -> Subtract;
Map[eqo /. {yh -> #[[1]], yd -> #[[2]], xh -> #[[3]]} &,ptsN] // MinMax
(*{-7.65864*10^-10, 4.83533*10^-11}*)

How to refine ContourPlot3D result?

2 Answers2