# How to obtain adaptive sampling as in Plot function?

Adaptive sampling in Plot function can capture the oscillation of a function with very few points. How can I get a similar sequence of point pairs without using Plot. You can somehow dig into Graphic object to get the points as following:

In:= Plot[Sin[x], {x, 0, 6*Pi}, PlotPoints -> 2][[1, 1, 3, 2, 1]]

Out= {{0.0000188496, 0.0000188496}, {0.294543,
0.290302}, {0.589066, 0.555585}, {0.88359, 0.773021}, {1.17811,
0.923886}, {1.47264, 0.995186}, {1.76716, 0.980782}, {2.06168,
0.881914}, {2.35621, 0.707097}, {2.65073, 0.471385}, {2.94526,
0.195078}, {3.23978, -0.0980295}, {3.5343, -0.382694}, {3.82883, \
-0.634402}, {4.12335, -0.831476}, {4.41787, -0.956943}, {4.7124, \
-1.}, {5.00692, -0.956938}, {5.30145, -0.831465}, {5.59597, \
-0.634387}, {5.89049, -0.382677}, {6.18502, -0.0980107}, {6.47954,
0.195096}, {6.77406, 0.471401}, {7.06859, 0.70711}, {7.36311,
0.881923}, {7.65764, 0.980786}, {7.95216, 0.995184}, {8.24668,
0.923879}, {8.54121, 0.773009}, {8.83573, 0.555569}, {9.13025,
0.290284}, {9.42478,
3.67394*10^-16}, {9.7193, -0.290284}, {10.0138, -0.555569}, \
{10.3083, -0.773009}, {10.6029, -0.923879}, {10.8974, -0.995184}, \
{11.1919, -0.980786}, {11.4864, -0.881923}, {11.781, -0.70711}, \
{12.0755, -0.471401}, {12.37, -0.195096}, {12.6645,
0.0980107}, {12.9591, 0.382677}, {13.2536, 0.634387}, {13.5481,
0.831465}, {13.8426, 0.956938}, {14.1372, 1.}, {14.4317,
0.956943}, {14.7262, 0.831476}, {15.0207, 0.634402}, {15.3153,
0.382694}, {15.6098,
0.0980295}, {15.9043, -0.195078}, {16.1988, -0.471385}, {16.4933, \
-0.707097}, {16.7879, -0.881914}, {17.0824, -0.980782}, {17.3769, \
-0.995186}, {17.6714, -0.923886}, {17.966, -0.773021}, {18.2605, \
-0.555585}, {18.555, -0.290302}, {18.8495, -0.0000188496}}


I think there should be a formal way to do this. I searched adaptive sampling in the documents, nothing interesting pops out.

## Conclusions

@Michael E2 provides an extensive answer for this question and similar questions. The ultimate solution for this problem is FunctionInterpolation, which can adjust accuracy of the result by changing PrecisionGoal and AccuracyGoal. The option MaxRecursion should be added, if high accuracy result is needed. It should also be noted MaxRecursion in this function does not has the limit of 15 as in Plot function.

Meanwhile, the solution in Transform an InterpolatingFunction is another very interesting way to solve this problem as well.

If the adaptive sampling in Plot function is needed, take the code from @george2079's answer.

• The key is what goal is guiding the adapting. With Plot, it is the angle between consecutive line segments of the plot, which is entirely local. With NIntegrate and NDSolve, it is the estimated error in the result; NIntegrate has both global and local adaptive strategies. – Michael E2 May 8 '15 at 14:13
• @MichaelE2 Why do you always know the internal implementation of those functions? – Kattern May 8 '15 at 14:18
• Aw, shucks. I've used it for 25 years, read the docs, attended conferences, read/participated on this site. Related: mathematica.stackexchange.com/a/59892 and its comments; mathematica.stackexchange.com/a/67642 – Michael E2 May 8 '15 at 14:34
• In fact, one neat trick for numerically evaluating an integral in Mathematica is to let Plot[] do the adaptive sampling and then use the output for the trapezoidal rule. – J. M.'s technical difficulties May 8 '15 at 18:39
• @MichaelE2 Here is why I am interesting in this problem. I am using a software, which require me to input a tabular data for calculation. Then the software will use the given value and do linear interpolation for its own calculation. What I want to do is to use minimum points of input to get a good approximation of the original function. This is why I am thinking about the adaptive sampling in Plot function. – Kattern May 9 '15 at 13:57

[Edit notice: I discovered a workaround to make the options PrecisionGoal and AccuracyGoal work.]

### FunctionInterpolation

I have had a prejudice against FunctionInterpolation. It does not really seem fully implemented. I've even asked people at Wolfram, who didn't seem to use it and even recommended using NDSolve. But it does do adaptive sampling. It takes PrecisionGoal, AccuracyGoal, MaxRecursion and InterpolationOrder options, as well as a few others. None are mentioned in the documentation, but one can find them with Options:

