Why is there no function for complex integrals?

classic Classic list List threaded Threaded
6 messages Options
Reply | Threaded
Open this post in threaded view
|

Why is there no function for complex integrals?

BGreen
I needed to evaluate a complex integral, but Octave only handles real integrals, so I wrote a very simple function to split a complex function into its real and imaginary parts to be integrated separately.

function cmplx_int_val = complex_integral(cmplx_fx_handle,lowerlim,upperlim)
fx_re = @(x) real(cmplx_fx_handle(x));
fx_im = @(x) imag(cmplx_fx_handle(x));
int_val_re = integral(fx_re,lowerlim,upperlim);
int_val_im = integral(fx_im,lowerlim,upperlim);
cmplx_int_val = int_val_re + i*int_val_im;
end

There is nothing novel about this, which is why I'm surprised I didn't read in the documentation about an existing function that serves the same purpose. Is there some subtle numerical reason why writing a complex integral this way may be a bad idea and produce unreliable results, or something like that?

- Brett Green


Reply | Threaded
Open this post in threaded view
|

Re: Why is there no function for complex integrals?

nrjank
Administrator
On Fri, Jun 19, 2020 at 1:59 PM Brett Green <[hidden email]> wrote:
I needed to evaluate a complex integral, but Octave only handles real integrals


I've used octave to do complex integration quite often.  Where were you looking for info on complex integration that you didn't come across the help for quadgk or other integral functions? Perhaps there's something we can do to improve the documentation  

see:


specifically, see the functions quadgk and integral

 


Reply | Threaded
Open this post in threaded view
|

Re: Why is there no function for complex integrals?

BGreen

On Sat, Jun 20, 2020 at 3:01 PM Nicholas Jankowski <[hidden email]> wrote:
On Fri, Jun 19, 2020 at 1:59 PM Brett Green <[hidden email]> wrote:
I needed to evaluate a complex integral, but Octave only handles real integrals


I've used octave to do complex integration quite often.  Where were you looking for info on complex integration that you didn't come across the help for quadgk or other integral functions? Perhaps there's something we can do to improve the documentation  

see:


specifically, see the functions quadgk and integral

This is strange. Reading the documentation, it sounds as if Octave should handle this MWE without any problems and without a specially-constructed function.

>> 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 = complex_integral(fx,-1,2)
cmplx_int_val = -24.0000 +  9.0000i

Restating complex_integral for completeness:

function cmplx_int_val = complex_integral(cmplx_fx_handle,lowerlim,upperlim)
fx_re = @(x) real(cmplx_fx_handle(x));
fx_im = @(x) imag(cmplx_fx_handle(x));
int_val_re = integral(fx_re,lowerlim,upperlim);
int_val_im = integral(fx_im,lowerlim,upperlim);
cmplx_int_val = int_val_re + i*int_val_im;
end

I notice that in the documentation, of the quadrature algorithms, only quadgk is described as handling complex integrals. Indeed, if I use this there is no problem.

>> cmplx_int_val = quadgk(fx,-1,2)
cmplx_int_val = -24.0000 +  9.0000i

So apparently the issue is that the wrapper function integral is choosing the wrong quadrature algorithm.

Similarly, the documentation states

integral

A compatibility wrapper function that will choose between quadv and quadgk depending on the integrand and options chosen.

...but the error in my first block of code shows that integral is using quadcc, contradicting the documentation which states integral chooses between quadv or quadgk.

I have not yet filed a bug report, in case I am misunderstanding something, but I will if you agree that this is a bug. It seems to me that both the function and the documentation have small problems regarding the quadrature algorithm chosen.


Reply | Threaded
Open this post in threaded view
|

Re: Why is there no function for complex integrals?

nrjank
Administrator
On Sat, Jun 20, 2020 at 3:36 PM Brett Green <[hidden email]> wrote:

On Sat, Jun 20, 2020 at 3:01 PM Nicholas Jankowski <[hidden email]> wrote:
On Fri, Jun 19, 2020 at 1:59 PM Brett Green <[hidden email]> wrote:
I needed to evaluate a complex integral, but Octave only handles real integrals


I've used octave to do complex integration quite often.  Where were you looking for info on complex integration that you didn't come across the help for quadgk or other integral functions? Perhaps there's something we can do to improve the documentation  

see:


specifically, see the functions quadgk and integral

This is strange. Reading the documentation, it sounds as if Octave should handle this MWE without any problems and without a specially-constructed function.

>> 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 = complex_integral(fx,-1,2)
cmplx_int_val = -24.0000 +  9.0000i

Restating complex_integral for completeness:

function cmplx_int_val = complex_integral(cmplx_fx_handle,lowerlim,upperlim)
fx_re = @(x) real(cmplx_fx_handle(x));
fx_im = @(x) imag(cmplx_fx_handle(x));
int_val_re = integral(fx_re,lowerlim,upperlim);
int_val_im = integral(fx_im,lowerlim,upperlim);
cmplx_int_val = int_val_re + i*int_val_im;
end

I notice that in the documentation, of the quadrature algorithms, only quadgk is described as handling complex integrals. Indeed, if I use this there is no problem.

>> cmplx_int_val = quadgk(fx,-1,2)
cmplx_int_val = -24.0000 +  9.0000i

So apparently the issue is that the wrapper function integral is choosing the wrong quadrature algorithm.

Similarly, the documentation states

integral

A compatibility wrapper function that will choose between quadv and quadgk depending on the integrand and options chosen.

...but the error in my first block of code shows that integral is using quadcc, contradicting the documentation which states integral chooses between quadv or quadgk.

I have not yet filed a bug report, in case I am misunderstanding something, but I will if you agree that this is a bug. It seems to me that both the function and the documentation have small problems regarding the quadrature algorithm chosen.

I worked on that quickie wrapper a little while ago. my guess is you're not wrong. if it's not using quadgk, it's choosing quadv.  quadv may use quadcc under the hood. I haven't looked at quadv in a while.  Might need to add an 'iscomplex' check.


Reply | Threaded
Open this post in threaded view
|

Re: Why is there no function for complex integrals?

nrjank
Administrator
In reply to this post by BGreen
Similarly, the documentation states

integral

A compatibility wrapper function that will choose between quadv and quadgk depending on the integrand and options chosen.


oh, and that's apparently a doc typo.  further down if you read the help for integral it states:

integral is a wrapper for quadcc (general scalar integrands), quadgk (integrals with specified integration paths), and quadv (array-valued integrands) that is intended to provide MATLAB compatibility. 

I'll look back through and try a few things and open a bug once i refresh my memory with what's going on.
 


Reply | Threaded
Open this post in threaded view
|

Re: Why is there no function for complex integrals?

BGreen

On Sat, Jun 20, 2020 at 4:05 PM Nicholas Jankowski <[hidden email]> wrote:
Similarly, the documentation states
integral

A compatibility wrapper function that will choose between quadv and quadgk depending on the integrand and options chosen.


oh, and that's apparently a doc typo.  further down if you read the help for integral it states:

integral is a wrapper for quadcc (general scalar integrands), quadgk (integrals with specified integration paths), and quadv (array-valued integrands) that is intended to provide MATLAB compatibility. 

I'll look back through and try a few things and open a bug once i refresh my memory with what's going on.
 

I'd missed that later on. O.K., I'll leave the bug report to you. Thank you!