Mercurial > octave-nkf

comparison scripts/specfun/legendre.m @ 12893:72ffa81a68d4 stable

Find changesets by keywords (author, files, the commit message), revision number or hash, or revset expression.
legendre.m: Allow ND-array inputs (Bug #33526). * legendre.m: Allow ND-array inputs (Bug #33526).
author Marco Caliari <marco.caliari@univr.it>
date Wed, 27 Jul 2011 11:21:21 -0700
parents c792872f8942
children 72c96de7a403
comparison
equal deleted inserted replaced
12892:67bf9b30f3f9 12893:72ffa81a68d4
164 164
165 if (!isscalar (n) || n < 0 || n != fix (n)) 165 if (!isscalar (n) || n < 0 || n != fix (n))
166 error ("legendre: N must be a non-negative scalar integer"); 166 error ("legendre: N must be a non-negative scalar integer");
167 endif 167 endif
168 168
169 if (!isvector (x) || !isreal (x) || any (x < -1 | x > 1)) 169 if (!isreal (x) || any (x(:) < -1 | x(:) > 1))
170 error ("legendre: X must be real-valued vector in the range -1 <= X <= 1"); 170 error ("legendre: X must be real-valued vector in the range -1 <= X <= 1");
171 endif 171 endif
172 172
173 if (nargin == 3) 173 if (nargin == 3)
174 normalization = lower (normalization); 174 normalization = lower (normalization);
185 scale = 1; 185 scale = 1;
186 otherwise 186 otherwise
187 error ('legendre: expecting NORMALIZATION option to be "norm", "sch", or "unnorm"'); 187 error ('legendre: expecting NORMALIZATION option to be "norm", "sch", or "unnorm"');
188 endswitch 188 endswitch
189 189
190 if (rows (x) != 1) 190 scale = scale * ones (size (x));
191 x = x';
192 endif
193 scale = scale * ones (1, numel (x));
194 191
195 ## Based on the recurrence relation below 192 ## Based on the recurrence relation below
196 ## m m m 193 ## m m m
197 ## (n-m+1) * P (x) = (2*n+1)*x*P (x) - (n+1)*P (x) 194 ## (n-m+1) * P (x) = (2*n+1)*x*P (x) - (n+1)*P (x)
198 ## n+1 n n-1 195 ## n+1 n n-1
199 ## http://en.wikipedia.org/wiki/Associated_Legendre_function 196 ## http://en.wikipedia.org/wiki/Associated_Legendre_function
200 197
201 overflow = false; 198 overflow = false;
199 retval = zeros([n+1, size(x)]);
202 for m = 1:n 200 for m = 1:n
203 lpm1 = scale; 201 lpm1 = scale;
204 lpm2 = (2*m-1) .* x .* scale; 202 lpm2 = (2*m-1) .* x .* scale;
205 lpm3 = lpm2; 203 lpm3 = lpm2;
206 for k = m+1:n 204 for k = m+1:n
215 || any (abs (lpm3) > realmax)) 213 || any (abs (lpm3) > realmax))
216 overflow = true; 214 overflow = true;
217 endif 215 endif
218 endif 216 endif
219 endfor 217 endfor
220 retval(m,:) = lpm3; 218 retval(m,:) = lpm3(:);
221 if (strcmp (normalization, "unnorm")) 219 if (strcmp (normalization, "unnorm"))
222 scale = -scale * (2*m-1); 220 scale = -scale * (2*m-1);
223 else 221 else
224 ## normalization == "sch" or normalization == "norm" 222 ## normalization == "sch" or normalization == "norm"
225 scale = scale / sqrt ((n-m+1)*(n+m))*(2*m-1); 223 scale = scale / sqrt ((n-m+1)*(n+m))*(2*m-1);
226 endif 224 endif
227 scale = scale .* sqrt(1-x.^2); 225 scale = scale .* sqrt(1-x.^2);
228 endfor 226 endfor
229 227
230 retval(n+1,:) = scale; 228 retval(n+1,:) = scale(:);
229
230 if (isvector (x))
231 ## vector case is special
232 retval = reshape (retval, n + 1, length (x));
233 endif
231 234
232 if (strcmp (normalization, "sch")) 235 if (strcmp (normalization, "sch"))
233 retval(1,:) = retval(1,:) / sqrt (2); 236 retval(1,:) = retval(1,:) / sqrt (2);
234 endif 237 endif
235 238
237 warning ("legendre: overflow - results may be unstable for high orders"); 240 warning ("legendre: overflow - results may be unstable for high orders");
238 warned_overflow = true; 241 warned_overflow = true;
239 endif 242 endif
240 243
241 endfunction 244 endfunction
245
242 246
243 %!test 247 %!test
244 %! result = legendre (3, [-1.0 -0.9 -0.8]); 248 %! result = legendre (3, [-1.0 -0.9 -0.8]);
245 %! expected = [ 249 %! expected = [
246 %! -1.00000 -0.47250 -0.08000 250 %! -1.00000 -0.47250 -0.08000
276 %! assert (result(end) < -1.7976e308 && all (isfinite (result(1:end-1)))); 280 %! assert (result(end) < -1.7976e308 && all (isfinite (result(1:end-1))));
277 281
278 %!test 282 %!test
279 %! result = legendre (150, 0); 283 %! result = legendre (150, 0);
280 %! ## This agrees with Matlab's result. 284 %! ## This agrees with Matlab's result.
281 %! assert (result(end), 3.7532741115719e+306, 0.0000000000001e+306) 285 %! assert (result(end), 3.7532741115719e+306, 0.0000000000001e+306);
282 286
283 %!test 287 %!test
284 %! result = legendre (0, 0:0.1:1); 288 %! result = legendre (0, 0:0.1:1);
285 %! assert (result, full(ones(1,11))) 289 %! assert (result, full(ones(1,11)));
290
291 %!test
292 %! result = legendre (3, [-1,0,1;1,0,-1]);
293 %! ## Test matrix input
294 %! expected(:,:,1) = [-1,1;0,0;0,0;0,0];
295 %! expected(:,:,2) = [0,0;1.5,1.5;0,0;-15,-15];
296 %! expected(:,:,3) = [1,-1;0,0;0,0;0,0];
297 %! assert (result, expected);
298
299 %!test
300 %! result = legendre (3, [-1,0,1;1,0,-1]');
301 %! expected(:,:,1) = [-1,0,1;0,1.5,0;0,0,0;0,-15,0];
302 %! expected(:,:,2) = [1,0,-1;0,1.5,0;0,0,0;0,-15,0];
303 %! assert (result, expected);
286 304
287 %% Check correct invocation 305 %% Check correct invocation
288 %!error legendre (); 306 %!error legendre ();
289 %!error legendre (1); 307 %!error legendre (1);
290 %!error legendre (1,2,3,4); 308 %!error legendre (1,2,3,4);