[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 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.
Plot, it is the angle between consecutive line segments of the plot, which is entirely local. WithNIntegrateandNDSolve, it is the estimated error in the result;NIntegratehas both global and local adaptive strategies. $\endgroup$Plot[]do the adaptive sampling and then use the output for the trapezoidal rule. $\endgroup$Plotfunction. $\endgroup$