Mercurial > forge

comparison main/splines/csape.m @ 0:6b33357c7561 octave-forge

Find changesets by keywords (author, files, the commit message), revision number or hash, or revset expression.
Initial revision
author pkienzle
date Wed, 10 Oct 2001 19:54:49 +0000
parents
children bac8128dc91a
comparison
equal deleted inserted replaced
-1:000000000000 0:6b33357c7561
1 ## Copyright (C) 2000,2001 Kai Habel
2 ##
3 ## This program is free software; you can redistribute it and/or modify
4 ## it under the terms of the GNU General Public License as published by
5 ## the Free Software Foundation; either version 2 of the License, or
6 ## (at your option) any later version.
7 ##
8 ## This program is distributed in the hope that it will be useful,
9 ## but WITHOUT ANY WARRANTY; without even the implied warranty of
10 ## MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
11 ## GNU General Public License for more details.
12 ##
13 ## You should have received a copy of the GNU General Public License
14 ## along with this program; if not, write to the Free Software
15 ## Foundation, Inc., 59 Temple Place, Suite 330, Boston, MA 02111-1307 USA
16
17 ## -*- texinfo -*-
18 ## @deftypefn {Function File} {@var{pp} = } csape (@var{x}, @var{y}, @var{cond}, @var{valc})
19 ## cubic spline interpolation with various end conditions.
20 ## creates the pp-form of the cubic spline.
21 ##
22 ## the following end conditions as given in @var{cond} are possible.
23 ## @table @asis
24 ## @item 'complete'
25 ## match slopes at first and last point as given in @var{valc}
26 ## @item 'not-a-knot'
27 ## third derivatives are continuous at the second and second last point
28 ## @item 'periodic'
29 ## match first and second derivative of first and last point
30 ## @item 'second'
31 ## match second derivative at first and last point as given in @var{valc}
32 ## @item 'variational'
33 ## set second derivative at first and last point to zero (natural cubic spline)
34 ## @end table
35 ##
36 ## @seealso{ppval, spline}
37 ## @end deftypefn
38
39 ## Author: Kai Habel <kai.habel@gmx.de>
40 ## Date: 23. nov 2000
41 ## Algorithms taken from G. Engeln-Muellges, F. Uhlig:
42 ## "Numerical Algorithms with C", Springer, 1996
43
44 ## Paul Kienzle, 19. feb 2001, csape supports now matrix y value
45
46 function pp = csape (x, y, cond, valc)
47
48 x = x(:);
49 n = length(x);
50
51 transpose = (columns(y) == n);
52 if (transpose) y = y'; endif
53
54 a = y;
55 b = c = zeros (size (y));
56 h = diff (x);
57 idx = ones (columns(y),1);
58
59 if (nargin < 3 || strcmp(cond,"complete"))
60 # specified first derivative at end point
61 if (nargin < 4)
62 valc = [0, 0];
63 endif
64
65 dg = 2 * (h(1:n - 2) .+ h(2:n - 1));
66 dg(1) = dg(1) - 0.5 * h(1);
67 dg(n - 2) = dg(n-2) - 0.5 * h(n - 1);
68
69 e = h(2:n - 2);
70
71 g = 3 * diff (a(2:n,:)) ./ h(2:n - 1,idx)\
72 - 3 * diff (a(1:n - 1,:)) ./ h(1:n - 2,idx);
73 g(1,:) = 3 * (a(3,:) - a(2,:)) / h(2) \
74 - 3 / 2 * (3 * (a(2,:) - a(1,:)) / h(1) - valc(1));
75 g(n - 2,:) = 3 / 2 * (3 * (a(n,:) - a(n - 1,:)) / h(n - 1) - valc(2))\
76 - 3 * (a(n - 1,:) - a(n - 2,:)) / h(n - 2);
77
78 c(2:n - 1,:) = trisolve(dg,e,g);
79 c(1,:) = (3 / h(1) * (a(2,:) - a(1,:)) - 3 * valc(1)
80 - c(2,:) * h(1)) / (2 * h(1));
81 c(n,:) = - (3 / h(n - 1) * (a(n,:) - a(n - 1,:)) - 3 * valc(2)
82 + c(n - 1,:) * h(n - 1)) / (2 * h(n - 1));
83 b(1:n - 1,:) = diff (a) ./ h(1:n - 1, idx)\
84 - h(1:n - 1,idx) / 3 .* (c(2:n,:) + 2 * c(1:n - 1,:));
85 d = diff (c) ./ (3 * h(1:n - 1, idx));
86
87 elseif (strcmp(cond,"variational") || strcmp(cond,"second"))
88
89 if ((nargin < 4) || strcmp(cond,"variational"))
90 ## set second derivatives at end points to zero
91 valc = [0, 0];
92 endif
93
94 c(1,:) = valc(1) / 2;
95 c(n,:) = valc(2) / 2;
96
97 g = 3 * diff (a(2:n,:)) ./ h(2:n - 1, idx)\
98 - 3 * diff (a(1:n - 1,:)) ./ h(1:n - 2, idx);
99
100 g(1,:) = g(1,:) - h(1) * c(1,:);
101 g(n - 2,:) = g(n-2,:) - h(n - 1) * c(n,:);
102
103 dg = 2 * (h(1:n - 2) .+ h(2:n - 1));
104 e = h(2:n - 2);
105
106 c(2:n - 1,:) = trisolve (dg,e,g);
107 b(1:n - 1,:) = diff (a) ./ h(1:n - 1,idx)\
108 - h(1:n - 1,idx) / 3 .* (c(2:n,:) + 2 * c(1:n - 1,:));
109 d = diff (c) ./ (3 * h(1:n - 1, idx));
110
111 elseif (strcmp(cond,"periodic"))
112
113 h = [h; h(1)];
114
115 ## XXX FIXME XXX --- the following gives a smoother periodic transition:
116 ## a(n,:) = a(1,:) = ( a(n,:) + a(1,:) ) / 2;
117 a(n,:) = a(1,:);
118
119 tmp = diff (shift ([a; a(2,:)], -1));
120 g = 3 * tmp(1:n - 1,:) ./ h(2:n,idx)\
121 - 3 * diff (a) ./ h(1:n - 1,idx);
122
123 if (n > 3)
124 dg = 2 * (h(1:n - 1) .+ h(2:n));
125 e = h(2:n - 1);
126 c(2:n,idx) = trisolve(dg,e,g,h(1),h(1));
127 elseif (n == 3)
128 A = [2 * (h(1) + h(2)), (h(1) + h(2));
129 (h(1) + h(2)), 2 * (h(1) + h(2))];
130 c(2:n,idx) = A \ g;
131 endif
132
133 c(1,:) = c(n,:);
134 b = diff (a) ./ h(1:n - 1,idx)\
135 - h(1:n - 1,idx) / 3 .* (c(2:n,:) + 2 * c(1:n - 1,:));
136 b(n,:) = b(1,:);
137 d = diff (c) ./ (3 * h(1:n - 1, idx));
138 d(n,:) = d(1,:);
139
140 elseif (strcmp(cond,"not-a-knot"))
141
142 if (n > 4)
143
144 dg = 2 * (h(1:n - 2) .+ h(2:n - 1));
145 dg(1) = dg(1) - h(1);
146 dg(n - 2) = dg(n-2) - h(n - 1);
147
148 ldg = udg = h(2:n - 2);
149 udg(1) = udg(1) - h(1);
150 ldg(n - 3) = ldg(n-3) - h(n - 1);
151
152 elseif (n == 4)
153
154 dg = [h(1) + 2 * h(2), 2 * h(2) + h(3)];
155 ldg = h(2) - h(3);
156 udg = h(2) - h(1);
157
158 endif
159 g = zeros(n - 2,columns(y));
160 g(1,:) = 3 / (h(1) + h(2)) * (a(3,:) - a(2,:)\
161 - h(2) / h(1) * (a(2,:) - a(1,:)));
162 if (n > 4)
163 g(2:n - 3,:) = 3 * diff (a(3:n - 1,:)) ./ h(2:n - 3,idx)\
164 - 3 * diff (a(2:n - 2,:)) ./ h(1:n - 4,idx);
165 endif
166 g(n - 2,:) = 3 / (h(n - 1) + h(n - 2)) *\
167 (h(n - 2) / h(n - 1) * (a(n,:) - a(n - 1,:)) -\
168 (a(n - 1,:) - a(n - 2,:)));
169 c(2:n - 1,:) = trisolve(ldg,dg,udg,g);
170 c(1,:) = c(2,:) + h(1) / h(2) * (c(2,:) - c(3,:));
171 c(n,:) = c(n - 1,:) + h(n - 1) / h(n - 2) * (c(n - 1,:) - c(n - 2,:));
172 b = diff (a) ./ h(1:n - 1, idx)\
173 - h(1:n - 1, idx) / 3 .* (c(2:n,:) + 2 * c(1:n - 1,:));
174 d = diff (c) ./ (3 * h(1:n - 1, idx));
175
176 else
177 msg = sprintf("unknown end condition: %s",cond);
178 error (msg);
179 endif
180
181 d = d(1:n-1,:); c=c(1:n-1,:); b=b(1:n-1,:); a=a(1:n-1,:);
182 coeffs = [d(:), c(:), b(:), a(:)];
183 pp = mkpp (x, coeffs);
184
185 endfunction
186
187
188 %!shared x,y,cond
189 %! x = linspace(0,2*pi,15)'; y = sin(x);
190
191 %!assert (ppval(csape(x,y),x), y, 10*eps);
192 %!assert (ppval(csape(x,y),x'), y', 10*eps);
193 %!assert (ppval(csape(x',y'),x'), y', 10*eps);
194 %!assert (ppval(csape(x',y'),x), y, 10*eps);
195 %!assert (ppval(csape(x,[y,y]),x), \
196 %! [ppval(csape(x,y),x),ppval(csape(x,y),x)], 10*eps)
197
198 %!test cond='complete';
199 %!assert (ppval(csape(x,y,cond),x), y, 10*eps);
200 %!assert (ppval(csape(x,y,cond),x'), y', 10*eps);
201 %!assert (ppval(csape(x',y',cond),x'), y', 10*eps);
202 %!assert (ppval(csape(x',y',cond),x), y, 10*eps);
203 %!assert (ppval(csape(x,[y,y],cond),x), \
204 %! [ppval(csape(x,y,cond),x),ppval(csape(x,y,cond),x)], 10*eps)
205
206 %!test cond='variational';
207 %!assert (ppval(csape(x,y,cond),x), y, 10*eps);
208 %!assert (ppval(csape(x,y,cond),x'), y', 10*eps);
209 %!assert (ppval(csape(x',y',cond),x'), y', 10*eps);
210 %!assert (ppval(csape(x',y',cond),x), y, 10*eps);
211 %!assert (ppval(csape(x,[y,y],cond),x), \
212 %! [ppval(csape(x,y,cond),x),ppval(csape(x,y,cond),x)], 10*eps)
213
214 %!test cond='second';
215 %!assert (ppval(csape(x,y,cond),x), y, 10*eps);
216 %!assert (ppval(csape(x,y,cond),x'), y', 10*eps);
217 %!assert (ppval(csape(x',y',cond),x'), y', 10*eps);
218 %!assert (ppval(csape(x',y',cond),x), y, 10*eps);
219 %!assert (ppval(csape(x,[y,y],cond),x), \
220 %! [ppval(csape(x,y,cond),x),ppval(csape(x,y,cond),x)], 10*eps)
221
222 %!test cond='periodic';
223 %!assert (ppval(csape(x,y,cond),x), y, 10*eps);
224 %!assert (ppval(csape(x,y,cond),x'), y', 10*eps);
225 %!assert (ppval(csape(x',y',cond),x'), y', 10*eps);
226 %!assert (ppval(csape(x',y',cond),x), y, 10*eps);
227 %!assert (ppval(csape(x,[y,y],cond),x), \
228 %! [ppval(csape(x,y,cond),x),ppval(csape(x,y,cond),x)], 10*eps)
229
230 %!test cond='not-a-knot';
231 %!assert (ppval(csape(x,y,cond),x), y, 10*eps);
232 %!assert (ppval(csape(x,y,cond),x'), y', 10*eps);
233 %!assert (ppval(csape(x',y',cond),x'), y', 10*eps);
234 %!assert (ppval(csape(x',y',cond),x), y, 10*eps);
235 %!assert (ppval(csape(x,[y,y],cond),x), \
236 %! [ppval(csape(x,y,cond),x),ppval(csape(x,y,cond),x)], 10*eps)