Interpolating data with a step

Introduction and background

While Interpolation complains, NDSolve is able to construct such interpolations:

yifn = NDSolveValue[{y'[t] == y[t], y[0] == 0, WhenEvent[t > 1, y[t] -> 1]},
   y, {t, 0, 2}];

But if we use Interpolation to reinterpolate the solution, we get a repeated-point error again:

Interpolation[Transpose[{yifn["Grid"], yifn["ValuesOnGrid"]}]] // Short

Interpolation::inddp: The point 1.` in dimension 1 is duplicated. >>

(*  Interpolation[{{{0.}, 0.},..., {{2.}, 2.71828}}]  *)

It's not mistaken -- as if it could be:

Min@Differences@yifn["Grid"]
(*  0.  *)

Plot[yifn[t], {t, 0, 2}]

Okay, so constructing an InterpolatingFunction with repeated abcissae is possible. How can it be done? Are there hidden options to Interpolation or is there some other way to do it?

Some remarks:

1. The NDSolve solution also contains derivative information that is lost in the reinterpolation, but including it does not solve the issue of repeated abscissae.

2. Albert Retey's answer to How to splice together several instances of InterpolatingFunction?, which uses Piecewise to switch between interpolating functions, may be the best solution in most cases, since it makes explicit the points of discontinuity or non-smoothness. But it is not what NDSolve does. Ian Hinck's modification to István Zachar's answer could be made more accurate if two InterpolatingFunctions could be joined without deleting the duplicate endpoints.

Answer

The key is setting a flag in a bitfield, which I will illustrate with an example. But that is far from the complete story, which I will explain further below. It would be nice if there were a way to do it with built-in methods.

I have a rather klunky way to construct an InterpolatingFunction directly. The information in What's inside InterpolatingFunction[{{1., 4.}}, <>]? is incomplete; further, it changes slightly with each new version of InterpolatingFunction.

If we examine the difference between two interpolating functions produced by NDSolve, one with a discontinuity and one without, we see that the bitfield yifn[[2,2]] is different. Setting the fourth low-order bit (2^4) in the bitfield (of 2nd element of second argument) permits repeated abscissae.

We need n+1 data values to construct an order n interpolation on each interval $[x_{i-1},x_{i}]$. For instance for an order 3 interpolation (by cubic polynomials), four data points are needed to determine the cubic polynomial. For an interior interval $[x_{i-1},x_{i}]$ of an ordinary interpolation of points, the points used are $$(x_{i-2},y_{i-2}),(x_{i-1},y_{i-1}),(x_{i},y_{i}),(x_{i+1},y_{i+1})$$ as shown in this answer adapted from the documentation . This cannot be done if $x_{i-2}=x_{i-1}$. Another way to construct the interpolating polynomial is to use both function and derivative values. For instance, with $$(x_{i-1},y_{i-1},y'_{i-1}),(x_{i},y_{i},y'_{i})$$ the four values $$y_{i-1},y'_{i-1},y_{i},y'_{i}$$ determine a cubic polynomial. With such data, it does not matter if $x_{i-2}=x_{i-1}$. This is the approach of NDSolve. When there are not enough derivative values at the endpoints of an interval, points and values from adjacent intervals are used. It is important, therefore, that there be enough derivative information at the duplicated points.

Example

Interpolation[{{{0.}, 1., 0.}, {{1.}, 0.8, 1.}, {{1.}, 2., 0.}, {{2.}, 1., -2.}}]

Interpolation::inddp: The point 1.` in dimension 1 is duplicated. >>

(*  Interpolation[{{{0.}, 1., 0.}, {{1.}, 0.8, 1.}, {{1.}, 2., 0.}, {{2.}, 1., -2.}}]  *)

murf = InterpolatingFunction @@ {
  {{0.`, 2.`}}, {5, 16 + 5, 1, {4}, {4}, 0, 0, 0, 0, Automatic, {}, {}, False},
   {{0.`, 1.`, 1.`, 2.`}},
   {Developer`PackedArrayForm,
    {0, 2, 4, 6, 8}, {1.`, 0.`, 0.8`, 1.`, 2.`, 0.`, 1.`, -2.`}},
   {Automatic}}
(*  InterpolatingFunction[{{0., 2.}}, <>] *)

murf[1]  (* evaluates left-side formula *)
Plot[murf[x], {x, 0, 2}]
Plot[murf[x], {x, 0, 2},
 Exclusions ->
  Thread[x == Pick[Most@First@murf["Coordinates"],
     Differences@First@murf["Coordinates"], 0.]]]

(*  0.8  *)

General-purpose function

Note that we call InterpolatingFunction directly. It seems to assume that the input has been checked (by Interpolation). Feeding it bad data can cause crashes. Caveat: I haven't checked all the things I didn't think of. I did put some constraints on the input to discontInterpolation that are more restrictive than necessary. For instance, we require the derivative values given at each point to be the to the same order. Not knowing how the internals work exactly, it was difficult to construct an accurate check of the requirements that would be minimally satisfactory to InterpolationFunction. And since Mathematica's way of letting me know I'd got it wrong was often to crash, I decided it would be best to mimic NDSolve output, which seems reliable.

ClearAll[discontInterpolation];
discontInterpolation::data =
  "Input data not of the form {{{x1}, y1,...},...}.";
discontInterpolation::darray =
  "The number of derivative values should be the same for each input.";
discontInterpolation::data =
  "Derivative order should be at least the interpolation order minus one.";
Options[discontInterpolation] =
  Join[Options[Interpolation], {"ExtrapolationHandler" -> Automatic}];

Begin["discontInterpolation`"];

(* parts of InterpolationFunction that we use *)
$domain = Sequence[1];              (* interpolation domain {{xmin, xmax}} *)
$xmin = Sequence[1, 1];
$xmax = Sequence[1, 2];
$bitfield = Sequence[2, 2];         (* flags indicating properties *)
$derivativeorder = Sequence[2, 3];  (* max derivative order *)
$ngrid = Sequence[2, 4];            (* number of points in each coordinate grid *)
$iorder = Sequence[2, 5];           (* interpolation order *)
$derivative = Sequence[2, 6];       (* derivative to evaluate: 0->f[x], 1->f'[x],... *)
$periodic = Sequence[2, 7];
$grid = Sequence[3];                (* abscissae of interpolation grid *)
$itype = Sequence[4, 1];            (* Developer`PackedArrayForm - this varies by type *)
$fnoffsets = Sequence[4, 2];        (* offsets in function/derivative array *)
$fnvalues = Sequence[4, 3];         (* interleaved function/derivative values *)

(* bit field positions - inferred, perhaps mistaken *)
$fullArrayBit = 1; (* whether the data (f, f',...) is a full array (not ragged) $*)
$repeatedBit = 4;  (* whether repeated abscissae are permitted $*)

checkArgs[data_, iorder_, dorder_] :=
 Block[{pattern},
  If[! MatchQ[data, {List[{_?NumericQ}, __?NumericQ] ..}],
   Message[discontInterpolation::data]; Return[False]];
  pattern = N@First[data] /. _Real -> Blank[];
  If[! MatchQ[data, {pattern ..}],
   Message[discontInterpolation::darray]; Return[False]];
  If[dorder < iorder - 2,
   Message[discontInterpolation::data]; Return[False]];
  True
  ];
discontInterpolation[data_List, opts : OptionsPattern[]] :=
  Block[{iorder, dorder},
   iorder = OptionValue[InterpolationOrder] /. Automatic -> 3;
   dorder = Length[First@data] - 2;
   With[{if0 = Interpolation[{0., 1.}, InterpolationOrder -> 1, opts]},
     ReplacePart[
       if0 /. InterpolatingFunction -> Inactive[InterpolatingFunction],
       {{$domain} -> {System`MinMax[data[[All, 1]]]},        (* Domain $*)
        {$bitfield} -> BitSet[if0[[$bitfield]], $repeatedBit],  (* Bitfield $*)
        {$derivativeorder} -> dorder,                        (* Max derivative order $*)
        {$ngrid} -> {Length[data]},                          (* Number of abscissae $*)
        {$iorder} -> {iorder + 1},                           (* Interpolation order $*)
        {$grid} -> {Flatten@ data[[All, 1]]},                (* Input grid coordinates $*)
        {$fnoffsets} ->             (* offsets of function values on input grid $*)
         Table[i, {i, 0, (dorder + 1) Length[data], dorder + 1}],
        {$fnvalues} -> Join @@ data[[All, 2 ;;]] (* function and derivative values $*)
        }] // Activate
     ] /; checkArgs[data, iorder, dorder]
   ];

