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