Mercurial > hg > octave-lyh
diff scripts/polynomial/spline.m @ 5838:376e02b2ce70
[project @ 2006-06-01 20:23:53 by jwe]
author | jwe |
---|---|
date | Thu, 01 Jun 2006 20:23:54 +0000 |
parents | 55404f3b0da1 |
children | 437f9086b967 |
line wrap: on
line diff
--- a/scripts/polynomial/spline.m +++ b/scripts/polynomial/spline.m @@ -112,88 +112,93 @@ b = c = zeros (size (a)); h = diff (x); - idx = ones (columns (a),1); + idx = ones (columns (a), 1); if (complete) if (n == 3) dg = 1.5 * h(1) - 0.5 * h(2); - c(2:n - 1,:) = 1/dg(1); + c(2:n-1,:) = 1/dg(1); else - dg = 2 * (h(1:n - 2) .+ h(2:n - 1)); + 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); + dg(n-2) = dg(n-2) - 0.5 * h(n-1); - e = h(2:n - 2); + e = h(2:n-2); size(a) size(h) n - 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) \ + 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); + 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; + 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)); - 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)); + - 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)); + 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)); 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,:))); + 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) - 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(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 = 2 * (h(1:n-2) .+ h(2:n-1)); dg(1) = dg(1) - h(1); - dg(n - 2) = dg(n-2) - h(n - 1); + dg(n-2) = dg(n-2) - h(n-1); - ldg = udg = h(2: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(:)]],[-1,0,1],n-2,n-2) \ g; + 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) 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; + 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); + 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)); + 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,:); c=c(1:n-1,:); b=b(1:n-1,:); a=a(1:n-1,:); + d = d(1:n-1,:); + c = c(1:n-1,:); + b = b(1:n-1,:); + a = a(1:n-1,:); coeffs = [d(:), c(:), b(:), a(:)]; ret = mkpp (x, coeffs, szy(1:end-1));