1## Copyright (C) 2000, 2001 Kai Habel 2## 3## This program is free software; you can redistribute it and/or modify it under 4## the terms of the GNU General Public License as published by the Free Software 5## Foundation; either version 3 of the License, or (at your option) any later 6## version. 7## 8## This program is distributed in the hope that it will be useful, but WITHOUT 9## ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or 10## FITNESS FOR A PARTICULAR PURPOSE. See the GNU General Public License for more 11## details. 12## 13## You should have received a copy of the GNU General Public License along with 14## this program; if not, see <http://www.gnu.org/licenses/>. 15 16## -*- texinfo -*- 17## @deftypefn {Function File} {@var{pp} = } csape (@var{x}, @var{y}, @var{cond}, @var{valc}) 18## cubic spline interpolation with various end conditions. 19## creates the pp-form of the cubic spline. 20## 21## @var{x} should be @var{n} by 1, @var{y} should be @var{n} by @var{m}, 22## @var{valc} should be 2 by @var{m} or 2 by 1 23## 24## The following end conditions as given in @var{cond} are possible: 25## @table @asis 26## @item 'complete' 27## match slopes at first and last point as given in @var{valc} 28## (default; if @var{valc} is not given, the slopes matched are those 29## of the cubic polynomials that interpolate the first and last four points) 30## @item 'not-a-knot' 31## third derivatives are continuous at the second and second last point 32## @item 'periodic' 33## match first and second derivative of first and last point 34## @item 'second' 35## match second derivative at first and last point as given in @var{valc} 36## @item 'variational' 37## set second derivative at first and last point to zero (natural cubic spline) 38## @end table 39## 40## @seealso{ppval, spline} 41## @end deftypefn 42 43## Author: Kai Habel <kai.habel@gmx.de> 44## Date: 23 Nov 2000 45## Algorithms taken from G. Engeln-Muellges, F. Uhlig: 46## "Numerical Algorithms with C", Springer, 1996 47 48function pp = csape (x, y, cond, valc) 49 50 x = x(:); 51 n = length(x); 52 if (n < 3) 53 error ("csape requires at least 3 points"); 54 endif 55 56 ## Check the size and shape of y 57 ndy = ndims (y); 58 szy = size (y); 59 if (ndy == 2 && (szy(1) == n || szy(2) == n)) 60 if (szy(2) == n) 61 a = y.'; 62 else 63 a = y; 64 szy = fliplr (szy); 65 endif 66 else 67 a = shiftdim (reshape (y, [prod(szy(1:end-1)), szy(end)]), 1); 68 endif 69 m = size (a, 2); 70 71 if exist("valc", "var") && size(valc) == [1 2] 72 valc = valc'; 73 endif 74 75 b = c = zeros (n, m); 76 h = diff (x); 77 idx = ones (columns(a),1); 78 79 if (nargin < 3 || strcmp(cond,"complete")) 80 81 # set first derivative at end points 82 if (nargin < 4) 83 valc = zeros(2, m); 84 n_use = min(n, 4); 85 for i = 1:m 86 valc(1, i) = polyval(polyder(polyfit(x(1:n_use), a(1:n_use, i), n_use-1)), x(1)); 87 valc(2, i) = polyval(polyder(polyfit(x(end-n_use+1:end), a(end-n_use+1:end, i), n_use-1)), x(end)); 88 endfor 89 endif 90 91 if (n == 3) 92 warning ("off", "Octave:broadcast", "local"); 93 e = 2 * [h(1) h(1:(end-1))+h(2:end) h(end)]; 94 A = spdiags([[h; 0], e(:), [0; h]], [-1,0,1], n, n); 95 d = diff(a) ./ h; #uses broadcasting if columns(a) > 1 96 g = 3 * diff(d); 97 A(1, 1) = 2 * h(1); 98 A(1, 2) = h(1); 99 A(n, n) = 2 * h(end); 100 A(end, end-1) = h(end); 101 g = [3*(d(1, :) - valc(1, :)); g; 3*(valc(2, :) - d(end, :))]; 102 c = A \ g; 103 else 104 dg = 2 * (h(1:n - 2) .+ h(2:n - 1)); 105 dg(1) = dg(1) - 0.5 * h(1); 106 dg(n - 2) = dg(n-2) - 0.5 * h(n - 1); 107 108 e = h(2:n - 2); 109 110 g = 3 * diff (a(2:n,:)) ./ h(2:n - 1,idx)... 111 - 3 * diff (a(1:n - 1,:)) ./ h(1:n - 2,idx); 112 g(1,:) = 3 * (a(3,:) - a(2,:)) / h(2) ... 113 - 3 / 2 * (3 * (a(2,:) - a(1,:)) / h(1) - valc(1,:)); 114 g(n - 2,:) = 3 / 2 * (3 * (a(n,:) - a(n - 1,:)) / h(n - 1) - valc(2,:))... 115 - 3 * (a(n - 1,:) - a(n - 2,:)) / h(n - 2); 116 117 c(2:n - 1,:) = spdiags([[e(:);0],dg,[0;e(:)]],[-1,0,1],n-2,n-2) \ g; 118 119 c(1,:) = (3 / h(1) * (a(2,:) - a(1,:)) - 3 * valc(1, :) - c(2,:) * h(1)) / (2 * h(1)); 120 c(n,:) = - (3 / h(n - 1) * (a(n,:) - a(n - 1,:)) - 3 * valc(2, :) + c(n - 1,:) * h(n - 1)) / (2 * h(n - 1)); 121 end 122 b(1:n - 1,:) = diff (a) ./ h(1:n - 1, idx)... 123 - h(1:n - 1,idx) / 3 .* (c(2:n,:) + 2 * c(1:n - 1,:)); 124 d = diff (c) ./ (3 * h(1:n - 1, idx)); 125 126 elseif (strcmp(cond,"variational") || strcmp(cond,"second")) 127 128 if ((nargin < 4) || strcmp(cond,"variational")) 129 ## set second derivatives at end points to zero 130 valc = zeros (2, 1); 131 endif 132 133 c(1,:) = valc(1, :) / 2; 134 c(n,:) = valc(2, :) / 2; 135 136 g = 3 * diff (a(2:n,:)) ./ h(2:n - 1, idx)... 137 - 3 * diff (a(1:n - 1,:)) ./ h(1:n - 2, idx); 138 139 g(1,:) = g(1,:) - h(1) * c(1,:); 140 g(n - 2,:) = g(n-2,:) - h(n - 1) * c(n,:); 141 142 if( n == 3) 143 dg = 2 * h(1); 144 c(2:n - 1,:) = g / dg; 145 else 146 dg = 2 * (h(1:n - 2) .+ h(2:n - 1)); 147 e = h(2:n - 2); 148 c(2:n - 1,:) = spdiags([[e(:);0],dg,[0;e(:)]],[-1,0,1],n-2,n-2) \ g; 149 end 150 151 b(1:n - 1,:) = diff (a) ./ h(1:n - 1,idx)... 152 - h(1:n - 1,idx) / 3 .* (c(2:n,:) + 2 * c(1:n - 1,:)); 153 d = diff (c) ./ (3 * h(1:n - 1, idx)); 154 155 elseif (strcmp(cond,"periodic")) 156 157 h = [h; h(1)]; 158 159 ## XXX FIXME XXX --- the following gives a smoother periodic transition: 160 ## a(n,:) = a(1,:) = ( a(n,:) + a(1,:) ) / 2; 161 a(n,:) = a(1,:); 162 163 tmp = diff (shift ([a; a(2,:)], -1)); 164 g = 3 * tmp(1:n - 1,:) ./ h(2:n,idx)... 165 - 3 * diff (a) ./ h(1:n - 1,idx); 166 167 if (n > 3) 168 dg = 2 * (h(1:n - 1) .+ h(2:n)); 169 e = h(2:n - 1); 170 171 ## Use Sherman-Morrison formula to extend the solution 172 ## to the cyclic system. See Numerical Recipes in C, pp 73-75 173 gamma = - dg(1); 174 dg(1) -= gamma; 175 dg(end) -= h(1) * h(1) / gamma; 176 z = spdiags([[e(:);0],dg,[0;e(:)]],[-1,0,1],n-1,n-1) \ ... 177 [[gamma; zeros(n-3,1); h(1)],g]; 178 fact = (z(1,2:end) + h(1) * z(end,2:end) / gamma) / ... 179 (1.0 + z(1,1) + h(1) * z(end,1) / gamma); 180 181 c(2:n,:) = z(:,2:end) - z(:,1) * fact; 182 endif 183 184 c(1,:) = c(n,:); 185 b = diff (a) ./ h(1:n - 1,idx)... 186 - h(1:n - 1,idx) / 3 .* (c(2:n,:) + 2 * c(1:n - 1,:)); 187 b(n,:) = b(1,:); 188 d = diff (c) ./ (3 * h(1:n - 1, idx)); 189 d(n,:) = d(1,:); 190 191 elseif (strcmp(cond,"not-a-knot")) 192 193 g = zeros(n - 2,columns(a)); 194 g(1,:) = 3 / (h(1) + h(2)) * (a(3,:) - a(2,:)... 195 - h(2) / h(1) * (a(2,:) - a(1,:))); 196 g(n - 2,:) = 3 / (h(n - 1) + h(n - 2)) *... 197 (h(n - 2) / h(n - 1) * (a(n,:) - a(n - 1,:)) -... 198 (a(n - 1,:) - a(n - 2,:))); 199 200 if (n > 4) 201 202 g(2:n - 3,:) = 3 * diff (a(3:n - 1,:)) ./ h(3:n - 2,idx)... 203 - 3 * diff (a(2:n - 2,:)) ./ h(2:n - 3,idx); 204 205 dg = 2 * (h(1:n - 2) .+ h(2:n - 1)); 206 dg(1) = dg(1) - h(1); 207 dg(n - 2) = dg(n-2) - h(n - 1); 208 209 ldg = udg = h(2:n - 2); 210 udg(1) = udg(1) - h(1); 211 ldg(n - 3) = ldg(n-3) - h(n - 1); 212 c(2:n - 1,:) = spdiags([[ldg(:);0],dg,[0;udg(:)]],[-1,0,1],n-2,n-2) \ g; 213 214 elseif (n == 4) 215 216 dg = [h(1) + 2 * h(2); 2 * h(2) + h(3)]; 217 ldg = h(2) - h(3); 218 udg = h(2) - h(1); 219 c(2:n - 1,:) = spdiags([[ldg(:);0],dg,[0;udg(:)]],[-1,0,1],n-2,n-2) \ g; 220 221 else # n == 3 222 223 dg= [h(1) + 2 * h(2)]; 224 c(2:n - 1,:) = g/dg(1); 225 226 endif 227 228 c(1,:) = c(2,:) + h(1) / h(2) * (c(2,:) - c(3,:)); 229 c(n,:) = c(n - 1,:) + h(n - 1) / h(n - 2) * (c(n - 1,:) - c(n - 2,:)); 230 b = diff (a) ./ h(1:n - 1, idx)... 231 - h(1:n - 1, idx) / 3 .* (c(2:n,:) + 2 * c(1:n - 1,:)); 232 d = diff (c) ./ (3 * h(1:n - 1, idx)); 233 234 else 235 msg = sprintf("unknown end condition: %s",cond); 236 error (msg); 237 endif 238 239 d = d(1:n-1,:); c=c(1:n-1,:); b=b(1:n-1,:); a=a(1:n-1,:); 240 pp = mkpp (x, cat (2, d'(:), c'(:), b'(:), a'(:)), szy(1:end-1)); 241 242endfunction 243 244 245%!shared x,x2,y,cond,pp,pp1,pp2,h,valc 246%! x = linspace(0,2*pi,5); y = sin(x); x2 = linspace(0,2*pi,16); 247 248%!assert (ppval(csape(x,y),x), y, 10*eps); 249%!assert (ppval(csape(x,y),x'), y', 10*eps); 250%!assert (ppval(csape(x',y'),x'), y', 10*eps); 251%!assert (ppval(csape(x',y'),x), y, 10*eps); 252%!assert (ppval(csape(x,[y;y]),x), ... 253%! [ppval(csape(x,y),x);ppval(csape(x,y),x)], 10*eps) 254%!assert (ppval(csape(x,[y;y]),x2), ... 255%! [ppval(csape(x,y),x2);ppval(csape(x,y),x2)], 10*eps) 256 257%!test cond='complete'; 258%!assert (ppval(csape(x,y,cond),x), y, 10*eps); 259%!assert (ppval(csape(x,y,cond),x'), y', 10*eps); 260%!assert (ppval(csape(x',y',cond),x'), y', 10*eps); 261%!assert (ppval(csape(x',y',cond),x), y, 10*eps); 262%!assert (ppval(csape(x,[y;y],cond),x), ... 263%! [ppval(csape(x,y,cond),x);ppval(csape(x,y,cond),x)], 10*eps) 264%!assert (ppval(csape(x,[y;y],cond),x2), ... 265%! [ppval(csape(x,y,cond),x2);ppval(csape(x,y,cond),x2)], 10*eps) 266 267%!test cond='variational'; 268%!assert (ppval(csape(x,y,cond),x), y, 10*eps); 269%!assert (ppval(csape(x,y,cond),x'), y', 10*eps); 270%!assert (ppval(csape(x',y',cond),x'), y', 10*eps); 271%!assert (ppval(csape(x',y',cond),x), y, 10*eps); 272%!assert (ppval(csape(x,[y;y],cond),x), ... 273%! [ppval(csape(x,y,cond),x);ppval(csape(x,y,cond),x)], 10*eps) 274%!assert (ppval(csape(x,[y;y],cond),x2), ... 275%! [ppval(csape(x,y,cond),x2);ppval(csape(x,y,cond),x2)], 10*eps) 276 277%!test cond='second'; 278%!assert (ppval(csape(x,y,cond),x), y, 10*eps); 279%!assert (ppval(csape(x,y,cond),x'), y', 10*eps); 280%!assert (ppval(csape(x',y',cond),x'), y', 10*eps); 281%!assert (ppval(csape(x',y',cond),x), y, 10*eps); 282%!assert (ppval(csape(x,[y;y],cond),x), ... 283%! [ppval(csape(x,y,cond),x);ppval(csape(x,y,cond),x)], 10*eps) 284%!assert (ppval(csape(x,[y;y],cond),x2), ... 285%! [ppval(csape(x,y,cond),x2);ppval(csape(x,y,cond),x2)], 10*eps) 286 287%!test cond='periodic'; 288%!assert (ppval(csape(x,y,cond),x), y, 10*eps); 289%!assert (ppval(csape(x,y,cond),x'), y', 10*eps); 290%!assert (ppval(csape(x',y',cond),x'), y', 10*eps); 291%!assert (ppval(csape(x',y',cond),x), y, 10*eps); 292%!assert (ppval(csape(x,[y;y],cond),x), ... 293%! [ppval(csape(x,y,cond),x);ppval(csape(x,y,cond),x)], 10*eps) 294%!assert (ppval(csape(x,[y;y],cond),x2), ... 295%! [ppval(csape(x,y,cond),x2);ppval(csape(x,y,cond),x2)], 10*eps) 296 297%!test cond='not-a-knot'; 298%!assert (ppval(csape(x,y,cond),x), y, 10*eps); 299%!assert (ppval(csape(x,y,cond),x'), y', 10*eps); 300%!assert (ppval(csape(x',y',cond),x'), y', 10*eps); 301%!assert (ppval(csape(x',y',cond),x), y, 10*eps); 302%!assert (ppval(csape(x,[y;y],cond),x), ... 303%! [ppval(csape(x,y,cond),x);ppval(csape(x,y,cond),x)], 10*eps) 304%!assert (ppval(csape(x,[y;y],cond),x2), ... 305%! [ppval(csape(x,y,cond),x2);ppval(csape(x,y,cond),x2)], 10*eps) 306%!assert (ppval(csape(x(1:4),y(1:4),cond),x(1:4)), y(1:4), 10*eps); 307%!assert (ppval(csape(x(1:4)',y(1:4)',cond),x(1:4)), y(1:4), 10*eps); 308 309%!test cond='complete'; 310%! h = pi/6; n = 3; x = linspace(0,(n-1)*h,n)'; y = sin(x); valc = cos([x(1) x(end)]); pp = csape(x, y, cond, valc); 311%!assert (ppval(csape(x,[y y],cond, valc),x)', ... 312%! repmat(ppval(pp, x), [1 2]), 10*eps) 313%!assert (polyval(pp.coefs(1, :), [0 h])', y(1:2), 10*eps) 314%!assert (polyval(pp.coefs(2, :), [0 h])', y(2:3), 10*eps) 315%!assert (polyval([3*pp.coefs(1, 1) 2*pp.coefs(1, 2) pp.coefs(1, 3)], 0), valc(1), 10*eps) 316%!assert (polyval([3*pp.coefs(2, 1) 2*pp.coefs(2, 2) pp.coefs(2, 3)], h), valc(2), 10*eps) 317%!assert (polyval([3*pp.coefs(1, 1) 2*pp.coefs(1, 2) pp.coefs(1, 3)], h), polyval([3*pp.coefs(2, 1) 2*pp.coefs(2, 2) pp.coefs(2, 3)], 0), 10*eps) 318%!assert (polyval([6*pp.coefs(1, 1) 2*pp.coefs(1, 2)], h), polyval([6*pp.coefs(2, 1) 2*pp.coefs(2, 2)], 0), 10*eps) 319%! y = cos(x); valc = -sin([x(1) x(end)]); pp1 = csape(x, y, cond, valc); 320%!assert (ppval(csape(x,[y y],cond, valc),x)', ... 321%! repmat(ppval(pp1, x), [1 2]), 10*eps) 322%!assert (polyval(pp1.coefs(1, :), [0 h])', y(1:2), 10*eps) 323%!assert (polyval(pp1.coefs(2, :), [0 h])', y(2:3), 10*eps) 324%!assert (polyval([3*pp1.coefs(1, 1) 2*pp1.coefs(1, 2) pp1.coefs(1, 3)], 0), valc(1), 10*eps) 325%!assert (polyval([3*pp1.coefs(2, 1) 2*pp1.coefs(2, 2) pp1.coefs(2, 3)], h), valc(2), 10*eps) 326%!assert (polyval([3*pp1.coefs(1, 1) 2*pp1.coefs(1, 2) pp1.coefs(1, 3)], h), polyval([3*pp1.coefs(2, 1) 2*pp1.coefs(2, 2) pp1.coefs(2, 3)], 0), 10*eps) 327%!assert (polyval([6*pp1.coefs(1, 1) 2*pp1.coefs(1, 2)], h), polyval([6*pp1.coefs(2, 1) 2*pp1.coefs(2, 2)], 0), 10*eps) 328%! y = [sin(x) cos(x)]; valc = [cos([x(1); x(end)]) -sin([x(1); x(end)])]; pp2 = csape(x, y, cond, valc); 329%!assert (pp2.coefs([1 3], :), pp.coefs) 330%!assert (pp2.coefs([2 4], :), pp1.coefs) 331