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

 
 

bug #47738: expint is inaccurate for certain imaginary inputs

Submitted by:  Colin Macdonald <cbm>
Submitted on:  Tue 19 Apr 2016 09:18:35 PM UTC  
 
Category: Octave FunctionSeverity: 3 - Normal
Priority: 5 - NormalItem Group: Inaccurate Result
Status: FixedAssigned to: None
Originator Name: Open/Closed: Closed
Release: devOperating System: Any

Add a New Comment(Rich Markup)
   

You are not logged in

Please log in, so followups can be emailed to you.

 

(Jump to the original submission 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>
Project 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>
Project 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>
Project 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>
Project 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>
Project 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>
Project 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>
Project 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 file #39052 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 file #39054 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>
Project 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>
Project Administrator
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

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.

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>
Project 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>
Project 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.

Matlab:

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 File(s):
   
   
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
  • -unavailable- added by m_ginesi (Updated the item)
  • -unavailable- added by caliari (Updated the item)
  • -unavailable- added by dastew (Posted a comment)
  • -unavailable- added by mtmiller (Posted a comment)
  • -unavailable- added by rik5 (Posted a comment)
  • -unavailable- added by cbm (Submitted the item)
  •  

    Do you think this task is very important?
    If so, you can click here to add your encouragement to it.
    This task has 0 encouragements so far.

    Only project members can vote.

     

    Please enter the title of George Orwell's famous dystopian book (it's a date):

     

     

    Follow 15 latest changes.

    Date Changed By Updated Field Previous Value => Replaced By
    Sun 22 Jan 2017 04:37:18 PM UTCrik5StatusPatch Submitted=>Fixed
      Open/ClosedOpen=>Closed
    Tue 10 Jan 2017 06:40:22 PM UTCm_ginesiAttached File-=>Added patch_expint.txt, #39428
    Tue 27 Dec 2016 07:42:53 PM UTCrik5StatusConfirmed=>Patch Submitted
    Wed 23 Nov 2016 08:24:16 AM UTCm_ginesiAttached 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
    Mon 02 May 2016 07:38:31 AM UTCcaliariAttached File-=>Added expint.ref, #37040
      Attached File-=>Added exponential_integral.m, #37041
      Attached File-=>Added testexpint.m, #37042
    Wed 20 Apr 2016 06:22:57 PM UTCmtmillerStatusNone=>Confirmed
    Wed 20 Apr 2016 04:21:14 PM UTCrik5Summaryexpint is inaccurate for imaginary inputs=>expint is inaccurate for certain imaginary inputs
    Tue 19 Apr 2016 09:19:29 PM UTCcbmAttached File-=>Added expint_err_on_imag_axis.png, #36965
      Attached File-=>Added expint_test.m, #36966

    Back to the top


    Powered by Savane 3.1-cleanup1