2
$\begingroup$

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?

$\endgroup$
  • $\begingroup$ If you are willing to try different interpolation methods then there may be something here of relevance. $\endgroup$ – Daniel Lichtblau Apr 4 '13 at 20:09
  • $\begingroup$ @DanielLichtblau That's some really neat ways to do interpolation, thanks! $\endgroup$ – ssch Apr 4 '13 at 20:23
1
$\begingroup$

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 *)
$\endgroup$

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service, privacy policy and cookie policy

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