Mathematica 如何创建 InterpolatingFunction 对象?
Mathematica 如何创建 InterpolatingFunction 对象?示例:
test1 = FunctionInterpolation[Sin[x],{x,0,2*Pi}]
test1 的 FullForm 很长,但主要是 x 值, 相应的 y 值。然而,插值不是线性的 (因为我没有设置 InterpolationOrder -> 1)
我知道 Mathematica 使用三次样条,部分原因是默认 InterpolationOrder 为 3,还因为:
Plot[D[test1[t],t,t,t,t] /. t->x, {x,0,2*Pi}]
表明四阶导数一致为 0。
那么,Mathematica 如何计算这个三次样条呢?
我的目标是在 Perl 中使用 FunctionInterpolation 对象。
编辑:谢谢你,萨沙!这正是我想要的,有一点点 毛刺。下面是我尝试重新实现 Hermite 插值 很容易转换为 Perl 的方式(也可以在 https://github.com/barrycarter /bcapps/blob/master/bc-approx-sun-ra-dec.m#L234)。
问题:最后 3 个图的值虽小但非零。我不能 告诉我是否错误地实现了 Hermite,或者这只是一个数字 毛刺。
(* the Hermite <h>(not Hermione)</h> polynomials *)
h00[t_] = (1+2*t)*(1-t)^2
h10[t_] = t*(1-t)^2
h01[t_] = t^2*(3-2*t)
h11[t_] = t^2*(t-1)
(*
This confirms my understanding of InterpolatingFunction by calculating
the value in a different, Perl-friendly, way; this probably does NOT
work for all InterpolatingFunction's, just the ones I'm using here.
f = interpolating function, t = value to evaluate at
*)
altintfuncalc[f_, t_] := Module[
{xvals, yvals, xint, tisin, tpos, m0, m1, p0, p1},
(* figure out x values *)
xvals = Flatten[f[[3]]];
(* and corresponding y values *)
yvals = Flatten[f[[4,3]]];
(* and size of each x interval; there are many other ways to do this *)
(* <h>almost all of which are better than this?</h> *)
xint = (xvals[[-1]]-xvals[[1]])/(Length[xvals]-1);
(* for efficiency, all vars above this point should be cached *)
(* which interval is t in?; interval i = x[[i]],x[[i+1]] *)
tisin = Min[Max[Ceiling[(t-xvals[[1]])/xint],1],Length[xvals]-1];
(* and the y values for this interval, using Hermite convention *)
p0 = yvals[[tisin]];
p1 = yvals[[tisin+1]];
(* what is t's position in this interval? *)
tpos = (t-xvals[[tisin]])/xint;
(* what are the slopes for the intervals immediately before/after this one? *)
(* we are assuming interval length of 1, so we do NOT divide by int *)
m0 = p0-yvals[[tisin-1]];
m1 = yvals[[tisin+2]]-p1;
(* return the Hermite approximation *)
(* <h>Whoever wrote the wp article was thinking of w00t</h> *)
h00[tpos]*p0 + h10[tpos]*m0 + h01[tpos]*p1 + h11[tpos]*m1
]
(* test cases *)
f1 = FunctionInterpolation[Sin[x],{x,0,2*Pi}]
f2 = FunctionInterpolation[x^2,{x,0,10}]
f3 = FunctionInterpolation[Exp[x],{x,0,10}]
Plot[{altintfuncalc[f1,t] - f1[t]},{t,0,2*Pi}]
Plot[{altintfuncalc[f2,t] - f2[t]},{t,0,10}]
Plot[{altintfuncalc[f3,t] - f3[t]},{t,0,10}]
How does Mathematica create an InterpolatingFunction object? Example:
test1 = FunctionInterpolation[Sin[x],{x,0,2*Pi}]
The FullForm of test1 is long, but is primarily x values with the
corresponding y values. However, the interpolation is not linear
(since I didn't set InterpolationOrder -> 1)
I know Mathematica uses cubic splines, in part because the default
InterpolationOrder is 3, but also because:
Plot[D[test1[t],t,t,t,t] /. t->x, {x,0,2*Pi}]
shows the 4th derivative is uniformly 0.
So, how does Mathematica compute this cubic spline?
My goal is to use a FunctionInterpolation object in Perl.
EDIT: Thank you, Sasha! That did exactly what I wanted, with a minor
glitch. Below is my attempt to reimplement Hermite interpolation in a
way that's easy to convert to Perl (also available at
https://github.com/barrycarter/bcapps/blob/master/bc-approx-sun-ra-dec.m#L234).
The problem: the last 3 plots have small, but nonzero values. I can't
tell if I implemented Hermite wrong, or this is just a numerical
glitch.
(* the Hermite <h>(not Hermione)</h> polynomials *)
h00[t_] = (1+2*t)*(1-t)^2
h10[t_] = t*(1-t)^2
h01[t_] = t^2*(3-2*t)
h11[t_] = t^2*(t-1)
(*
This confirms my understanding of InterpolatingFunction by calculating
the value in a different, Perl-friendly, way; this probably does NOT
work for all InterpolatingFunction's, just the ones I'm using here.
f = interpolating function, t = value to evaluate at
*)
altintfuncalc[f_, t_] := Module[
{xvals, yvals, xint, tisin, tpos, m0, m1, p0, p1},
(* figure out x values *)
xvals = Flatten[f[[3]]];
(* and corresponding y values *)
yvals = Flatten[f[[4,3]]];
(* and size of each x interval; there are many other ways to do this *)
(* <h>almost all of which are better than this?</h> *)
xint = (xvals[[-1]]-xvals[[1]])/(Length[xvals]-1);
(* for efficiency, all vars above this point should be cached *)
(* which interval is t in?; interval i = x[[i]],x[[i+1]] *)
tisin = Min[Max[Ceiling[(t-xvals[[1]])/xint],1],Length[xvals]-1];
(* and the y values for this interval, using Hermite convention *)
p0 = yvals[[tisin]];
p1 = yvals[[tisin+1]];
(* what is t's position in this interval? *)
tpos = (t-xvals[[tisin]])/xint;
(* what are the slopes for the intervals immediately before/after this one? *)
(* we are assuming interval length of 1, so we do NOT divide by int *)
m0 = p0-yvals[[tisin-1]];
m1 = yvals[[tisin+2]]-p1;
(* return the Hermite approximation *)
(* <h>Whoever wrote the wp article was thinking of w00t</h> *)
h00[tpos]*p0 + h10[tpos]*m0 + h01[tpos]*p1 + h11[tpos]*m1
]
(* test cases *)
f1 = FunctionInterpolation[Sin[x],{x,0,2*Pi}]
f2 = FunctionInterpolation[x^2,{x,0,10}]
f3 = FunctionInterpolation[Exp[x],{x,0,10}]
Plot[{altintfuncalc[f1,t] - f1[t]},{t,0,2*Pi}]
Plot[{altintfuncalc[f2,t] - f2[t]},{t,0,10}]
Plot[{altintfuncalc[f3,t] - f3[t]},{t,0,10}]
一般来说,它使用分段Hermite三次插值。不过,我不确定节点的选择。看起来它们是在整个区间内统一选择的。我确信假设平滑函数所要求的精度有间隔下限的结果,但我没有详细信息。
Generally it uses piecewise Hermite cubic interpolation. I am not sure about the choice of nodes, though. It seems they are chosen uniformly across the interval. I am sure there are results for the lower bounds of intervals for a requested precision assuming smooth function, but I do not have details.