octave_packages/specfun-1.1.0/laguerre.m

   1 ## Copyright (C) 2008 Eric Chassande-Mottin
   2 ##
   3 ## This program is free software; you can redistribute it and/or modify
   4 ## it under the terms of the GNU General Public License as published by
   5 ## the Free Software Foundation; either version 3 of the License, or
   6 ## (at your option) any later version.
   7 ##
   8 ## This program is distributed in the hope that it will be useful,
   9 ## but WITHOUT ANY WARRANTY; without even the implied warranty of
  10 ## MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
  11 ## GNU General Public License for more details.
  12 ##
  13 ## You should have received a copy of the GNU General Public License
  14 ## along with this program; If not, see <http://www.gnu.org/licenses/>.
  15
  16 ## -*- texinfo -*-
  17 ## @deftypefn {Function File} {@var{y} = } laguerre (@var{x},@var{n})
  18 ## @deftypefnx {Function File} {[@var{y} @var{p}]= } laguerre (@var{x},@var{n})
  19 ##
  20 ## Compute the value of the Laguerre polynomial of order @var{n} for each element of @var{x}
  21 ##
  22 ## @end deftypefn
  23
  24 function [y,p]=laguerre(x,n)
  25
  26   if (nargin != 2)
  27     print_usage;
  28   elseif (n < 0 || !isscalar (n))
  29     error("second argument 'n' must be a positive integer");
  30   endif
  31
  32   p0=1;
  33   p1=[-1 1];
  34
  35   if (n==0)
  36     p=p0;
  37   elseif (n==1)
  38     p=p1;
  39   elseif (n > 1)
  40     % recursive calculation of the polynomial coefficients
  41     for k=2:n
  42       p=zeros(1,k+1);
  43       p(1) = -p1(1)/k;
  44       p(2) = ((2*k-1)*p1(1)-p1(2))/k;
  45       if (k > 2)
  46         p(3:k) = ((2*k-1)*p1(2:k-1)-p1(3:k)-(k-1)*p0(1:k-2))/k;
  47       endif
  48       p(k+1) = ((2*k-1)*p1(k)-(k-1)*p0(k-1))/k;
  49       p0=p1;
  50       p1=p;
  51     endfor
  52   endif
  53
  54   y=polyval(p,x);
  55
  56 endfunction
  57
  58 %!test
  59 %! x=rand;
  60 %! y1=laguerre(x,0);
  61 %! p0=[1];
  62 %! y2=polyval(p0,x);
  63 %! assert(y1-y2,0,eps);
  64
  65 %!test
  66 %! x=rand;
  67 %! y1=laguerre(x,1);
  68 %! p1=[-1 1];
  69 %! y2=polyval(p1,x);
  70 %! assert(y1-y2,0,eps);
  71
  72 %!test
  73 %! x=rand;
  74 %! y1=laguerre(x,2);
  75 %! p2=[.5 -2 1];
  76 %! y2=polyval(p2,x);
  77 %! assert(y1-y2,0,eps);
  78
  79 %!test
  80 %! x=rand;
  81 %! y1=laguerre(x,3);
  82 %! p3=[-1/6 9/6 -18/6 1];
  83 %! y2=polyval(p3,x);
  84 %! assert(y1-y2,0,eps);
  85
  86 %!test
  87 %! x=rand;
  88 %! y1=laguerre(x,4);
  89 %! p4=[1/24 -16/24 72/24 -96/24 1];
  90 %! y2=polyval(p4,x);
  91 %! assert(y1-y2,0,eps);