Mathematica Stack Exchange is a question and answer site for users of Mathematica. Join them; it only takes a minute:

Sign up
Here's how it works:
  1. Anybody can ask a question
  2. Anybody can answer
  3. The best answers are voted up and rise to the top

I want to find roots of a function that is very slow to calculate (the function itself involves a bunch of FindMaxima) so I interpolate the function in a region where I think it likely the the root is, since I have very little certainty where the root is I have to interpolate an absurdly large region.

Is there a way to have InterpolatingFunction remember the original function and expand its domain when evaluation is attempted outside it?

share|improve this question
If you are willing to try different interpolation methods then there may be something here of relevance. – Daniel Lichtblau Apr 4 '13 at 20:09
@DanielLichtblau That's some really neat ways to do interpolation, thanks! – ssch Apr 4 '13 at 20:23

Here is a way that sadly it has a little bit of overhead making it slower to evaluate. When argument is outside domain it interpolates towards that direction and creates a new InterpolatingFunction by merging together with the old one.

(*   f: function to interpolate
    x0: Starting position
  size: minimum interval to interpolate each time
 order: InterpolatingOrder 
AutoInterpolatingFunction[f_, x0_, size_, opts:OptionsPattern[FunctionInterpolation]] :=
  wrap["if"] = 
    f[\[FormalX]], {\[FormalX], x0 - size, x0 + size}, 
  wrap["f"] = f;
  wrap["xmin"] = x0 - size;
  wrap["xmax"] = x0 + size;
  (* Argument in domain *)
  wrap[x_?NumericQ] := wrap["if"][x] /; wrap["xmin"] <= x <= wrap["xmax"];
  (* Interpolate a bit further *)
  wrap[x_?NumericQ] :=
    If[x > wrap["xmax"],
     newf = FunctionInterpolation[
       {\[FormalX], wrap["xmax"], Max[wrap["xmax"] + size, x]},
     wrap["if"] = 
      Interpolation[{Join[wrap["if"]["Grid"], Rest@newf["Grid"]],
     wrap["xmax"] = wrap["if"]["Domain"][[1, 2]]
     newf = 
       wrap["f"][\[FormalX]], {\[FormalX], 
        Min[x, wrap["xmin"] - size], wrap["xmin"]},
     wrap["if"] = 
      Interpolation[{Join[newf["Grid"], Rest@wrap["if"]["Grid"]],
     wrap["xmin"] = wrap["if"]["Domain"][[1, 1]]

A simple speed test:

aip = AutoInterpolatingFunction[Sin[#] Cos[#] &, 0, 1];
ip = FunctionInterpolation[Sin[x] Cos[x], {x, 0, 1}];
AbsoluteTiming[Do[ip[0.5], {100000}]]
AbsoluteTiming[Do[aip[0.5], {100000}]]
(* 1.7s vs 0.7s *)
share|improve this answer

Your Answer


By posting your answer, you agree to the privacy policy and terms of service.

Not the answer you're looking for? Browse other questions tagged or ask your own question.