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.
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
*)
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}},<>]*)