MATLAB: Convergence of Laguerre Function

Question

Hi, I need calculate the lagauerre function defined by :

fi(x)=exp(-x/2)* Li(x) for various order i when Li(x) is the Laguerre polynomials ( http://mathworld.wolfram.com/LaguerrePolynomial.html ) given by the following function:

function L = Funlaguerre(n,x)
sum= 0;
for i=0:n
sum = sum + ((i^(-1)* (factorial(n)/(factoria(l(n-i) * factorial(i) * factorial(i))).*x^k));
end
L= sum;
end

According to the theory , Laguerre function must converge to 0 when the order of the Laguerre polynomial is high. I calculate the laguerre function based on the following parameters

Tf= 1e-007;
M=101;
delta_t=Tf/100;
T=0:delta_t:(M-1)*delta_t;
S=10^9; % scaling Factor
lag=60; % order of Laguerre polynomial
for n=0:lag
    for t=1:M
F(n+1,t)=exp(-s*T(t)/2)* Funlaguerre (n,s*T(t));
    end
end

I obtain correct results for the order : [1…30] . But, when the order becomes more than 30, the laguerre function diverges as seen in the attached curves. i don't understand why the function diverges? best regards

%

Answer 1

Hamdi, Your trouble is entirely due to round-off errors whose effect becomes large for large order Laguerre polynomials, together with the fact that you didn't let x get large enough to see that your function was actually converging to zero. The multiplicative factor exp(-x/2) guarantees that it must eventually converge to zero. However, the higher order polynomials will require larger values of x before their convergence to zero becomes evident.

There is a method of computation of Laguerre polynomials using recursion that will greatly reduce the noise due to round-off. It is:

 L(n,x) = ((2*n-1-x)*L(n-1,x)-(n-1)*L(n-2,x))/n;

See: http://mathworld.wolfram.com/LaguerrePolynomial.html

Here is code I used to generate your function using this recursion:

 clear
 n = 100; % <-- Select the order
 N = 8192; % <-- Choose the number of points to plot
 X = linspace(0,500,N); % <-- Choose the range of x
 F = zeros(1,N);
 for ix = 1:N
   x = X(ix);
   t0 = 1; % <-- The zero-th order polynomial
   t1 = 1-x; % <-- The first order polynomial
   for k = 2:n
     t = ((2*k-1-x)*t1-(k-1)*t0)/k; % <-- The recursion formula
     t0 = t1; t1 = t;
   end
   F(ix) = t1*exp(-x/2);
 end
 plot(X,F,'y-')