1

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

Michael E2
  • 235,386
  • 17
  • 334
  • 747
Kater Tatianow
  • 71
  • 3
  • 1
    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. – flinty Jul 22 '20 at 23:57
  • 1
    Take a look at the answers to this question. – Rohit Namjoshi Jul 23 '20 at 00:43
  • 2
    Related: https://mathematica.stackexchange.com/questions/82912/how-to-obtain-adaptive-sampling-as-in-plot-function – Michael E2 Jul 23 '20 at 02:05
  • What's your goal? Plot is not particularly good at picking numerical points. – Michael E2 Jul 23 '20 at 02:24

2 Answers2

2

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]

enter image description here

You can get the function values with:

approx@"ValuesOnGrid"

(*  {0.`, 0.0002, ..., -0.0338323, -4.9266110^-16}  )
Michael E2
  • 235,386
  • 17
  • 334
  • 747
0

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]]}]

enter image description here

The sampling density increases with x:

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

enter image description here

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}

This answer provides more detail.

Sjoerd Smit
  • 23,370
  • 46
  • 75