help-octave
[Top][All Lists]
Advanced

[Date Prev][Date Next][Thread Prev][Thread Next][Date Index][Thread Index]

## Bessel functions

 From: Eyal Doron Subject: Bessel functions Date: Mon, 2 Oct 1995 16:52:37 +0100 (MET)

Hi,
Ok, here they are. The following .m file calculates the Bessel and
Newmann functions for real orders and complex arguments. It returns
all orders from {\nu} to \nu in increments of one, where \nu is the
maximal requested order. Also works for negative \nu, and can return
orders from -\nu to \nu .

This .m file uses various variations on the subject of forward-backward
recurrences + normalization. It seems to work to about 12 digits at least,
usually to 14-15 digits, and is fully vectorized.

I have taken care of various sections in the complex plane, but I'm not
sure that all of them are OK, so use at your own risk. Also, it should
be possible to extend it to complex orders, but I'm too lazy.

Enjoy!

Eyal Doron

-----------------------jybess.m------------------------
function [Jn,Yn] = jybess(n,x,ToDo);
% FUNCTION [Jn,Yn] = jybess(n,x,ToDo):
%
% returns the Bessel and Newmann functions for complex arguments and real
% order. Returns a size [length(x), |n|+1] or [length(x), 2|n|+1] matrix.
%
% ToDo - Optional command string. 'Y' returns only Y_n(x) instead of
%        J_n(x), 'L' returns only the highest computed order.
% n    - Maximal order, real scalar.
% x    - Real or complex argument vector.
% ToDo - Optional string of command characters:
%        'S' - returns the orders -n to n rather than 0 to n. (also ToDo=1).
%        'Y' - if nargout<2, return Y_n(x) instead of J_n(x).
%        'L' - return only the highest computed order.
% -----------------------------------------
% Jn   - Bessel functions of the first kind.
% Yn   - Bessel functions of the second kind (Newmann functions), optional.

% Algorithm:
% ----------
% 1) Calculate J_n by backwards recursion from a calculated starting point.
%    Minimal order is \nu=|n|-floor(|n|).
% 2) Normalization: Normalize J_\nu(x) by:
%    a) for integer n and |\Im(x)|<log(10), use
%                         1 + 2 sum_{k=1}^\infty J_{2k}(x) = 1
%    b) for integer n and |\Im(x)|>log(10), use
%                         1 + 2 sum_{k=1}^\infty (-)^k J_{2k}(x) = cos(x)
%    c) for half-integer n, use explicit formula for the spherical Bessels.
%    d) for real n and |\Im(x)|<5, use
%       sum_{k=0}^\infty (\nu+2k)\Gamma(\nu+k)/k! J_{\nu+2k}(x) = (x/2)^\nu
%    e) for real n and |\Im(x)|>5, use
%       e^{i\phi\nu} \sum_{k=0}^\infty [r/2*(1-e^{2i\phi})]^k/k! J_{\nu+k}(r)
%                   = J_\nu(r e^i\phi)
%       This requires a recursive call to jybess to obtain the J_{\nu+k}(r).
% 3) If Y_n, or J_n for negative non-integer n, are required, calculate
%    the Y_n(x) as follows:
% 4) Evaluate Y_\nu:
%    a) for nu=0, use
%      2/\pi(\log(x/2)+\gamma)J_0(x) - 4/\pi\sum_{k=1}^\infty (-)^k J_{2k}(x)/k
%                  = Y_0(x)
%    b) for nu=0.5, use explicit formula for the spherical Bessels.
%    c) otherwise, do the following:
%        i) Calculate J_{1-\nu}(x) and J_{2-\nu}(x) using a recursive call.
%       ii) Calculate J_{-\nu}(x) using one backwards recursion step.
%      iii) Use [J_\nu(x)\cos(\nu\pi)-J_{-\nu}(x)]/\sin(\nu\pi) = Y_\nu(x)
% 5) Obtain the rest of the Y_{\nu+k}(x) using the Wronskian relation
%    J_{\nu+1}(x)Y_\nu(x) - J_\nu(x)Y_{\nu+1}(x) = 2/(\pi x)
%    This turns out to be much more stable for complex arguments than the
%    more conventional forward recurrence technique.
% 6) If J for negative orders is required:
%    a) for nu=0, use J_n(x) = (-)^n J_{-n}(x)
%    b) otherwise, use 4.c.iii
% 7) If Y for negative orders is required:
%    a) for nu=0, use Y_n(x) = (-)^n Y_{-n}(x)
%    b) otherwise, use 4.c.iii

