home *** CD-ROM | disk | FTP | other *** search
/ Liren Large Software Subsidy 10 / 10.iso / l / l460 / 2.ddi / SPECFUN.DI$ / BESSELA.M < prev    next >
Encoding:
Text File  |  1993-03-07  |  3.5 KB  |  114 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. %
  13. %    See also BESSEL, BESSELJ, BESSELI.
  14. %
  15. %    For some values of the arguments, this computation is
  16. %    severely contaminated by roundoff error.  Consequently
  17. %    [J,digits] = bessela(nu,x) returns an estimate of the
  18. %    number of correct significant digits in the computed
  19. %    result.  digits is the log10 of the estimated relative error,
  20. %    so digits = 14 or 15 corresponds to nearly full accuracy
  21. %    in IEEE or VAX arithmetic, while digits = 1 or 2 indicates
  22. %    nearly useless results.  Any negative value of digits is
  23. %    replaced by zero, the corresponding J set to NaN and a
  24. %    division by zero warning message is generated.
  25. %    If either nu or x is less than 50, digits will be at least 8.
  26. %    In the (nu,x) plane, the region of least accuracy is near
  27. %    the line nu = x, so small values of nu and large values of x,
  28. %    or vice versa, give the most accurate results.
  29. %    Here some representative values of digits:
  30. %
  31. %                                      nu
  32. %                                  25      75
  33. %                              --------------
  34. %                     x   25   |  11.8    14.4
  35. %                         75   |  14.5     1.9
  36.  
  37. %    Reference: Abramowitz and Stegun, Handbook of Mathematical
  38. %      Functions, National Bureau of Standards, Applied Math.
  39. %      Series no. 55, formulas 9.1.10 and 9.2.5.
  40. %    C.B. Moler, 7-13-87, 4-20-90.
  41. %    Copyright (c) 1984-93 by The MathWorks, Inc.
  42.  
  43. % Temporarily replace x=0 with x=1
  44. xz = abs(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 = abs(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)./sqrt(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.  
  108. if min(digits <= 12) & (nargout < 2)
  109.    disp(' ');
  110.    disp('Warning: Results from BESSEL may be inaccurate.')
  111.    disp('Use BESSELA with two output arguments.')
  112.    disp('[result,digits] = bessela(nu,x).')
  113. end
  114.