Mercurial > octave
changeset 11196:d17cb8a1271d
spline.m: Fix algorithm for input with 3 elements (bug #31098)
| author | Marco Caliari <marco.caliari@univr.it> |
|---|---|
| date | Sat, 06 Nov 2010 18:43:06 +0100 |
| parents | 8f67fe9dd64e |
| children | 836427db633b |
| files | scripts/ChangeLog scripts/polynomial/spline.m |
| diffstat | 2 files changed, 49 insertions(+), 32 deletions(-) [+] |
line wrap: on
line diff
--- a/scripts/ChangeLog Sat Nov 06 09:08:47 2010 -0400 +++ b/scripts/ChangeLog Sat Nov 06 18:43:06 2010 +0100 @@ -1,3 +1,8 @@ +2010-11-06 Marco Caliari <marco.caliari@univr.it> + + * polynomial/spline.m: Fit a parabola for input with 3 + elements (bug #31098). + 2010-11-04 Rik <octave@nomad.inbox5.com> * plot/__fltk_ginput__.m: Use semicolons to prevent internal
--- a/scripts/polynomial/spline.m Sat Nov 06 09:08:47 2010 -0400 +++ b/scripts/polynomial/spline.m Sat Nov 06 18:43:06 2010 +0100 @@ -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)]; + + else - 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); + 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 > 4) - 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(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); - 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(:)]], + 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(:)]], + [-1, 0, 1], n-2, n-2) \ g; - 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,:);