# How to splice together several instances of InterpolatingFunction?

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].

Update: information below updated with values from version 10.0.0

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", "ElementMesh", "Evaluate", "Grid",
"InterpolationMethod", "InterpolationOrder", "MethodInformation", "Methods",
"OutputDimensions", "Periodicity", "PlottableQ", "Properties", "QuantityUnits",
"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]@ElementMesh[] returns the element mesh if one is present.

• 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]@InterpolationMethod[] returns the method used for interpolation.

• 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]@OutputDimensions[] returns the output dimensions of the interpolating function.

• InterpolatingFunction[domain, data]@Periodicity[] returns whether the interpolating function is periodic in the respective dimensions.

• InterpolatingFunction[domain, data]@PlottableQ[] returns whether the interpolating function is plottable or not.

• InterpolatingFunction[domain, data]@Properties gives the list of possible properties.

• InterpolatingFunction[domain, data]@QuantityUnits[] returns the quantity units associated with abscissa and ordinates.

• 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:

Print /@ f1 /@ {"Coordinates", "DerivativeOrder", "Domain", "ElementMesh", Evaluate[], "Grid",
"InterpolationMethod", "InterpolationOrder", "Methods", "OutputDimensions",
"Periodicity", "PlottableQ", "Properties", "QuantityUnits", "ValuesOnGrid"};

{{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}}

None

{{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}}

Hermite

{3}

{Coordinates,DerivativeOrder,Domain,ElementMesh,Evaluate,Grid,InterpolationMethod,InterpolationOrder,MethodInformation,Methods,OutputDimensions,Periodicity,PlottableQ,Properties,QuantityUnits,ValuesOnGrid}

{}

{False}

True

{Properties}

{None,None}

{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}

• I've seen f1["MethodInformation"@#] & ~Scan~ f1["Methods"] before. The Call of Cthulhu, chapter 3. Feb 5, 2013 at 5:32
• @belisarius you can add it to the list of my scary coding. Feb 5, 2013 at 5:36
• In version 9 there is additionally the "ElementMesh" method, no doubt courtesy of @ruebenko. It's funny... if this were any other language then your comment about this being little-known would probably be taken to mean that it was hidden away in an obscure corner of the documentation, rather than being completely undocumented and known about only through the grapevine... Feb 8, 2013 at 0:35
• @Oleksandr Sometimes it's good to be "Mr. Wizard" rather than "Mr. Architect" -- using Mathematica can be like programming with incantations. :o) Feb 8, 2013 at 4:58
• After some trial and error I guess I got to understand a bit of this really scary coding... f1["MethodInformation"["Domain"]] returns some information and is the basic structure that is repeated by that command. Also, both Mathematica's documentation and this site have no mention on how to use "MethodInformation" that I could find.
– ivbc
Jun 23, 2016 at 0:04

## Chaining extrapolation handlers

We can chain together the extrapolation handlers. It will overwrite any existing extrapolation handler except in the last interpolating function; however, that seems consistent with the goal of splicing together interpolating functions.

