normalized ALF (Assotiated Legendre Function)

[kruvalig]

In matlab enxist function legendre, and it can compute fully normalized  
assotiated function, and normalized by Schmidt. Does exist anything the  
same in octave.

I now that legendre-function in octave can compute un-normalized function.

And i interested in fully normalized ALF

                 (2-delta(0,m))(2n+1)(n-m)!
P_nm(t) = sqrt(----------------------------) P_nm(t)
                         (n+m)!

_______________________________________________
Help-octave mailing list
Help-octave@...
https://www.cae.wisc.edu/mailman/listinfo/help-octave

[Ben Abbott]

kruvalig wrote:

In matlab enxist function legendre, and it can compute fully normalized  
assotiated function, and normalized by Schmidt. Does exist anything the  
same in octave.

I've posted a patch and changelog to the maintainers list. Attached to this
email is the modified script.

legendre.m

[kruvalig]

Thank's a lot for quick answer. I very suprised of it.

On Sun, 10 Feb 2008 01:37:10 +0300, Ben Abbott <bpabbott@...> wrote:

> I've posted a patch and changelog to the maintainers list. Attached to  
> this
> email is the modified script.
> http://www.nabble.com/file/p15391002/legendre.m legendre.m

Can any one test output of this file and output matlab file to high degree  
(>70).

   result_mlab=legendre(80,[-1:0.1:1], "norm");
   result_octave=legendre(80,[-1:0.1:1], "norm");
   max(max(abs(result_mlab-result_octave)))

Here ftp://77.108.213.161/result_octave_80 is a file, generated in octave  
by code:

   result_octave=legendre(80,[-1:0.1:1], "norm");
   save -v7 result_octave_80

I think this code
   scale = (-1).^m .* sqrt ((n+0.5) .* factorial (n-m) ./ factorial (n+m));
whill give bad answer in high order of legendre.
And we have to write: To what order it compute in right.

_______________________________________________
Help-octave mailing list
Help-octave@...
https://www.cae.wisc.edu/mailman/listinfo/help-octave

[Ben Abbott]

On Feb 11, 2008, at 3:01 AM, kruvalig wrote:

> Thank's a lot for quick answer. I very suprised of it.
>
> On Sun, 10 Feb 2008 01:37:10 +0300, Ben Abbott <bpabbott@...>  
> wrote:
>
>> I've posted a patch and changelog to the maintainers list. Attached  
>> to
>> this
>> email is the modified script.
>> http://www.nabble.com/file/p15391002/legendre.m legendre.m
>
> Can any one test output of this file and output matlab file to high  
> degree
> (>70).
>
>   result_mlab=legendre(80,[-1:0.1:1], "norm");
>   result_octave=legendre(80,[-1:0.1:1], "norm");
>   max(max(abs(result_mlab-result_octave)))

        warning: legendre is unstable for higher orders
        ans =  1.4343e+16

This is what I get.

> Here ftp://77.108.213.161/result_octave_80 is a file, generated in  
> octave
> by code:
>
>   result_octave=legendre(80,[-1:0.1:1], "norm");
>   save -v7 result_octave_80
>
> I think this code
>   scale = (-1).^m .* sqrt ((n+0.5) .* factorial (n-m) ./ factorial (n
> +m));
> whill give bad answer in high order of legendre.
> And we have to write: To what order it compute in right.

I'm not so sure a problem lies with the normalization. When I compare  
results for legendre(80, [-1:0.1:1]), I get

        warning: legendre is unstable for higher orders
        ans =  5.4154e+126

I haven't looked into the details as how this function works yet. I'll  
do that today.

In any event, can you (someone else?) propose a change? ... or a  
reference for how it should be done?

Ben

_______________________________________________
Help-octave mailing list
Help-octave@...
https://www.cae.wisc.edu/mailman/listinfo/help-octave

[Marco Caliari-4]

Hi.

Please, try the enclosed legendre.m. It is based on the recurrence
relation found in
http://en.wikipedia.org/wiki/Associated_Legendre_function. I think it
should be much more stable.

Best regards,

Marco
function L = legendre(n,x,varargin)
% Based on the first recurrence relation found in
% http://en.wikipedia.org/wiki/Associated_Legendre_function
x = x(:).';
L = zeros(n+1,length(x));
if (nargin == 2)
  normalization = "unnorm";
