Code covered by the BSD License  

Highlights from
Bessel Polynomial Linearization Coefficients

image thumbnail

Bessel Polynomial Linearization Coefficients

by

 

Linearization coefficients for Bessel polynomials, and certain rational polynomial spectra.

bpqcc(n,a)
function [c,k] = bpqcc(n,a)
%BPQCC   Bessel polynomial connection coefficients (n,0).
%
%   Special case of expansion in Bessel polynomials ([1], eq. 7)
% 
%       q_n(a u) = SUM_k c_k^{(n)}(a) q_k(u)
%
%   k = 0 to n where the q_k(u) are Bessel polynomials of order k.
%   Connection coefficients c_k^{(n)}(a) are equal to the special case of
%   linearization coefficients beta_k^{(n,0)}(a).
%
%   [c,n] = bpqcc(n,a)  Column vector c of n+1 Bessel polynomial connection
%   coefficients, and column vector k of indexes 0 to n, describing an
%   expansion of q_n(a u) into a weighted sum of q_k(u).
%
%   Example:
%    n = 3; a = 0.3; [c,k] = bpqcc(n,a);
%    bp = bpqlc(n,a); [b,kk] = bp.bkvec(n,0);
%    format short g; disp([k c kk b]);
%         0          0.7            0          0.7
%         1       0.1974            1       0.1974
%         2       0.0756            2       0.0756
%         3        0.027            3        0.027
%
% Algorithm: Equation from Theorem 2.1 of [1].
%
% References:
% [1] Christian Berg and Christophe Vignat, Linearization coefficients of
% Bessel polynomials, arXiv:math/0506458v1 [math.PR] 

% 2013-06-01 Robert Dickson

ik = 0:n-1;
fac = a.^ik.*exp(gammaln(n+1)-gammaln(n-ik+1)-gammaln(ik+1)- ...
    gammaln(2*n+1)+gammaln(2*(n-ik)+1)+gammaln(2*ik+1));

c = [zeros(n,1); a^n];

for kk = ik
    r = 1:min(n-kk,kk+1);
    v = (1-a).^r.*exp(gammaln(n+2)-gammaln(n-kk+r+1)-gammaln(kk+2-r)+ ... 
    	gammaln(n-kk)-gammaln(n-kk-r+1)-gammaln(r)); 
    c(kk+1) = fac(kk+1)*sum(v);
end

k = (0:n)';

end

Contact us