# HG changeset patch # User Rik # Date 1378249266 25200 # Node ID 9aca7020c89f3ce854d37e7872fa407061a14281 # Parent 6e1a3b8fc3124fe82117b64833cbdafaf84839e6 expint.m: Overhaul function. Better accuracy and faster performance. Improved documentation and added %!tests. * scripts/specfun/expint.m: Improved accuracy of calculations by adjusting quad() tolerance options. Speeded up performance by breaking out of series summation as soon as possible. Used persistent variables to avoid re-calculating expensive intermediate values. Added %!tests verified against arbitrary precision Mathematica calculations. diff -r 6e1a3b8fc312 -r 9aca7020c89f scripts/specfun/expint.m --- a/scripts/specfun/expint.m Tue Sep 03 23:35:55 2013 +0100 +++ b/scripts/specfun/expint.m Tue Sep 03 16:01:06 2013 -0700 @@ -20,24 +20,54 @@ ## -*- texinfo -*- ## @deftypefn {Function File} {} expint (@var{x}) -## Compute the exponential integral, +## Compute the exponential integral: ## @tex ## $$ -## E_1 (x) = \int_x^\infty {e^{-t} \over t} dt. +## {\rm E_1} (x) = \int_x^\infty {e^{-t} \over t} dt ## $$ ## @end tex ## @ifnottex ## ## @example ## @group -## infinity -## / -## expint (x) = | exp (-t)/t dt -## / -## x +## infinity +## / +## E_1 (x) = | exp (-t)/t dt +## / +## x ## @end group ## @end example +## @end ifnottex ## +## Note: For compatibility, this functions uses the @sc{matlab} definition +## of the exponential integral. Most other sources refer to this particular +## value as @math{E_1 (x)}, and the exponential integral is +## @tex +## $$ +## {\rm Ei} (x) = - \int_{-x}^\infty {e^{-t} \over t} dt. +## $$ +## @end tex +## @ifnottex +## +## @example +## @group +## infinity +## / +## Ei (x) = - | exp (-t)/t dt +## / +## -x +## @end group +## @end example +## @end ifnottex +## +## The two definititions are related, for positive real values of @var{x}, by +## @tex +## $ +## E_1 (-x) = -{\rm Ei} (x) - i\pi. +## $ +## @end tex +## @ifnottex +## @w{@code{E_1 (-x) = -Ei (x) - i*pi}}. ## @end ifnottex ## @end deftypefn @@ -47,98 +77,85 @@ print_usage (); endif - y = expint_E1 (x); - -endfunction + y = x; # Copy over all values, including NaNs -## -*- texinfo -*- -## @deftypefn {Function File} {@var{y} =} expint_E1 (@var{x}) -## Compute the exponential integral, -## @verbatim -## infinity -## / -## expint(x) = | exp(-t)/t dt -## / -## x -## @end verbatim -## @end deftypefn - -function y = expint_E1 (x) + if (isreal (x)) + idx = (x >= 0); + y(idx) = -expint_Ei (-x(idx)); - if (nargin != 1) - print_usage (); - endif - - y = x; - - idx = (imag (x) > 0 & imag (x) != 0); - y(idx) = -expint_Ei (-y(idx)) - i.*pi; + idx = (x < 0); + y(idx) = -expint_Ei (-x(idx)) - i*pi; + else + idx = (imag (x) > 0); + y(idx) = -expint_Ei (-x(idx)) - i*pi; - idx = (imag (x) < 0 & imag (x) != 0); - y(idx) = -expint_Ei (-y(idx)) + i.*pi; + idx = (imag (x) < 0); + y(idx) = -expint_Ei (-x(idx)) + i*pi; - idx = (real (x) >= 0 & imag (x) == 0); - y(idx) = -expint_Ei (-y(idx)); + isreal_idx = (imag (x) == 0); + idx = (isreal_idx & real (x) >= 0); + y(idx) = -expint_Ei (-x(idx)); - idx = (real (x) < 0 & imag (x) == 0); - y(idx) = -expint_Ei (-y(idx)) - i.*pi; + idx = (isreal_idx & real (x) < 0); + y(idx) = -expint_Ei (-x(idx)) - i*pi; + endif endfunction ## -*- texinfo -*- ## @deftypefn {Function File} {@var{y} =} expint_Ei (@var{x}) -## Compute the exponential integral, +## Compute the exponential integral: ## @verbatim -## infinity -## / -## expint_Ei(x) = - | exp(t)/t dt -## / -## -x +## infinity +## / +## expint_Ei (x) = - | exp(-t)/t dt +## / +## -x ## @end verbatim ## @end deftypefn function y = expint_Ei (x) - if (nargin != 1) - print_usage (); - endif - y = zeros (size (x)); F = @(x) exp (-x)./x; - s = prod (size (x)); - for t = 1:s; - if (x(t) < 0 && imag (x(t)) == 0) - y(t) = -quad (F, -x(t), Inf); + for t = 1:numel (x) + xt = x(t); + if (xt < 0 && imag (xt) == 0) + ## Direct integration for most real inputs + y(t) = -quad (F, -xt, Inf, [0, 1e-10]); + elseif (xt > 2 && imag (xt) == 0) + persistent Ei_2 = 4.954234356001890; + y(t) = Ei_2 - quad (F, -xt, -2); + elseif (abs (xt) < 10) + ## Series Expansion for real (range [0,2]) or complex inputs (r < 10) + k = 1; + do + term = xt^k / (k*factorial (k)); + y(t) += term; + until (abs (term) < eps (abs (y(t))) / 2 || k++ >= 100) + y(t) = 0.57721566490153286 + log (xt) + y(t); else - if (abs (x(t)) > 2 && imag (x(t)) == 0) - y(t) = expint_Ei (2) - quad (F, -x(t), -2); + ## FIXME: This expansion is accurate to only 1e-13 at the beginning + ## near 10+i, although it becomes more accurate as the magnitude + ## of xt grows. + if (imag (xt) <= 0) + persistent a1 = 4.03640; + persistent a2 = 1.15198; + persistent b1 = 5.03637; + persistent b2 = 4.19160; + y(t) = -(xt^2 - a1*xt + a2) ... + / ((xt^2 - b1*xt + b2) * (-xt) * exp (-xt)) ... + - i*pi; else - if (abs (x(t)) >= 10) - if (imag (x(t)) <= 0) - a1 = 4.03640; - a2 = 1.15198; - b1 = 5.03637; - b2 = 4.19160; - y(t) = -(x(t).^2 - a1.*x(t) + a2) ... - ./ ((x(t).^2 - b1.*x(t) + b2) .* (-x(t)) .* exp (-x(t))) ... - - i.*pi; - else - y(t) = conj (expint_Ei (conj (x(t)))); - endif; - ## Serie Expansion - else - for k = 1:100; - y(t) = y(t) + x(t).^k ./ (k.*factorial (k)); - endfor - y(t) = 0.577215664901532860606512090082402431 + log (x(t)) + y(t); - endif - endif + y(t) = conj (expint_Ei (conj (xt))); + endif; endif endfor endfunction -%% Test against A&S Table 5.1 + +## Test against A&S Table 5.1 %!test %! x = [5:5:50]'/100; %! gamma = 0.5772156649; @@ -198,7 +215,52 @@ %! y = expint (x); %! assert (y, y_exp, 1e-9); -%% Test input validation +## Series expansion (-2 < x < 0) +## Expected values from Mathematica +%!test +%! x = [-0.1; -0.5; -1; -1.5; -2]; +%! y_exp = [ 1.6228128139692767 - i*pi; +%! -0.45421990486317358 - i*pi; +%! -1.8951178163559368 - i*pi; +%! -3.3012854491297978 - i*pi; +%! -4.9542343560018902 - i*pi]; +%! y = expint (x); +%! assert (y, y_exp, eps (real (y_exp))); + +## (x < -2, x real) +%!test +%! x = [-2.5; -3; -10;-15; -25]; +%! y_exp = [-7.0737658945786007 - i*pi; +%! -9.9338325706254165 - i*pi; +%! -2492.2289762418777 - i*pi; +%! -234955.85249076830 - i*pi; +%! -3.0059509065255486e9 - i*pi]; +%! y = expint (x); +%! assert (y, y_exp, 8*eps (real (y_exp))); + +## Complex values +%!test +%! x = [i; -1-i; 10-i; 10+i]; +%! y_exp = [-0.33740392290096813 - i*0.62471325642771360; +%! -1.7646259855638540 + i*0.75382280207927082; +%! 1.90746381979783120e-6 + i*3.67354374003294739e-6; +%! 1.90746381979783120e-6 - i*3.67354374003294739e-6]; +%! y = expint (x); +%! assert (y, y_exp, 1e-12); + +## Exceptional values (-Inf, Inf, NaN, 0, 0.37250741078) +%!test +%! x = [-Inf; Inf; NaN; 0; -0.3725074107813668]; +%! y_exp = [-Inf - i*pi; +%! -Inf; # should be 0; +%! NaN; +%! Inf; +%! 0 - i*pi]; +%! y = expint (x); +%! assert (y, y_exp, 5*eps); + +## Test input validation %!error expint () %!error expint (1,2) +