三次样条插值

nshukwrd

1. function output = spline(x,y,xx)
   
2. %三次样条插值
   
3. %yi=spline(x,y,xi)等价于YI=interp(x,y,xi,'spline'),
   
4. %     根据数据(x,y)给出在xi的线性插值结果yi..
   
5. %     使用"非扭结"端点条件, 即强迫第一﹑二端多项式三次项系数相同,
   
6. %     最后一段和倒数第二段三次项系数相同.
   
7. %pp=spline(x,y)返回样条插值的分段多项式(pp形式).
   
8. %[breaks,coefs]=unmkpp(pp)将pp形式展开,其中breaks为结点,coefs为各段多项式系数.
   
9. %yi=ppval(pp,xi),pp形式在xi的函数值.
   
10. %例 考虑数据
    
11. %          x  | 1  2  4  5
    
12. %         ---|-------------
    
13. %          y  | 1  3  4  2
    
14. %      clear;close;
    
15. %      x=[1 2 4 5];y=[1 3 4 2];
    
16. %      p=spline(x,y);
    
17. %      xi=1:0.1:5;yi=ppval(p,xi);
    
18. %      plot(x,y,'o',xi,yi,'k');
    
19. %      title('not-a-knot SPLINE');
    
20. %      [b,c]=unmkpp(p)
    
21. %另一个例子见下列英文部分
    
22. %
    
23. %SPLINE Cubic spline data interpolation.
    
24. %   YY = SPLINE(X,Y,XX) uses cubic spline interpolation
    
25. %   to find a vector YY corresponding to XX.  X and Y are the
    
26. %   given data vectors and XX is the new abscissa vector.
    
27. %
    
28. %   The ordinates Y may be vector-valued, in which case Y(:,j) is
    
29. %   the j-th ordinate.
    
30. %
    
31. %   PP = SPLINE(X,Y) returns the pp-form of the cubic spline
    
32. %   interpolant instead, for later use with  PPVAL, etc.
    
33. %
    
34. %   Ordinarily, the not-a-knot end conditions are used. However, if
    
35. %   Y contains two more ordinates than X has entries, then the first
    
36. %   and last ordinate in Y are used as the endslopes for the cubic spline.
    
37. %
    
38. %   Here's an example that generates a coarse sine curve, then
    
39. %   samples the spline over a finer mesh:
    
40. %
    
41. %       x = 0:10;  y = sin(x);
    
42. %       xx = 0:.25:10;
    
43. %       yy = spline(x,y,xx);
    
44. %       plot(x,y,'o',xx,yy)
    
45. %
    
46. %   Here is an example that features a vector-valued spline, along with complete
    
47. %   spline interpolation, i.e., fitting to given end slopes (instead of using the
    
48. %   not-a-knot end condition); it uses SPLINE to generate a circle:
    
49. %
    
50. %       circle = spline( 0:4, [0 1 0 -1 0 1 0; pi/2 0 1 0 -1 0 pi/2] );
    
51. %       xx = 0:.1:4; cc = ppval(circle, xx); plot(cc(1,:), cc(2,:)), axis equal
    
52. %
    
53. %   See also INTERP1, PPVAL, SPLINES (The Spline Toolbox).
    
54. 
    
55. %   Carl de Boor 7-2-86
    
56. %   Revised 11-24-87 JNL, 6-16-92 CBM, 10-14-97 CB.
    
57. %   Copyright (c) 1984-98 by The MathWorks, Inc.
    
58. %   $Revision: 5.11 $  $Date: 1997/12/03 19:22:33 $
    
59. 
    
60. % Generate the cubic spline interpolant in pp form, depending on the
    
61. % number of data points (and usually using the not-a-knot end condition).
    
62. 
    
63. output=[];
    
64. n=length(x);
    
65. if n<2, error('There should be at least two data points.'), end
    
66. 
    
67. if any(diff(x)<0), [x,ind]=sort(x); else, ind=1:n; end
    
