# How to add an interpolating point to InterpolatingFunction?

Suppose I have an InterpolatingFunction object defined on $[\epsilon,1]$; let's call it f. What is the best way to construct a function g that coincides with f on $[\epsilon,1]$ and satisfies g[0]==a (g should be an InterpolatingFunction, too)? Small perturbations near $\epsilon$ are allowed.

Any generalizations are welcome as well.

-
I have a difficulty in choosing appropriate tags here. Is [functions] enough? –  faleichik Jan 28 '12 at 20:11
I am not sure if this is possible at all unless you use linear interpolation. Higher order interpolation uses more than one point to compute function values. –  Szabolcs Jan 28 '12 at 20:12
In fact nonlinear extensions are possible, but I've added a note that it is not very necessary to strictly reproduce f near epsilon. –  faleichik Jan 28 '12 at 20:17
What I meant when I said not sure if it's possible was that I'm not sure if Mathematica's InterpolatingFunction makes this technically possible. If you are willing to accept a different type of Mathematica object (not InterpolatingFunction), then it is definitely possible. –  Szabolcs Jan 28 '12 at 20:19
I wonder which of these solutions will be able to handle the result of Interpolation[{{{0.1}, 0.11, 0.73, 0.92}, {{0.3}, 0.91}, {{0.4}, 0.95}, {{0.6}, 0.49, 0.87}, {{0.9}, 0.95, 0.25}}]... –  Ｊ. Ｍ. Jan 28 '12 at 23:51

When I told you that this was not possible, I was wrong.

My understanding is that you have points $(x_1, y_1), (x_2, y_2), ..., (x_n, y_n)$ through which you construct an interpolating function $f$. Now you need to add another point $(x_0, y_0)$, and construct a new interpolating function $f^*$ for which it is true that $f^*(x) = f(x)$ for all $x \in [x_1, x_n]$. I thought it was not possible to keep the function value unchanged in the interval $[x_1, x_n]$ when using higher value interpolation, but this is not the case. See below:

Let's make an interpolation function from cosine values between 0.1 and 1.0:

ifun = Interpolation@Table[{x, Cos[x]}, {x, .1, 1, .1}]


It looks like this:

The trick is that when we add an extra point at $x=0$, we need to keep all derivatives unchanged in $x = 0.1$ up to the order of interpolation.

You can get the order of interpolation like this:

ifun["InterpolationOrder"]

(* ==> 3 *)


Let's get the derivative values in the first point:

derivs = Table[
Derivative[i][ifun][ifun["Grid"][[1, 1]]],
{i, 0, First@ifun["InterpolationOrder"] - 1}]

(* ==> {0.995004, -0.0995897, -1.00396} *)


And inject them back into the function, while adding a new value $f(0) = 2$:

ifun2 = Interpolation@Join[
{{{0}, 2},
{ifun["Grid"][[1]], Sequence @@ derivs}},
]


Notice that the function is unchanged for all values greater than 0.1:

Plot[{ifun2[x], ifun[x]}, {x, 0, 1}, PlotRange -> All]


If you are wondering where I got this special API to InterpolatingFunction where we do things like ifun["Grid"]: I simply looked into the DifferentialEquationsInterpolatingFunctionAnatomy package that the other answers used.

-

You will have to make assumptions for the derivatives of your function at zero, unless your interpolating order is 1 (which is what @Szabolcs said in the comments). There is a package caled DifferentialEquationsInterpolatingFunctionAnatomy, which should help (it does, but is not enough). The way to really reconstruct your data is outlined in my answer to this question. I will reproduce the main function from that answer here (see the answer for details):

Clear[reconstructInterpolatingFunction];
reconstructInterpolatingFunction[intf_InterpolatingFunction] :=
With[{data = intf[[4, 3]],
step = Subtract @@ Reverse[ Take[intf[[4, 2]], 2]],
order =
DeveloperFromPackedArray@
InterpolatingFunctionInterpolationOrder[intf],
grid = InterpolatingFunctionGrid[intf]
},
Interpolation[
InterpolationOrder -> order
]
];


You will have to modify it, by changing MapThread[Prepend, {Partition[data, step], grid}] to Prepend[MapThread[Prepend, {Partition[data, step], grid}],{a,f'[0],...}], where f'[0] is the value of first derivative of f at zero, etc.

-
yes, yours is more complete as it more accurately reconstructs the original interpolation function in all details. I'll delete my answer as this is better. –  acl Jan 28 '12 at 20:36
@acl Why don't you undelete your answer, since it is simpler, and the procedure is anyways rather ad hoc, unless values for derivatives at zero are known precisely. You could just make a remark on the limitations. –  Leonid Shifrin Jan 28 '12 at 20:39
Ahhh ... I just notice you posted the same thing I did a full 30 minutes before ... –  Szabolcs Jan 28 '12 at 20:58
@Szabolcs I had an advantage: I did it before, and merely reproduced part of my past answer :). You also made a nice visual analysis which I did not. The point that I made but you did not however is that InterpolatingFunction object does keep all the values for derivatives, so there is no need to compute them. The price to pay is that I rely on undocumented internal details of the format in which InterpolatingFunction keeps its data. This would be unnecessary if the designers of the InterpolatingFunctionAnatomy  API did not overlook that info for derivatives is not accessible. –  Leonid Shifrin Jan 28 '12 at 21:10
Have a look at intf["Methods"]; Furthermore, if the input to the function specifies g[a] and g'[a] you could get the derivatives of intf, by using D[intf[x],x] and then extract the data from that interpolation function append the g'[a] and re-interpolate. –  user21 Jan 28 '12 at 21:20