normalized ALF (Assotiated Legendre Function)
Marco Caliari
marco.caliari at univr.it
Wed Feb 13 08:58:25 CST 2008
I discovered that the warning is misleading. Consider
legendre(151,0)
The last entry is -Inf and the script "should" returns a warning (you
should replace scale == inf with abs(scale) == inf). But, there is no
instability: the value -Inf is "right", in the sense that it is too large
in magnitude to be represented. Now consider
legendre(151,eps)
Again a warning (due to the check on a and b), but, again, no instability.
Moreover, Matlab has no warnings and it is less robust. Consider
legendre(151,-0.9)
in Matlab and with the script now enclosed. To conclude, I would suggest
to leave the warnings out.
Marco
>
> On Feb 13, 2008, at 8:16 AM, Marco Caliari wrote:
>
>>> That line was incorporated from your version. Should it be something else?
>>
>> Something like that to be added:
>>
>> if (max(abs(scale)) == Inf)
>> underflow = true;
>> endif
>
> ok.
>
> How about the attached?
>
> We still need some tests added for higher orders. I'll get around to that
> latter today.
>
> Ben
>
>
-------------- next part --------------
## Copyright (C) 2000, 2006, 2007 Kai Habel
## Copyright (C) 2008 Marco Caliari
##
## This file is part of Octave.
##
## Octave is free software; you can redistribute it and/or modify it
## under the terms of the GNU General Public License as published by
## the Free Software Foundation; either version 3 of the License, or (at
## your option) any later version.
##
## Octave is distributed in the hope that it will be useful, but
## WITHOUT ANY WARRANTY; without even the implied warranty of
## MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the GNU
## General Public License for more details.
##
## You should have received a copy of the GNU General Public License
## along with Octave; see the file COPYING. If not, see
## <http://www.gnu.org/licenses/>.
## -*- texinfo -*-
## @deftypefn {Function File} {@var{L} =} legendre (@var{n}, @var{X})
## @deftypefnx {Function File} {@var{L} =} legendre (@var{n}, @var{X}, "unnorm")
##
## Legendre Function of degree n and order m
## where all values for m = 0.. at var{n} are returned.
## @var{n} must be a scalar in the range [0..255].
## The return value has one dimension more than @var{x}.
##
## @example
## The Legendre Function of degree n and order m
##
## @group
## m m 2 m/2 d^m
## P(x) = (-1) * (1-x ) * ---- P (x)
## n dx^m n
## @end group
##
## with:
## Legendre polynomial of degree n
##
## @group
## 1 d^n 2 n
## P (x) = ------ [----(x - 1) ]
## n 2^n n! dx^n
## @end group
##
## legendre(3,[-1.0 -0.9 -0.8]) returns the matrix
##
## @group
## x | -1.0 | -0.9 | -0.8
## ------------------------------------
## m=0 | -1.00000 | -0.47250 | -0.08000
## m=1 | 0.00000 | -1.99420 | -1.98000
## m=2 | 0.00000 | -2.56500 | -4.32000
## m=3 | 0.00000 | -1.24229 | -3.24000
## @end group
## @end example
##
## @deftypefnx {Function File} {@var{L} =} legendre (@var{n}, @var{X}, "sch")
##
## Computes the Schmidt semi-normalized associated Legendre function.
## The Schmidt semi-normalized associated Legendre function is related
## to the unnormalized Legendre functions by
##
## @example
## For Legendre functions of degree n and order 0
##
## @group
## 0 0
## SP (x) = P (x)
## n n
## @end group
##
## For Legendre functions of degree n and order m
##
## @group
## m m m 2(n-m)! 0.5
## SP (x) = P (x) * (-1) * [-------]
## n n (n+m)!
## @end group
## @end example
##
## @deftypefnx {Function File} {@var{L} =} legendre (@var{n}, @var{X}, "norm")
##
## Computes the fully normalized associated Legendre function.
## The fully normalized associated Legendre function is related
## to the unnormalized Legendre functions by
##
## @example
## For Legendre functions of degree n and order m
##
## @group
## m m m (n+0.5)(n-m)! 0.5
## NP (x) = P (x) * (-1) * [-------------]
## n n (n+m)!
## @end group
## @end example
##
## @end deftypefn
## Author: Marco Caliari <marco.caliari at univr.it>
function L = legendre (n, x, normalization)
if (nargin < 2 || nargin > 3)
print_usage ();
endif
if nargin == 3
normalization = lower (normalization);
else
normalization = "unnorm";
endif
if (! isscalar (n) || n < 0 || n != fix (n))
error ("n must be a integer between 0 and 255]");
endif
if (! isvector (x) || any (x < -1 || x > 1))
error ("x must be vector in range -1 <= x <= 1");
endif
switch normalization
case "norm"
scale1 = sqrt (n+0.5);
case "sch"
scale1 = sqrt (2);
case "unnorm"
scale1 = 1;
otherwise
print_usage ();
endswitch
scale2 = x .* scale1;
## 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
underflow = false;
for m = 1:n
LPM = zeros(n+2-m,length(x));
LPM(1,:) = scale1;
LPM(2,:) = (2*m-1) .* scale2;
for l = m+1:n
a = (2*l-1) .* x .* LPM(l-m+1,:);
b = (l+m-2) .* LPM(l-m,:);
LPM(l-m+2,:) = ((2*l-1) .* x .* LPM(l-m+1,:)...
-(l+m-2) .* LPM(l-m,:))/(l-m+1);
endfor
L(m,:) = LPM(n+2-m,:);
if (strcmp (normalization, "unnorm"))
scale1 = -scale1 * (2*m-1);
scale2 = -scale2 * (2*m-1);
else
## normalization == "sch" or normalization == "norm"
scale1 = scale1 / sqrt ((n-m+1)*(n+m))*(2*m-1);
scale2 = scale2 / sqrt ((n-m+1)*(n+m))*(2*m-1);
endif
scale1 = scale1 .* sqrt(1-x.^2);
scale2 = scale2 .* sqrt(1-x.^2);
endfor
L(n+1,:) = scale1;
if (strcmp (normalization, "sch"))
L(1,:) = L(1,:) / sqrt (2);
endif
endfunction
%!test
%! result=legendre(3,[-1.0 -0.9 -0.8]);
%! expected = [
%! -1.00000 -0.47250 -0.08000
%! 0.00000 -1.99420 -1.98000
%! 0.00000 -2.56500 -4.32000
%! 0.00000 -1.24229 -3.24000
%! ];
%! assert(result,expected,1e-5);
%!test
%! result=legendre(3,[-1.0 -0.9 -0.8], "sch");
%! expected = [
%! -1.00000 -0.47250 -0.08000
%! 0.00000 0.81413 0.80833
%! -0.00000 -0.33114 -0.55771
%! 0.00000 0.06547 0.17076
%! ];
%! assert(result,expected,1e-5);
%!test
%! result=legendre(3,[-1.0 -0.9 -0.8], "norm");
%! expected = [
%! -1.87083 -0.88397 -0.14967
%! 0.00000 1.07699 1.06932
%! -0.00000 -0.43806 -0.73778
%! 0.00000 0.08661 0.22590
%! ];
%! assert(result,expected,1e-5);
More information about the Help-octave
mailing list