if nargin==0
help jybess
return
end
if nargin<2
error('Not enough input parameters!');
end
if max(max(size(n)))>1 | any(imag(n)~=0)
error('Max order must be a real scalar!')
end
sym=0; Return_Y=(nargout==2); Return_Last=0;
if nargin>2
sym=any(ToDo==1) | any(ToDo=='s') | any(ToDo=='S');
Return_Y=Return_Y | any(ToDo=='y' | ToDo=='Y');
Return_Last=any(ToDo=='l' | ToDo=='L');
end
orig_n=n;
n_is_negative=(n<0); n=abs(n);
nu=n-floor(n); n=floor(n);
if nargin<3, sym=0; end
sym=(sym~=0);
x=x(:).'; xlen=length(x);
calc_Y=(Return_Y | sym | n_is_negative);
Acc_Jnorm=(nu~=0.5);
Acc_Ynorm=calc_Y & nu==0;

z = x==0; x = x + z;            % Temporarily replace x=0 with x=1

tiny = 16^(-250);

%                             Starting index for backwards recurrence.
c = [ 0.9507    1.4208   14.1850
0.9507    1.4208   14.1850
0.9507    1.4208   14.1850
0.7629    1.4222   13.9554
0.7369    1.4289   13.1756
0.7674    1.4311   12.4523
0.8216    1.4354   11.2121
0.8624    1.4397    9.9718
0.8798    1.4410    8.9217
0.9129    1.4360    7.8214
0.9438    1.5387    6.5014
0.9609    1.5216    5.7256
0.9693    1.5377    5.3565
0.9823    1.5220    4.5659
0.9934    1.5049    3.7902
0.9985    1.4831    3.2100
1.0006    1.4474    3.0239
0.9989    1.4137    2.8604
0.9959    1.3777    2.7760
1.0005    1.3500    2.3099]';
j = 1+min(n,19);
m = c(1,j).*max(3,j) + c(2,j).*(max(1,abs(x))-1) + ...
c(3,j)./(1-log(min(1,abs(x))));
m=max([m; n+10+0*m]);
m = 4*ceil(m/4);

%                                       Prevent underflow
logtiny=log(tiny); logx=log(abs(x)/2)+1;
II=1:xlen;
while ~isempty(II)
est=-log(2*pi*m(II))/2+m(II).*(logx(II)-log(m(II)))<logtiny;
II=II(find(est)); m(II)=m(II)-4;
end;
mm = max(m(:));

% Normalization summation coefficients

if nu==0
Nrm=[1; 2*ones(ceil(mm/2),1)];
UseCos=find(abs(imag(x))>log(10));
if ~isempty(UseCos)                  % Use norm. to cos(z) instead of 1.
Nrm=Nrm(:,ones(xlen,1));
[r,c]=size(Nrm);
Nrm(2:2:r,UseCos)=-Nrm(2:2:r,UseCos);
end
elseif nu~=0.5               % nu=0.5 doesn't use normalizations
k=(1:ceil(mm/2)).';
Nrm=cumprod([nu*gamma(nu); (nu+2*k).*(nu+k-1)./k./(nu+2*k-2)]);
Nrm=Nrm(:,ones(xlen,1));
end

%                           Backwards recursion for the Jn
Jn=zeros(n+1,xlen);
k = mm;
bkp1 = 0*x;
bk = tiny*(m==k);
if Acc_Jnorm, t = Nrm(k/2+1,:).*bk; end
if Acc_Ynorm, y = 2*bk/k; end
sig = 1;
for k = mm-1:-1:0
bkp2 = bkp1;
bkp1 = bk;
bk = 2*(k+1+nu)*bkp1./x - bkp2 + tiny*(m==k);
if k<=n
Jn(k+1,:)=bk;
end
if (floor(k/2)*2-k==0)
if Acc_Jnorm, t = t + Nrm(k/2+1,:).*bk; end
if Acc_Ynorm & k>0
sig = -sig;
y = y + sig*2*bk/k;
end
end
end
clear Nrm bkp2 sig

