octave_packages/m/specfun/legendre.m

   1 ## Copyright (C) 2000-2012 Kai Habel
   2 ## Copyright (C) 2008 Marco Caliari
   3 ##
   4 ## This file is part of Octave.
   5 ##
   6 ## Octave is free software; you can redistribute it and/or modify it
   7 ## under the terms of the GNU General Public License as published by
   8 ## the Free Software Foundation; either version 3 of the License, or (at
   9 ## your option) any later version.
  10 ##
  11 ## Octave is distributed in the hope that it will be useful, but
  12 ## WITHOUT ANY WARRANTY; without even the implied warranty of
  13 ## MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the GNU
  14 ## General Public License for more details.
  15 ##
  16 ## You should have received a copy of the GNU General Public License
  17 ## along with Octave; see the file COPYING.  If not, see
  18 ## <http://www.gnu.org/licenses/>.
  19
  20 ## -*- texinfo -*-
  21 ## @deftypefn  {Function File} {@var{l} =} legendre (@var{n}, @var{x})
  22 ## @deftypefnx {Function File} {@var{l} =} legendre (@var{n}, @var{x}, @var{normalization})
  23 ## Compute the Legendre function of degree @var{n} and order
  24 ## @var{m} = 0 @dots{} N@.  The optional argument, @var{normalization},
  25 ## may be one of @code{"unnorm"}, @code{"sch"}, or @code{"norm"}.
  26 ## The default is @code{"unnorm"}.  The value of @var{n} must be a
  27 ## non-negative scalar integer.
  28 ##
  29 ## If the optional argument @var{normalization} is missing or is
  30 ## @code{"unnorm"}, compute the Legendre function of degree @var{n} and
  31 ## order @var{m} and return all values for @var{m} = 0 @dots{} @var{n}.
  32 ## The return value has one dimension more than @var{x}.
  33 ##
  34 ## The Legendre Function of degree @var{n} and order @var{m}:
  35 ##
  36 ## @tex
  37 ## $$
  38 ## P^m_n(x) = (-1)^m (1-x^2)^{m/2}{d^m\over {dx^m}}P_n (x)
  39 ## $$
  40 ## @end tex
  41 ## @ifnottex
  42 ##
  43 ## @example
  44 ## @group
  45 ##  m        m       2  m/2   d^m
  46 ## P(x) = (-1) * (1-x  )    * ----  P(x)
  47 ##  n                         dx^m   n
  48 ## @end group
  49 ## @end example
  50 ##
  51 ## @end ifnottex
  52 ##
  53 ## @noindent
  54 ## with Legendre polynomial of degree @var{n}:
  55 ##
  56 ## @tex
  57 ## $$
  58 ## P(x) = {1\over{2^n n!}}\biggl({d^n\over{dx^n}}(x^2 - 1)^n\biggr)
  59 ## $$
  60 ## @end tex
  61 ## @ifnottex
  62 ##
  63 ## @example
  64 ## @group
  65 ##           1    d^n   2    n
  66 ## P(x) = ------ [----(x - 1) ]
  67 ##  n     2^n n!  dx^n
  68 ## @end group
  69 ## @end example
  70 ##
  71 ## @end ifnottex
  72 ##
  73 ## @noindent
  74 ## @code{legendre (3, [-1.0, -0.9, -0.8])} returns the matrix:
  75 ##
  76 ## @example
  77 ## @group
  78 ##  x  |   -1.0   |   -0.9   |   -0.8
  79 ## ------------------------------------
  80 ## m=0 | -1.00000 | -0.47250 | -0.08000
  81 ## m=1 |  0.00000 | -1.99420 | -1.98000
  82 ## m=2 |  0.00000 | -2.56500 | -4.32000
  83 ## m=3 |  0.00000 | -1.24229 | -3.24000
  84 ## @end group
  85 ## @end example
  86 ##
  87 ## If the optional argument @code{normalization} is @code{"sch"},
  88 ## compute the Schmidt semi-normalized associated Legendre function.
  89 ## The Schmidt semi-normalized associated Legendre function is related
  90 ## to the unnormalized Legendre functions by the following:
  91 ##
  92 ## For Legendre functions of degree n and order 0:
  93 ##
  94 ## @tex
  95 ## $$
  96 ## SP^0_n (x) = P^0_n (x)
  97 ## $$
  98 ## @end tex
  99 ## @ifnottex
 100 ##
 101 ## @example
 102 ## @group
 103 ##   0      0
 104 ## SP(x) = P(x)
 105 ##   n      n
 106 ## @end group
 107 ## @end example
 108 ##
 109 ## @end ifnottex
 110 ##
 111 ## For Legendre functions of degree n and order m:
 112 ##
 113 ## @tex
 114 ## $$
 115 ## SP^m_n (x) = P^m_n (x)(-1)^m\biggl({2(n-m)!\over{(n+m)!}}\biggl)^{0.5}
 116 ## $$
 117 ## @end tex
 118 ## @ifnottex
 119 ##
 120 ## @example
 121 ## @group
 122 ##   m      m         m    2(n-m)! 0.5
 123 ## SP(x) = P(x) * (-1)  * [-------]
 124 ##   n      n              (n+m)!
 125 ## @end group
 126 ## @end example
 127 ##
 128 ## @end ifnottex
 129 ##
 130 ## If the optional argument @var{normalization} is @code{"norm"},
 131 ## compute the fully normalized associated Legendre function.
 132 ## The fully normalized associated Legendre function is related
 133 ## to the unnormalized Legendre functions by the following:
 134 ##
 135 ## For Legendre functions of degree @var{n} and order @var{m}
 136 ##
 137 ## @tex
 138 ## $$
 139 ## NP^m_n (x) = P^m_n (x)(-1)^m\biggl({(n+0.5)(n-m)!\over{(n+m)!}}\biggl)^{0.5}
 140 ## $$
 141 ## @end tex
 142 ## @ifnottex
 143 ##
 144 ## @example
 145 ## @group
 146 ##   m      m         m    (n+0.5)(n-m)! 0.5
 147 ## NP(x) = P(x) * (-1)  * [-------------]
 148 ##   n      n                  (n+m)!
 149 ## @end group
 150 ## @end example
 151 ##
 152 ## @end ifnottex
 153 ## @end deftypefn
 154
 155 ## Author: Marco Caliari <marco.caliari@univr.it>
 156
 157 function retval = legendre (n, x, normalization)
 158
 159   persistent warned_overflow = false;
 160
 161   if (nargin < 2 || nargin > 3)
 162     print_usage ();
 163   endif
 164
 165   if (!isscalar (n) || n < 0 || n != fix (n))
 166     error ("legendre: N must be a non-negative scalar integer");
 167   endif
 168
 169   if (!isreal (x) || any (x(:) < -1 | x(:) > 1))
 170     error ("legendre: X must be real-valued vector in the range -1 <= X <= 1");
 171   endif
 172
 173   if (nargin == 3)
 174     normalization = lower (normalization);
 175   else
 176     normalization = "unnorm";
 177   endif
 178
 179   switch (normalization)
 180     case "norm"
 181       scale = sqrt (n+0.5);
 182     case "sch"
 183       scale = sqrt (2);
 184     case "unnorm"
 185       scale = 1;
 186     otherwise
 187       error ('legendre: expecting NORMALIZATION option to be "norm", "sch", or "unnorm"');
 188   endswitch
 189
 190   scale = scale * ones (size (x));
 191
 192   ## Based on the recurrence relation below
 193   ##            m                 m              m
 194   ## (n-m+1) * P (x) = (2*n+1)*x*P (x)  - (n+1)*P (x)
 195   ##            n+1               n              n-1
 196   ## http://en.wikipedia.org/wiki/Associated_Legendre_function
 197
 198   overflow = false;
 199   retval = zeros([n+1, size(x)]);
 200   for m = 1:n
 201     lpm1 = scale;
 202     lpm2 = (2*m-1) .* x .* scale;
 203     lpm3 = lpm2;
 204     for k = m+1:n
 205       lpm3a = (2*k-1) .* x .* lpm2;
 206       lpm3b = (k+m-2) .* lpm1;
 207       lpm3 = (lpm3a - lpm3b)/(k-m+1);
 208       lpm1 = lpm2;
 209       lpm2 = lpm3;
 210       if (! warned_overflow)
 211         if (any (abs (lpm3a) > realmax)
 212             || any (abs (lpm3b) > realmax)
 213             || any (abs (lpm3)  > realmax))
 214           overflow = true;
 215         endif
 216       endif
 217     endfor
 218     retval(m,:) = lpm3(:);
 219     if (strcmp (normalization, "unnorm"))
 220       scale = -scale * (2*m-1);
 221     else
 222       ## normalization == "sch" or normalization == "norm"
 223       scale = scale / sqrt ((n-m+1)*(n+m))*(2*m-1);
 224     endif
 225     scale = scale .* sqrt(1-x.^2);
 226   endfor
 227
 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
 234
 235   if (strcmp (normalization, "sch"))
 236     retval(1,:) = retval(1,:) / sqrt (2);
 237   endif
 238
 239   if (overflow && ! warned_overflow)
 240     warning ("legendre: overflow - results may be unstable for high orders");
 241     warned_overflow = true;
 242   endif
 243
 244 endfunction
 245
 246
 247 %!test
 248 %! result = legendre (3, [-1.0 -0.9 -0.8]);
 249 %! expected = [
 250 %!    -1.00000  -0.47250  -0.08000
 251 %!     0.00000  -1.99420  -1.98000
 252 %!     0.00000  -2.56500  -4.32000
 253 %!     0.00000  -1.24229  -3.24000
 254 %! ];
 255 %! assert (result, expected, 1e-5);
 256
 257 %!test
 258 %! result = legendre (3, [-1.0 -0.9 -0.8], "sch");
 259 %! expected = [
 260 %!    -1.00000  -0.47250  -0.08000
 261 %!     0.00000   0.81413   0.80833
 262 %!    -0.00000  -0.33114  -0.55771
 263 %!     0.00000   0.06547   0.17076
 264 %! ];
 265 %! assert (result, expected, 1e-5);
 266
 267 %!test
 268 %! result = legendre (3, [-1.0 -0.9 -0.8], "norm");
 269 %! expected = [
 270 %!    -1.87083  -0.88397  -0.14967
 271 %!     0.00000   1.07699   1.06932
 272 %!    -0.00000  -0.43806  -0.73778
 273 %!     0.00000   0.08661   0.22590
 274 %! ];
 275 %! assert (result, expected, 1e-5);
 276
 277 %!test
 278 %! result = legendre (151, 0);
 279 %! ## Don't compare to "-Inf" since it would fail on 64 bit systems.
 280 %! assert (result(end) < -1.7976e308 && all (isfinite (result(1:end-1))));
 281
 282 %!test
 283 %! result = legendre (150, 0);
 284 %! ## This agrees with Matlab's result.
 285 %! assert (result(end), 3.7532741115719e+306, 0.0000000000001e+306);
 286
 287 %!test
 288 %! result = legendre (0, 0:0.1:1);
 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);
 304
 305 %% Check correct invocation
 306 %!error legendre ();
 307 %!error legendre (1);
 308 %!error legendre (1,2,3,4);
 309 %!error legendre ([1, 2], [-1, 0, 1]);
 310 %!error legendre (-1, [-1, 0, 1]);
 311 %!error legendre (1.1, [-1, 0, 1]);
 312 %!error legendre (1, [-1+i, 0, 1]);
 313 %!error legendre (1, [-2, 0, 1]);
 314 %!error legendre (1, [-1, 0, 2]);
 315 %!error legendre (1, [-1, 0, 1], "badnorm");