# HG changeset patch # User Marco Caliari # Date 1287420512 -7200 spline.m: Fix algorithm for input with 3 elements (bug #31098) diff -r bd53bcd30d03 scripts/ChangeLog --- a/scripts/ChangeLog Fri Sep 10 20:57:15 2010 +0200 +++ b/scripts/ChangeLog Mon Oct 18 18:48:32 2010 +0200 @@ -0,0 +1,5 @@ +2010-10-18 Marco Caliari + + * polynomial/spline.m: Fit a parabola for input with 3 + elements (bug #31098). + diff -r bd53bcd30d03 scripts/polynomial/spline.m --- a/scripts/polynomial/spline.m Fri Sep 10 20:57:15 2010 +0200 +++ b/scripts/polynomial/spline.m Mon Oct 18 18:48:32 2010 +0200 @@ -144,48 +144,60 @@ else - g = zeros (n-2, columns (a)); - g(1,:) = 3 / (h(1) + h(2)) ... - * (a(3,:) - a(2,:) - h(2) / h(1) * (a(2,:) - a(1,:))); - g(n-2,:) = 3 / (h(n-1) + h(n-2)) ... - * (h(n-2) / h(n-1) * (a(n,:) - a(n-1,:)) - (a(n-1,:) - a(n-2,:))); + if (n == 3) - if (n > 4) + n = 2; + c = (a(1,:) - a(3,:)) / ((x(3) - x(1)) * (x(2) - x(3))) ... + + (a(2,:) - a(1,:)) / ((x(2) - x(1)) * (x(2) - x(3))); + b = (a(2,:) - a(1,:)) * (x(3) - x(1)) ... + / ((x(2) - x(1)) * (x(3) - x(2))) ... + + (a(1,:) - a(3,:)) * (x(2) - x(1)) ... + / ((x(3) - x(1)) * (x(3) - x(2))); + a = a(1,:); + d = zeros (size (a)); + x = [min(x), max(x)]; - g(2:n - 3,:) = 3 * diff (a(3:n-1,:)) ./ h(3:n-2,idx) ... - - 3 * diff (a(2:n-2,:)) ./ h(2:n - 3,idx); + else - dg = 2 * (h(1:n-2) .+ h(2:n-1)); - dg(1) = dg(1) - h(1); - dg(n-2) = dg(n-2) - h(n-1); + g = zeros (n-2, columns (a)); + g(1,:) = 3 / (h(1) + h(2)) ... + * (a(3,:) - a(2,:) - h(2) / h(1) * (a(2,:) - a(1,:))); + g(n-2,:) = 3 / (h(n-1) + h(n-2)) ... + * (h(n-2) / h(n-1) * (a(n,:) - a(n-1,:)) - (a(n-1,:) - a(n-2,:))); - ldg = udg = h(2:n-2); - udg(1) = udg(1) - h(1); - ldg(n - 3) = ldg(n-3) - h(n-1); - c(2:n-1,:) = spdiags ([[ldg(:); 0], dg, [0; udg(:)]], + if (n > 4) + + g(2:n - 3,:) = 3 * diff (a(3:n-1,:)) ./ h(3:n-2,idx) ... + - 3 * diff (a(2:n-2,:)) ./ h(2:n - 3,idx); + + dg = 2 * (h(1:n-2) .+ h(2:n-1)); + dg(1) = dg(1) - h(1); + dg(n-2) = dg(n-2) - h(n-1); + + ldg = udg = h(2:n-2); + udg(1) = udg(1) - h(1); + ldg(n - 3) = ldg(n-3) - h(n-1); + c(2:n-1,:) = spdiags ([[ldg(:); 0], dg, [0; udg(:)]], [-1, 0, 1], n-2, n-2) \ g; - elseif (n == 4) + elseif (n == 4) - dg = [h(1) + 2 * h(2); 2 * h(2) + h(3)]; - ldg = h(2) - h(3); - udg = h(2) - h(1); - c(2:n-1,:) = spdiags ([[ldg(:);0], dg, [0; udg(:)]], + dg = [h(1) + 2 * h(2); 2 * h(2) + h(3)]; + ldg = h(2) - h(3); + udg = h(2) - h(1); + c(2:n-1,:) = spdiags ([[ldg(:);0], dg, [0; udg(:)]], [-1, 0, 1], n-2, n-2) \ g; - - else # n == 3 - - dg = h(1) + 2 * h(2); - c(2:n-1,:) = g/dg(1); + + endif + + c(1,:) = c(2,:) + h(1) / h(2) * (c(2,:) - c(3,:)); + c(n,:) = c(n-1,:) + h(n-1) / h(n-2) * (c(n-1,:) - c(n-2,:)); + b = diff (a) ./ h(1:n-1, idx) ... + - h(1:n-1, idx) / 3 .* (c(2:n,:) + 2 * c(1:n-1,:)); + d = diff (c) ./ (3 * h(1:n-1, idx)); endif - c(1,:) = c(2,:) + h(1) / h(2) * (c(2,:) - c(3,:)); - c(n,:) = c(n-1,:) + h(n-1) / h(n-2) * (c(n-1,:) - c(n-2,:)); - b = diff (a) ./ h(1:n-1, idx) ... - - h(1:n-1, idx) / 3 .* (c(2:n,:) + 2 * c(1:n-1,:)); - d = diff (c) ./ (3 * h(1:n-1, idx)); - endif d = d(1:n-1,:);