if nu==0.5               % Use explicit formulae for the spherical Bessels
si=sin(x);
II=(abs(si)>0.1);
F1=find(II); F2=find(1-II);
nrm=ones(1,xlen);
SpFact=sqrt(2*x/pi);
if ~isempty(F1)
nrm(F1)=SpFact(F1).*(si(F1)./x(F1))./Jn(1,F1);
end
if ~isempty(F2)
nrm(F2)=SpFact(F2).*(si(F2)./x(F2).^2-cos(x(F2))./x(F2))./bkp1(F2);
end
elseif nu==0
nrm=ones(size(x));
if ~isempty(UseCos)
nrm(UseCos)=cos(x(UseCos));
end
nrm=nrm./t;
else
II=(abs(imag(x))<5);
F1=find(II); F2=find(1-II);
nrm=ones(1,xlen);
if ~isempty(F1)
nrm(F1)=(x(F1)/2).^nu ./t(F1);
end
if ~isempty(F2)               % Use mult. Theorem to calc J_0
r=abs(x(F2)); ph=angle(x(F2)); rn=max(ceil(r));
MaxOrder=ceil(18+1.25*rn);
Jtmp=jybess(nu+MaxOrder,r);
k=(0:MaxOrder).';
Kg=cumsum([0;log(1:MaxOrder).']);          % log(k!)
v=i*ph-i*pi/2+log(imag(x(F2)));

nrm(F2)=exp(i*nu*ph).*sum(exp(k*v-Kg(:,ones(length(F2),1))).*Jtmp.')./Jn(1,F2);
end
end
Jn=Jn.*nrm(ones(n+1,1),:);           % Normalizing condition.

%                Restore results for x = 0; j0(0) = 1, jn(0) = 0 for n > 0.
if any(z)
Ii=find(z); LI=max(size(Ii));
Jn(:,Ii)=[ones(1,LI);zeros(n,LI)];
end
Jn=Jn.';

if calc_Y                                    % Also Yn
%                         First, get Y_0(x)
if n==0
J1 = bkp1.*nrm;
else
J1=Jn(:,2).';
end
if nu==0                          % Use dependence on sums of J
gamma = 0.57721566490153286;
lx=log(x/2);
II=find(imag(x)==0 & real(x)<0);
if ~isempty(II)
lx(II)=log(abs(x(II))/2);
end
y = 2/pi*(lx + gamma).*bk - 4/pi*y;
Y0 = y.*nrm;
elseif nu==0.5                    % explicit formula
Y0=-cos(x)./x.*SpFact;
else                              % Use linear relations between J and Y
Jtmp=jybess(2-nu,x);
Jminus=2*(1-nu)./x.'.*Jtmp(:,1)-Jtmp(:,2);
Y0=(Jn(:,1).'*cos(pi*nu)-Jminus.')/sin(pi*nu);
end

Yn=[Y0.', zeros(xlen,n)];
for l=1:n                      % Iterate using Wronskian relation
Yn(:,l+1)=(Jn(:,l+1).*Yn(:,l)-2/pi./x)./Jn(:,l);
end

if any(z), Yn(find(z),:)=-inf*ones(size(Yn(find(z),:))); end
end

if nu==0
if n_is_negative                                 % negative order
Jn(:,2:2:n+1)=-Jn(:,2:2:n+1);
if Return_Y
Yn(:,2:2:n+1)=-Yn(:,2:2:n+1);
end
end
if sym & n>0                                     % symmetric output
Rev=n:-1:1; Rev=2*(floor(Rev/2)*2==Rev)-1;
Jn=[Jn(:,n+1:-1:2).*Rev(ones(xlen,1),:) Jn];
if Return_Y
Yn=[Yn(:,n+1:-1:2).*Rev(ones(xlen,1),:) Yn];
end
end
else
if n_is_negative | sym
ord=nu+(0:n); co=cos(pi*ord); so=sin(pi*ord);
Jminus=Jn.*co(ones(xlen,1),:)-Yn.*so(ones(xlen,1),:);
if Return_Y
ord=nu+(0:n);co=cos(pi*ord); so=sin(-pi*ord);
Yminus=(Jminus.*co(ones(xlen,1),:)-Jn)./so(ones(xlen,1),:);
end
end
if sym
Jn=[fliplr(Jminus) Jn];
if n_is_negative, Jn=fliplr(Jn); end
if Return_Y
Yn=[fliplr(Yminus) Yn];
if n_is_negative, Yn=fliplr(Yn); end
end
elseif n_is_negative
Jn=Jminus;
if Return_Y, Yn=Yminus; end
end
end
if Return_Last
[r,c]=size(Jn);
Jn=Jn(:,c);
if calc_Y, Yn=Yn(:,c); end
end
if nargout<2 & Return_Y, Jn=Yn; end



reply via email to

 [Prev in Thread] Current Thread [Next in Thread]
• Bessel functions, Eyal Doron <=