octave_packages/splines-1.0.7/csape.m

   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, 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 ## the following end conditions as given in @var{cond} are possible.
  22 ## @table @asis
  23 ## @item 'complete'
  24 ##    match slopes at first and last point as given in @var{valc}
  25 ## @item 'not-a-knot'
  26 ##    third derivatives are continuous at the second and second last point
  27 ## @item 'periodic'
  28 ##    match first and second derivative of first and last point
  29 ## @item 'second'
  30 ##    match second derivative at first and last point as given in @var{valc}
  31 ## @item 'variational'
  32 ##    set second derivative at first and last point to zero (natural cubic spline)
  33 ## @end table
  34 ##
  35 ## @seealso{ppval, spline}
  36 ## @end deftypefn
  37
  38 ## Author:  Kai Habel <kai.habel@gmx.de>
  39 ## Date: 23. nov 2000
  40 ## Algorithms taken from G. Engeln-Muellges, F. Uhlig:
  41 ## "Numerical Algorithms with C", Springer, 1996
  42
  43 ## Paul Kienzle, 19. feb 2001,  csape supports now matrix y value
  44
  45 function pp = csape (x, y, cond, valc)
  46
  47   x = x(:);
  48   n = length(x);
  49   if (n < 3)
  50     error("csape requires at least 3 points");
  51   endif
  52
  53   ## Check the size and shape of y
  54   ndy = ndims (y);
  55   szy = size (y);
  56   if (ndy == 2 && (szy(1) == n || szy(2) == n))
  57     if (szy(2) == n)
  58       a = y.';
  59     else
  60       a = y;
  61       szy = fliplr (szy);
  62     endif
  63   else
  64     a = shiftdim (reshape (y, [prod(szy(1:end-1)), szy(end)]), 1);
  65   endif
  66
  67
  68   b = c = zeros (size (a));
  69   h = diff (x);
  70   idx = ones (columns(a),1);
  71
  72   if (nargin < 3 || strcmp(cond,"complete"))
  73     # specified first derivative at end point
  74     if (nargin < 4)
  75       valc = [0, 0];
  76     endif
  77
  78     if (n == 3)
  79       dg = 1.5 * h(1) - 0.5 * h(2);
  80       c(2:n - 1,:) = 1/dg(1);
  81     else
  82       dg = 2 * (h(1:n - 2) .+ h(2:n - 1));
  83       dg(1) = dg(1) - 0.5 * h(1);
  84       dg(n - 2) = dg(n-2) - 0.5 * h(n - 1);
  85
  86       e = h(2:n - 2);
  87
  88       g = 3 * diff (a(2:n,:)) ./ h(2:n - 1,idx)\
  89         - 3 * diff (a(1:n - 1,:)) ./ h(1:n - 2,idx);
  90       g(1,:) = 3 * (a(3,:) - a(2,:)) / h(2) \
  91           - 3 / 2 * (3 * (a(2,:) - a(1,:)) / h(1) - valc(1));
  92       g(n - 2,:) = 3 / 2 * (3 * (a(n,:) - a(n - 1,:)) / h(n - 1) - valc(2))\
  93           - 3 * (a(n - 1,:) - a(n - 2,:)) / h(n - 2);
  94
  95       c(2:n - 1,:) = spdiags([[e(:);0],dg,[0;e(:)]],[-1,0,1],n-2,n-2) \ g;
  96
  97     end
  98
  99     c(1,:) = (3 / h(1) * (a(2,:) - a(1,:)) - 3 * valc(1)
 100               - c(2,:) * h(1)) / (2 * h(1));
 101     c(n,:) = - (3 / h(n - 1) * (a(n,:) - a(n - 1,:)) - 3 * valc(2)
 102
 103                 + c(n - 1,:) * h(n - 1)) / (2 * h(n - 1));
 104     b(1:n - 1,:) = diff (a) ./ h(1:n - 1, idx)\
 105       - h(1:n - 1,idx) / 3 .* (c(2:n,:) + 2 * c(1:n - 1,:));
 106     d = diff (c) ./ (3 * h(1:n - 1, idx));
 107
 108   elseif (strcmp(cond,"variational") || strcmp(cond,"second"))
 109
 110     if ((nargin < 4) || strcmp(cond,"variational"))
 111       ## set second derivatives at end points to zero
 112       valc = [0, 0];
 113     endif
 114
 115     c(1,:) = valc(1) / 2;
 116     c(n,:) = valc(2) / 2;
 117
 118     g = 3 * diff (a(2:n,:)) ./ h(2:n - 1, idx)\
 119       - 3 * diff (a(1:n - 1,:)) ./ h(1:n - 2, idx);
 120
 121     g(1,:) = g(1,:) - h(1) * c(1,:);
 122     g(n - 2,:) = g(n-2,:) - h(n - 1) * c(n,:);
 123
 124     if( n == 3)
 125       dg = 2 * h(1);
 126       c(2:n - 1,:) = g / dg;
 127     else
 128       dg = 2 * (h(1:n - 2) .+ h(2:n - 1));
 129       e = h(2:n - 2);
 130       c(2:n - 1,:) = spdiags([[e(:);0],dg,[0;e(:)]],[-1,0,1],n-2,n-2) \ g;
 131     end
 132
 133     b(1:n - 1,:) = diff (a) ./ h(1:n - 1,idx)\
 134       - h(1:n - 1,idx) / 3 .* (c(2:n,:) + 2 * c(1:n - 1,:));
 135     d = diff (c) ./ (3 * h(1:n - 1, idx));
 136
 137   elseif (strcmp(cond,"periodic"))
 138
 139     h = [h; h(1)];
 140
 141     ## XXX FIXME XXX --- the following gives a smoother periodic transition:
 142     ##    a(n,:) = a(1,:) = ( a(n,:) + a(1,:) ) / 2;
 143     a(n,:) = a(1,:);
 144
 145     tmp = diff (shift ([a; a(2,:)], -1));
 146     g = 3 * tmp(1:n - 1,:) ./ h(2:n,idx)\
 147       - 3 * diff (a) ./ h(1:n - 1,idx);
 148
 149     if (n > 3)
 150       dg = 2 * (h(1:n - 1) .+ h(2:n));
 151       e = h(2:n - 1);
 152
 153       ## Use Sherman-Morrison formula to extend the solution
 154       ## to the cyclic system. See Numerical Recipes in C, pp 73-75
 155       gamma = - dg(1);
 156       dg(1) -=  gamma;
 157       dg(end) -= h(1) * h(1) / gamma;
 158       z = spdiags([[e(:);0],dg,[0;e(:)]],[-1,0,1],n-1,n-1) \ ...
 159           [[gamma; zeros(n-3,1); h(1)],g];
 160       fact = (z(1,2:end) + h(1) * z(end,2:end) / gamma) / ...
 161           (1.0 + z(1,1) + h(1) * z(end,1) / gamma);
 162
 163       c(2:n,idx) = z(:,2:end) - z(:,1) * fact;
 164     endif
 165
 166     c(1,:) = c(n,:);
 167     b = diff (a) ./ h(1:n - 1,idx)\
 168       - h(1:n - 1,idx) / 3 .* (c(2:n,:) + 2 * c(1:n - 1,:));
 169     b(n,:) = b(1,:);
 170     d = diff (c) ./ (3 * h(1:n - 1, idx));
 171     d(n,:) = d(1,:);
 172
 173   elseif (strcmp(cond,"not-a-knot"))
 174
 175     g = zeros(n - 2,columns(a));
 176     g(1,:) = 3 / (h(1) + h(2)) * (a(3,:) - a(2,:)\
 177           - h(2) / h(1) * (a(2,:) - a(1,:)));
 178     g(n - 2,:) = 3 / (h(n - 1) + h(n - 2)) *\
 179         (h(n - 2) / h(n - 1) * (a(n,:) - a(n - 1,:)) -\
 180          (a(n - 1,:) - a(n - 2,:)));
 181
 182     if (n > 4)
 183
 184       g(2:n - 3,:) = 3 * diff (a(3:n - 1,:)) ./ h(3:n - 2,idx)\
 185         - 3 * diff (a(2:n - 2,:)) ./ h(2:n - 3,idx);
 186
 187       dg = 2 * (h(1:n - 2) .+ h(2:n - 1));
 188       dg(1) = dg(1) - h(1);
 189       dg(n - 2) = dg(n-2) - h(n - 1);
 190
 191       ldg = udg = h(2:n - 2);
 192       udg(1) = udg(1) - h(1);
 193       ldg(n - 3) = ldg(n-3) - h(n - 1);
 194       c(2:n - 1,:) = spdiags([[ldg(:);0],dg,[0;udg(:)]],[-1,0,1],n-2,n-2) \ g;
 195
 196     elseif (n == 4)
 197
 198       dg = [h(1) + 2 * h(2), 2 * h(2) + h(3)];
 199       ldg = h(2) - h(3);
 200       udg = h(2) - h(1);
 201       c(2:n - 1,:) = spdiags([[ldg(:);0],dg,[0;udg(:)]],[-1,0,1],n-2,n-2) \ g;
 202
 203     else # n == 3
 204
 205       dg= [h(1) + 2 * h(2)];
 206       c(2:n - 1,:) = g/dg(1);
 207
 208     endif
 209
 210     c(1,:) = c(2,:) + h(1) / h(2) * (c(2,:) - c(3,:));
 211     c(n,:) = c(n - 1,:) + h(n - 1) / h(n - 2) * (c(n - 1,:) - c(n - 2,:));
 212     b = diff (a) ./ h(1:n - 1, idx)\
 213       - h(1:n - 1, idx) / 3 .* (c(2:n,:) + 2 * c(1:n - 1,:));
 214     d = diff (c) ./ (3 * h(1:n - 1, idx));
 215
 216   else
 217     msg = sprintf("unknown end condition: %s",cond);
 218     error (msg);
 219   endif
 220
 221   d = d(1:n-1,:); c=c(1:n-1,:); b=b(1:n-1,:); a=a(1:n-1,:);
 222   pp = mkpp (x, cat (2, d'(:), c'(:), b'(:), a'(:)), szy(1:end-1));
 223
 224 endfunction
 225
 226
 227 %!shared x,y,cond
 228 %! x = linspace(0,2*pi,15); y = sin(x);
 229
 230 %!assert (ppval(csape(x,y),x), y, 10*eps);
 231 %!assert (ppval(csape(x,y),x'), y', 10*eps);
 232 %!assert (ppval(csape(x',y'),x'), y', 10*eps);
 233 %!assert (ppval(csape(x',y'),x), y, 10*eps);
 234 %!assert (ppval(csape(x,[y;y]),x), \
 235 %!        [ppval(csape(x,y),x);ppval(csape(x,y),x)], 10*eps)
 236
 237 %!test cond='complete';
 238 %!assert (ppval(csape(x,y,cond),x), y, 10*eps);
 239 %!assert (ppval(csape(x,y,cond),x'), y', 10*eps);
 240 %!assert (ppval(csape(x',y',cond),x'), y', 10*eps);
 241 %!assert (ppval(csape(x',y',cond),x), y, 10*eps);
 242 %!assert (ppval(csape(x,[y;y],cond),x), \
 243 %!        [ppval(csape(x,y,cond),x);ppval(csape(x,y,cond),x)], 10*eps)
 244
 245 %!test cond='variational';
 246 %!assert (ppval(csape(x,y,cond),x), y, 10*eps);
 247 %!assert (ppval(csape(x,y,cond),x'), y', 10*eps);
 248 %!assert (ppval(csape(x',y',cond),x'), y', 10*eps);
 249 %!assert (ppval(csape(x',y',cond),x), y, 10*eps);
 250 %!assert (ppval(csape(x,[y;y],cond),x), \
 251 %!        [ppval(csape(x,y,cond),x);ppval(csape(x,y,cond),x)], 10*eps)
 252
 253 %!test cond='second';
 254 %!assert (ppval(csape(x,y,cond),x), y, 10*eps);
 255 %!assert (ppval(csape(x,y,cond),x'), y', 10*eps);
 256 %!assert (ppval(csape(x',y',cond),x'), y', 10*eps);
 257 %!assert (ppval(csape(x',y',cond),x), y, 10*eps);
 258 %!assert (ppval(csape(x,[y;y],cond),x), \
 259 %!        [ppval(csape(x,y,cond),x);ppval(csape(x,y,cond),x)], 10*eps)
 260
 261 %!test cond='periodic';
 262 %!assert (ppval(csape(x,y,cond),x), y, 10*eps);
 263 %!assert (ppval(csape(x,y,cond),x'), y', 10*eps);
 264 %!assert (ppval(csape(x',y',cond),x'), y', 10*eps);
 265 %!assert (ppval(csape(x',y',cond),x), y, 10*eps);
 266 %!assert (ppval(csape(x,[y;y],cond),x), \
 267 %!        [ppval(csape(x,y,cond),x);ppval(csape(x,y,cond),x)], 10*eps)
 268
 269 %!test cond='not-a-knot';
 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',cond),x'), y', 10*eps);
 273 %!assert (ppval(csape(x',y',cond),x), y, 10*eps);
 274 %!assert (ppval(csape(x,[y;y],cond),x), \
 275 %!        [ppval(csape(x,y,cond),x);ppval(csape(x,y,cond),x)], 10*eps)