How to extract all data from InterpolatingFunction to create identical?

Question

In the comments on the previous question, it was proposed to extract all the data from the InterpolatingFunction generated by Mathematica 5.2, and then create another one in Mathematica 8. The proposed approach is to use the functions defined in the DifferentialEquations`InterpolatingFunctionAnatomy` package to extract data from the InterpolatingFunction . Trying naively

 Needs["DifferentialEquations`InterpolatingFunctionAnatomy`"] ifun = First[ x /. NDSolve[{x'[t] == Exp[x[t]] - x[t], x[0] == 1}, x, {t, 0, 10}]]; data = Transpose@ {InterpolatingFunctionGrid[ifun], InterpolatingFunctionValuesOnGrid[ifun]}; interpolationOrder = Developer` FromPackedArray@ InterpolatingFunctionInterpolationOrder[ifun]; ifun2 = Interpolation[data, InterpolationOrder -> interpolationOrder]; Table[ifun[x] - ifun2[x], {x, 0, 0.5160191740198963`, .1}] 

I get a significant difference between the original function and the reconstructed:

 {0., 2.13061*10^-7, 2.05087*10^-7, 2.46198*10^-7, 6.68937*10^-7, 1.5624*10^-7} 

A look at the InputForm these functions shows that they are not identical. Is it possible to restore an InterpolatingFunction by extracting all the data from it and calling Interpolation on the extracted data?

Answer 1

EDIT

Here is a general solution, in code:

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

See below for more details. However, note that the above code uses some details of the InterpolatingFunction object that may be version specific, since the DifferentialEquations`InterpolatingFunctionAnatomy` API does not seem to allow full restoration of the original object when values ​​for derived functions are important.

End edit

---

It looks like NDSolve contains derivative information when building an InterpolatingFunction , which makes sense. In your case, this will be the first derivative, since your equation is first order. But this information is lost when we use the functions from the DifferentialEquations` InterpolatingFunctionAnatomy` package. The way to get this is to access the original InterpolatingFunction object directly. Here is a simple example:

 In[156]:= ifun=First[x/.NDSolve[{x'[t]==2x[t],x[0]==1},x,{t,0,0.1}]]; In[157]:= ifun[[4,3]] Out[157]= {1.,2.,1.00002,2.00004,1.00004,2.00008,1.00457,2.00913,1.00911,2.01823, 1.01368,2.02736,1.02787,2.05573,1.04225,2.0845,1.05684,2.11368,1.07163,2.14326, 1.09328,2.18655,1.11536,2.23073,1.13789,2.27579,1.16088,2.32176,1.18433,2.36867, 1.20272,2.40545,1.2214,2.44281} 

This suggests that each value is followed by the value of the derivative at this grid point. So the way to create a new object is something like this:

 ifun5 = Interpolation[ MapThread[Prepend, {Partition[ifun[[4, 3]], 2], InterpolatingFunctionGrid[ifun]}], InterpolationOrder -> (interpolationOrder)] 

It uses an extended form of Interpolation , where you can also specify values ​​for derivatives. This passes our test:

 In[161]:= Table[(1-ifun5[x]/ifun[x]),{x,0,0.1,.01}] Out[161]= {0.,0.,0.,0.,0.,0.,0.,0.,0.,0.,0.} 

The way to determine to which derivative we have information in this InterpolatingFunction is to look at this part:

 In[176]:= ifun[[4,2]] Out[176]= {0,2,4,6,8,10,12,14,16,18,20,22,24,26,28,30,32,34} 

In this case, the step is 2, so we have the value plus the first derivative. For, say, a second-order equation, the step will be 3, and you will need Partition[...,3] . So, you determine the order by going through the step in this part.

Now, the real thing:

 In[162]:= ifun=First[x/.NDSolve[{x'[t]==Exp[x[t]]-x[t],x[0]==1},x,{t,0,10}]]; interpolationOrder=Developer` FromPackedArray@InterpolatingFunctionInterpolationOrder [ifun]; ifunnew = Interpolation[MapThread[Prepend, {Partition[ifun[[4,3]],2],InterpolatingFunctionGrid[ifun]}], InterpolationOrder->(interpolationOrder)]; Table[(1-ifunnew[x]/ifun[x]),{x,0,0.5,.1}] During evaluation of In[162]:= NDSolve::ndsz: At t == 0.5160191740198969`, step size is effectively zero; singularity or stiff system suspected. >> Out[165]= {0.,0.,0.,1.11022*10^-16,0.,0.}