# Function Interpolation with Automatic / Algorithmic Values Mesh

I have an underlying function f(x,y,z) that is computationally intensive, but is smooth and continuous. I'm needing to find the function values along a line in xyz. Currently, I'm calculating f at discrete steps and I'm wondering if there is a way to get Mathematica to automatically chose the step size and work out an interpolating function based on some inputs like how accurate I need it to be and what the initial step size should be. I'm thinking it would kind of work like the mesh function does in plotting by evaluating functions more often in areas of higher complexity.

Does anyone have any ideas where to start?

-
Have you tried FunctionInterpolation? It takes many of the same options as Interpolation. I find I often have to use them since automation fails for all but simple functions however. It might be helpful to give your function even if it is complicated. – Andy Ross Mar 31 '12 at 4:11
Please see this question and the answer. I believe this is a duplicate, although I'd wait to close it as such. Let me know if that helped. – R. M. Mar 31 '12 at 4:15
This question may be useful too. – Sjoerd C. de Vries Mar 31 '12 at 6:48
Based on its documentation, FunctionInterpolation definitely should do what Chris is asking for. However, as Andy mentions, it's very brittle. See for example this question. Still, playing around with its options is probably less work than reverse-engineering the Plot function. – Jess Riedel May 16 at 22:00

As far as I know, there is no built-in function to do this. However, what you can do is a heavy abuse of Part to extract the points from a Plot object:

g = Plot[Sin[x], {x, 0, 2 Pi}]


The InputForm, i.e. how Mathematica sees this picture internally, looks like like this:

Graphics[{{{}, {}, {Hue[0.67, 0.6, 0.6], Line[{{1.2*^-7, 1.2*^-7}, [long list of points]} ... (options etc)


You can now extract these points by using Part, i.e. [[ ]]:

points = g[[1, 1, 3, 2, 1]]

{{0., 0.}, {0.01, 0.01}, {0.02, 0.02}, {0.03, 0.03}, {0.06, 0.06}, ...


This can now be used for a ListPlot to visualize the points extracted:

ListPlot[points]


The density of points and the recursion can be set usual when plotting, i.e. PlotPoints and MaxRecursion.

Remarks:

1. I've rounded the values above heavily so it doesn't blow up the answer. Usually the numbers are much more precise.
2. Bear in mind that this does extractions from Plot that are not supposed to be done. For this reason, you should be especially careful setting plotting parameters and things like these, as I could imagine that some of them may change the structure of the Graphics object you're extracting the data from.
-
Here are some issues I see. First, Plot will give numbers to machine precision and only machine precision. Second, Plot cannot be used to find values for a function f(x,y,z) assuming this is a plot in 4 dimensions. Finally, any attempt to extract parts from the input form of an object can change from version to version so don't expect it to work if you upgrade. All that being said, I do this myself all the time +1 :) – Andy Ross Mar 31 '12 at 4:22
I agree, it's a pretty hacky workaround. I really wish they included the subdivision algorithm. – David Mar 31 '12 at 4:27
@Andy How about evaluating the function afresh using the sampling points (I.e. x values) from the plot to overcome machine precision? – JxB Mar 31 '12 at 14:32
@JxB I had commented with some code to do just that but it took up so much space I deleted it :) I don't know but it might be worthwhile to set WorkingPrecision on Plot before hand (even though the results will need to be coerced via SetPrecision) – Andy Ross Mar 31 '12 at 16:16
@AndyRoss an $\mathbb{R}^3 \to \mathbb{R}$ function can be plotted using Plot, provided you parameterize it. The OP wants this along a line in $\mathbb{R}^3$, so that's covered. Also, you can increase the WorkingPrecision of Plot to go beyond machine precision. – rcollyer Apr 1 '12 at 2:53

Whenever I want to do cheap adaptive sampling along a function, I used to do what David did back in old versions of Mathematica. Nowadays, I proceed like so:

pts = Cases[Normal[Plot[Sin[x], {x, 0, 2 π}, Mesh -> All]], Point[pt_] :> pt, ∞];


ListPlot[pts] should yield an image similar to the one in David's answer. A similar procedure can be done for parametrically-defined plane curves (via ParametricPlot[]) and space curves (via ParametricPlot3D[]).

-