Mercurial > octave-nkf

comparison scripts/polynomial/spline.m @ 10549:95c3e38098bf

Find changesets by keywords (author, files, the commit message), revision number or hash, or revset expression.
Untabify .m scripts
author Rik <code@nomad.inbox5.com>
date Fri, 23 Apr 2010 11:28:50 -0700
parents 31900e17b5f5
children c66a4657d764
comparison
equal deleted inserted replaced
10548:479536c5bb10 10549:95c3e38098bf
123 - 3 / 2 * (3 * (a(2,:) - a(1,:)) / h(1) - dfs); 123 - 3 / 2 * (3 * (a(2,:) - a(1,:)) / h(1) - dfs);
124 g(n-2,:) = 3 / 2 * (3 * (a(n,:) - a(n-1,:)) / h(n-1) - dfe) ... 124 g(n-2,:) = 3 / 2 * (3 * (a(n,:) - a(n-1,:)) / h(n-1) - dfe) ...
125 - 3 * (a(n-1,:) - a(n-2,:)) / h(n-2); 125 - 3 * (a(n-1,:) - a(n-2,:)) / h(n-2);
126 126
127 c(2:n-1,:) = spdiags ([[e(:); 0], dg, [0; e(:)]], 127 c(2:n-1,:) = spdiags ([[e(:); 0], dg, [0; e(:)]],
128 [-1, 0, 1], n-2, n-2) \ g; 128 [-1, 0, 1], n-2, n-2) \ g;
129 endif 129 endif
130 130
131 c(1,:) = (3 / h(1) * (a(2,:) - a(1,:)) - 3 * dfs 131 c(1,:) = (3 / h(1) * (a(2,:) - a(1,:)) - 3 * dfs
132 - c(2,:) * h(1)) / (2 * h(1)); 132 - c(2,:) * h(1)) / (2 * h(1));
133 c(n,:) = - (3 / h(n-1) * (a(n,:) - a(n-1,:)) - 3 * dfe 133 c(n,:) = - (3 / h(n-1) * (a(n,:) - a(n-1,:)) - 3 * dfe
134 + c(n-1,:) * h(n-1)) / (2 * h(n-1)); 134 + c(n-1,:) * h(n-1)) / (2 * h(n-1));
135 b(1:n-1,:) = diff (a) ./ h(1:n-1, idx) ... 135 b(1:n-1,:) = diff (a) ./ h(1:n-1, idx) ...
136 - h(1:n-1,idx) / 3 .* (c(2:n,:) + 2 * c(1:n-1,:)); 136 - h(1:n-1,idx) / 3 .* (c(2:n,:) + 2 * c(1:n-1,:));
137 d = diff (c) ./ (3 * h(1:n-1, idx)); 137 d = diff (c) ./ (3 * h(1:n-1, idx));
138 138
139 else 139 else
140 140
141 g = zeros (n-2, columns (a)); 141 g = zeros (n-2, columns (a));
142 g(1,:) = 3 / (h(1) + h(2)) ... 142 g(1,:) = 3 / (h(1) + h(2)) ...
143 * (a(3,:) - a(2,:) - h(2) / h(1) * (a(2,:) - a(1,:))); 143 * (a(3,:) - a(2,:) - h(2) / h(1) * (a(2,:) - a(1,:)));
144 g(n-2,:) = 3 / (h(n-1) + h(n-2)) ... 144 g(n-2,:) = 3 / (h(n-1) + h(n-2)) ...
145 * (h(n-2) / h(n-1) * (a(n,:) - a(n-1,:)) - (a(n-1,:) - a(n-2,:))); 145 * (h(n-2) / h(n-1) * (a(n,:) - a(n-1,:)) - (a(n-1,:) - a(n-2,:)));
146 146
147 if (n > 4) 147 if (n > 4)
148 148
149 g(2:n - 3,:) = 3 * diff (a(3:n-1,:)) ./ h(3:n-2,idx) ... 149 g(2:n - 3,:) = 3 * diff (a(3:n-1,:)) ./ h(3:n-2,idx) ...
150 - 3 * diff (a(2:n-2,:)) ./ h(2:n - 3,idx); 150 - 3 * diff (a(2:n-2,:)) ./ h(2:n - 3,idx);
155 155
156 ldg = udg = h(2:n-2); 156 ldg = udg = h(2:n-2);
157 udg(1) = udg(1) - h(1); 157 udg(1) = udg(1) - h(1);
158 ldg(n - 3) = ldg(n-3) - h(n-1); 158 ldg(n - 3) = ldg(n-3) - h(n-1);
159 c(2:n-1,:) = spdiags ([[ldg(:); 0], dg, [0; udg(:)]], 159 c(2:n-1,:) = spdiags ([[ldg(:); 0], dg, [0; udg(:)]],
160 [-1, 0, 1], n-2, n-2) \ g; 160 [-1, 0, 1], n-2, n-2) \ g;
161 161
162 elseif (n == 4) 162 elseif (n == 4)
163 163
164 dg = [h(1) + 2 * h(2); 2 * h(2) + h(3)]; 164 dg = [h(1) + 2 * h(2); 2 * h(2) + h(3)];
165 ldg = h(2) - h(3); 165 ldg = h(2) - h(3);
166 udg = h(2) - h(1); 166 udg = h(2) - h(1);
167 c(2:n-1,:) = spdiags ([[ldg(:);0], dg, [0; udg(:)]], 167 c(2:n-1,:) = spdiags ([[ldg(:);0], dg, [0; udg(:)]],
168 [-1, 0, 1], n-2, n-2) \ g; 168 [-1, 0, 1], n-2, n-2) \ g;
169 169
170 else # n == 3 170 else # n == 3
171 171
172 dg = h(1) + 2 * h(2); 172 dg = h(1) + 2 * h(2);
173 c(2:n-1,:) = g/dg(1); 173 c(2:n-1,:) = g/dg(1);
174 174
175 endif 175 endif
176 176
177 c(1,:) = c(2,:) + h(1) / h(2) * (c(2,:) - c(3,:)); 177 c(1,:) = c(2,:) + h(1) / h(2) * (c(2,:) - c(3,:));
178 c(n,:) = c(n-1,:) + h(n-1) / h(n-2) * (c(n-1,:) - c(n-2,:)); 178 c(n,:) = c(n-1,:) + h(n-1) / h(n-2) * (c(n-1,:) - c(n-2,:));
179 b = diff (a) ./ h(1:n-1, idx) ... 179 b = diff (a) ./ h(1:n-1, idx) ...
180 - h(1:n-1, idx) / 3 .* (c(2:n,:) + 2 * c(1:n-1,:)); 180 - h(1:n-1, idx) / 3 .* (c(2:n,:) + 2 * c(1:n-1,:));
181 d = diff (c) ./ (3 * h(1:n-1, idx)); 181 d = diff (c) ./ (3 * h(1:n-1, idx));
182 182
183 endif 183 endif
184 184
185 d = d(1:n-1,:); 185 d = d(1:n-1,:);