Mercurial > octave-nkf
changeset 8255:3f1199ad212f
legendre.m: Warn once on under/overflow.
| author | Ben Abbott <bpabbott@mac.com> |
|---|---|
| date | Tue, 21 Oct 2008 21:17:29 -0400 |
| parents | 8c29549c66de |
| children | dba0037e6602 |
| files | scripts/ChangeLog scripts/specfun/legendre.m |
| diffstat | 2 files changed, 25 insertions(+), 2 deletions(-) [+] |
line wrap: on
line diff
--- a/scripts/ChangeLog Tue Oct 21 21:07:05 2008 -0400 +++ b/scripts/ChangeLog Tue Oct 21 21:17:29 2008 -0400 @@ -1,5 +1,7 @@ 2008-10-21 Ben Abbott <bpabbott@mac.com> + * specfun/legendre.m: Warn once on under/overflow. + * plot/clf.m: Improve Matlab compatibility. 2008-10-21 John W. Eaton <jwe@octave.org>
--- a/scripts/specfun/legendre.m Tue Oct 21 21:07:05 2008 -0400 +++ b/scripts/specfun/legendre.m Tue Oct 21 21:17:29 2008 -0400 @@ -111,6 +111,8 @@ function retval = legendre (n, x, normalization) + persistent warned_overflow = false; + if (nargin < 2 || nargin > 3) print_usage (); endif @@ -140,20 +142,32 @@ error ("legendre: expecting normalization option to be \"norm\", \"sch\", or \"unnorm\""); endswitch + scale = scale * ones (1, numel (x)); + ## Based on the recurrence relation below ## m m m ## (n-m+1) * P (x) = (2*n+1)*x*P (x) - (n+1)*P (x) ## n+1 n n-1 ## http://en.wikipedia.org/wiki/Associated_Legendre_function + overflow = false; for m = 1:n lpm1 = scale; lpm2 = (2*m-1) .* x .* scale; lpm3 = lpm2; for k = m+1:n - lpm3 = ((2*k-1) .* x .* lpm2 - (k+m-2) .* lpm1)/(k-m+1); + lpm3a = (2*k-1) .* x .* lpm2; + lpm3b = (k+m-2) .* lpm1; + lpm3 = (lpm3a - lpm3b)/(k-m+1); lpm1 = lpm2; lpm2 = lpm3; + if (! warned_overflow) + if (any (abs (lpm3a) > realmax) + || any (abs (lpm3b) > realmax) + || any (abs (lpm3) > realmax)) + overflow = true; + endif + endif endfor retval(m,:) = lpm3; if (strcmp (normalization, "unnorm")) @@ -171,6 +185,11 @@ retval(1,:) = retval(1,:) / sqrt (2); endif + if (overflow && ! warned_overflow) + warning ("legendre: overflow - results may be unstable for high orders."); + warned_overflow = true; + endif + endfunction %!test @@ -213,4 +232,6 @@ %! ## This agrees with Matlab's result. %! assert (result(end), 3.7532741115719e+306, 0.0000000000001e+306) - +%!test +%! result = legendre (0, 0:0.1:1); +%! assert (result, ones(1,11))