# Evaluate function at "good" points used by Plot in numerical problems

Plot function with option Mesh->All shows how mathematica evaluates function to make it most optimal for plotting.

I'd like to evaluate some physical function in those points - in other words, I'd like to have more points where function behaves aggressively and only a few where it is constant. Is there any way to do it? Will it cost much computational time? (Plotting of functions is much slower than computing them, though I don't know the reason)

In the end instead of evaluating f[x] for x in Subdivide[0,1,n], I'm seeking to evaluate them in points selected by Mathematica algorithm

• From the docs: Plot initially evaluates f at a number of equally spaced sample points specified by PlotPoints. Then it uses an adaptive algorithm to choose additional sample points, subdividing a given interval at most MaxRecursion times . It would be interesting to know what the condition is to bisect a subdivision further; is it based on midpoint error, derivative, curvature etc. ? It's probably proprietary to Wolfram. You can always extract the curve from the plot by clicking on it though. As for computation, if you evaluate it at those exact points it should take about the same time. Jul 22, 2020 at 23:57
• Take a look at the answers to this question. Jul 23, 2020 at 0:43
• Jul 23, 2020 at 2:05
• What's your goal? Plot is not particularly good at picking numerical points. Jul 23, 2020 at 2:24

You can integrate a DAE for your function with NDSolve. I used a low-order integration rule to get dense sampling when the second derivative is large in magnitude. I used a low PrecisionGoal so that the number of points would be low, which helps the visualization below show where the sampling is denser. You can change the precision as desired.

approx = NDSolveValue[{
x'[t] == 1, x[0] == 0,       (* dummy DE *)
y[t] == Sin[3 t] - Sin[t]},  (* function to integrate *)
y, {t, 0, 2 Pi},
Method -> {"IDA", "MaxDifferenceOrder" -> 1},
PrecisionGoal -> 2, AccuracyGoal -> 3];

Plot[Sin[3 t] - Sin[t], {t, 0, 2 Pi},
Mesh -> {Flatten@approx@"Grid"}, MeshStyle -> Red]


You can get the function values with:

approx@"ValuesOnGrid"

(*  {0., 0.0002, ..., -0.0338323, -4.92661*10^-16}  *)


In principle you should be able to do this with FunctionInterpolation, but it's not super well documented.

As an example, the function Sin[x^2] becomes progressively more curved and therefore needs progressively denser sampling:

int = FunctionInterpolation[Sin[x^2], {x, 0, 10}, MaxRecursion -> 20, InterpolationOrder -> 2]


You can inspect the internals of the interpolation function by evaluating:

List @@ int


You'll notice that the x-values are stored in the 3rd argument and the y-values in the 4th. You can plot them like this:

ListPlot[Transpose @ {int[[3, 1]], int[[4, All, 1]]}]


The sampling density increases with x:

Histogram[int[[3, 1]], 50]


FunctionInterpolation has a number of options that control how it samples the function. You may need to tinker with these to get a good result (especially the MaxRecursion option often needs increasing). You can find them by evaluating:

Options[FunctionInterpolation]
`

{AccuracyGoal -> Automatic, InterpolationOrder -> 3, InterpolationPoints -> 11, InterpolationPrecision -> Automatic, MaxRecursion -> 6, PrecisionGoal -> Automatic}