[Octave-bug-tracker] [bug #47738] expint is inaccurate for imaginary inputs

[Rik]

Follow-up Comment #3, bug #47738 (project octave):

This is clearly about the series expansion in expint.m being inaccurate.  The
relevant code from expint.m is


    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
      ## 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;


Your value of 10i is exactly the worst input for this function according to
the FIXME note.  You can get a sense of this by using a value that is close,
but goes through a different code path.  See below.


exp = 0.0454564330044554 +   0.0875512674239774i;
format long
x = 10i - eps(10)*i
x =  0.000000000000000 + 9.999999999999998i
obs = expint (x)
obs =  0.0454564330044520 + 0.0875512674239927i
err = exp - obs
err =  3.38618022510673e-15 - 1.52655665885959e-14i
abs (err)
ans =    1.56366153558805e-14


Clearly that works.

If you look at the code for the expansion the constants a1,a2,b1,b2 are only
accurate to 6 significant digits.  It is therefore no surprise that the answer
is also only acurate to approximately 1e-6.  I don't know where this expansion
came from, but finding more accurate versions of the coefficients might be
enough to make this work.  Maybe some internet searching will yield the answer
quickly.


    _______________________________________________________

Reply to this item at:

  <http://savannah.gnu.org/bugs/?47738>

_______________________________________________
  Message sent via/by Savannah
  http://savannah.gnu.org/