Options[FunctionInterpolation]
(*
{InterpolationOrder -> 3, InterpolationPrecision -> Automatic,
AccuracyGoal -> Automatic, PrecisionGoal -> Automatic,
InterpolationPoints -> 11, MaxRecursion -> 6}
*)


I'm not sure what InterpolationPrecision precision does; it does not seem to be the same as WorkingPrecision in other numerical solvers. The option InterpolationPoints controls the initial sampling; I say controls, because the actual initial sampling is not always the same. (For instance, when set to 11, 12, or 13, the initial sampling consists of 13 points in each case.)

Using the PrecisionGoal and AccuracyGoal with any settings other than Automatic results in errors, when the domain is given in terms of exact numbers. If the domain is given in terms of approximate real numbers (machine or arbitrary precision), then the options and FunctionInterpolation work. The default settings for PrecisionGoal and AccuracyGoal appear to be 6 or thereabouts. More importantly, it adapts its sampling to the InterpolationOrder and seems to do its job in meeting the precision+accuracy goal. (The goal is that the absolute error be less than 10^-acc + 10^-prec * Abs[f[x]].)

You can get the points from an InterpolatingFunction like this:

points = Transpose[{Flatten[if["Grid"]], if["ValuesOnGrid"]}]


Linear interpolation.

In all the examples, the gold curve shows the precision+accuracy goal. When the blue plot is below it, the precision+accuracy goal is met.

An objective function that oscillates increasingly rapidly:

obj[x_] = 10 + Sin[x^2]
if = FunctionInterpolation[obj[x], {x, 0., 10.}, MaxRecursion -> 15,
InterpolationOrder -> 1];
"Points" -> Length[if["Grid"]]
With[{prec = 6, acc = 6},
LogPlot[   (* error plot *)
Evaluate[Flatten[{Abs[if[x] - obj[x]], 10^-acc + 10^-prec*Abs[obj[x]]}]],
{x, 0., 10.}, PlotPoints -> 1000]
]

(*  "Points" -> 11634  *) An objective function that grows larger:

obj[x_] = Exp[x + Sin[4 x]];
if = FunctionInterpolation[obj[x], {x, 0., 10.}, MaxRecursion -> 15,
InterpolationOrder -> 1];
"Points" -> Length[if["Grid"]]
With[{prec = 6, acc = 6},
LogPlot[   (* error plot *)
Evaluate[Flatten[{Abs[if[x] - obj[x]], 10^-acc + 10^-prec*Abs[obj[x]]}]],
{x, 0., 10.}, PlotPoints -> 1000]
]

(*  "Points" -> 17603  *) The first example with PrecisionGoal -> 4, AccuracyGoal -> 4:

obj[x_] = 10 + Sin[x^2]
With[{prec = 4, acc = 4},
if = FunctionInterpolation[obj[x], {x, 0., 10.}, MaxRecursion -> 15,
InterpolationOrder -> 1, PrecisionGoal -> prec, AccuracyGoal -> acc];
Print["Points" -> Length[if["Grid"]]];
LogPlot[   (* error plot *)
Evaluate[Flatten[{Abs[if[x] - obj[x]], 10^-acc + 10^-prec*Abs[obj[x]]}]],
{x, 0., 10.}, PlotPoints -> 1000]
]
(*  Points->1191  *) Cubic interpolation.

InterpolationOrder -> 3 is the default:

obj[x_] = 10 + Sin[x^2]
if = FunctionInterpolation[obj[x], {x, 0., 10.}, MaxRecursion -> 15,
InterpolationOrder -> 3];
"Points" -> Length[if["Grid"]]
With[{prec = 6, acc = 6},
LogPlot[   (* error plot *)
Evaluate[Flatten[{Abs[if[x] - obj[x]], 10^-acc + 10^-prec*Abs[obj[x]]}]],
{x, 0., 10.}, PlotPoints -> 1000]
]

(*  "Points" -> 856  *) ### Other choices - Summary

A main advantage to these over FunctionInterpolation is being able to change the precision/accuracy goals of the interpolation. To me, there now seems no great advantage to the methods below; indeed, the interpolation produced by FunctionInterpolation seems excellent -- I'll have to rethink my prejudice. FunctionInterpolation is a bit slow by comparison with Plot and NDSolve. On the other hand, FunctionInterpolation allows higher settings of MaxRecursion, whereas in Plot, it is limited to 15 -- okay, not a severe restriction probably. Overall, I would say FunctionInterpolation is superior to the methods below.

