home *** CD-ROM | disk | FTP | other *** search
/ Liren Large Software Subsidy 10 / 10.iso / l / l450 / 1.ddi / MATPAK.1 / arc / BESSELA.M < prev    next >
Encoding:
Text File  |  1989-11-22  |  3.7 KB  |  107 lines
  1. function [J,digits] = bessela(nu,x)
  2. %BESSELA  Fractional order Bessel functions and accuracy estimate.
  3. %       Warning: May be inaccurate for large arguments.  See below.
  4. %       J = bessela(nu,X) = Bessel function, J-sub-nu (X).
  5. %       Bessel functions are solutions to Bessel's differential
  6. %       equation of order nu:
  7. %                2                    2    2
  8. %               x * y'' +  x * y' + (x - nu ) * y = 0
  9. %
  10. %       The function is evaluated for every point in the array X.
  11. %       The order, nu, must be a nonnegative scalar.
  12. %       See also BESSEL, BESSELN and BESSELH.
  13. %
  14. %       For some values of the arguments, this computation is
  15. %       severely contaminated by roundoff error.  Consequently
  16. %       [J,digits] = bessela(nu,x) returns an estimate of the
  17. %       number of correct significant digits in the computed
  18. %       result.  digits is the log10 of the estimated relative error,
  19. %       so digits = 14 or 15 corresponds to nearly full accuracy
  20. %       in IEEE or VAX arithmetic, while digits = 1 or 2 indicates
  21. %       nearly useless results.  Any negative value of digits is
  22. %       replaced by zero, the corresponding J set to NaN and a
  23. %       division by zero warning message is generated.
  24. %       If either nu or x is less than 50, digits will be at least 8.
  25. %       In the (nu,x) plane, the region of least accuracy is near
  26. %       the line nu = x, so small values of nu and large values of x,
  27. %       or vice versa, give the most accurate results.
  28. %       Here some representative values of digits:
  29. %
  30. %                                       nu
  31. %                                    25     75
  32. %                                 --------------
  33. %                        x   25   |  11.8    14.4
  34. %                            75   |  14.5     1.9
  35.  
  36. %       C.B. Moler, 7-13-87
  37. %       Copyright (c) 1987 by the MathWorks, Inc.
  38. %
  39. %       Reference: Abramowitz and Stegun, Handbook of Mathematical
  40. %         Functions, National Bureau of Standards, Applied Math.
  41. %         Series #55, formulas 9.1.10 and 9.2.5.
  42.  
  43.    % Temporarily replace x=0 with x=1
  44.    xz = x==0;
  45.    x = x + xz;
  46.  
  47.    % Asymptotic series
  48.    KMIN = 12; % Try the series for at least KMIN/2 terms.
  49.    z = (-1)./(8*x).^2;
  50.    t = 0*x;
  51.    q = t;
  52.    t = t + 1;
  53.    p = t;
  54.    r = abs(t);
  55.    mu = 2*nu;
  56.    k = 0;
  57.    l = -1;
  58.    oka = real(x) > 0.9*nu;
  59.    while any(any(abs(t) > eps*max(abs(p),abs(q)) & oka))
  60.       s = t;
  61.       k = k+1;
  62.       l = l+2;
  63.       t = (mu-l)*(mu+l)/k*t;
  64.       q = q + t;
  65.       k = k+1;
  66.       l = l+2;
  67.       t = (mu-l)*(mu+l)/k*t.*z;
  68.       p = p+t;
  69.       r = max(r,abs(t));
  70.       if k >= 2*KMIN, oka = oka .* (abs(t) <= 5*abs(s)); end
  71.    end
  72.    q = q./(8*x);
  73.    y = x - (2*nu+1)*pi/4;
  74.    Ja = sqrt((2)./(pi*x)) .* (p.*cos(y) - q.*sin(y));
  75.    da = max(oka.*log10(max(abs(p),abs(q))./(eps*abs(r))),0);
  76.    oka = oka.*(da > 0);
  77.  
  78.    % If necessary, use power series expansion.
  79.    % Requires the gamma function from another M file
  80.    okp = da < 12;
  81.    if any(any(okp))
  82.       t = okp .* (x/2).^nu / gamma(nu+1);
  83.       s = t;
  84.       r = abs(t);
  85.       z = -(x/2).^2;
  86.       k = 0;
  87.       while any(any(abs(t) > eps*abs(s) & okp))
  88.          k = k + 1;
  89.          t = z.*t/(k*(nu+k));
  90.          s = s + t;
  91.          r = max(r,abs(t));
  92.          okp = okp .* (eps*abs(r) <= abs(s));
  93.       end
  94.       Jp = s;
  95.       dp = max(okp.*log10((1-okp+abs(s))./(1-okp+eps*abs(r))),0);
  96.       okp = okp.*(dp > 0).*(dp>da);
  97.       oka = oka.*(da>=dp);
  98.       digits = oka.*da + okp.*dp;
  99.       J = oka.*Ja + okp.*Jp + (0)./digits;
  100.    else
  101.       J = Ja;
  102.       digits = da;
  103.    end
  104.  
  105.    % Restore results for x = 0; J0(0) = 1, Jnu(0) = 0 for nu > 0.
  106.    J = (1-xz).*J + xz*(nu==0);
  107.