这篇文章是由我的同事Cosmin Ionita。

的“makima”立方插值法最近介绍了MATLAB®R2017b释放的新选项interp1,interp2,interp3,interpn,griddedInterpolant。它的实现并不是用户可见的;因此,我们已经收到我们的用户询问关于这个新立方方法的细节。

在下面,我们解决我们用户的询问回答这些问题:

简而言之,“makima”是短暂的修改Akima分段立方埃尔米特插值。它代表了MATLAB-specific修改Akima的导数公式和具有以下关键属性:

内容

Akima分段立方埃尔米特插值

每个区间[美元x_i ~间{i + 1})的一组输入数据的节点x美元和值v美元美元,分段立方埃尔米特插值不仅发现一个三次多项式插入给定的数据值v_i美元和美元v_ {i + 1}间隔的节点x_i美元和美元间的{i + 1} $,但也有特定的衍生品d_i美元和美元d_ {i + 1} x_i美元和美元间的{i + 1} $。更多细节我们参考插值章在克里夫与MATLAB数值计算教科书。

三次埃尔米特插值的关键是d_i美元衍生品的选择。

在他1970年的纸Hiroshi Akima引入导数公式避免过度当地起伏:

h . Akima”,新的插值方法和基于当地的光滑曲线拟合程序”,JACM诉节,p。589 - 602年,1970年。

让美元\ delta_i = \离开(v_ {i + 1} v_i \右)/ \离开(间{i + 1} -x_i \右)是间隔的斜率(美元x_i ~间{i + 1})美元。Akima x_i美元的导数的定义是:

$ $ d_i = \压裂{| \ delta_ {i + 1} - \ delta_i | \ delta_张{}+ | \ delta_张{}- \ delta_{我2}| \ delta_i} {| \ delta_ {i + 1} - \ delta_i | + | \ delta_张{}- \ delta_{我2}|},$ $

和代表一个斜坡之间的加权平均$ \ delta_张{}$和$ \ delta_i的间隔(美元间张{}~ x_i) $和$ [x_i ~间{i + 1}):美元

$ $ d_i = \压裂{w_1} {w_1 + w_2} \ delta_张{}+ \压裂{w_2} {w_1 + w_2} \ delta_i \,, \四\四w_1 = | \ delta_ {i + 1} - \ delta_i |, \四\四w_2 = | \ delta_张{}- \ delta_{我2}| $ $

注意Akima的导数x_i本地计算的美元五分美元间{我2}$,$间张{}$,x_i美元,美元间{i + 1},美元和美元间的{我+ 2}$。为结束点x_1 x_n美元,美元需要山坡上美元\ delta_ {1}, \ delta_0 $和$ \ delta_n \ delta_ {n + 1} $。因为这些斜坡输入数据不可用,Akima提议使用二次外推计算他们美元\ delta_0 = 2 \ delta_1 - \ delta_2美元,美元\ delta_ {1} = 2 \ delta_0 - \ delta_1 $和$ \ delta_n = 2 \ delta_ {n} - \ delta_ {2} $, $ \ delta_ {n + 1} = 2 \ delta_n - \ delta_ {n} $。

但是等等!

MATLAB已经有两个立方埃尔米特插值方法(见克里夫的博客样条函数和Pchips):

Akima的公式比较如何样条的和“pchip”吗?

Akima的公式发现之间的一个中间立场样条的和“pchip”:

这是一个代表性的例子比较这三个立方埃尔米特interpolants:

x1 = [1 2 3 4 5 5.5 9.5 7 8 9 10];v1 = [0 0 0 0.5 0.4 - 1.2 1.2 - 0.1 0.3 - 0.6);xq1 = 0.75:0.05:10.25;compareCubicPlots (x1, v1, xq1,假的,“不”)

修改Akima插值,“makima”

当然Akima的公式产生好结果。然而,我们还没有准备完全采用。

边缘的边坡

当较低的斜坡上是相等的,美元\ delta_{我2}= \ delta_张{},美元和上层的斜坡上是相等的,$ \ delta_i = \ delta_ {i + 1} $,分子和分母成为0和Akima的公式的回报南。Akima自己意识到了这个问题并提出的平均的上下斜坡边缘案例:$ d_i = \离开(\ delta_张{}+ \ delta_i \右)/ 2美元。

让我们试试这个修复以下数据集的时间间隔(3 ~ 4)美元和美元(4 ~ 5)平等斜坡\ delta_3 = \ delta_4美元= 1和区间[5 ~ 6)美元和美元[6 ~ 7)也有相等的斜坡\ delta_5美元= \ delta_6 = 0美元:

