Take the 2-minute tour ×
Mathematica Stack Exchange is a question and answer site for users of Mathematica. It's 100% free, no registration required.

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
add comment

1 Answer

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]] :=
 Module[{
   wrap,
   if},
  wrap["if"] = 
   FunctionInterpolation[
    f[\[FormalX]], {\[FormalX], x0 - size, x0 + size}, 
    opts];
  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] :=
   Module[{newf},
    If[x > wrap["xmax"],
     newf = FunctionInterpolation[
       wrap["f"][\[FormalX]],
       {\[FormalX], wrap["xmax"], Max[wrap["xmax"] + size, x]},
       opts];
     wrap["if"] = 
      Interpolation[{Join[wrap["if"]["Grid"], Rest@newf["Grid"]],
         Join[wrap["if"]["ValuesOnGrid"], 
          Rest@newf["ValuesOnGrid"]]}\[Transpose]];
     wrap["xmax"] = wrap["if"]["Domain"][[1, 2]]
     ,
     newf = 
      FunctionInterpolation[
       wrap["f"][\[FormalX]], {\[FormalX], 
        Min[x, wrap["xmin"] - size], wrap["xmin"]},
       opts];
     wrap["if"] = 
      Interpolation[{Join[newf["Grid"], Rest@wrap["if"]["Grid"]],
         Join[newf["ValuesOnGrid"], 
          Rest@wrap["if"]["ValuesOnGrid"]]}\[Transpose]];
     wrap["xmin"] = wrap["if"]["Domain"][[1, 1]]
     ];
    wrap["if"][x]];
  wrap
  ]

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
add comment

Your Answer

 
discard

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.