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));