End[];

OP's example

We have to supplement data with (at least first) derivative information. In the toy example, this easy.

data = {{1, 1}, {2, 2}, {3, 3}, {3, 4}, {4, 5}, {5, 6}};
ddata = data /. {x_?NumericQ, y_?NumericQ} :> {{x}, y, 1};
op = discontInterpolation[ddata, InterpolationOrder -> 1]
(*  InterpolatingFunction[{{1, 5}}, <>]  *)

Plot[op[x], {x, 1, 5}]

I suppose you can construct a Piecewise expression.

data = {{0, 1}, {1, 2}, {2, 1}, {3, 3}, {3, 4}, {4, 5}, {5, 6},
        {6, 5}, {7, 6}, {7, 4}, {8, 3}, {9, 3}, {10, 5}};

ifs = Interpolation[#, InterpolationOrder -> 1] & /@ 
   Split[data, First[#1 - #2] != 0 &];

piece[if_, t_: t] := {if[t], #1 < t < #2 & @@ First@if["Domain"]};
pw = Piecewise[piece /@ ifs, Indeterminate];

excl = Times @@ (t - Union@Flatten[#["Domain"] & /@ ifs][[2 ;; -2]]);

Plot[pw, {t, 0, 10}, PlotPoints -> {35, {3, 7}}, 
 Exclusions -> excl == 0]

I made a random choice about boundary conditions, that is, what the value should be at the jumps (Indeterminate), but this can be changed.

You could parametrize, e.g.

data = {{0, 1}, {1, 2}, {2, 1}, {3, 3}, {3, 4}, {4, 5}, {5, 6}, {6, 
    5}, {7, 6}, {7, 4}, {8, 3}, {9, 3}, {10, 5}};
{x, y} = Transpose[data];
xf = ListInterpolation[x, {0, 1}, InterpolationOrder -> 1];
yf = ListInterpolation[y, {0, 1}, InterpolationOrder -> 1];
ParametricPlot[{xf[t], yf[t]}, {t, 0, 1}, Frame -> True, 
 AxesOrigin -> {0, 0}, Epilog -> {Red, PointSize[0.02], Point[data]}]