|
5824
|
1 ## Copyright (C) 2000,2001 Kai Habel |
|
|
2 ## Copyright (C) 2006 David Bateman |
|
|
3 ## |
|
|
4 ## This file is part of Octave. |
|
|
5 ## |
|
|
6 ## Octave is free software; you can redistribute it and/or modify it |
|
|
7 ## under the terms of the GNU General Public License as published by |
|
|
8 ## the Free Software Foundation; either version 2, or (at your option) |
|
|
9 ## any later version. |
|
|
10 ## |
|
|
11 ## Octave is distributed in the hope that it will be useful, but |
|
|
12 ## WITHOUT ANY WARRANTY; without even the implied warranty of |
|
|
13 ## MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the GNU |
|
|
14 ## General Public License for more details. |
|
|
15 ## |
|
|
16 ## You should have received a copy of the GNU General Public License |
|
|
17 ## along with Octave; see the file COPYING. If not, write to the Free |
|
|
18 ## Software Foundation, Inc., 51 Franklin Street, Fifth Floor, Boston, MA |
|
|
19 ## 02110-1301, USA. |
|
|
20 |
|
|
21 ## -*- texinfo -*- |
|
|
22 ## @deftypefn {Function File} {@var{pp} = } spline (@var{x}, @var{y}) |
|
|
23 ## @deftypefnx {Function File} {@var{yi} = } spline (@var{x}, @var{y}, @var{xi}) |
|
|
24 ## |
|
|
25 ## Returns the cubic spline interpolation of @var{y} at the point |
|
5837
|
26 ## @var{x}. Called with two arguments the piece-wise polynomial @var{pp} |
|
5824
|
27 ## that may later be used with @code{ppval} to evaluate the polynomial |
|
|
28 ## at specific points. |
|
|
29 ## |
|
|
30 ## The variable @var{x} must be a vector of length @var{n}, and @var{y} |
|
5837
|
31 ## can be either a vector or array. In the case where @var{y} is a |
|
5824
|
32 ## vector, it can have a length of either @var{n} or @code{@var{n} + 2}. |
|
|
33 ## If the length of @var{y} is @var{n}, then the 'not-a-knot' end |
|
|
34 ## condition is used. If the length of @var{y} is @code{@var{n} + 2}, |
|
|
35 ## then the first and last values of the vector @var{y} are the first |
|
|
36 ## derivative of the cubic spline at the end-points. |
|
|
37 ## |
|
|
38 ## If @var{y} is an array, then the size of @var{y} must have the form |
|
|
39 ## @iftex |
|
|
40 ## @tex |
|
|
41 ## $$[s_1, s_2, \cdots, s_k, n]$$ |
|
|
42 ## @end tex |
|
|
43 ## @end iftex |
|
|
44 ## @ifinfo |
|
|
45 ## @code{[@var{s1}, @var{s2}, @dots{}, @var{sk}, @var{n}]} |
|
|
46 ## @end ifinfo |
|
|
47 ## or |
|
|
48 ## @iftex |
|
|
49 ## @tex |
|
|
50 ## $$[s_1, s_2, \cdots, s_k, n + 2]$$. |
|
|
51 ## @end tex |
|
|
52 ## @end iftex |
|
|
53 ## @ifinfo |
|
|
54 ## @code{[@var{s1}, @var{s2}, @dots{}, @var{sk}, @var{n} + 2]}. |
|
|
55 ## @end ifinfo |
|
|
56 ## The array is then reshaped internally to a matrix where to leading |
|
|
57 ## dimension is given by |
|
|
58 ## @iftex |
|
|
59 ## @tex |
|
|
60 ## $$s_1 s_2 \cdots s_k$$ |
|
|
61 ## @end tex |
|
|
62 ## @end iftex |
|
|
63 ## @ifinfo |
|
|
64 ## @code{@var{s1} * @var{s2} * @dots{} * @var{sk}} |
|
|
65 ## @end ifinfo |
|
|
66 ## and each row this matrix is then treated seperately. Note that this |
|
|
67 ## is exactly the opposite treatment than @code{interp1} and is done |
|
|
68 ## for compatiability. |
|
|
69 ## |
|
|
70 ## Called with a third input argument, @code{spline} evaluates the |
|
|
71 ## piece-wise spline at the points @var{xi}. There is an equivalence |
|
|
72 ## between @code{ppval (spline (@var{x}, @var{y}), @var{xi})} and |
|
|
73 ## @code{spline (@var{x}, @var{y}, @var{xi})}. |
|
|
74 ## |
|
|
75 ## @seealso{ppval, mkpp, unmkpp} |
|
|
76 ## @end deftypefn |
|
|
77 |
|
|
78 ## This code is based on csape.m from octave-forge, but has been |
|
|
79 ## modified to use the sparse solver code in octave that itself allows |
|
|
80 ## special casing of tri-diagonal matrices, modified for NDArrays and |
|
|
81 ## for the treatment of vectors y 2 elements longer than x as complete |
|
|
82 ## splines. |
|
|
83 |
|
|
84 function ret = spline (x, y, xi) |
|
|
85 |
|
|
86 x = x(:); |
|
|
87 n = length (x); |
|
|
88 if (n < 3) |
|
|
89 error ("spline: requires at least 3 points"); |
|
|
90 endif |
|
|
91 |
|
|
92 ## Check the size and shape of y |
|
|
93 ndy = ndims (y); |
|
|
94 szy = size (y); |
|
|
95 if (ndy == 2 && (szy(1) == 1 || szy(2) == 1)) |
|
|
96 if (szy(1) == 1) |
|
6014
|
97 a = y.'; |
|
5824
|
98 else |
|
|
99 a = y; |
|
|
100 szy = fliplr (szy); |
|
|
101 endif |
|
|
102 else |
|
6014
|
103 a = reshape (y, [prod(szy(1:end-1)), szy(end)]).'; |
|
5824
|
104 endif |
|
|
105 complete = false; |
|
|
106 if (size (a, 1) == n + 2) |
|
|
107 complete = true; |
|
|
108 dfs = a(1,:); |
|
|
109 dfe = a(end,:); |
|
|
110 a = a(2:end-1,:); |
|
|
111 endif |
|
|
112 |
|
|
113 b = c = zeros (size (a)); |
|
|
114 h = diff (x); |
|
5838
|
115 idx = ones (columns (a), 1); |
|
5824
|
116 |
|
|
117 if (complete) |
|
|
118 |
|
|
119 if (n == 3) |
|
|
120 dg = 1.5 * h(1) - 0.5 * h(2); |
|
5838
|
121 c(2:n-1,:) = 1/dg(1); |
|
5824
|
122 else |
|
5838
|
123 dg = 2 * (h(1:n-2) .+ h(2:n-1)); |
|
5824
|
124 dg(1) = dg(1) - 0.5 * h(1); |
|
5838
|
125 dg(n-2) = dg(n-2) - 0.5 * h(n-1); |
|
5824
|
126 |
|
5838
|
127 e = h(2:n-2); |
|
5824
|
128 |
|
|
129 size(a) |
|
|
130 size(h) |
|
|
131 n |
|
|
132 |
|
5838
|
133 g = 3 * diff (a(2:n,:)) ./ h(2:n-1,idx) ... |
|
|
134 - 3 * diff (a(1:n-1,:)) ./ h(1:n-2,idx); |
|
|
135 g(1,:) = 3 * (a(3,:) - a(2,:)) / h(2) ... |
|
5824
|
136 - 3 / 2 * (3 * (a(2,:) - a(1,:)) / h(1) - dfs); |
|
5838
|
137 g(n-2,:) = 3 / 2 * (3 * (a(n,:) - a(n-1,:)) / h(n-1) - dfe) ... |
|
|
138 - 3 * (a(n-1,:) - a(n-2,:)) / h(n-2); |
|
5824
|
139 |
|
5838
|
140 c(2:n-1,:) = spdiags ([[e(:); 0], dg, [0; e(:)]], |
|
|
141 [-1, 0, 1], n-2, n-2) \ g; |
|
5824
|
142 endif |
|
|
143 |
|
|
144 c(1,:) = (3 / h(1) * (a(2,:) - a(1,:)) - 3 * dfs |
|
5838
|
145 - c(2,:) * h(1)) / (2 * h(1)); |
|
|
146 c(n,:) = - (3 / h(n-1) * (a(n,:) - a(n-1,:)) - 3 * dfe |
|
|
147 + c(n-1,:) * h(n-1)) / (2 * h(n-1)); |
|
|
148 b(1:n-1,:) = diff (a) ./ h(1:n-1, idx) ... |
|
|
149 - h(1:n-1,idx) / 3 .* (c(2:n,:) + 2 * c(1:n-1,:)); |
|
|
150 d = diff (c) ./ (3 * h(1:n-1, idx)); |
|
5824
|
151 |
|
|
152 else |
|
|
153 |
|
5838
|
154 g = zeros (n-2, columns (a)); |
|
|
155 g(1,:) = 3 / (h(1) + h(2)) ... |
|
|
156 * (a(3,:) - a(2,:) - h(2) / h(1) * (a(2,:) - a(1,:))); |
|
|
157 g(n-2,:) = 3 / (h(n-1) + h(n-2)) ... |
|
|
158 * (h(n-2) / h(n-1) * (a(n,:) - a(n-1,:)) - (a(n-1,:) - a(n-2,:))); |
|
5824
|
159 |
|
|
160 if (n > 4) |
|
|
161 |
|
5838
|
162 g(2:n - 3,:) = 3 * diff (a(3:n-1,:)) ./ h(3:n-2,idx) ... |
|
|
163 - 3 * diff (a(2:n-2,:)) ./ h(2:n - 3,idx); |
|
5824
|
164 |
|
5838
|
165 dg = 2 * (h(1:n-2) .+ h(2:n-1)); |
|
5824
|
166 dg(1) = dg(1) - h(1); |
|
5838
|
167 dg(n-2) = dg(n-2) - h(n-1); |
|
5824
|
168 |
|
5838
|
169 ldg = udg = h(2:n-2); |
|
5824
|
170 udg(1) = udg(1) - h(1); |
|
5838
|
171 ldg(n - 3) = ldg(n-3) - h(n-1); |
|
|
172 c(2:n-1,:) = spdiags ([[ldg(:); 0], dg, [0; udg(:)]], |
|
|
173 [-1, 0, 1], n-2, n-2) \ g; |
|
5824
|
174 |
|
|
175 elseif (n == 4) |
|
|
176 |
|
6248
|
177 dg = [h(1) + 2 * h(2); 2 * h(2) + h(3)]; |
|
5824
|
178 ldg = h(2) - h(3); |
|
|
179 udg = h(2) - h(1); |
|
5838
|
180 c(2:n-1,:) = spdiags ([[ldg(:);0], dg, [0; udg(:)]], |
|
|
181 [-1, 0, 1], n-2, n-2) \ g; |
|
5824
|
182 |
|
|
183 else # n == 3 |
|
|
184 |
|
5838
|
185 dg = h(1) + 2 * h(2); |
|
|
186 c(2:n-1,:) = g/dg(1); |
|
5824
|
187 |
|
|
188 endif |
|
|
189 |
|
|
190 c(1,:) = c(2,:) + h(1) / h(2) * (c(2,:) - c(3,:)); |
|
5838
|
191 c(n,:) = c(n-1,:) + h(n-1) / h(n-2) * (c(n-1,:) - c(n-2,:)); |
|
|
192 b = diff (a) ./ h(1:n-1, idx) ... |
|
|
193 - h(1:n-1, idx) / 3 .* (c(2:n,:) + 2 * c(1:n-1,:)); |
|
|
194 d = diff (c) ./ (3 * h(1:n-1, idx)); |
|
5824
|
195 |
|
|
196 endif |
|
|
197 |
|
5838
|
198 d = d(1:n-1,:); |
|
|
199 c = c(1:n-1,:); |
|
|
200 b = b(1:n-1,:); |
|
|
201 a = a(1:n-1,:); |
|
5824
|
202 coeffs = [d(:), c(:), b(:), a(:)]; |
|
|
203 ret = mkpp (x, coeffs, szy(1:end-1)); |
|
|
204 |
|
|
205 if (nargin == 3) |
|
|
206 ret = ppval (ret, xi); |
|
|
207 endif |
|
|
208 |
|
|
209 endfunction |
|
|
210 |
|
|
211 %!demo |
|
|
212 %! x = 0:10; y = sin(x); |
|
|
213 %! xspline = 0:0.1:10; yspline = spline(x,y,xspline); |
|
|
214 %! title("spline fit to points from sin(x)"); |
|
6702
|
215 %! plot(xspline,sin(xspline),"r",xspline,yspline,"g-",x,y,"b+"); |
|
|
216 %! legend("original","interpolation","interpolation points"); |
|
5824
|
217 %! %-------------------------------------------------------- |
|
|
218 %! % confirm that interpolated function matches the original |
|
|
219 |
|
6686
|
220 %!shared x,y,abserr |
|
|
221 %! x = [0:10]; y = sin(x); abserr = 1e-14; |
|
|
222 %!assert (spline(x,y,x), y, abserr); |
|
|
223 %!assert (spline(x,y,x'), y', abserr); |
|
|
224 %!assert (spline(x',y',x'), y', abserr); |
|
|
225 %!assert (spline(x',y',x), y, abserr); |
|
5824
|
226 %!assert (isempty(spline(x',y',[]))); |
|
|
227 %!assert (isempty(spline(x,y,[]))); |
|
6686
|
228 %!assert (spline(x,[y;y],x), [spline(x,y,x);spline(x,y,x)],abserr) |
|
6014
|
229 %! y = cos(x) + i*sin(x); |
|
6686
|
230 %!assert (spline(x,y,x), y, abserr) |
|
|
231 %!assert (real(spline(x,y,x)), real(y), abserr); |
|
|
232 %!assert (real(spline(x,y,x.')), real(y).', abserr); |
|
|
233 %!assert (real(spline(x.',y.',x.')), real(y).', abserr); |
|
|
234 %!assert (real(spline(x.',y,x)), real(y), abserr); |
|
|
235 %!assert (imag(spline(x,y,x)), imag(y), abserr); |
|
|
236 %!assert (imag(spline(x,y,x.')), imag(y).', abserr); |
|
|
237 %!assert (imag(spline(x.',y.',x.')), imag(y).', abserr); |
|
|
238 %!assert (imag(spline(x.',y,x)), imag(y), abserr); |