x2 = 1:8;v2 = [1 1 1 0 1 1 1 1];xq2 = 0.75:0.05:8.25;compareCubicPlots (x2, v2, xq2,假的,“本身”)ylim ([-1.2 - 1.2])

在这种情况下,Akima南导数x_5 = 5美元被替换为平均斜率d_5美元= \离开(\ delta_4 + \ delta_5 \右)/ 2 = 0.5美元。

但是我们还没有准备好!

我们应该从Akima平均公式的公式?如果斜坡几乎等于多少?

让我们改变上面的例子中,$ \ delta_5 = \ varepsilon $和$ \ delta_6 = - \ varepsilon通过设置v_6 = 1 + \ varepsilon美元美元,而美元\ varepsilon = 2 ^{-52}代表美元机ε。我们现在可以比较Akima的完整公式与平均公式d_5两美元数据集不同只有\ varepsilon:美元

% v2eps:新的数据集,只有不同的eps v2 (6)v2eps = v2;v2eps (6) = v2eps(6) +每股收益;情节(x2, v2,“柯”,“线宽”2,“MarkerSize”10“DisplayName的”,“输入数据v2”)举行在情节(xq2 akima (x2, v2, xq2),“颜色”(1 2/3 0),“线宽”2,…“DisplayName的”,“v2(6) = 1使用平均公式”)情节(xq2 akima (x2, v2eps, xq2),“-”。,“颜色”(1 2/3 0),“线宽”2,…“DisplayName的”,“v2 (6) = 1 + eps使用Akima”年代公式)举行从传奇(“位置”,“本身”)标题(“Akima interpolant(几乎)等于边坡的)ylim ([-1.2 - 1.2])

平均公式(实线)返回d_5 = 0.5美元而Akima导数的公式(虚线)返回d_5 = 1美元。因此,有一个多余的不小区别x_5 = 5美元左右两Akima interpolants,即使两个基础数据集不同只有\ varepsilon v_6 $美元美元。

要求(1)

我们不希望Akima interpolant切换到一个不同的边界情况的公式。

过度如果数据是恒定的连续超过两个节点

前面的例子是战略选择也透露的另一个属性Akima interpolant:如果数据是恒定的连续超过两个节点,像v_5 = v_6 = v_7 = 1美元在上一个示例中,然后Akima interpolant可能产生超调,即interpolant的波动(5 ~ 6)间隔超过美元。

这种特殊类型的过度在许多工程应用中是不可取的。

同样的问题保持低于出现在[2 ~ 3)间隔美元以上。

要求(2)

我们不希望Akima interpolant产生过度(或低于)当数据是恒定的连续超过两个节点。

修改Akima公式

在MATLAB中,“makima”三次埃尔米特插值地址需求(1)和(2)上面列出。

消除过度,避免边界情况的分子和分母都等于0,我们修改Akima的导数公式通过调整权重w_1 w_2美元美元的斜坡\ delta_美元张{}和$ \ delta_i:美元

$ $ d_i = \压裂{w_1} {w_1 + w_2} \ delta_张{}+ \压裂{w_2} {w_1 + w_2} \ delta_i \,, \四\四w_1 = | \ delta_ {i + 1} - \ delta_i | + | \ delta_ {i + 1} + \ delta_i | / 2 \四\四w_2 = | \ delta_张{}- \ delta_{我2}| + | \ delta_张{}+ \ delta_{我2}| / 2。$ $

这个修改后的公式表示“makima”导数公式中使用MATLAB:

让我们试一试“makima”在上面的公式过度的例子:

compareCubicPlots (x2, v2, xq2,真的,“本身”)ylim ([-1.2 - 1.2])

的确,“makima”不产生超调,如果数据是恒定的超过两个节点(v_5 = v_6 = v_7 = 1美元以上)。

但是这是什么意思的起伏,我们看到在我们的第一个例子吗?

compareCubicPlots (x1, v1, xq1,真的,“不”)

请注意,“makima”与Akima密切遵循结果的公式。事实上,结果是如此相似,很难分辨他们的阴谋。

因此,“makima”仍然保留Akima可取的属性之间的一个中间立场样条的和“pchip”由此产生的波动。

概括一天电网

但是我们还没有完成!

Akima公式和我们的修改“makima”公式有另一种可取的属性:他们推广到高维天网格数据。在1974年他Akima注意到这个属性纸。

h . Akima“二元插值的方法和基于局部表面光滑拟合程序”,ACM通讯,诉丹麦队,p。18 - 20,1974年。

