wolfram-mathematica - 如何从InterpolatingFunction提取所有数据以创建相同的数据？

在对previous question的评论中，建议从Mathematica 5.2生成的InterpolatingFunction中提取所有数据，然后在Mathematica 8中创建另一个数据。建议的方法是使用DifferentialEquations`InterpolatingFunctionAnatomy`包中定义的函数从InterpolatingFunction中提取数据。天真地尝试

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}]

我在原始函数和重构函数之间有了明显的区别:

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

查看这些函数的InputForm可以发现它们并不相同。是否可以通过从中提取所有数据并在提取的数据上调用InterpolatingFunction来重建Interpolation？

最佳答案

编辑

这是一个通用的解决方案，使用代码:

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
     ]
   ];

请参阅下面的说明。但是请注意，上述代码确实依赖于InterpolatingFunction对象的某些细节，这些细节可能是特定于版本的，因为当函数导数的值很重要时，显然DifferentialEquations`InterpolatingFunctionAnatomy`的API似乎无法完全重构原始对象。

结束编辑

似乎在构造NDSolve时InterpolatingFunction包含有关派生的信息，这很有意义。在您的情况下，这将包括一阶导数，因为您的方程是一阶的。但是，当我们使用DifferentialEquations`InterpolatingFunctionAnatomy`包中的函数时，这些信息会丢失。获得它的方法是直接访问初始InterpolatingFunction对象。这是一个简单的示例:

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}

这表明，在此网格点处，每个值后面都跟随着导数的值。因此，构造新对象的方法如下所示:

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

这使用了Interpolation的扩展形式，其中还可以指定导数的值。这通过了我们的测试:

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.}

确定给定InterpolatingFunction中的信息最多取决于哪一个导数的方法是查看以下部分:

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}

在这种情况下，步骤为2，因此我们有一个值加上一阶导数。例如，对于一个二阶方程，步长将为3，并且您将需要Partition[...,3]。因此，您可以通过执行本部分中的步骤来确定顺序。

现在，真正的事情是:

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.}