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


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

JoinInterpolatingFunction[intervals_List, flist_List] := 
  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;
   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].

share|improve this question

4 Answers 4

up vote 10 down vote accepted

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}];
{"Coordinates", "DerivativeOrder", "Domain", "ElementMesh", "Evaluate", "Grid",
"InterpolationMethod", "InterpolationOrder", "MethodInformation", "Methods",
"OutputDimensions", "Periodicity", "PlottableQ", "Properties", "QuantityUnits",

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"};













share|improve this answer
I've seen f1["MethodInformation"@#] & ~Scan~ f1["Methods"] before. The Call of Cthulhu, chapter 3. –  belisarius Feb 5 '13 at 5:32
@belisarius you can add it to the list of my scary coding. –  Mr.Wizard Feb 5 '13 at 5:36
Damn! I thought I was the craziest around. –  belisarius Feb 5 '13 at 5:39
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... –  Oleksandr R. Feb 8 '13 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) –  Mr.Wizard Feb 8 '13 at 4:58

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)

share|improve this answer
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 '13 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... –  Albert Retey Feb 5 '13 at 18:57
@ruebenko: Is that typo in your comment to be interpreted as a Freudian slip? –  Albert Retey Mar 13 '13 at 13:27
Oh, not, it's to be interpolated ;-) –  user21 Mar 13 '13 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:


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

Mathematica graphics

share|improve this answer
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. –  Oleksandr R. Feb 8 '13 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... :) –  István Zachar Feb 8 '13 at 8:18

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