Integrating complex functions

[jondea]

In the README you recommend integrating complex functions as vectors with PAIRED error. Would you be averse to having a simple Julia wrapper which accepts functions which return complex values? As in does it go against the library philosophy or is it a matter of time? If it is the latter I would be happy to contribute.

[stevengj]

I don't want a whole other set of cubature functions. I guess we could have a wrapper

complexfunc(f) = (x,v) -> let c = f(x); v[1] = real(c); v[2] = imag(c); end
vec2complex(v) = Complex(v[1],v[2])

So that they you could do e.g. vec2complex(hquadrature(2, complexfunc(f), ...)[1]).

[jondea]

I suppose there is no way to do a clever dispatch based on what f returns (or even the type of arguments it takes).

[stevengj]

No, unfortunately, Julia doesn't have arrow types.

[jondea]

That's very neat, I have never seen the concept expressed like that, thank you. Also on an off topic note, I must thank you for your online maths resources, which I have found useful on several occasions (namely your PML tutorial and differentiation by FFT).

[dazsmith]

Here's a slight generalisation of the wrapper which also allows for complex intervals of integration:

complexfunc(f) = (x,v) -> let c = f(x); v[1] = real(c); v[2] = imag(c); end;
vec2complex(v) = Complex(v[1],v[2]);
hquadratureCOMPLEX(f,a,b) = vec2complex(hquadrature(2, complexfunc(t -> f(a+t*(b-a))*(b-a)),0,1,abstol=1e-8)[1]);

It must be possible to improve this so that the named argument abstol can be passed to hquadratureCOMPLEX rather than set in the wrapper itself, but I'm new to Julia. I'll update if I figure it out.

[stevengj]

Ultimately, the right way to solve this is to rewrite in pure Julia.

[stevengj]

See https://github.com/stevengj/HCubature.jl for a pure-Julia re-implementation of hcubature that handles complex integrands, etc.