Suppose I have some data with a step in it:
data = {{1, 1}, {2, 2}, {3, 3}, {3, 4}, {4, 5}, {5, 6}};
Interpolation will complain about this and not give you an interpolation function:
Interpolation[data, InterpolationOrder -> 1]
Interpolation::inddp: The point 3 in dimension 1 is duplicated.
Is there any way to easily avoid this? Can I get Mathematica to give me Piecewise instead, or jiggle my points?
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.
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?
The NDSolve solution also contains derivative information that is lost in the reinterpolation, but including it does not solve the issue of repeated abscissae.
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.
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 *)
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 specified 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 InterpolatingFunction. 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 InterpolatingFunction 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 *)
$extrapBit = 0; (* whether to warn about extrapolation *)
$fullArrayBit = 1; (* whether the data (f, f',...) is a full array (not ragged) $*)
$machPrecBit = 2; (* whether the data (f,f',...) is MachinePrecision *)
$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[];
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}]
InterpolatingFunction), it should also be fully manipulatable and/or constructible (or at least documented).
$\endgroup$
Oct 2, 2020 at 21:46
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]}]
Your data do not pass the 'vertical line' test and thus do not comprise a function. This can be seen by the fact that an input of 3 has two different outputs.
If you replace {3,3} by {2.9, 3} it will work fine. You may wish to use a different workaround but you cannot use the original data without modification.
data = {{1, 1}, {2, 2}, {2.9, 3}, {3, 4}, {4, 5}, {5, 6}};
f = Interpolation[data, InterpolationOrder -> 1]
Plot[f[x], {x, 1, 5}]
In this specific case you might presume your (corrected data) function is a piecewise linear function, in which case ListPlot would be appropriate:
ListPlot[data, Joined -> True]
(* same output as above *)
But normally it's better to use Interpolation.
If you are certain that your data are a list of ordered pairs of numbers, you may use the following to test whether they comprise a function:
functionQ[dat_] := Length[dat] == Length@GatherBy[dat, First]
functionQ[{{1, 1}, {2, 2}, {3, 3}, {3, 4}, {4, 5}, {5, 6}}]
functionQ[{{1, 1}, {2, 2}, {2.9, 3}, {3, 4}, {4, 5}, {5, 6}}]
False
True
You might also check whether there are any violations of the vertical line test.
verticalLineFails[dat_] := Cases[GatherBy[dat, First], x_ /; Length[x] > 1]
verticalLineFails[{{1, 1}, {2, 2}, {3, 3}, {3, 4}, {4, 5}, {5, 6}}]
{{{3, 3}, {3, 4}}}
