Add a useful package (from Source forge) for octave

[CreaPhase.git] / octave_packages / specfun-1.1.0 / laguerre.m

diff --git a/octave_packages/specfun-1.1.0/laguerre.m b/octave_packages/specfun-1.1.0/laguerre.m
new file mode 100644 (file)
index 0000000..a6b5f36
--- /dev/null
+++ b/octave_packages/specfun-1.1.0/laguerre.m
@@ -0,0 +1,91 @@
+## Copyright (C) 2008 Eric Chassande-Mottin
+##
+## This program 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.
+##
+## This program 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 this program; If not, see <http://www.gnu.org/licenses/>.
+
+## -*- texinfo -*-
+## @deftypefn {Function File} {@var{y} = } laguerre (@var{x},@var{n})
+## @deftypefnx {Function File} {[@var{y} @var{p}]= } laguerre (@var{x},@var{n})
+##
+## Compute the value of the Laguerre polynomial of order @var{n} for each element of @var{x}
+##
+## @end deftypefn
+
+function [y,p]=laguerre(x,n)
+
+  if (nargin != 2)
+    print_usage;
+  elseif (n < 0 || !isscalar (n))
+    error("second argument 'n' must be a positive integer");
+  endif
+
+  p0=1;
+  p1=[-1 1];
+
+  if (n==0)
+    p=p0;
+  elseif (n==1)
+    p=p1;
+  elseif (n > 1)
+    % recursive calculation of the polynomial coefficients
+    for k=2:n
+      p=zeros(1,k+1);
+      p(1) = -p1(1)/k;
+      p(2) = ((2*k-1)*p1(1)-p1(2))/k;
+      if (k > 2)
+        p(3:k) = ((2*k-1)*p1(2:k-1)-p1(3:k)-(k-1)*p0(1:k-2))/k;
+      endif
+      p(k+1) = ((2*k-1)*p1(k)-(k-1)*p0(k-1))/k;
+      p0=p1;
+      p1=p;
+    endfor
+  endif
+
+  y=polyval(p,x);
+
+endfunction
+
+%!test
+%! x=rand;
+%! y1=laguerre(x,0); 
+%! p0=[1]; 
+%! y2=polyval(p0,x);
+%! assert(y1-y2,0,eps);
+
+%!test
+%! x=rand;
+%! y1=laguerre(x,1); 
+%! p1=[-1 1]; 
+%! y2=polyval(p1,x);
+%! assert(y1-y2,0,eps);
+
+%!test
+%! x=rand;
+%! y1=laguerre(x,2); 
+%! p2=[.5 -2 1];
+%! y2=polyval(p2,x);
+%! assert(y1-y2,0,eps);
+
+%!test
+%! x=rand;
+%! y1=laguerre(x,3); 
+%! p3=[-1/6 9/6 -18/6 1];
+%! y2=polyval(p3,x);
+%! assert(y1-y2,0,eps);
+
+%!test
+%! x=rand;
+%! y1=laguerre(x,4); 
+%! p4=[1/24 -16/24 72/24 -96/24 1];
+%! y2=polyval(p4,x);
+%! assert(y1-y2,0,eps);