分段三次hermite插值

久病成欢

% 这是MATLAB里面的pchip.m文件，这里把它的注释改写成汉语，主要是想弄清楚它是怎么计算在节点处的导数的。

1. function v = pchip(x,y,xx)
   
2. %输入：n个插值节点的纵坐标向量x；横坐标向量y；插值点xx。
   
3. %输出：分段三次Hermite插值结果。
   
4. 
   
5. %   PCHIP  Piecewise Cubic Hermite Interpolating Polynomial.
   
6. %   PP = PCHIP(X,Y)为X处的值Y提供了一种特定的保形分段三次厄尔米特插值（shape-preserving piecewise cubic Hermite interpolant）
   
7. %   的分段多项式形式，在后面的PPVAL和样条功能UNMKPP（spline utility UNMKPP）将用到这个函数。
   
8. %   X必须是个向量。
   
9. %   如果Y是个向量，则Y的第j个元素Y(j)被取为和X的第j个元素X(j)匹配的值，因此Y和X的长度必须一样。
   
10. %   如果Y是一个矩阵，或者N维数组，则Y(:,...,:,j)被取为和X(j)相匹配的值，因此Y的最后一维必须等于length(X).
    
11. %   YY = PCHIP(X,Y,XX)和YY = PPVAL(PCHIP(X,Y),XX)是一样的,因此在YY中给出了在XX处的插值。
    
12. %   PCHIP插值函数p(x)满足:
    
13. %   在每个子区间X(k) <= x <= X(k+1),p(x)都是三阶Hermite插值多项式（给定插值点和两个端点的斜率）。
    
14. %   因而，p(x) interpolates Y,也就是说，p(X(j)) = Y(:,j),并且一阶导数Dp(x)是连续的，但是
    
15. %   二阶导数D^2p(x)可能不是连续的;在X(j)处可能会出现跳跃.
    
16. %   在X(j)处的斜率的选取方法，确保了p(x)是"shape preserving"和"respects monotonicity"的，
    
17. %   这意味着，在那些数据是单调的区间里，p(x)也是单调的；在那些数据是局部极值（local extremum）的点，
    
18. %   p(x)也取局部极值。
    
19. %
    
20. %  PCHIP与SPLINE的对比:
    
21. %   SPLINE提供的函数s(x)的构建方法和PCHIP里面的函数p(x)完全相同，只不过在X(j)处的斜率的选择方法不一样，
    
22. %   SPLINE函数的s(x)在X(j)的二阶导数D^2s(x)也是连续的，这导致了如下结果：
    
23. %   SPLINE更加光滑，也就是说，D^2s(x)是连续的。
    
24. %   如果数据是一个光滑函数的值，则SPLINE更加精确。
    
25. %   如果数据不是光滑的，则PCHIP没有overshoots，也不太震荡（less oscillation）。
    
26. %   PCHIP建立的难度较小（is less expensive to set up）.
    
27. %   这两种函数估计的难度是一样的。
    
28. 
    
29. %   样条比pchip光滑，样条的两阶导数连续，而pchip一阶导数连续。不连续的两阶导数隐含着不连续的曲率。人的眼睛可以检测出图形上曲率的不连续。另一方面，pchip是保形状的，而样条不一定保形状。
    
30. 
    
31. %
    
32. %   例子:
    
33. %     x = -3:3;
    
34. %     y = [-1 -1 -1 0 1 1 1];
    
35. %     t = -3:.01:3;
    
36. %     plot(x,y,'o',t,[pchip(x,y,t); spline(x,y,t)])
    
37. %     legend('data','pchip','spline',4)
    
38. %
    
39. %   Class support for inputs x, y, xx:
    
40. %      float: double, single
    
41. %
    
42. %   还可参见INTERP1, SPLINE, PPVAL, UNMKPP.
    
43. 
    
44. % 参考文献:
    
45. %   F. N. Fritsch and R. E. Carlson, "Monotone Piecewise Cubic
    
46. %   Interpolation", SIAM J. Numerical Analysis 17, 1980, 238-246.
    
47. %   David Kahaner, Cleve Moler and Stephen Nash, Numerical Methods
    
48. %   and Software, Prentice Hall, 1988.
    
49. %
    
50. %   Copyright 1984-2004 The MathWorks, Inc.
    
51. %   $Revision: 1.7.4.4 [        DISCUZ_CODE_2        ]nbsp; $Date: 2004/03/02 21:47:53 $
    
52. 
    
53. % 检验数据的可接受性，如果不可接受，则对其进行适当的调整
    
54. [x,y,sizey] = chckxy(x,y);          %chckxy返回三个变量：x,y,和sizey。但是不知道chckxy是什么意思。
    
