4

When trying to recreate an InterpolationFunction produced by NDSolve I faced very strange problem with InterpolationOrder option of Interpolation. Consider the following InterpolationFunction (an example function from the Documentation):

ifun = First[
  x /. NDSolve[{x'[t] == Exp[x[t]] - x[t], x[0] == 1}, x, {t, 0, 10}]]

Now let us to try to reconstruct it. Here is the data:

Needs["DifferentialEquations`InterpolatingFunctionAnatomy`"]
data = Transpose@{InterpolatingFunctionGrid[ifun], 
   InterpolatingFunctionValuesOnGrid[ifun]};

And here is InterpolationOrder:

interpolationOrder = InterpolatingFunctionInterpolationOrder[ifun]
(*=> {3}*)

Now we try to construct the InterpolatingFunction:

Interpolation[data, InterpolationOrder -> interpolationOrder];

and get error Message:

Interpolation::inord: Value of option InterpolationOrder -> {3} should be a non-negative machine-sized integer or a list of integers with length equal to the number of dimensions, 1. >>

But if we specify InterpolationOrder by hands, it is OK:

Interpolation[data, InterpolationOrder -> {3}]
(*=> InterpolatingFunction[{{0.,0.516019}},<>]*)

Can anyone explain why InterpolationOrder -> interpolationOrder does not work while InterpolationOrder -> {3} works although interpolationOrder must be replaced with {3} BEFORE calling Interpolation according to the standard evaluation sequence?

P.S. The problem occurs in Mathematica 7.0.1 and 8.0.1 but not in Mathematica 5.2.

UPDATE

I have found one workaround for this bug:

Interpolation[data, 
 ToExpression@ToString[InterpolationOrder -> interpolationOrder]]

works as expected.

It seems that expressions generated by evaluation of Rule[InterpolationOrder,interpolationOrder] and Rule[InterpolationOrder,{3}] has different internal structure in spite of the fact that they are identical:

ByteCount // Attributes
ByteCount[InterpolationOrder -> interpolationOrder]
ByteCount[InterpolationOrder -> {3}]
Order[InterpolationOrder -> interpolationOrder, 
 InterpolationOrder -> {3}]

(*=> 
{Protected}
192
112    
0
*)
Alexey Popkov
  • 9,355
  • 4
  • 42
  • 93

2 Answers2

5

It seems that I have found the reason for this behavior. It is because InterpolatingFunctionInterpolationOrder function returns PackedArray:

Developer`PackedArrayQ@InterpolatingFunctionInterpolationOrder[ifun]
(*=> True*)

We can convert {3} into PackedArray ourselves:

Interpolation[data, 
  InterpolationOrder -> Developer`ToPackedArray@{3}];

(*=> gives the error Message*)

So the reason is that Interpolate does not support PackedArray as a value for InterpolationOrder option. The workaround is to unpack it manually:

Interpolation[data, 
 InterpolationOrder -> Developer`FromPackedArray@interpolationOrder]
(*=> InterpolatingFunction[{{0.,0.516019}},<>]*)
Alexey Popkov
  • 9,355
  • 4
  • 42
  • 93
2

Very strange behaviour indeed. Something like

a = {3};
Interpolation[data, InterpolationOrder -> a]

works fine, and both ??interpolationOrder and OwnValues[interpolationOrder] seem to indicate that interpolationOrder is just equal to {3}. Even weirder is that this does seem to work

interpolationOrder = 2 InterpolatingFunctionInterpolationOrder[ifun]/2
Interpolation[data, InterpolationOrder -> interpolationOrder]
Heike
  • 24,102
  • 2
  • 31
  • 45