How to splice together several instances of InterpolatingFunction?

Question

I have a set of InterpolatingFunction returned by NDSolve which are valid over different (but overall continuous) domains. How do I splice them together into one single InterpolatingFunction over all the domains?

Piecewise seems to promising, but I can't manage to return the piecewise function from another function then use it later the same way as InterpolatingFunction.

I guess there is also the brute force way of generating a grid of points using the original set of InterpolatingFunction then interpolating points again, but that's very elaborate and CPU-consuming, not to mention potentially inaccurate if the interpolation grid is not chosen properly.

Thoughts?

---

Thanks to the answer from Mr. Wizard, this is the solution I ended up using:

JoinInterpolatingFunction[intervals_List, flist_List] := 
 Module[{getGrid},
  getGrid[f_InterpolatingFunction, min_?NumericQ, 
    max_?NumericQ] := {{min, f[min]}}~
     Join~(Transpose@{f["Grid"] // Flatten, f["ValuesOnGrid"]} // 
       Select[#, (min < #[[1]] < max) &] & )~Join~{{max, f[max]}} // N;
  Interpolation[
   Table[getGrid[flist[[i]], intervals[[i]], intervals[[i + 1]]], {i, 
       Length@flist}] // Flatten[#, 1] & // 
    DeleteDuplicates[#, (#1[[1]] == #2[[1]]) &] &, 
   InterpolationOrder -> 2]]

JoinInterpolatingFunction[{I1,I2,..,In},{func1,func2,...func(n-1)}] gives an InterpolatingFunction that takes values of func1 between [I1,I2], func2 between (I2,I3], ... func(n-1) between (I(n-1),In].

Answer 1

I expect that if the InterpolationOrder is the same between functions it should be possible to merge them into one. If not Piecewise may be the best you can do.

This is an incomplete answer but hopefully a useful signpost that may lead you to a solution. You can get the constituent parts (or at least their related forms) using the little-known "Methods" syntax, which is akin to the "Properties" of SparseArray if you have seen that before.

Here is a list of the "Methods":

f1 = Interpolation @ Table[{i, Sin[i]}, {i, 0, Pi, 0.1}];
f1["Methods"]

{"Coordinates", "DerivativeOrder", "Domain", "Evaluate", "Grid", "InterpolationOrder",
  "MethodInformation", "Methods", "Properties", "ValuesOnGrid"}

Here are the internal usage messages:

f1["MethodInformation"@#] & ~Scan~ f1["Methods"]

InterpolatingFunction[domain, data]@Coordinates[] returns the grid coordinates in each dimension.

InterpolatingFunction[domain, data]@DerivativeOrder[] returns what derivative of the interpolated function will be computed upon evaluation.

InterpolatingFunction[domain, data]@Domain[] returns the domain of the InterpolatingFunction.

InterpolatingFunction[domain, data]@Evaluate[arg] evaluates the InterpolatingFunction at the argument arg.

InterpolatingFunction[domain, data]@Grid[] gives the grid of points where the interpolated data is defined.

InterpolatingFunction[domain, data]@InterpolationOrder[] returns the degree of polynomials used for computing interpolated values.

InterpolatingFunction[domain, data]@MethodInformation[method] gives information about a particular method.

InterpolatingFunction[domain, data]@Methods[pat] gives the list of methods matching the string pattern pat.

InterpolatingFunction[domain, data]@Propertiesgives the list of possible properties.

InterpolatingFunction[domain, data]@ValuesOnGrid[] gives the function values at each grid point. In some cases, this may be faster than evaluating at each of the grid points.

Here is the actual output when applying these "Methods" to the example InterpolatingFunction above:

f1 /@ {"Coordinates", "DerivativeOrder", "Domain", "Grid", 
  "InterpolationOrder", "Properties", "ValuesOnGrid"} // Column

{{0., 0.1, 0.2, 0.3, 0.4, 0.5, 0.6, 0.7, 0.8, 0.9, 1., 1.1, 1.2, 1.3, 1.4, 1.5,
  1.6, 1.7, 1.8, 1.9, 2., 2.1, 2.2, 2.3, 2.4, 2.5, 2.6, 2.7, 2.8, 2.9, 3., 3.1}}

0

{{0., 3.1}}

{{0.},{0.1},{0.2},{0.3},{0.4},{0.5},{0.6},{0.7},{0.8},{0.9},{1.},{1.1},{1.2},
 {1.3},{1.4},{1.5},{1.6},{1.7},{1.8},{1.9},{2.},{2.1},{2.2},{2.3},{2.4},{2.5},
 {2.6},{2.7},{2.8},{2.9},{3.},{3.1}}

{3}

{0., 0.0998334, 0.198669, 0.29552, 0.389418, 0.479426, 0.564642, 0.644218,
 0.717356, 0.783327, 0.841471, 0.891207, 0.932039, 0.963558, 0.98545, 0.997495,
 0.999574, 0.991665, 0.973848, 0.9463, 0.909297, 0.863209, 0.808496, 0.745705,
 0.675463, 0.598472, 0.515501, 0.42738, 0.334988, 0.239249, 0.14112, 0.0415807}

Answer 2

It could well be that one of the other suggestions will lead you to what you'll be using in the end. I think you should still know about the most straightforward way to create a combination of interpolating functions using Piecewise and a pure function:

ipf1 = Interpolation[Table[{x, Sin[x]}, {x, 0, 1, 0.1}]]
ipf2 = Interpolation[Table[{x, Sin[x]}, {x, 1, Pi, 0.1}]]
ipfCombined = Function[Piecewise[{{ipf1[#], # <= 1}, {ipf2[#], # > 1}}]]

the result can almost everywhere be used just like an InterpolatingFunction:

Plot[ipfCombined[x], {x, 0, Pi}]
Integrate[ipfCombined[x], {x, 0, Pi}]

(if you want to show a continuous plot you can add the option Exclusions -> None)

Answer 3

If you explicitly want to join the list of coordinates for e.g. two univariate InterpolatingFunctions (i.e. not opting for Piecewise), you can use the InterpolatingFunctionAnatomy package, that allows extraction of coordinates and grid points of ip.functions:

Needs["DifferentialEquations`InterpolatingFunctionAnatomy`"];

Options[joinInterpolation] = Options[Interpolation];
joinInterpolation[int : {__InterpolatingFunction}, opts___] := Module[{coord, vals},
   coord = Join @@ (First@InterpolatingFunctionCoordinates@# & /@ int);
   vals = Join @@ (InterpolatingFunctionValuesOnGrid /@ int);
   Interpolation[MapThread[List, {coord, vals}], opts]
   ];

Test the function:

ipf1 = Interpolation[Table[{x, Sin[x]}, {x, 0, Pi, 0.1}]];
ipf2 = Interpolation[Table[{x, Sin[x]}, {x, Pi, 2 Pi, 0.1}]];
joined = joinInterpolation[{ipf1, ipf2}];

Row@{Plot[ipf1[t], {t, 0, Pi}],
     Plot[ipf2[t], {t, Pi, 2 Pi}],
     Plot[joined[t], {t, 0, 2 Pi}]
    }