68. 
    
69. x=x(:); dx = diff(x);
    
70. if all(dx)==0, error('The data abscissae should be distinct.'), end
    
71. 
    
72. [yd,yn] = size(y); % if Y happens to be a column matrix, change it to
    
73.                    % the expected row matrix.
    
74. if yn==1, yn=yd; y=reshape(y,1,yn); yd=1; end
    
75. 
    
76. if yn==n
    
77.    notaknot = 1;
    
78. elseif yn==n+2
    
79.    notaknot = 0; endslopes = y(:,[1 n+2]).'; y(:,[1 n+2])=[];
    
80. else
    
81.    error('Abscissa and ordinate vector should be of the same length.')
    
82. end
    
83. 
    
84. yi=y(:,ind).'; dd = ones(1,yd);
    
85. dx = diff(x); divdif = diff(yi)./dx(:,dd);
    
86. if n==2
    
87.    if notaknot, % the interpolant is a straight line
    
88.       pp=mkpp(x.',[divdif.' yi(1,:).'],yd);
    
89.    else         % the interpolant is the cubic Hermite polynomial
    
90.       divdif2 = diff([endslopes(1,:);divdif;endslopes(2,:)])./dx([1 1],dd);
    
91.       pp = mkpp(x,...
    
92.       [(diff(divdif2)./dx(1,dd)).' ([2 -1]*divdif2).' ...
    
93.                                            endslopes(1,:).' yi(1,:).'],yd);
    
94.    end
    
95. elseif n==3¬aknot, % the interpolant is a parabola
    
96.    yi(2:3,:)=divdif;
    
97.    yi(3,:)=diff(divdif)/(x(3)-x(1));
    
98.    yi(2,:)=yi(2,:)-yi(3,:)*dx(1);
    
99.    pp = mkpp([x(1),x(3)],yi([3 2 1],:).',yd);
    
100. else % set up the sparse, tridiagonal, linear system for the slopes at  X .
     
101.    b=zeros(n,yd);
     
102.    b(2:n-1,:)=3*(dx(2:n-1,dd).*divdif(1:n-2,:)+dx(1:n-2,dd).*divdif(2:n-1,:));
     
103.    if notaknot
     
104.       x31=x(3)-x(1);xn=x(n)-x(n-2);
     
105.       b(1,:)=((dx(1)+2*x31)*dx(2)*divdif(1,:)+dx(1)^2*divdif(2,:))/x31;
     
106.       b(n,:)=...
     
107.       (dx(n-1)^2*divdif(n-2,:)+(2*xn+dx(n-1))*dx(n-2)*divdif(n-1,:))/xn;
     
108.    else
     
109.       x31 = 0; xn = 0; b([1 n],:) = dx([1 n-2],dd).*endslopes;
     
110.    end
     
111.    c = spdiags([ [dx(2:n-1);xn;0] ...
     
112.         [dx(2);2*[dx(2:n-1)+dx(1:n-2)];dx(n-2)] ...
     
113.         [0;x31;dx(1:n-2)] ],[-1 0 1],n,n);
     
114. 
     
115.    % sparse linear equation solution for the slopes
     
116.    mmdflag = spparms('autommd');
     
117.    spparms('autommd',0);
     
118.    s=c\b;
     
119.    spparms('autommd',mmdflag);
     
120.    % convert to pp form
     
121.    c4=(s(1:n-1,:)+s(2:n,:)-2*divdif(1:n-1,:))./dx(:,dd);
     
122.    c3=(divdif(1:n-1,:)-s(1:n-1,:))./dx(:,dd) - c4;
     
123.    pp=mkpp(x.',...
     
124.      reshape([(c4./dx(:,dd)).' c3.' s(1:n-1,:).' yi(1:n-1,:).'],(n-1)*yd,4),yd);
     
125. end
     
126. if nargin==2
     
127.    output=pp;
     
128. else
     
129.    output=ppval(pp,xx);
     
130. end
     

复制代码