55. n = length(x);                      %n为向量x的长度。也就是后面要用的节点数目。
    
56. h = diff(x);                        %diff表示把向量x的相邻元素相减。得到h=[X(2)-X(1) X(3)-X(2) ... X(n)-X(n-1)]
    
57. m = prod(sizey);                    %
    
58. 
    
59. %确保插值点是实数
    
60. if nargin==3 && any(~isreal(reshape(xx,numel(xx),1)))
    
61.   error('MATLAB:pchip:ComplexInterpPts',...
    
62.         'The interpolation points should be real.')
    
63. end
    
64. 
    
65. %计算斜率
    
66. del = diff(y,1,2)./repmat(h,m,1);
    
67. % diff(y,n,dim)是在标量dim指定的维度上进行n次差分。如果阶数n等于或超过第dim维的长度，则diff返回一个空的数组。
    
68. % 例如y=[1 3 4 6 9 10 8 12];则diff(y,1,2)=[2 1 2 3 1 -2 4];
    
69. % repmat(h,m,1)把矩阵h在纵向方面复制m次。二者相除就是一阶差商。
    
70. slopes = zeros(size(y));                                % 设定一个全是0的向量，准备存放斜率数值。
    
71. for r = 1:m
    
72.       slopes(r,:) = pchipslopes(x,y(r,:),del(r,:));     %调用函数见下。
    
73. end
    
74. 
    
75. % 对上述值x,y,和斜率计算分段三次Hermite插值
    
76. v = pwch(x,y,slopes,h,del); v.dim = sizey;
    
77. if nargin == 3   % if values are wanted instead, provide them
    
78.    v = ppval(v,xx);
    
79. end

复制代码

% 下面是计算节点处的斜率的函数pchipslopes------------------------------------------

1. function d = pchipslopes(x,y,del)
   
2. %PCHIPSLOPES Derivative values for shape-preserving Piecewise Cubic Hermite Interpolation.
   
3. % d = pchipslopes(x,y,del)计算一阶导数d(k)=P'(x(k)).
   
4. 
   
5. % 特殊情况：n=2,此时使用线性插值.
   
6.    n = length(x);
   
7.    if n==2
   
8.       d = repmat(del(1),size(y));
   
9.       return
   
10.    end
    
11. 
    
12. %  内点（interior points）处的斜率.
    
13. %  如果第k个节点处的左右差商del(k-1)和del(k)符号相同，则设定d(k)等于二者的加权平均。
    
14. %  如果第k个节点处的左右差商del(k-1)和del(k)符号相反，或者其中一个为0，则设定d(k)=0.
    
15.    d = zeros(size(y));
    
16.    if isreal(del)                                          %如果del是实数。
    
17.       k = find(sign(del(1:n-2)).*sign(del(2:n-1)) > 0);    %则把其左右差商同号的那个序号赋值给k.
    
18.    else
    
19.       k = find(~(del(1:n-2) == 0 & del(2:n-1) == 0));
    
20.    end
    
21.    h = diff(x);
    
22.    hs = h(k)+h(k+1);
    
23.    w1 = (h(k)+hs)./(3*hs);
    
24.    w2 = (hs+h(k+1))./(3*hs);
    
25.    dmax = max(abs(del(k)), abs(del(k+1)));
    
26.    dmin = min(abs(del(k)), abs(del(k+1)));
    
27.    d(k+1) = dmin./conj(w1.*(del(k)./dmax) + w2.*(del(k+1)./dmax));
    
28. %函数congj(a)返回数组a的每个元素的共轭复数组成的数组。
    
29. 
    
30.   
    
31. %  区间端点处的斜率（end points）.
    
32. %  Set d(1) and d(n) via non-centered, shape-preserving three-point formulae.
    
33. 
    
34.    d(1) = ((2*h(1)+h(2))*del(1) - h(1)*del(2))/(h(1)+h(2));
    
35.    if isreal(d) && (sign(d(1)) ~= sign(del(1)))
    
36.       d(1) = 0;
    
37.    elseif (sign(del(1)) ~= sign(del(2))) && (abs(d(1)) > abs(3*del(1)))
    
38.       d(1) = 3*del(1);
    
39.    end
    
40.    d(n) = ((2*h(n-1)+h(n-2))*del(n-1) - h(n-1)*del(n-2))/(h(n-1)+h(n-2));
    
41.    if isreal(d) && (sign(d(n)) ~= sign(del(n-1)))
    
42.       d(n) = 0;
    
43.    elseif (sign(del(n-1)) ~= sign(del(n-2))) && (abs(d(n)) > abs(3*del(n-1)))
    
44.       d(n) = 3*del(n-1);
    
45.    end

复制代码

happy

谢谢分享相关代码