bugGNU Octave - Bugs: bug #58636, Integral fails to call quadgk for...

 
 

bug #58636: Integral fails to call quadgk for complex integration

Submitter:  Nicholas Jankowski <nrjank>
Submitted:  Sat 20 Jun 2020 08:19:36 PM UTC
   
 
Category:  Octave Function Severity:  3 - Normal
Priority:  5 - Normal Item Group:  Unexpected Error or Warning
Status:  Fixed Assigned to:  None
Originator Name:  Nicholas Jankowski 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 13 Sep 2020 06:30:54 PM UTC, comment #6: 

I reviewed the patch, made a few changes for Octave coding style, and checked it in under your name on the development branch here: http://hg.savannah.gnu.org/hgweb/octave/rev/4f1e8a79bd44.

Marking as fixed and closing report.

Rik <rik5>
Group administrator
Sun 21 Jun 2020 03:50:41 AM UTC, comment #5: 

ok. this version of integral has the limit checks, the feval(a,b) check, and a try/catch if those checks pass but the function still returns a complex answer. It currently works for the existing quadcc.  I added some complex input tests and it passes.

pending patch will tweak quadcc docstring to note it doesn't handle complex inputs or outputs, and once I can test that error messages split in quadcc, I can do that too and update integral to catch that particular error message.

doc fix in quad.txi and that should do it.

(file #49335)

Nicholas Jankowski <nrjank>
Group Member
Sun 21 Jun 2020 01:40:05 AM UTC, comment #4: 

looking at quadcc, there's really only one error message that needs to be caught.  there are two that involve errors for complex numbers, but one is for complex values in the singularities list, and integral doesn't include that detail so never passes it to quadcc. The other one is:


if (fvals.length () != 1 || ! fvals(0).is_real_matrix ())
  error ("quadcc: integrand F must return a single, real-valued vector");


at three places in the code.  I don't know if I want this to catch that and run quadgk for both. looking at the code I'm not sure what sort of input would trigger that that first length error. I tried a couple inputs but couldn't get it to pop up. Maybe if I figure out how to trigger that I could see if it needs to be separated out.

Separating them is easy enough into:


if (fvals.length () != 1)
  error ("quadcc: integrand F must return a single vector");

if (! fvals(0).is_real_matrix ())
  error ("quadcc: integrand F must return a real-valued vector");


So if it doesn't cause issues integral could try/catch based on the first one.  I also noticed there are no BISTs for that error, so I'd want to add one for each after the split, if I can find a fn for the latter.

will play with the try/catch in integral now.


Nicholas Jankowski <nrjank>
Group Member
Sat 20 Jun 2020 09:35:32 PM UTC, comment #3: 

need to see if there's a way to catch complex functions. something link

integral (exp (i*x), 0, 1) evaluates to real values at the real endpoints, so it'll sneak past the check. Also, not sure about whether there are other not-so-obvious complex functions. even something like


>> fx = @(x) complex(x)
fx =

@(x) complex (x)

>> quadgk(fx,0,1)
ans =  0.50000
>> quadcc(fx,0,1)
error: quadcc: integrand F must return a single, real-valued vector


needs to be caught.

Nicholas Jankowski <nrjank>
Group Member
Sat 20 Jun 2020 09:15:51 PM UTC, comment #2: 

and this version checks the value of the function at the limits to catch things like
integral (sqrt(x), -1,1)



(file #49333)

Nicholas Jankowski <nrjank>
Group Member
Sat 20 Jun 2020 08:54:28 PM UTC, comment #1: 

attached version of integral now works for complex limits using an iscomplex? check at the beginning. if there's an easy way to test the integrand, should be easy to add to that check and the whole function will work.  need to add some complex tests as well then this would be fixed.

(file #49332)

Nicholas Jankowski <nrjank>
Group Member
Sat 20 Jun 2020 08:19:36 PM UTC, original submission:  

As reported by Brett on the mailing list [1]:

Integral is a wrapper for quadcc, quadgk, and quadv, and defaults to quadcc where possible. quadgk is the only quadrature method built for complex integration and limits, including complex path integration, and should be called.  There is currently no check for complex functions or limits, so if waypoints are not specified the script chooses quadcc which issues an error, while both quadgk and quadv would have successfully computed the result.

Complex function case:

>> fx = @(x) (x+i.*3).^2;
>> cmplx_int_val = integral(fx,-1,2)
error: quadcc: integrand F must return a single, real-valued vector
error: called from
    integral at line 114 column 7
>> cmplx_int_val = quadcc(fx,-1,2)
error: quadcc: integrand F must return a single, real-valued vector
>> cmplx_int_val = quadgk(fx,-1,2)
cmplx_int_val = -24.0000 +  9.0000i
>> cmplx_int_val = quadv(fx,-1,2)
cmplx_int_val = -24 +  9i


Complex limits case:

>> fx = @(x) (x).^2;
>> cmplx_int_val = integral(fx,-i,i)
error: quadcc: lower limit of integration (A) must be a real scalar
error: called from
    integral at line 114 column 7
>> cmplx_int_val = quadcc(fx,-i,i)
error: quadcc: lower limit of integration (A) must be a real scalar
>> cmplx_int_val = quadgk(fx,-i,i)
cmplx_int_val =  0.00000 - 0.66667i
>> cmplx_int_val = quadv(fx,-i,i)
cmplx_int_val = -0.00000 - 0.66667i


the isreal check on the limits should be easy enough.  the check on the function may be a bit more complicated, and I'm not sure the best way to check that.

it would be inefficient, but lacking something better could just if (catch specific complex error message from quadcc) --> kick it to quadgk. 

better suggestions welcome

(also, need to fix the docs for any change including fixing the preamble to the one function integration section that doesn't mention integral calls quadcc)

[1] https://octave.1599824.n4.nabble.com/Why-is-there-no-function-for-complex-integrals-tp4697602.html

Nicholas Jankowski <nrjank>
Group Member

 

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

Attach Files:
   
   
Comment:
   

Attached Files
file #49335:  integral.m added by nrjank (11KiB - text/plain)
file #49333:  integral.m added by nrjank (11KiB - text/plain - works for cmplex limits and functions where f(a) or f(b) are complex)
file #49332:  integral.m added by nrjank (10KiB - text/plain - integral function now working for complex limits, still not for cmplx integrand)

 

Depends on the following items: None found

Items that depend on this one: None found

 

Carbon-Copy List
  • -email is unavailable- added by rik5 (Posted a comment)
  • -email is unavailable- added by nrjank (Submitted the item)
  • -email is unavailable- added by nrjank (reported on mailing list)
  •  

    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 7 latest changes.

    Date Changed by Updated Field Previous Value => Replaced by
    2020-09-13 rik5 StatusNone Fixed
        Open/ClosedOpen Closed
        Release5.2.0 dev
    2020-06-21 nrjank Attached File- Added integral.m, #49335
    2020-06-20 nrjank Attached File- Added integral.m, #49333
    2020-06-20 nrjank Attached File- Added integral.m, #49332
    2020-06-20 nrjank Carbon-Copy- Added -email is unavailable-

    Back to the top

    Powered by Savane 3.13-d3ae.
    Corresponding source code