octave_packages/m/polynomial/spline.m

   1 ## Copyright (C) 2000-2012 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 3 of the License, or (at
   9 ## your option) 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, see
  18 ## <http://www.gnu.org/licenses/>.
  19
  20 ## -*- texinfo -*-
  21 ## @deftypefn  {Function File} {@var{pp} =} spline (@var{x}, @var{y})
  22 ## @deftypefnx {Function File} {@var{yi} =} spline (@var{x}, @var{y}, @var{xi})
  23 ## Return the cubic spline interpolant of points @var{x} and @var{y}.
  24 ##
  25 ## When called with two arguments, return the piecewise polynomial @var{pp}
  26 ## that may be used with @code{ppval} to evaluate the polynomial at specific
  27 ## points.  When called with a third input argument, @code{spline} evaluates
  28 ## the spline at the points @var{xi}.  The third calling form @code{spline
  29 ## (@var{x}, @var{y}, @var{xi})} is equivalent to @code{ppval (spline
  30 ## (@var{x}, @var{y}), @var{xi})}.
  31 ##
  32 ## The variable @var{x} must be a vector of length @var{n}.  @var{y} can be
  33 ## either a vector or array.  If @var{y} is a vector it must have a length of
  34 ## either @var{n} or @code{@var{n} + 2}.  If the length of @var{y} is
  35 ## @var{n}, then the "not-a-knot" end condition is used.  If the length of
  36 ## @var{y} is @code{@var{n} + 2}, then the first and last values of the
  37 ## vector @var{y} are the values of the first derivative of the cubic spline
  38 ## at the endpoints.
  39 ##
  40 ## If @var{y} is an array, then the size of @var{y} must have the form
  41 ## @tex
  42 ## $$[s_1, s_2, \cdots, s_k, n]$$
  43 ## @end tex
  44 ## @ifnottex
  45 ## @code{[@var{s1}, @var{s2}, @dots{}, @var{sk}, @var{n}]}
  46 ## @end ifnottex
  47 ## or
  48 ## @tex
  49 ## $$[s_1, s_2, \cdots, s_k, n + 2].$$
  50 ## @end tex
  51 ## @ifnottex
  52 ## @code{[@var{s1}, @var{s2}, @dots{}, @var{sk}, @var{n} + 2]}.
  53 ## @end ifnottex
  54 ## The array is reshaped internally to a matrix where the leading
  55 ## dimension is given by
  56 ## @tex
  57 ## $$s_1 s_2 \cdots s_k$$
  58 ## @end tex
  59 ## @ifnottex
  60 ## @code{@var{s1} * @var{s2} * @dots{} * @var{sk}}
  61 ## @end ifnottex
  62 ## and each row of this matrix is then treated separately.  Note that this
  63 ## is exactly opposite to @code{interp1} but is done for @sc{matlab}
  64 ## compatibility.
  65 ##
  66 ## @seealso{pchip, ppval, mkpp, unmkpp}
  67 ## @end deftypefn
  68
  69 ## This code is based on csape.m from octave-forge, but has been
  70 ## modified to use the sparse solver code in octave that itself allows
  71 ## special casing of tri-diagonal matrices, modified for NDArrays and
  72 ## for the treatment of vectors y 2 elements longer than x as complete
  73 ## splines.
  74
  75 function ret = spline (x, y, xi)
  76
  77   x = x(:);
  78   n = length (x);
  79   if (n < 2)
  80     error ("spline: requires at least 2 points");
  81   endif
  82
  83   ## Check the size and shape of y
  84   ndy = ndims (y);
  85   szy = size (y);
  86   if (ndy == 2 && (szy(1) == n || szy(2) == n))
  87     if (szy(2) == n)
  88       a = y.';
  89     else
  90       a = y;
  91       szy = fliplr (szy);
  92     endif
  93   else
  94     a = shiftdim (reshape (y, [prod(szy(1:end-1)), szy(end)]), 1);
  95   endif
  96
  97   for k = (1:columns (a))(any (isnan (a)))
  98     ok = ! isnan (a(:,k));
  99     a(!ok,k) = spline (x(ok), a(ok,k), x(!ok));
 100   endfor
 101
 102   complete = false;
 103   if (size (a, 1) == n + 2)
 104     complete = true;
 105     dfs = a(1,:);
 106     dfe = a(end,:);
 107     a = a(2:end-1,:);
 108   endif
 109
 110   if (~issorted (x))
 111     [x, idx] = sort(x);
 112     a = a(idx,:);
 113   endif
 114
 115   b = c = zeros (size (a));
 116   h = diff (x);
 117   idx = ones (columns (a), 1);
 118
 119   if (complete)
 120
 121     if (n == 2)
 122       d = (dfs + dfe) / (x(2) - x(1)) ^ 2 + ...
 123           2 * (a(1,:) - a(2,:)) / (x(2) - x(1)) ^ 3;
 124       c = (-2 * dfs - dfe) / (x(2) - x(1)) - ...
 125           3 * (a(1,:) - a(2,:)) / (x(2) - x(1)) ^ 2;
 126       b = dfs;
 127       a = a(1,:);
 128
 129       d = d(1:n-1,:);
 130       c = c(1:n-1,:);
 131       b = b(1:n-1,:);
 132       a = a(1:n-1,:);
 133     else
 134       if (n == 3)
 135         dg = 1.5 * h(1) - 0.5 * h(2);
 136         c(2:n-1,:) = 1/dg(1);
 137       else
 138         dg = 2 * (h(1:n-2) .+ h(2:n-1));
 139         dg(1) = dg(1) - 0.5 * h(1);
 140         dg(n-2) = dg(n-2) - 0.5 * h(n-1);
 141
 142         e = h(2:n-2);
 143
 144         g = 3 * diff (a(2:n,:)) ./ h(2:n-1,idx) ...
 145           - 3 * diff (a(1:n-1,:)) ./ h(1:n-2,idx);
 146         g(1,:) = 3 * (a(3,:) - a(2,:)) / h(2) ...
 147           - 3 / 2 * (3 * (a(2,:) - a(1,:)) / h(1) - dfs);
 148         g(n-2,:) = 3 / 2 * (3 * (a(n,:) - a(n-1,:)) / h(n-1) - dfe) ...
 149           - 3 * (a(n-1,:) - a(n-2,:)) / h(n-2);
 150
 151         c(2:n-1,:) = spdiags ([[e(:); 0], dg, [0; e(:)]],
 152                               [-1, 0, 1], n-2, n-2) \ g;
 153       endif
 154
 155       c(1,:) = (3 / h(1) * (a(2,:) - a(1,:)) - 3 * dfs
 156              - c(2,:) * h(1)) / (2 * h(1));
 157       c(n,:) = - (3 / h(n-1) * (a(n,:) - a(n-1,:)) - 3 * dfe
 158              + c(n-1,:) * h(n-1)) / (2 * h(n-1));
 159       b(1:n-1,:) = diff (a) ./ h(1:n-1, idx) ...
 160         - h(1:n-1,idx) / 3 .* (c(2:n,:) + 2 * c(1:n-1,:));
 161       d = diff (c) ./ (3 * h(1:n-1, idx));
 162
 163       d = d(1:n-1,:);
 164       c = c(1:n-1,:);
 165       b = b(1:n-1,:);
 166       a = a(1:n-1,:);
 167     endif
 168   else
 169
 170     if (n == 2)
 171       b = (a(2,:) - a(1,:)) / (x(2) - x(1));
 172       a = a(1,:);
 173       d = [];
 174       c = [];
 175       b = b(1:n-1,:);
 176       a = a(1:n-1,:);
 177     elseif (n == 3)
 178
 179       n = 2;
 180       c = (a(1,:) - a(3,:)) / ((x(3) - x(1)) * (x(2) - x(3))) ...
 181           + (a(2,:) - a(1,:)) / ((x(2) - x(1)) * (x(2) - x(3)));
 182       b = (a(2,:) - a(1,:)) * (x(3) - x(1)) ...
 183           / ((x(2) - x(1)) * (x(3) - x(2))) ...
 184           + (a(1,:) - a(3,:)) * (x(2) - x(1)) ...
 185           / ((x(3) - x(1)) * (x(3) - x(2)));
 186       a = a(1,:);
 187       d = [];
 188       x = [min(x), max(x)];
 189
 190       c = c(1:n-1,:);
 191       b = b(1:n-1,:);
 192       a = a(1:n-1,:);
 193     else
 194
 195       g = zeros (n-2, columns (a));
 196       g(1,:) = 3 / (h(1) + h(2)) ...
 197           * (a(3,:) - a(2,:) - h(2) / h(1) * (a(2,:) - a(1,:)));
 198       g(n-2,:) = 3 / (h(n-1) + h(n-2)) ...
 199           * (h(n-2) / h(n-1) * (a(n,:) - a(n-1,:)) - (a(n-1,:) - a(n-2,:)));
 200
 201       if (n > 4)
 202
 203         g(2:n - 3,:) = 3 * diff (a(3:n-1,:)) ./ h(3:n-2,idx) ...
 204             - 3 * diff (a(2:n-2,:)) ./ h(2:n - 3,idx);
 205
 206         dg = 2 * (h(1:n-2) .+ h(2:n-1));
 207         dg(1) = dg(1) - h(1);
 208         dg(n-2) = dg(n-2) - h(n-1);
 209
 210         ldg = udg = h(2:n-2);
 211         udg(1) = udg(1) - h(1);
 212         ldg(n - 3) = ldg(n-3) - h(n-1);
 213         c(2:n-1,:) = spdiags ([[ldg(:); 0], dg, [0; udg(:)]],
 214                               [-1, 0, 1], n-2, n-2) \ g;
 215
 216       elseif (n == 4)
 217
 218         dg = [h(1) + 2 * h(2); 2 * h(2) + h(3)];
 219         ldg = h(2) - h(3);
 220         udg = h(2) - h(1);
 221         c(2:n-1,:) = spdiags ([[ldg(:);0], dg, [0; udg(:)]],
 222                               [-1, 0, 1], n-2, n-2) \ g;
 223
 224       endif
 225
 226       c(1,:) = c(2,:) + h(1) / h(2) * (c(2,:) - c(3,:));
 227       c(n,:) = c(n-1,:) + h(n-1) / h(n-2) * (c(n-1,:) - c(n-2,:));
 228       b = diff (a) ./ h(1:n-1, idx) ...
 229           - h(1:n-1, idx) / 3 .* (c(2:n,:) + 2 * c(1:n-1,:));
 230       d = diff (c) ./ (3 * h(1:n-1, idx));
 231
 232       d = d(1:n-1,:);d = d.'(:);
 233       c = c(1:n-1,:);c = c.'(:);
 234       b = b(1:n-1,:);b = b.'(:);
 235       a = a(1:n-1,:);a = a.'(:);
 236     endif
 237
 238   endif
 239   ret = mkpp (x, cat (2, d, c, b, a), szy(1:end-1));
 240
 241   if (nargin == 3)
 242     ret = ppval (ret, xi);
 243   endif
 244
 245 endfunction
 246
 247 %!demo
 248 %! x = 0:10; y = sin(x);
 249 %! xspline = 0:0.1:10; yspline = spline(x,y,xspline);
 250 %! title("spline fit to points from sin(x)");
 251 %! plot(xspline,sin(xspline),"r",xspline,yspline,"g-",x,y,"b+");
 252 %! legend("original","interpolation","interpolation points");
 253 %! %--------------------------------------------------------
 254 %! % confirm that interpolated function matches the original
 255
 256 %!shared x,y,abserr
 257 %! x = [0:10]; y = sin(x); abserr = 1e-14;
 258 %!assert (spline(x,y,x), y, abserr);
 259 %!assert (spline(x,y,x'), y', abserr);
 260 %!assert (spline(x',y',x'), y', abserr);
 261 %!assert (spline(x',y',x), y, abserr);
 262 %!assert (isempty(spline(x',y',[])));
 263 %!assert (isempty(spline(x,y,[])));
 264 %!assert (spline(x,[y;y],x), [spline(x,y,x);spline(x,y,x)],abserr)
 265 %!assert (spline(x,[y;y],x'), [spline(x,y,x);spline(x,y,x)],abserr)
 266 %!assert (spline(x',[y;y],x), [spline(x,y,x);spline(x,y,x)],abserr)
 267 %!assert (spline(x',[y;y],x'), [spline(x,y,x);spline(x,y,x)],abserr)
 268 %! y = cos(x) + i*sin(x);
 269 %!assert (spline(x,y,x), y, abserr)
 270 %!assert (real(spline(x,y,x)), real(y), abserr);
 271 %!assert (real(spline(x,y,x.')), real(y).', abserr);
 272 %!assert (real(spline(x.',y.',x.')), real(y).', abserr);
 273 %!assert (real(spline(x.',y,x)), real(y), abserr);
 274 %!assert (imag(spline(x,y,x)), imag(y), abserr);
 275 %!assert (imag(spline(x,y,x.')), imag(y).', abserr);
 276 %!assert (imag(spline(x.',y.',x.')), imag(y).', abserr);
 277 %!assert (imag(spline(x.',y,x)), imag(y), abserr);
 278 %!test
 279 %! xnan = 5;
 280 %! y(x==xnan) = NaN;
 281 %! ok = ! isnan (y);
 282 %! assert (spline (x, y, x(ok)), y(ok), abserr);
 283 %!test
 284 %! ok = ! isnan (y);
 285 %! assert (! isnan (spline (x, y, x(!ok))));
 286 %!test
 287 %! x = [1,2];
 288 %! y = [1,4];
 289 %! assert (spline (x,y,x), [1,4], abserr);
 290 %!test
 291 %! x = [2,1];
 292 %! y = [1,4];
 293 %! assert (spline (x,y,x), [1,4], abserr);
 294 %!test
 295 %! x = [1,2];
 296 %! y = [1,2,3,4];
 297 %! pp = spline (x,y);
 298 %! [x,P] = unmkpp (pp);
 299 %! assert (norm (P-[3,-3,1,2]), 0, abserr);
 300 %!test
 301 %! x = [2,1];
 302 %! y = [1,2,3,4];
 303 %! pp = spline (x,y);
 304 %! [x,P] = unmkpp (pp);
 305 %! assert (norm (P-[7,-9,1,3]), 0, abserr);