else
  normalization = varargin{1};
end
for m = 1:n+1
  if (strcmp(normalization,"norm"))
    scale = sqrt((n+0.5)/prod(n-m+2:n+m-1));
  elseif (strcmp(normalization,"sch"))
    scale = sqrt(2/prod(n-m+2:n+m-1));
  else
    scale = (-1).^(m-1);
  end
  LP = ones(n+1,length(x))*scale;
  LP(m,:) = prod(1:2:2*(m-1)-1)*LP(m,:).*sqrt(1-x.^2).^(m-1);
  if (m == n+1)
    L(m,:) = LP(m,:);
    break
  end  
  LP(m+1,:) = (2*m-1).*x.*LP(m,:);
  for l = m+1:n
    LP(l+1,:) = ((2*(l-1)+1).*x.*LP(l,:)-(l-1+m-1).*LP(l-1,:))./(l-m+1);
  end
  L(m,:) = LP(n+1,:);
end
if (strcmp(normalization,"sch"))
  L(1,:) = L(1,:)/sqrt(2/prod(n+1:n));
end
_______________________________________________
Help-octave mailing list
Help-octave@...
https://www.cae.wisc.edu/mailman/listinfo/help-octave

[Ben Abbott]

On Feb 11, 2008, at 10:20 AM, Marco Caliari wrote:

Please, try the enclosed legendre.m. It is based on the recurrence relation found in http://en.wikipedia.org/wiki/Associated_Legendre_function. I think it should be much more stable.

I compared Matlab's script and Octave's and your proposal.

  result_matlab = legendre (80, [-1:0.1:1]);

result_octave0 respects the original and result_octave1 respects your algorithm.

I calculated an error with respect to the matlab version (I'm not sure Matlab's is to be trusted as correct in all cases).

d0 = abs (result_matlab - result_octave0) / abs (result_matlab); d1 = abs (result_matlab - result_octave1) / abs (result_matlab); er1 = max (d1, [], 2); er0 = max (d0, [], 2);

[er0(:), er1(:)] produces the following

ans =

   2.3424e-119   2.0877e-144

   1.1666e-120   2.8601e-144

   3.1132e-119   1.3575e-140

   8.1518e-118   1.0167e-139

   2.1137e-116   8.8310e-137

   5.4094e-115   1.6560e-135

   1.3778e-113   5.7097e-133

   3.5365e-112   1.3705e-131

   9.1806e-111   3.3573e-129

   2.5238e-109   1.3295e-127

   8.9744e-108   1.8397e-125

   3.8727e-106   1.2779e-123

   1.9273e-104   9.8135e-122

   1.0123e-102   9.4790e-120

   5.4782e-101   4.6351e-118

    2.9220e-99   6.1185e-116

    1.5443e-97   1.6927e-114

    7.9263e-96   5.0573e-112

    3.9915e-94   4.8575e-111

    1.9686e-92   2.5731e-108

    9.4747e-91   1.0210e-106

    4.4378e-89   1.2955e-104

    2.0255e-87   9.2168e-103

    8.9000e-86   5.0162e-101

    3.7992e-84    7.4455e-99

    1.5630e-82    6.0867e-98

    6.3553e-81    4.0858e-95

    2.4881e-79    1.5141e-93

    9.3534e-78    1.6912e-91

    3.4021e-76    1.8648e-89

    1.2025e-74    5.4167e-88

    4.1090e-73    1.0303e-85

    1.3681e-71    4.2479e-84

    4.2235e-70    3.6317e-82

    1.2357e-68    4.2982e-80

    3.4224e-67    7.9412e-79

    9.0728e-66    2.0485e-76

    2.1559e-64    1.1802e-74

    4.2186e-63    5.7797e-73

    4.8458e-62    7.6818e-71

    1.0928e-61    4.0895e-69

    4.4609e-59    2.5363e-67

    3.6901e-57    2.3634e-65

    2.5598e-55    3.8612e-64

    8.8479e-54    8.5079e-62

    3.4999e-52    7.1329e-60

    1.0974e-50    1.5563e-58

    3.8537e-49    2.3469e-56

    8.7346e-48    2.4322e-54

    2.2104e-46    6.8924e-53

    7.5648e-45    4.6632e-51

    1.6279e-43    6.4798e-49

    3.0916e-42    2.7505e-47

    7.0616e-42    7.3142e-46

    8.1111e-40    7.6535e-44

    7.9301e-39    6.1412e-42

    5.7348e-37    3.2699e-40

    6.0936e-36    1.4796e-38

    4.3591e-34    9.0850e-37

    1.2343e-32    5.8704e-35

    3.4634e-31    3.9178e-33

    1.0903e-29    2.0041e-31

    8.0223e-29    9.2888e-30

    2.4566e-27    3.6487e-28

    4.2821e-26    1.7018e-26

    1.9256e-25    6.5910e-25

    8.1766e-25    2.4639e-23

    4.2551e-23    6.8962e-22

    1.1459e-22    1.7888e-20

    3.2583e-20    3.9508e-19

    1.8758e-18    6.5279e-18

    3.1305e-17    1.1644e-16

    9.8748e-16    1.9183e-15

    1.5795e-16    3.3544e-14

    2.9368e-13    4.7596e-13

    2.8870e-12    7.9144e-12

    5.9192e-11    9.2655e-11

    6.6729e-10    9.5730e-10

    4.7224e-09    7.8467e-09

    2.2427e-08    3.8634e-08

    5.3147e-08    1.0190e-07