Plot -- Precision control is a little tricky, but can be adequately managed with PlotPoints, MaxRecursion, and Method -> {Refinement -> {ControlValue -> bound}} (MaxBend is now deprecated), where bound is a bound on the angle in radians (see e.g. How does Plot work?). Plot tends to oversample unnecessarily in some neighborhoods with high settings of MaxRecursion, sampling well beyond the meeting the maximum bend bound between line segments. FunctionInterpolation produces fewer sample points meeting the precision and accuracy goal of 6.

NDSolve -- Good precision control. Usually reasonably fast. Seems to produce fewer sample points than Plot, more than FunctionInterpolation. See my answer, Transform an InterpolatingFunction, for two basic methods for interpolating a function with NDSolve.

NIntegrate -- Worst. I can't really imagine recommending it in any situation, unless the interpolation is to be integrated. While you can control the precision of the resulting integral, the precision of an interpolation is different (except when integrating it). The option setting Method -> {"GlobalAdaptive", Method -> {"TrapezoidalRule", "SymbolicProcessing" -> 0, "RombergQuadrature" -> False}} is the best you can do for approximating linear interpolation (that I found). Sometimes it manages to get good precision and sampling efficiency. It usually is much slower than any other method, and for what you get, it hardly seems worth it.

Manual global recursive subdivision -- Not too hard to code. If the function tamely oscillates, you can get good control over precision and good sampling efficiency.

Manual local adaptive step-size iteration -- By which I mean, choosing a step size at each step to meet the precision/accuracy goals (for linear interpolation, the error estimate is based on the second derivative). For tame functions, this can produce the best sampling faster than the other methods.

FWIW, I have code and data to back this up. Just ask. I've toyed with this problem on and off in my spare time over the last few weeks. It seems to come up in one way or another every now then on this site. What I have, though, has ballooned to the point that I'm embarrassed.

• A really cool answer, takes time to read. I am curious how can I get the method of InterpolatingFunction, like ["Grid"] in Transpose[{Flatten[if["Grid"]], if["ValuesOnGrid"]}]. I do not find them in the documents. – Kattern Jun 4 '15 at 3:02
• @kattern Thanks! Here are some references on InterpolatingFunction: mathematica.stackexchange.com/questions/28337/…, reference.wolfram.com/language/tutorial/…, mathematica.stackexchange.com/questions/19042/… – Michael E2 Jun 4 '15 at 3:09
• The answer is a bit confusing about whether the FunctionInterpolation can control the precision/accuracy goals. With a little test, I think Using the precision and accuracy options results in errors, so they're useless. seems accurate. – Kattern Jun 4 '15 at 7:12
• @kattern Thanks for your comment. For whatever reason, it made me think to try something that hadn't occurred to me. I seem to have found a way to make FunctionInterpolation work properly with the options. So you now can set the precision/accuracy goals and get an interpolation that is that precise/accurate (as defined in Mathematica). Let me know if you have trouble. – Michael E2 Jun 4 '15 at 12:43
• Additional information on undocumented syntax of FunctionInterpolation can be found in this thread. – Alexey Popkov May 16 '16 at 22:16

I know you said "without plot", but why reinvent the wheel?

  adaptivepoints=Sort@Last@Last@
Reap[Plot[f = Sin[x], {x, 0, 100},
EvaluationMonitor :> Sow[{f, x}]]];


You can similarly use NIntegrate, which might be more robust in terms of warning you if the adaptive scheme is not converging.

   g[x_?NumericQ] := Sin[x];
Sort@Last@Last@Reap@
NIntegrate[f = g[x], {x, 0, 1000},
EvaluationMonitor :> Sow[{x, f}]];


note for NIntegrate you should force numeric evaluation, otherwise it might employ some quadrature scheme that minimizes the number of eval points.

• OK, that seems a solution. – Kattern May 8 '15 at 13:55
• Why is NIntegrate more robust? It looks like they both controlled by MaxRecursion. – Kattern May 8 '15 at 14:02
• NIntegrate is (much) more robust in the criteria it uses to confirm convergence. Plot just looks for the result to be smooth, NIntegate looks for successive recursions of the integral to converge. – george2079 May 8 '15 at 14:32
• Thanks to share the trick on disable quadrature in NIntegrate. Theoretically, you also can do quadrature on numeric functions, how does this trick work? – Kattern May 8 '15 at 15:03
• @Kattern From the symbolic form of the integrand NIntegrate can detect some its features like periodicity in order to minimize the number of evaluation points. AFAIK it will apply symbolic until you switch it off with Method->{Automatic, "SymbolicPreprocessing"->None} (instead of Automatic may be explicit method specification) or by using the "black box" method (_?NumericQ). Both ways do not disable the quadrature scheme. – Alexey Popkov Jun 10 '15 at 4:58