We can find the position of the extrapolation handler this way (see also What's inside InterpolatingFunction[{{1., 4.}}, <>]? for more on the structure of InterpolatingFunction):

Block[{f = Unique["ExtrapolationHandler"]},
First@Position[Interpolation[Range[4], "ExtrapolationHandler" -> {f}], f]]
(*  {2, 10}  *)


Then we can fold together the interpolating function thus:

With[{extrapHandlerPos = Block[{f = Unique["ExtrapolationHandler"]},
First@Position[Interpolation[Range[4], "ExtrapolationHandler" -> {f}], f]]},
interpolationJoin[ifns__] :=
Fold[ReplacePart[#2, extrapHandlerPos -> #1] &, Reverse@Flatten[{ifns}]]]


Test case:

ifns = Table[
Interpolation@Table[{x + 2 Pi i, Sin[i^2 x]}, {x, -2 Pi, 0, 2 Pi/60}], {i, 3}];
if = interpolationJoin[ifns];

Plot[if[x], {x, 0, 6 Pi}]


## V10: Joining multivariable interpolating functions

Using Piecewise with a Cuboid region for each domain, we can piece together functions of any number of variables (as long as the number of variables is the same).

Clear[interpolationJoin];
interpolationJoin[ifns__, vars_] /;
Apply[Equal, Length[#["Domain"]] & /@ Flatten[{ifns}]] :=
Piecewise @@
{{# @@ Flatten[{vars}], Flatten[{vars}] ∈ Cuboid @@ Transpose[#["Domain"]]} & /@
Flatten[{ifns}]}


Test case:

ifns2d = {NDSolveValue[
{Laplacian[u[x, y], {x, y}] == 1, DirichletCondition[u[x, y] == 0, True]},
u, {x, y} ∈ Disk[],
"ExtrapolationHandler" -> {Indeterminate &, "WarningMessage" -> False}],
NDSolveValue[
{Laplacian[u[x, y], {x, y}] == 1, DirichletCondition[u[x, y] == 0, True]},
u, {x, y} ∈ Cuboid[{1, -1}, {2, 1}],
"ExtrapolationHandler" -> {Indeterminate &, "WarningMessage" -> False}]};

if2 = interpolationJoin[ifns2d, {x, y}];

Plot3D[if2, {x, -1, 2}, {y, -1, 1}]


Note: For a more sophisticated approach, one could test for an ElementMesh domain in each interpolating function and use that instead of a cuboid, when an ElementMesh is present.

• Just a note: the resulting IF will have the domain of the first input IF and not the combined domain of all IFs. Sep 30, 2020 at 14:09
• @IstvánZachar By domain, I take it you mean what is returned by IF["Domain"], not the input values for which you get a well-defined (numeric) return value. Sep 30, 2020 at 16:16
• Yes, exactly. Your extrapolation-chaining method is elegant and fast, and it results in a well-formed IF. However, while Plot[if[x], {x, 0, combinedX}] nicely plots over the combined domain, ListLinePlot[if] does not and it may confuse unassuming users. In my use case I needed self-contained IFs that hold the domain info in a (uniformly) queriable fashion. Do you think that your method can be extended to modify in-place the domain and possibly other properties of the IF? (I have tried and failed.) Oct 1, 2020 at 8:23
• @IstvánZachar I think this method is a more promising one to adapt. The IFs being joined would have to have identical methods. There would be a problem if the domains intersected in an interval or were disjoint (instead of just a point). If there was a gap, then one might be able to interpolate Undefined or Indeterminate, but would have to convert the IFs from PackedArrayForm if they came from NDSolve. Since the linked method tends to crash the kernel on bad input, I didn't pursue it. But I think it's possible. Oct 1, 2020 at 13:45

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)

• I think this is somewhat preferable, for exmaple this also works when the underlying grid (e.g. in 2D) is the same. It will even work when the underlying grid is unstructured.
– user21
Feb 5, 2013 at 18:41
• I think this is the most general approach and should work almost always. Whenever it works there is no reason to do something more complicated. There might be cases where one wants to do something more complicated -- e.g. for performance improvement. I think those cases are rare, though... Feb 5, 2013 at 18:57
• @ruebenko: Is that typo in your comment to be interpreted as a Freudian slip? Mar 13, 2013 at 13:27
• Oh, not, it's to be interpolated ;-)
– user21
Mar 13, 2013 at 20:25

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["DifferentialEquationsInterpolatingFunctionAnatomy"];

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


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}]
}


• Just a small note: the information one can get with the functions in this package is the same as that available by using the undocumented properties listed by @Mr.Wizard above (this is apparent on looking at the package source code). While being documented is certainly an advantage, IMO these package functions are quite impractical due to their extremely verbose names. Feb 8, 2013 at 0:24
• @OleksandrR. I wasn't exactly sure that it does the same thing, thanks for the confirmation! Agreed on the long names, but of course we all know that Mathematica is notorious on using impractically long names... :) Feb 8, 2013 at 8:18

You could use NDSolve to piece together the interpolating functions. For example:

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

if = NDSolveValue[
{y'[x] == Piecewise[{{ipf1'[x], x<1}}, ipf2'[x]], y[0] == ipf1[0]},
y,
{x, 0, 3.1}
];
if //OutputForm


InterpolatingFunction[{{0., 3.1}}, <>]

Visualization:

Plot[if[x], {x, 0, 3.1}]


Nothing much to add here...this is just a slight modification to István's answer where the call to Needs has been removed in favour of Mr. Wizard's undocumented trick. The calling syntax is slightly different, and I also call DeleteDuplitacesBy to avoid the Interpolation warning about multiple points at a single coordinate.

Options[JoinInterpolatingFunction] = Options[Interpolation];
JoinInterpolatingFunction[int__InterpolatingFunction, opts:OptionsPattern[]] := Module[{data},
data = {Join @@ (First[#["Coordinates"]] & /@ {int}), Join @@ (#["ValuesOnGrid"] & /@ {int})}\[Transpose];
data = DeleteDuplicatesBy[data,First];
Interpolation[data, opts]
];


Test the function:

ipf1 = Interpolation[Table[{x, Sin[x]}, {x, 0, Pi, Pi/16}]];
ipf2 = Interpolation[Table[{x, Sin[x]}, {x, Pi, 2 Pi, Pi/16}]];
joined = JoinInterpolatingFunction[ipf1, ipf2];

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