【求助】如何让InterpolationFunction现原形?
mathematica吧
全部回复
仅看楼主
level 8
flumer 楼主
想拟合一个曲面啊,出来带入到NIntegrate真是好慢啊,想让它现个原形算快点,要是多项式拟合就好了,但是点很多1500多个,InterpolationPolynomial用不出。纠结
2013年08月01日 11点08分 1
吧务
level 15
……啥叫现原形啊?
2013年08月01日 12点08分 2
就是,它的函数表达式
2013年08月01日 12点08分
回复 flumer :……你对InterpolatingFunction的理解有点问题,这会儿没啥时间了,你先仔细看看帮助吧。
2013年08月01日 13点08分
level 8
flumer 楼主
这问题好纠结,自己都想删掉算了
2013年08月01日 12点08分 3
吧务
level 15
……奇怪,13年的9月我怎么忘了继续回这帖了。不过那时节确实也没多少可说的,因为:
1. InterpolatingFunction 里确实没有存储与插值函数对应的符号多项式,所以我们没法简单地从中把所需的多项式提取出来。
2. InterpolatingFunction 本身已经是优化过的数值函数,改用符号式的插值多项式大概也占不了多少便宜(编译一下或许有效果?)。在NIntegrate那边想想办法可能还实际一点,比如给NIntegrate添加Method -> {Automatic, SymbolicProcessing -> 0}选项,或是直接改用梯形规则或是高斯积分规则。(当然了,这些方法对于一些性质不太好的积分效果不佳,但是既然你的被积函数里含了插值函数,那你的被积函数性质应该是比较好的?)
最后,需要补充的是,InterpolatingFunction 在版本11(16年才发布嗯)引入了一个"GetPolynomial"方法,这个方法可以从插值方法为"Hermite"的
一元函数里面提取多项式,Carl Woll(WRI员工,Wolfram|Alpha首席开发者)还为此写了一个函数:
InterpolationToPiecewise[if_, x_] := Module[{main, default, grid}, grid = if["Grid"]; Piecewise[ {if @ "GetPolynomial"[
#, x-#
], x < First @ #}& /@ grid[[2 ;; -2]], if @ "GetPolynomial"[
#, x-#
]& @ grid[[-1]] ]] /; if["InterpolationMethod"] == "Hermite"
更多内容可参考《Extracting the function from InterpolatingFunction object》(SE编号 59944 )。
我们拿这个函数测试一下前面的猜测:
dat = Table[{x, Sin[x]}, {x, 0, 2 Pi, 2 Pi/25.}];
func = Interpolation[dat];
piece[x_] = InterpolationToPiecewise[func, x];
cf = Compile[x, #, RuntimeOptions -> EvaluateSymbolically -> False] &@
piece[x];
NIntegrate[x piece[x], {x, 0, 2 Pi}] // RepeatedTiming
(* {0.0655642, -6.28287} *)
NIntegrate[x func[x], {x, 0, 2 Pi}] // RepeatedTiming
(* {0.013448, -6.28287} *)
NIntegrate[x cf[x], {x, 0, 2 Pi}] // RepeatedTiming
(* {0.00544045, -6.28292} *)
NIntegrate[x func[x], {x, 0, 2 Pi},
Method -> {Automatic, SymbolicProcessing -> 0}] // RepeatedTiming
(* {0.000853251, -6.28292} *)
NIntegrate[x piece[x], {x, 0, 2 Pi},
Method -> {Automatic, SymbolicProcessing -> 0}] // RepeatedTiming
(* {0.000951117, -6.28292} *)
NIntegrate[x cf[x], {x, 0, 2 Pi},
Method -> {Automatic, SymbolicProcessing -> 0}] // RepeatedTiming
(* {0.000528332, -6.28292} *)
好嘛,未编译的分段多项式反而是最慢的。对符号分段多项式进行编译可以提速,但是改NIntegrate选项综合而言还是最优的(编程简单,而且提速显著)。当然了,一维的经验可能不能直接推广至二维,但我想上面这个例子还是有启发性的。
2025年03月29日 05点03分 5
吧务
level 15
此外,对于不规则的(unstructured)网格近年来也出现了一些相关讨论,不过这个话题就更难了(主要坑点是自变量的定位,或者说,分段函数的分段条件),一个相关的SE帖子是《2D Interpolation for symbolic function values on unstructured grid》(编号 311670 ),有兴趣的可以参考。
2025年03月29日 05点03分 6
1