Mercurial > hg > octave-lojdl
changeset 14575:2e23cd0a9e40
Fix complete spline with three points (bug #35739)
* spline.m: Correctly compute the coefficients of a complete
spline with three points. Add two tests to check for this bug.
author | Marco Caliari <marco.caliari@univr.it> |
---|---|
date | Wed, 11 Apr 2012 12:54:41 +0200 |
parents | 9546bb28648c |
children | 4dbb47d09219 |
files | scripts/polynomial/spline.m |
diffstat | 1 files changed, 20 insertions(+), 27 deletions(-) [+] |
line wrap: on
line diff
--- a/scripts/polynomial/spline.m +++ b/scripts/polynomial/spline.m @@ -131,31 +131,12 @@ b = b(1:n-1,:); a = a(1:n-1,:); else - if (n == 3) - dg = 1.5 * h(1) - 0.5 * h(2); - c(2:n-1,:) = 1/dg(1); - else - dg = 2 * (h(1:n-2) .+ h(2:n-1)); - dg(1) = dg(1) - 0.5 * h(1); - dg(n-2) = dg(n-2) - 0.5 * h(n-1); - - e = h(2:n-2); - - g = 3 * diff (a(2:n,:)) ./ h(2:n-1,idx) ... - - 3 * diff (a(1:n-1,:)) ./ h(1:n-2,idx); - g(1,:) = 3 * (a(3,:) - a(2,:)) / h(2) ... - - 3 / 2 * (3 * (a(2,:) - a(1,:)) / h(1) - dfs); - g(n-2,:) = 3 / 2 * (3 * (a(n,:) - a(n-1,:)) / h(n-1) - dfe) ... - - 3 * (a(n-1,:) - a(n-2,:)) / h(n-2); - - c(2:n-1,:) = spdiags ([[e(:); 0], dg, [0; e(:)]], - [-1, 0, 1], n-2, n-2) \ g; - endif - - c(1,:) = (3 / h(1) * (a(2,:) - a(1,:)) - 3 * dfs - - c(2,:) * h(1)) / (2 * h(1)); - c(n,:) = - (3 / h(n-1) * (a(n,:) - a(n-1,:)) - 3 * dfe - + c(n-1,:) * h(n-1)) / (2 * h(n-1)); + g(1,:) = (a(2,:) - a(1,:)) / h(1) - dfs; + g(2:n-1,:) = (a(3:n,:) - a(2:n-1,:)) ./ h(2:n-1) - ... + (a(2:n-1,:) - a(1:n-2,:)) ./ h(1:n-2); + g(n,:) = dfe - (a(n,:) - a(n-1,:)) / h(n-1); + c = spdiags([[h/6;0],[h(1)/3;(h(1:n-2)+h(2:n-1))/3;h(n-1)/3],[0;h/6]],... + [-1,0,1],n,n) \ (g / 2); b(1:n-1,:) = 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)); @@ -297,11 +278,23 @@ %! y = [1,2,3,4]; %! pp = spline (x,y); %! [x,P] = unmkpp (pp); -%! assert (norm (P-[3,-3,1,2]), 0, abserr); +%! assert (P, [3,-3,1,2], abserr); %!test %! x = [2,1]; %! y = [1,2,3,4]; %! pp = spline (x,y); %! [x,P] = unmkpp (pp); -%! assert (norm (P-[7,-9,1,3]), 0, abserr); +%! assert (P, [7,-9,1,3], abserr); +%!test +%! x = [0,1,2]; +%! y = [0,0,1,0,0]; +%! pp = spline (x,y); +%! [x,P] = unmkpp (pp); +%! assert (P, [-2,3,0,0;2,-3,0,1], abserr); +%!test +%! x = [0,1,2,3]; +%! y = [0,0,1,1,0,0]; +%! pp = spline (x,y); +%! [x,P] = unmkpp (pp); +%! assert (P, [-1,2,0,0;0,-1,1,1;1,-1,-1,1], abserr);