bugGNU Octave - Bugs: bug #47738, expint is inaccurate for certain...

 
 

bug #47738: expint is inaccurate for certain imaginary inputs

Submitter:  Colin Macdonald <cbm>
Submitted:  Tue 19 Apr 2016 09:18:35 PM UTC
   
 
Category:  Octave Function Severity:  3 - Normal
Priority:  5 - Normal Item Group:  Inaccurate Result
Status:  Fixed Assigned to:  None
Originator Name:  Open/Closed:  * Closed
Release:  * dev Operating System:  * Any
Fixed Release:  None Planned Release:  None
* Mandatory Fields

Add a New Comment Rich Markup
   

Jump to the original submission

Sun 22 Jan 2017 04:37:18 PM UTC, comment #18: 

I pushed the cset here (http://hg.savannah.gnu.org/hgweb/octave/rev/c56d30ea6cd4).  It took a while to review because there were actually a lot differences from the Octave coding guidelines.  Some examples for the next time to look : single quotes instead of double quotes, '~' instead of '!' for logical not operator, space between function name and opening parenthesis to differentiate between indexing and function calls, parentheses around the conditional in an if clause, and long lines greater than 80 characters.

Rik <rik5>
Group administrator
Fri 13 Jan 2017 08:49:30 PM UTC, comment #17: 

@Michele: Thanks for contributing.
@Marco: Thanks for reviewing.

I'll look it over this weekend and push it.

Rik <rik5>
Group administrator
Fri 13 Jan 2017 09:00:24 AM UTC, comment #16: 

The patch follows coding conventions now, as far as I know. Moreover the class of the input is preserved now. I think the patch is ready to be pushed.

Marco Caliari <caliari>
Group Member
Tue 10 Jan 2017 06:45:10 PM UTC, comment #15: 

Dear all,
in file #39428 there is my patch to expint. I think that all (or at least most of) the points to correct are OK now. Let me know if there is still something to correct.

Michele Ginesi <m_ginesi>
Tue 10 Jan 2017 06:40:22 PM UTC, comment #14: 

Dear all,


(file #39428)

Michele Ginesi <m_ginesi>
Tue 27 Dec 2016 07:42:53 PM UTC, comment #13: 

@Michele: Can you update your patch to use Octave coding conventions?  Guidelines are at http://wiki.octave.org/Contribution_guidelines.  The guidelines specifically for m-files are at http://wiki.octave.org/Octave_style_guide.

Rik <rik5>
Group administrator
Fri 16 Dec 2016 08:31:56 AM UTC, comment #12: 

@Rik: from the mathematical and algorithmic point of view, I think Michele made a good job. He should add a little bit of documentation and better follow the style of m-files in Octave. After that, maybe with the help of some active developer in chat, he could try to make a proper patch for the development branch of Octave.

Marco Caliari <caliari>
Group Member
Wed 14 Dec 2016 05:17:14 PM UTC, comment #11: 

@Marco: Can you comment on myexpint by Michele?  The accuracy and improvement in running time make it seem like a good idea to move the function onto the development branch of Octave.

Rik <rik5>
Group administrator
Wed 14 Dec 2016 02:26:04 PM UTC, comment #10: 

Dear Marco,
I've made some comparison and these are the results: the worst (relative) error for myexpint is 4.7825e-14, while for matlab is 1.0334e-13.
In the parts of the complex plane where Matlab make a better approximation, the worst error for Matlab is 2.8209e-15, for myexpint is 4.7825e-14.
Talking about the execution time, Octave's expint needs 5.69309 sec, Matlab's expint needs 0.515487 sec and myexpint needs 0.420296 sec.

Michele Ginesi <m_ginesi>
Mon 05 Dec 2016 08:25:13 AM UTC, comment #9: 

Dear Michele,

in file #39054, what do you mean by saying Matlab is better than Octave? How much better? Moreover, can you test the efficiency of your code by comparing it (tic toc) with the current code and the code in Matlab? The discussion in bug #47800 may help you in improving, if needed, your code.


Marco Caliari <caliari>
Group Member
Wed 23 Nov 2016 08:32:40 AM UTC, comment #8: 

Dear all,

I also made some experiment with this function. In file #39051 I used three strategies to evaluate expint: series expansion, asymptotic series expansion and continued fraction. The show which strategy make the less error (exact values were computed with the octave symbolic package in the complex square [-100,100]x[-100,100].
Then I wrote a function (file #39053) which divides the input and uses a different strategy depending on where the input is in the complex plane. In there is a comparison between myexpint.m and the Matlab expint.

Michele Ginesi <m_ginesi>
Mon 02 May 2016 07:38:31 AM UTC, comment #7: 

Dear all,

some time ago I made some experiments with expint. I implemented four formulas that I found on Abramowitz&Stegun and Numerical Recipes in Fortran77. The result was that the classical series expansion (Abramowitz 5.1.11) was the most accurate in the region around the segment [-40,0] in the complex plane. The continued fraction expansion (Abramowitz 5.1.22) or its modification (Numerical Recipes) was the most accurate elsewhere. In the computational algorithm proposed by Numerical Recipes, they use the continued fraction expansion and, if it fails, they use the series.

I do not remember the details now, but I attach my experiments. exponential_integral.m computes the values in the four ways, expint.ref contains reference values (Octave's binary format) and testexpint.m compare the results.

(file #37040, file #37041, file #37042)

Marco Caliari <caliari>
Group Member
Wed 27 Apr 2016 05:35:58 PM UTC, comment #6: 

SciPy uses the specfun Fortran library (unrelated to our specfun pkg) corresponding to the book "Computation of Special Functions" by Zhang and Jin.

Largest relative error I've seen in SciPy is 30*eps, which is still large but at least better.  I'll look at the licensing (SciPy is BSD-style so should be fine).  Or I can just reimplement the algorithm from their book.  But it has a 500 loop so might be better to compile it.

I need to look further into this anyway, because sinint and cosint are also poor (cosint uses same implementation as expint but sinint has its own much worse implementation which blows up for x>300: https://github.com/cbm755/octsympy/issues/416)

I haven't tested Boost (@Mike do you know what accuracy they claim?  Does it work for complex inputs?)

@Doug: it looks like the error could be as large as 2e-8, which looks poor.

Colin Macdonald <cbm>
Wed 20 Apr 2016 07:16:08 PM UTC, comment #5: 

What about using (5.1.56) from this page:
https://books.google.ca/books?id=KiPCAgAAQBAJ&pg=PA231&lpg=PA231&dq=4.03640;+1.15198;+5.03637;+4.19160;&source=bl&ots=kfpQohRmno&sig=Y_P8kcXMiLXcs7I8Ifxc8s4Hcf0&hl=en&sa=X&ved=0ahUKEwi1vaTQ9Z3MAhXptoMKHcl2BToQ6AEIJTAC#v=onepage&q=4.03640%3B%201.15198%3B%205.03637%3B%204.19160%3B&f=false

It has a better range and more digits.

Doug Stewart <dastew>
Wed 20 Apr 2016 06:22:57 PM UTC, comment #4: 

I found a reference to this expansion in Abramowitz & Stegun with exactly those terms to six significant figures (5.1.55). There are other papers out there on approximations of the exponential integral. Boost uses http://www.jstor.org/stable/2004541?seq=1#page_scan_tab_contents.

Mike Miller <mtmiller>
Group Member
Wed 20 Apr 2016 04:20:53 PM UTC, comment #3: 

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.

Rik <rik5>
Group administrator
Wed 20 Apr 2016 04:07:57 PM UTC, comment #2: 

The answer to the second question is that the '.m' file extension is not unique.  On my Kubuntu 15.10 install it thinks the .m files are Objective-C files by default.  Of course I just overrode that and made a local association.

Rik <rik5>
Group administrator
Tue 19 Apr 2016 09:20:52 PM UTC, comment #1: 

p.s.: I wonder why does savannah thinks my .m file is a 'application/vnd.wolfram.mathematica.package'

Colin Macdonald <cbm>
Tue 19 Apr 2016 09:18:35 PM UTC, original submission:  

For some imaginary inputs, the error can be as large as 10^-6.


>> A = expint(10i)
A =   0.0454579311577793 +   0.0875455463635251i
>> S = double(expint(sym(10i)))
S =   0.0454564330044554 +   0.0875512674239774i
>> err = A - S
err = 1.49815332391168e-06 - 5.72106045235166e-06i


Matlab:

>> M = expint(10i)
M =
         0.0454564330044556 +     0.087551267423977i
>> M - S
ans =
       1.80411241501588e-16 -  3.46944695195361e-16i


I'll attach a figure of the error along line in the complex plane.

Looking at the code shows it switches between two implementations at abs(x) == 10.  There is a FIXME there, but I think the inaccuracy is more serious than it suggests.

Colin Macdonald <cbm>

 

(Note: upload size limit is set to 16384 kB, after insertion of the required escape characters.)

Attach Files:
   
   
Comment:
   

Attached Files
file #39428:  patch_expint.txt added by m_ginesi (17KiB - text/plain)
file #39051:  different_expint.m added by m_ginesi (2KiB - text/x-objcsrc)
file #39052:  testdifferent.png added by m_ginesi (70KiB - image/png)
file #39053:  myexpint.m added by m_ginesi (5KiB - text/x-objcsrc)
file #39054:  octave_vs_matlab.png added by m_ginesi (55KiB - image/png)
file #37040:  expint.ref added by caliari (65KiB - application/octet-stream)
file #37041:  exponential_integral.m added by caliari (2KiB - text/x-objcsrc)
file #37042:  testexpint.m added by caliari (871B - text/x-objcsrc)
file #36965:  expint_err_on_imag_axis.png added by cbm (26KiB - image/png)
file #36966:  expint_test.m added by cbm (254B - application/vnd.wolfram.mathematica.package)

 

Depends on the following items: None found

Items that depend on this one: None found

 

Carbon-Copy List
  • -email is unavailable- added by m_ginesi (Updated the item)
  • -email is unavailable- added by caliari (Updated the item)
  • -email is unavailable- added by dastew (Posted a comment)
  • -email is unavailable- added by mtmiller (Posted a comment)
  • -email is unavailable- added by rik5 (Posted a comment)
  • -email is unavailable- added by cbm (Submitted the item)
  •  

    There are 0 votes so far. Votes easily highlight which items people would like to see resolved in priority, independently of the priority of the item set by tracker managers.

    Only group members can vote.

     

    Follow 15 latest changes.

    Date Changed by Updated Field Previous Value => Replaced by
    2017-01-22 rik5 StatusPatch Submitted Fixed
        Open/ClosedOpen Closed
    2017-01-10 m_ginesi Attached File- Added patch_expint.txt, #39428
    2016-12-27 rik5 StatusConfirmed Patch Submitted
    2016-11-23 m_ginesi Attached File- Added different_expint.m, #39051
        Attached File- Added testdifferent.png, #39052
        Attached File- Added myexpint.m, #39053
        Attached File- Added octave_vs_matlab.png, #39054
    2016-05-02 caliari Attached File- Added expint.ref, #37040
        Attached File- Added exponential_integral.m, #37041
        Attached File- Added testexpint.m, #37042
    2016-04-20 mtmiller StatusNone Confirmed
    2016-04-20 rik5 Summaryexpint is inaccurate for imaginary inputs expint is inaccurate for certain imaginary inputs
    2016-04-19 cbm Attached File- Added expint_err_on_imag_axis.png, #36965
        Attached File- Added expint_test.m, #36966

    Back to the top

    Powered by Savane 3.13-cf05.
    Corresponding source code