I'm inclined to agree that the recursion form should work better. I'm suspicious that Matlab's version is reliable for such high order legendre polynomials.

Anyone, is there a reliable method to verify the correct answers?

Ben

_______________________________________________
Help-octave mailing list
Help-octave@...
https://www.cae.wisc.edu/mailman/listinfo/help-octave

[Dmitri A. Sergatskov]

FWIW: Marko's script seems to agree with "legendre_Plm"
wrapper of GSL's function (gsl package from octave-forge).

Sincerely,

Dmitri.
--
_______________________________________________
Help-octave mailing list
Help-octave@...
https://www.cae.wisc.edu/mailman/listinfo/help-octave

[Marco Caliari-4]

> I compared Matlab's script and Octave's and your proposal.
>
> result_matlab = legendre (80, [-1:0.1:1]);
>
> result_octave0 respects the original and result_octave1 respects your
> algorithm.
>
> I calculated an error with respect to the matlab version (I'm not sure
> Matlab's is to be trusted as correct in all cases).
>
> d0 = abs (result_matlab - result_octave0) / abs (result_matlab);
> d1 = abs (result_matlab - result_octave1) / abs (result_matlab);
> er1 = max (d1, [], 2);
> er0 = max (d0, [], 2);
>
> [er0(:), er1(:)] produces the following

Dear Ben,

the results are even more impressive if you consider that
legendre0(80,[-1:0.1:1])(1,1) gives 6.7015e+14, whereas
legendre1(80,[-1:0.1:1])(1,1) gives 1.

> I'm inclined to agree that the recursion form should work better. I'm
> suspicious that Matlab's version is reliable for such high order legendre
> polynomials.
>
> Anyone, is there a reliable method to verify the correct answers?

There is a LegendreP function in Maple doing exactly the same. Can you
select few cases (a degree and a scalar value for x) for which the
difference between Matlab and my script is quite large? I will check
them with Maple.

Best regards,

Marco
_______________________________________________
Help-octave mailing list
Help-octave@...
https://www.cae.wisc.edu/mailman/listinfo/help-octave

[Ben Abbott]

On Feb 12, 2008, at 1:44 AM, Dmitri A. Sergatskov wrote:

FWIW: Marko's script seems to agree with "legendre_Plm"
wrapper of GSL's function (gsl package from octave-forge).

Sincerely,

Dmitri.
--

That is good news!

It would be nice to do a direct comparison for the example below.

result_matlab = legendre (80, [-1:0.1:1]);

Any chance anyone who knows c/c++ has the inclination to write a short program using gsl to calculate that example and post the results?

In addition, I noticed the comment below in the reference manual for gsl,

"The following functions compute the associated Legendre Polynomials P_l^m(x). Note that this function grows combinatorially with l and can overflow for l larger than about 150. There is no trouble for small m, but overflow occurs when m and l are both large. Rather than allow overflows, these functions refuse to calculate P_l^m(x) and return GSL_EOVRFLW when they can sense that l and m are too big."

It would be more proper to restrict the issuing of the warning for the cases where overflow actually occurs. For example, even this trivial calculation warns of overflow.

octave:73> legendre(0,0)

warning: legendre is unstable for higher orders

ans =  1

I'll take a look at modifying the script to catch overflows and issue the warning only when such occurs.

Ben

_______________________________________________
Help-octave mailing list
Help-octave@...
https://www.cae.wisc.edu/mailman/listinfo/help-octave

[Ben Abbott]

On Feb 12, 2008, at 3:06 AM, Marco Caliari wrote:

>> I compared Matlab's script and Octave's and your proposal.
>>
>> result_matlab = legendre (80, [-1:0.1:1]);
>>
>> result_octave0 respects the original and result_octave1 respects  
>> your algorithm.
>>
>> I calculated an error with respect to the matlab version (I'm not  
>> sure Matlab's is to be trusted as correct in all cases).
>>
>> d0 = abs (result_matlab - result_octave0) / abs (result_matlab);
>> d1 = abs (result_matlab - result_octave1) / abs (result_matlab);
>> er1 = max (d1, [], 2);
>> er0 = max (d0, [], 2);
>>
>> [er0(:), er1(:)] produces the following
>
> Dear Ben,
>
> the results are even more impressive if you consider that  
> legendre0(80,[-1:0.1:1])(1,1) gives 6.7015e+14, whereas legendre1(80,
> [-1:0.1:1])(1,1) gives 1.
>
>> I'm inclined to agree that the recursion form should work better.  
>> I'm suspicious that Matlab's version is reliable for such high  
>> order legendre polynomials.
>>
>> Anyone, is there a reliable method to verify the correct answers?
>
> There is a LegendreP function in Maple doing exactly the same. Can  
> you select few cases (a degree and a scalar value for x) for which  
> the difference between Matlab and my script is quite large? I will  
> check them with Maple.

I haven't used Maple in several years. Is it not possible to evaluate  
the example, legendre (80, [-1:0.1:1]),  and post the results (as an  
attachment)?

It can then read it as ...

        result_maple = load ("result_maple.txt");

and compared directly to octave and Matlab. I realize that we'll loose  
some precision, but if 16 digits are preserved, we should obtain  
relative errors on the order 10^(-15).

I notice Maxima also has a legendre function, legendre_p(n,m,x). Which  
uses `Abramowitz and Stegun, Handbook of Mathematical Functions' as  
the reference. I'll take a look at calculating the example there.

Ben
_______________________________________________
Help-octave mailing list
Help-octave@...
https://www.cae.wisc.edu/mailman/listinfo/help-octave

[Marco Caliari-4]

Dear all,

please consider the cleaner enclosed script.

I tried to run Ben's example in Maple for ten minutes, with no success. I
don't know how to give a vector in input to LegendreP.

Best regards,

Marco

On Tue, 12 Feb 2008, Ben Abbott wrote:

>
> On Feb 12, 2008, at 1:44 AM, Dmitri A. Sergatskov wrote:
>
>> FWIW: Marko's script seems to agree with "legendre_Plm"
>> wrapper of GSL's function (gsl package from octave-forge).
>>
>> Sincerely,
>>
>> Dmitri.
>> --
>
>
> That is good news!
>
> It would be nice to do a direct comparison for the example below.
>
> result_matlab = legendre (80, [-1:0.1:1]);
>
> Any chance anyone who knows c/c++ has the inclination to write a short
> program using gsl to calculate that example and post the results?
>
> In addition, I noticed the comment below in the reference manual for gsl,
>
> "The following functions compute the associated Legendre Polynomials
> P_l^m(x). Note that this function grows combinatorially with l and can
> overflow for l larger than about 150. There is no trouble for small m, but
> overflow occurs when m and l are both large. Rather than allow overflows,
> these functions refuse to calculate P_l^m(x) and return GSL_EOVRFLW when they
> can sense that l and m are too big."
>
> It would be more proper to restrict the issuing of the warning for the cases
> where overflow actually occurs. For example, even this trivial calculation
> warns of overflow.
>
> octave:73> legendre(0,0)
> warning: legendre is unstable for higher orders
> ans =  1
>
> I'll take a look at modifying the script to catch overflows and issue the
> warning only when such occurs.
>
> Ben

## 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..@var{n} are returned.
## @var{n} must be an integer.
## 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@...>
function L = legendre(n,x,varargin)
x = x(:).';
L = zeros(n+1,length(x));
if (nargin == 2)
  normalization = "unnorm";
else
  normalization = varargin{1};
end
if (strcmp(normalization,"norm"))
  scale = sqrt(n+0.5);
elseif (strcmp(normalization,"sch"))
  scale = sqrt(2);
elseif (strcmp(normalization,"unnorm"))
  scale = 1;
else
  error("Normalization option not recognized.")
end
% Based on the first recurrence relation found in
% http://en.wikipedia.org/wiki/Associated_Legendre_function
for m = 1:n
  LP = ones(n+1,length(x))*scale;
  LP(m,:) = prod(1:2:2*(m-1)-1)*LP(m,:).*sqrt(1-x.^2).^(m-1);
  if (m == n+1)
    L(m,:) = LP(m,:);
    break
  end  
  LP(m+1,:) = (2*m-1).*x.*LP(m,:);
  for l = m+1:n
    LP(l+1,:) = ((2*(l-1)+1).*x.*LP(l,:)-(l-1+m-1).*LP(l-1,:))./(l-m+1);
  end
  L(m,:) = LP(n+1,:);
  if (strcmp(normalization,"norm"))
    scale = scale/sqrt((n-m+1)*(n+m));
  elseif (strcmp(normalization,"sch"))
    scale = scale/sqrt((n-m+1)*(n+m));
  else
    scale = -scale;
  end
end
L(n+1,:) = prod(1:2:2*n-1)*scale.*sqrt(1-x.^2).^n;
if (strcmp(normalization,"sch"))
  L(1,:) = L(1,:)/sqrt(2);
end
_______________________________________________
Help-octave mailing list
Help-octave@...
https://www.cae.wisc.edu/mailman/listinfo/help-octave

[Marco Caliari-4]

Hi.

The normalized Lagrange functions should (almost) never give NaNs of Infs.
The enclosed script fixes a problem in the previous:

octave:1> legendreold(151,-0.9,"norm")(end-1:end)
ans =

   -3.3248e-53
           Inf
octave:2> legendre(151,-0.9,"norm")(end-1:end)
ans =

   -3.3248e-53
    9.2660e-55

Marco

> "The following functions compute the associated Legendre Polynomials
> P_l^m(x). Note that this function grows combinatorially with l and can
> overflow for l larger than about 150. There is no trouble for small m, but
> overflow occurs when m and l are both large. Rather than allow overflows,
> these functions refuse to calculate P_l^m(x) and return GSL_EOVRFLW when they
> can sense that l and m are too big."
## 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..@var{n} are returned.
## @var{n} must be an integer.
## 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@...>
function L = legendre(n,x,varargin)
x = x(:).';
L = zeros(n+1,length(x));
if (nargin == 2)
  normalization = "unnorm";
else
  normalization = varargin{1};
end
if (strcmp(normalization,"norm"))
  scale = sqrt(n+0.5);
elseif (strcmp(normalization,"sch"))
  scale = sqrt(2);
elseif (strcmp(normalization,"unnorm"))
  scale = 1;
else
  error("Normalization option not recognized.")
end
% Based on the first recurrence relation found in
% http://en.wikipedia.org/wiki/Associated_Legendre_function
for m = 1:n
  LP = ones(n+1,length(x));
  LP(m,:) = scale*LP(m,:).*sqrt(1-x.^2).^(m-1);
  LP(m+1,:) = (2*m-1).*x.*LP(m,:);
  for l = m+1:n
    LP(l+1,:) = ((2*l-1).*x.*LP(l,:)-(l+m-2).*LP(l-1,:))./(l-m+1);
  end
  L(m,:) = LP(n+1,:);
  if (strcmp(normalization,"unnorm"))
    scale = -scale*(2*m-1);
  else # normalization == "sch" or normalization == "norm"
    scale = scale/sqrt((n-m+1)*(n+m))*(2*m-1);
  end
end
L(n+1,:) = scale.*sqrt(1-x.^2).^n;
if (strcmp(normalization,"sch"))
  L(1,:) = L(1,:)/sqrt(2);
end
_______________________________________________
Help-octave mailing list
Help-octave@...
https://www.cae.wisc.edu/mailman/listinfo/help-octave

[Ben Abbott]

On Tuesday, February 12, 2008, at 10:02AM, "Marco Caliari" <marco.caliari@...> wrote:

>Hi.
>
>The normalized Lagrange functions should (almost) never give NaNs of Infs.
>The enclosed script fixes a problem in the previous:
>
>octave:1> legendreold(151,-0.9,"norm")(end-1:end)
>ans =
>
>   -3.3248e-53
>           Inf
>octave:2> legendre(151,-0.9,"norm")(end-1:end)
>ans =
>
>   -3.3248e-53
>    9.2660e-55
>
>Marco
>

I've tried to combine the various implementations. I still need to look into adding a warning in then event of an overflow.

We should also add more tests to include the current improvements.

Ben



_______________________________________________
Help-octave mailing list
Help-octave@...
https://www.cae.wisc.edu/mailman/listinfo/help-octave

[Ben Abbott]

On Feb 12, 2008, at 1:38 PM, Ben Abbott wrote:

> On Tuesday, February 12, 2008, at 10:02AM, "Marco Caliari" <marco.caliari@...
> > wrote:
>> Hi.
>>
>> The normalized Lagrange functions should (almost) never give NaNs  
>> of Infs.
>> The enclosed script fixes a problem in the previous:
>>
>> octave:1> legendreold(151,-0.9,"norm")(end-1:end)
>> ans =
>>
>>  -3.3248e-53
>>          Inf
>> octave:2> legendre(151,-0.9,"norm")(end-1:end)
>> ans =
>>
>>  -3.3248e-53
>>   9.2660e-55
>>
>> Marco
>>
>
> I've tried to combine the various implementations. I still need to  
> look into adding a warning in then event of an overflow.
>
> We should also add more tests to include the current improvements.
>
> Ben

FYI,

I asked Jason Riedy to run the example below for me with his patch/
workaround for xgelsd

        legendre  (80, [-1:0.1:1]);

I have been getting the warning below.

        warning: dgelsd: rank deficient 21x81 matrix, rank = 7

With the patch/workaround Jason did not get the warning.

Ben

_______________________________________________
Help-octave mailing list
Help-octave@...
https://www.cae.wisc.edu/mailman/listinfo/help-octave

[Ben Abbott]

I've added a check for underflow. This is my first time doing such, so  
I may be quite off on what is commonly accepted as proper.

The method I used detects an underflow for orders greater than 151.

        octave:104> result_octave1 = legendre  (152, [-1:0.1:1]);
        warning: legendre: results may be  unstable for high orders.

I've attached the script.

Ben






On Feb 12, 2008, at 1:38 PM, Ben Abbott wrote:

> On Tuesday, February 12, 2008, at 10:02AM, "Marco Caliari" <marco.caliari@...
> > wrote:
>> Hi.
>>
>> The normalized Lagrange functions should (almost) never give NaNs  
>> of Infs.
>> The enclosed script fixes a problem in the previous:
>>
>> octave:1> legendreold(151,-0.9,"norm")(end-1:end)
>> ans =
>>
>>  -3.3248e-53
>>          Inf
>> octave:2> legendre(151,-0.9,"norm")(end-1:end)
>> ans =
>>
>>  -3.3248e-53
>>   9.2660e-55
>>
>> Marco
>>
>
> I've tried to combine the various implementations. I still need to  
> look into adding a warning in then event of an overflow.
>
> We should also add more tests to include the current improvements.
>
> Ben
>
> <legendre.m>_______________________________________________
> Help-octave mailing list
> Help-octave@...
> https://www.cae.wisc.edu/mailman/listinfo/help-octave

_______________________________________________
Help-octave mailing list
Help-octave@...
https://www.cae.wisc.edu/mailman/listinfo/help-octave

[Marco Caliari-4]

Hi Ben.

This one still gives problems

legendre(151,0)

with an Inf. I would check the line

scale = -scale * (2*m-1);

and try to prevent that scale becomes Inf.
Moreover, I don't understand the limit 255 on n. With the normalizations,
there are no problems for n>255.

Best regards,

Marco

On Tue, 12 Feb 2008, Ben Abbott wrote:

> I've added a check for underflow. This is my first time doing such, so I may
> be quite off on what is commonly accepted as proper.
>
> The method I used detects an underflow for orders greater than 151.
>
> octave:104> result_octave1 = legendre  (152, [-1:0.1:1]);
> warning: legendre: results may be  unstable for high orders.
>
> I've attached the script.
>
> Ben
>

_______________________________________________
Help-octave mailing list
Help-octave@...
https://www.cae.wisc.edu/mailman/listinfo/help-octave

[Ben Abbott]

On Feb 13, 2008, at 6:47 AM, Marco Caliari wrote:

> Hi Ben.
>
> This one still gives problems
>
> legendre(151,0)
>
> with an Inf. I would check the line
>
> scale = -scale * (2*m-1);

That line was incorporated from your version. Should it be something  
else?

> and try to prevent that scale becomes Inf.
> Moreover, I don't understand the limit 255 on n. With the  
> normalizations, there are no problems for n>255.

Agreed.

I'll post an update in couple of hours.

Ben
_______________________________________________
Help-octave mailing list
Help-octave@...
https://www.cae.wisc.edu/mailman/listinfo/help-octave

[Marco Caliari-4]

> 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

>> and try to prevent that scale becomes Inf.
>> Moreover, I don't understand the limit 255 on n. With the normalizations,
>> there are no problems for n>255.
>
> Agreed.

The enclosed version has a further improvement:

# the warning when scale is Inf
octave:4> legendre(151,0)(end)
warning: legendre: results may be  unstable for high orders.
warning: legendre: results may be  unstable for high orders.
ans = -Inf
octave:5> legendreold(151,0)(end)
ans = -Inf

# NaN prevented in case of Inf*0
octave:8> legendre(151,-1)(end-1:end)
ans =

   -0
   -0
octave:9> legendreold(151,-1)(end-1:end)
ans =

     -0
    NaN

# Inf prevented in case of Inf*(small quantity)
octave:10> legendre(151,-0.9)(end-1:end)
ans =

   -8.1986e+254
   -3.9708e+254
octave:11> legendreold(151,-0.9)(end-1:end)
ans =

   -8.1986e+254
           -Inf

Best regards,

Marco
## 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..@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@...>

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 > 255 || 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"
      scale = sqrt (n+0.5);
    case "sch"
      scale = sqrt (2);
    case "unnorm"
      scale = 1;
    otherwise
      print_usage ();
  endswitch
  ## 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
    LP = ones (n+1, length (x));
    LP(m,:) = scale .* LP(m,:);
    LP(m+1,:) = (2*m-1) .* x .* LP(m,:);
    for l = m+1:n
      a = (2*l-1) .* x .* LP(l,:);
      b = (l+m-2) .* LP(l-1,:);
      if (any((a-b) == a & b != 0) || any((a-b) == -b & a != 0))
        a-b
        underflow = true;
      endif
      LP(l+1,:) = (a-b)./(l-m+1);
      ##LP(l+1,:) = ((2*l-1).*x.*LP(l,:)-(l+m-2).*LP(l-1,:))./(l-m+1);
    endfor
    L(m,:) = LP(n+1,:);
    if (strcmp (normalization, "unnorm"))
      scale = -scale * (2*m-1) .* sqrt (1-x.^2);
      if (max(abs(scale)) == Inf)
        underflow = true;
      endif
    else
      ## normalization == "sch" or normalization == "norm"
      scale = scale / sqrt ((n-m+1)*(n+m))*(2*m-1) .* sqrt (1-x.^2);
    endif
  endfor
  L(n+1,:) = scale;
  if (strcmp (normalization, "sch"))
    L(1,:) = L(1,:) / sqrt (2);
  endif
  if (underflow)
    warning ("legendre: results may be  unstable for high orders.");
  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);



_______________________________________________
Help-octave mailing list
Help-octave@...
https://www.cae.wisc.edu/mailman/listinfo/help-octave

[Ben Abbott]

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






_______________________________________________
Help-octave mailing list
Help-octave@...
https://www.cae.wisc.edu/mailman/listinfo/help-octave

[Marco Caliari-4]

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

## 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..@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@...>

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



_______________________________________________
Help-octave mailing list
Help-octave@...
https://www.cae.wisc.edu/mailman/listinfo/help-octave