例如,美元在二维网格插入\离开(x, y \右)美元,我们首先应用“makima”导数公式分别在x和y获得美元美元为每个网格节点两方向导数。然后,我们也计算cross-derivative xy美元为每个网格节点。第一cross-derivative公式计算二维划分差异对应二维网格数据和应用“makima”这些差异沿着x美元导数公式;然后,需要结果和应用导数公式以及y美元。衍生品和cross-derivatives然后插入两变量的系数表示二维interpolant立方埃尔米特多项式。

同样的程序适用于高维天网格与n \通用电气2美元和需要计算所有可能的cross-directions cross-derivatives。因此,“makima”不仅在万博1manbetx支持interp1,而且在interp2,interp3,interpn,griddedInterpolant。

的二维“makima”插值与二维样条的在生成网格数据的插值山峰功能:

(X1, Y1, V1 =山峰(5);[Xq1, Yq1] = meshgrid (3:0.1:3 3:0.1:3);Vqs1 = interp2 (X1, Y1, V1, Xq1, Yq1,样条的);冲浪(Xq1 Yq1 Vqs1)轴紧标题(”二维样条”)snapnow Vqm1 = interp2 (X1, Y1, V1, Xq1, Yq1,“makima”);冲浪(Xq1 Yq1 Vqm1)轴紧标题(“二维“makima””)

注意到小起伏(或摆动)生成的“makima”。

最后,让我们试一试“makima”上一个例子和一些二维峰数据有锐利的边缘和步骤:

V2 = 0 (10);V2 (2:5, 2:5) = 3/7;%的第一步V2(者者)= (1/4 1/5;1/5 1/4);%的中间步骤V2 (8:9 8:9) = 1/2;%最后一步V2 =翻转(V2, 2);[Xq2, Yq2] = meshgrid (1:0.2:10 1:0.2:10);Vqs2 = interp2 (V2, Xq2 Yq2,样条的);冲浪(Xq2 Yq2 Vqs2)轴紧标题(”二维样条”)snapnow Vqm2 = interp2 (V2, Xq2 Yq2,“makima”);冲浪(Xq2 Yq2 Vqm2)轴紧标题(“二维“makima””)

的“makima”结果小起伏(或摆动)样条。

总结

最后,我们总结的属性感兴趣的三个立方埃尔米特插值方法在MATLAB支持:万博1manbetx

摘要=表(…(“C2”;“使用所有点”;“Not-a-knot条件”;“是的”),…(“C1”;“使用4点”;“单边坡”;“不”),…(“C1”;“使用5点”;“二次推断”;“是的”),…“VariableNames”,…(“样条”“pchip”“makima”),…“RowNames”,…(“连续性”;“位置”;“终点”;“万博1manbetx支持一天”]);disp(总结)

花键pchip makima ______________________ _________________ _________________________连续性“C2”“C1”“C1”地区“使用所有点”“使用4点”“使用5点“端点”Not-a-knot条件”“单边坡”“二次推断”支持一天“是”“否”“是”万博1manbetx

代码

compareCubicPlots

函数compareCubicPlots (x, v, xq, showMakima legendLocation)%的阴谋pchip,花键,Akima, makima的插值结果vqp = pchip (x, v, xq);%一样vq = interp1 (x, v, xq, pchip)vq =花键(x, v, xq);%一样vq = interp1 (x, v, xq,样条)酒瓶= akima (x, v, xq);如果showMakima vqm = makima (x, v, xq);%一样vq = interp1 (x, v, xq, makima)结束情节(x, v,“柯”,“线宽”2,“MarkerSize”10“DisplayName的”,输入数据的)举行在情节(vqp xq,“线宽”4“DisplayName的”,“pchip”)情节(xq,矢量量化,“线宽”2,“DisplayName的”,“样条”)情节(xq,酒瓶,“线宽”2,“DisplayName的”,Akima”年代公式)如果showMakima情节(vqm xq,“——”,“颜色”(3/4 0),“线宽”2,…“DisplayName的”,“makima”)结束持有从xticks (x (1) 1: x(结束)+ 1)传说(“位置”legendLocation)标题(“三次埃尔米特插值)结束

akima

函数vq = akima (x, v, xq)% AKIMA AKIMA分段立方埃尔米特插值。%% vq = AKIMA (x, v, xq)篡改数据(x, v) xq查询点%使用Akima立方埃尔米特插值公式。%%的引用:%% h . Akima”,插值和平滑的曲线拟合的新方法%根据当地程序”,JACM诉节,p。589 - 602年,1970年。%% AKIMA与PCHIP与花键:%%——Akima立方公式是花键之间的中间地带,PCHIP:%比花键lower-amplitude摆动,但没有攻击性%减少PCHIP摆动。%——Akima立方公式和花键推广天网格。%%的例子:比较AKIMA MAKIMA,花键,PCHIP%% x = (1 2 3 4 5 5.5 9.5 7 8 9 10];% v = [0 0 0 0.5 0.4 1.2 1.2 0.1 0.3 - 0.6);% xq = 0.75:0.05:10.25;% vq =花键(x, v, xq);% vqp = pchip (x, v, xq);%酒瓶= akima (x, v, xq);% vqm = makima (x, v, xq);%%图%的阴谋(x, v,“柯”,“线宽”,2,MarkerSize, 10),等等%的阴谋(vqp xq,“线宽”,4)%的阴谋(vq xq, xq,酒瓶,xq, vqm,“——”,“线宽”,2)%传奇(“数据”,“pchip””,“样条“”、“Akima”年代公式”,“‘makima“”)%的标题(“MATLAB三次埃尔米特插值”)%%也看到MAKIMA,花键,PCHIP。% mailto: cosmin.ionita@mathworks.com断言(isvector (x) & & isvector (v) & &(元素个数(x) = =元素个数(v)))断言(元素个数(x) > 2) x = x (:)。';v = v (:)。';%在pchip形状一样。m和spline.m%计算Akima斜坡h = diff (x);δ= diff (v) / h;斜坡= akimaSlopes(δ);%计算分段立方埃尔米特interpolantvq = ppval (pwch (x, v,山坡,h, delta), xq);结束

akimaSlopes

函数s = akimaSlopes(δ)% Akima立方埃尔米特插值的导数值在网格节点x % Akima导数估计(我)需要四个有限%的差异五个网格节点对应x(我2:我+ 2)。%%的边界网格节点(1:2)和x (n - 1: n),附加有限的差异%将对应于(1:0)和x (n + 1: n + 2)通过使用以下%偏心差公式correspondin二次推断使用定义的二次多项式数据% x (1:3)% (Akima的论文2.3节):n =元素个数(δ)+ 1;%的网格节点xdelta_0 = 2 *(1) -δ(2);delta_m1 = 2 * delta_0 -δ(1);delta_n = 2 * (n - 1) -δ(n - 2);delta_n1 = 2 * delta_n -δ(n - 1);δ= [delta_m1 delta_0三角洲delta_n delta_n1];% Akima导数估计公式(方程式(1)摘要):%% h . Akima”,插值和平滑的曲线拟合的新方法%根据当地程序”,JACM诉节,p。589 - 602年,1970年。%% s (i) = (| d - d (i + 1)(我)| * d(张)+ | d(张)- d(我2)| * d(我))% / (| d (i + 1) - d (i) | + | d(张)- d(我2)|)重量= abs (diff(δ));weights1 =重量(1:n);% | d - d(我2)|(张)weights2 =重量(3:结束);% | d (i + 1) - d(我)|delta1 =δ(2:n + 1);% d(张)delta2 =δ(3:n + 2);% d(我)weights12 = weights1 + weights2;s = (weights2. / weights12)。* delta1 + (weights1. / weights12)。* delta2;%,以避免0/0,Akima提出平均划分不同d(张)%和d (i)的边界情况(我2)= d(张)和d (i) = d (i + 1):印第安纳州= weights1 = = 0 & weights2 = = 0;(印第安纳州)= (delta1(印第安纳州)+ delta2(印第安纳州))/ 2;结束

makima

函数vq = makima (x, v, xq)% MAKIMA修改Akima分段立方埃尔米特插值。%%我们修改Akima公式消除过度和脱靶%当数据是恒定的连续超过两个节点。%% vq = MAKIMA (x, v, xq)篡改数据(x, v) xq查询点%使用修改后的Akima立方埃尔米特插值公式。%%的引用:%% h . Akima”,插值和平滑的曲线拟合的新方法%根据当地程序”,JACM诉节,p。589 - 602年,1970年。%% MAKIMA与PCHIP与花键:%% - MAKIMA是花键之间的中间地带,PCHIP:%比花键lower-amplitude摆动,但没有攻击性%减少PCHIP摆动。%——MAKIMA和花键概括到一天网格。%%的例子:不超过和低于MAKIMA数据时%常数连续超过两个节点%% x = 1:7;% v = (1 1 1 0 1 1 1);% xq = 0.75:0.05:7.25;% vq =花键(x, v, xq);% vqp = pchip (x, v, xq);%酒瓶= akima (x, v, xq);% vqm = makima (x, v, xq);%%图%的阴谋(x, v,“柯”,“线宽”,2,MarkerSize, 10),等等%的阴谋(vqp xq,“线宽”,4)%的阴谋(vq xq, xq,酒瓶,xq, vqm,“线宽”,2)%的传说(“数据”,“pchip”’,”样条“”、“Akima”年代公式”,…%”makima””、“位置”、“SE”)%标题(“makima”没有过度(1:3)和x (7)”)%%的例子:比较MAKIMA AKIMA,花键,PCHIP%% x = (1 2 3 4 5 5.5 9.5 7 8 9 10];% v = [0 0 0 0.5 0.4 1.2 1.2 0.1 0.3 - 0.6);% xq = 0.75:0.05:10.25;% vq =花键(x, v, xq);% vqp = pchip (x, v, xq);%酒瓶= akima (x, v, xq);% vqm = makima (x, v, xq);%%图%的阴谋(x, v,“柯”,“线宽”,2,MarkerSize, 10),等等%的阴谋(vqp xq,“线宽”,4)%的阴谋(vq xq, xq,酒瓶,xq, vqm,“——”,“线宽”,2)%传奇(“数据”,“pchip””,“样条“”、“Akima”年代公式”,“‘makima“”)%的标题(“MATLAB三次埃尔米特插值”)%%也看到AKIMA,花键,PCHIP。% mailto: cosmin.ionita@mathworks.com断言(isvector (x) & & isvector (v) & &(元素个数(x) = =元素个数(v)))断言(元素个数(x) > 2) x = x (:)。';v = v (:)。';%在pchip形状一样。m和spline.m%计算修改Akima斜坡h = diff (x);δ= diff (v) / h;斜坡= makimaSlopes(δ);%计算分段立方埃尔米特interpolantvq = ppval (pwch (x, v,山坡,h, delta), xq);结束

makimaSlopes

函数s = makimaSlopes(δ)%的导数值修改Akima立方埃尔米特插值在网格节点x % Akima导数估计(我)需要四个有限%的差异五个网格节点对应x(我2:我+ 2)。%%的边界网格节点(1:2)和x (n - 1: n),附加有限的差异%将对应于(1:0)和x (n + 1: n + 2)通过使用以下%偏心差公式correspondin二次推断使用定义的二次多项式数据% x (1:3)% (Akima的论文2.3节):n =元素个数(δ)+ 1;%的网格节点xdelta_0 = 2 *(1) -δ(2);delta_m1 = 2 * delta_0 -δ(1);delta_n = 2 * (n - 1) -δ(n - 2);delta_n1 = 2 * delta_n -δ(n - 1);δ= [delta_m1 delta_0三角洲delta_n delta_n1];% Akima导数估计公式(方程式(1)摘要):%% h . Akima”,插值和平滑的曲线拟合的新方法%根据当地程序”,JACM诉节,p。589 - 602年,1970年。%% s (i) = (| d - d (i + 1)(我)| * d(张)+ | d(张)- d(我2)| * d(我))% / (| d (i + 1) - d (i) | + | d(张)- d(我2)|)%%消除过度和脱靶当数据是恒定的%比连续两个节点,在MATLAB的我们修改Akima makima%的公式通过添加额外的平均权重。% s (i) = ((| d - d (i + 1)(我)| + | d (i + 1) + d (i) | / 2) * d(张)+% (| d(张)- d(我2)| + | d(张)+ d(我2)| / 2)* d (i))% / ((| d - d (i + 1)(我)| + | d (i + 1) + d (i) | / 2) +% (| d(张)- d(我2)| + | d(张)+ d(我2)| / 2)重量= abs (diff(δ))+ abs (((1: end-1) +δ(2:结束))/ 2);weights1 =重量(1:n);% | d - d(我2)|(张)weights2 =重量(3:结束);% | d (i + 1) - d(我)|delta1 =δ(2:n + 1);% d(张)delta2 =δ(3:n + 2);% d(我)weights12 = weights1 + weights2;s = (weights2. / weights12)。* delta1 + (weights1. / weights12)。* delta2;%如果数据是恒定的连续超过4个节点,然后%分母为零和的公式产生不必要的南的结果。%这南替换为0。s (weights12 = = 0) = 0;结束

版权2019年MathWorks公司。
得到了MATLAB代码

发表与MATLAB®R2018b