[julia-users] Re: finding all roots of complex function
Say I would like to do the integral using quadgk. There are 2 caveats, however. Quoting from the Julia manual These quadrature rules work best for smooth functions within each interval, so if your function has a known discontinuity or other singularity, it is best to subdivide your interval to put the singularity at an endpoint. My function does have singularities (i.e. poles) inside of the integration contour, and it is not clear to me how I should split the interval. Am I missinterpreting this? The other issue is that quadgk seems to only accept function of one variable, wheres mine are 2-variable functions. On Thursday, July 17, 2014 7:53:38 PM UTC+2, Steven G. Johnson wrote: On Thursday, July 17, 2014 9:19:04 AM UTC-4, Andrei Berceanu wrote: I should perhaps mention that this is part of a bigger scheme, to first find all the poles of G(x,y)/F(x,y) and then use the residue theorem for solving a complex integral of the type integral( G(x,y)/F(x,y), (x,y)) Unless F(x,y) is very special (e.g. a polynomial), I suspect that it would be much faster to just do the integral. Since you have analytic functions, 1d/contour integration is very efficient (with an appropriate algorithm, e.g. the built-in quadgk function) unless you have poles lying very close to your integration contour (and even then it is not too bad with an adaptive scheme like quadgk.) In contrast, finding *all* the zeros of an arbitrary analytic function is hard, usually harder than integrating it unless you have a good idea of where the zeros are. In general, it's not practical to guarantee that you have found all the zeros unless you can restrict your search to some finite portion of the complex plane. For finding the roots of analytic functions inside a given contour, some of the best algorithms actually involve doing a sequence of integrals ( http://www.chebfun.org/examples/roots/ComplexRoots.html) that are just as hard as your original integral above. So, you might as well just do the integral to begin with.
[julia-users] Re: finding all roots of complex function
First, I think, Steven is talking about line integrals along a path, and
these integrals are one-dimensional.
And secondly, of course you should avoid poles along your path, by making
the radius big enough, or similar. These things are standard techniques in
complex analysis.
As an example, I append a function `line_integral` from my package
NumericalMath (not in METADATA) that does integrate along a path given by
`points`.
function line_integral(f::Function, points::Array{Complex{Float64},1}; tol =
1e-15)
local n::Int = length(points)
local Q::Complex = 0.0 + 0.0im
if n == 1; return Q; end
for i = 2:n
a = points[i-1]
b = points[i]
d = b - a
f1(t) = real(f(a + t*d))
f2(t) = imag(f(a + t*d))
Qre = quadgk(f1, 0.0, 1.0)[1]
Qim = quadgk(f2, 0.0, 1.0)[1]
Q += d * (Qre + Qim*1.0im)
end
return Q
end
On Friday, July 18, 2014 11:48:55 AM UTC+2, Andrei Berceanu wrote:
Say I would like to do the integral using quadgk. There are 2 caveats,
however.
Quoting from the Julia manual
These quadrature rules work best for smooth functions within each
interval, so if your function has a known discontinuity or other
singularity, it is best to subdivide your interval to put the singularity
at an endpoint.
My function does have singularities (i.e. poles) inside of the integration
contour, and it is not clear to me how I should split the interval. Am I
missinterpreting this?
The other issue is that quadgk seems to only accept function of one
variable, wheres mine are 2-variable functions.
On Thursday, July 17, 2014 7:53:38 PM UTC+2, Steven G. Johnson wrote:
On Thursday, July 17, 2014 9:19:04 AM UTC-4, Andrei Berceanu wrote:
I should perhaps mention that this is part of a bigger scheme, to first
find all the poles of G(x,y)/F(x,y) and then use the residue theorem for
solving a complex integral of the type
integral( G(x,y)/F(x,y), (x,y))
Unless F(x,y) is very special (e.g. a polynomial), I suspect that it
would be much faster to just do the integral. Since you have analytic
functions, 1d/contour integration is very efficient (with an appropriate
algorithm, e.g. the built-in quadgk function) unless you have poles lying
very close to your integration contour (and even then it is not too bad
with an adaptive scheme like quadgk.)
In contrast, finding *all* the zeros of an arbitrary analytic function is
hard, usually harder than integrating it unless you have a good idea of
where the zeros are. In general, it's not practical to guarantee that you
have found all the zeros unless you can restrict your search to some finite
portion of the complex plane. For finding the roots of analytic functions
inside a given contour, some of the best algorithms actually involve doing
a sequence of integrals (
http://www.chebfun.org/examples/roots/ComplexRoots.html) that are just
as hard as your original integral above. So, you might as well just do
the integral to begin with.
[julia-users] Re: finding all roots of complex function
On Friday, July 18, 2014 8:06:25 AM UTC-4, Hans W Borchers wrote:
First, I think, Steven is talking about line integrals along a path, and
these integrals are one-dimensional.
And secondly, of course you should avoid poles along your path, by making
the radius big enough, or similar. These things are standard techniques in
complex analysis.
Note also that if you have poles on your path, then you need more
information in order to define what your integral even means. e.g. do you
mean the limit as the pole approaches the path (or vice versa), and if so
you need to specify from what side you are taking the limit. Then you can
(for example) deform your integration contour away from the pole (in the
correct direction) as Hans suggests.
As an example, I append a function `line_integral` from my package
NumericalMath (not in METADATA) that does integrate along a path given by
`points`.
Hans, note that quadgk already supports line integration over straight
line-segments in the complex plane along any arbitrary sequence
points::Array{Complex} of points. Just do:
quadgk(f, points...)
This is better than doing the integral over each line segment separately,
because it is globally adaptive: in order to decide whether it needs to
subdivide any interval to improve accuracy, it will look at *all* the
intervals, not just the intervals within a single line segment.
Also, there is no need to integrate the real and imaginary parts
separately. quadgk can integrate complex-valued functions (and
vector-valued functions, and matrix-valued functions, and...) just fine.
By default, it measures the error in the usual norm sqrt(real^2 + imag^2),
but if you want to use min(|real|,|imag|) instead (which is effectively
what you are doing if you integrate them separately), then you can pass
that norm for the norm keyword argument.
--SGJ
[julia-users] Re: finding all roots of complex function
Okay, thanks, what was not so clear to me is how to use the three dots
notation to hand some points over to the integration function. For example,
to integrate 1/z along the unit circle around the pole at the origin:
julia x = linspace(0, 2*pi);
julia points = exp(x*1im);
julia quadgk(z - 1 ./ z, points...)
(-2.1630337481358435e-17 + 6.283185307179583im,5.957695483631124e-16)
That is the correct solution 2 \pi i to a very good accuracy.
On Friday, July 18, 2014 6:53:58 PM UTC+2, Steven G. Johnson wrote:
On Friday, July 18, 2014 8:06:25 AM UTC-4, Hans W Borchers wrote:
First, I think, Steven is talking about line integrals along a path, and
these integrals are one-dimensional.
And secondly, of course you should avoid poles along your path, by making
the radius big enough, or similar. These things are standard techniques in
complex analysis.
Note also that if you have poles on your path, then you need more
information in order to define what your integral even means. e.g. do you
mean the limit as the pole approaches the path (or vice versa), and if so
you need to specify from what side you are taking the limit. Then you can
(for example) deform your integration contour away from the pole (in the
correct direction) as Hans suggests.
As an example, I append a function `line_integral` from my package
NumericalMath (not in METADATA) that does integrate along a path given by
`points`.
Hans, note that quadgk already supports line integration over straight
line-segments in the complex plane along any arbitrary sequence
points::Array{Complex} of points. Just do:
quadgk(f, points...)
This is better than doing the integral over each line segment separately,
because it is globally adaptive: in order to decide whether it needs to
subdivide any interval to improve accuracy, it will look at *all* the
intervals, not just the intervals within a single line segment.
Also, there is no need to integrate the real and imaginary parts
separately. quadgk can integrate complex-valued functions (and
vector-valued functions, and matrix-valued functions, and...) just fine.
By default, it measures the error in the usual norm sqrt(real^2 + imag^2),
but if you want to use min(|real|,|imag|) instead (which is effectively
what you are doing if you integrate them separately), then you can pass
that norm for the norm keyword argument.
--SGJ
[julia-users] Re: finding all roots of complex function
On Friday, July 18, 2014 3:08:34 PM UTC-4, Hans W Borchers wrote: Okay, thanks, what was not so clear to me is how to use the three dots notation to hand some points over to the integration function. (The three dots are called splicing or splatting, see http://julia.readthedocs.org/en/latest/manual/functions/#varargs-functions)
[julia-users] Re: finding all roots of complex function
Also, there is the cis(x) function which does the same thing as exp(x*im). But I notice that it isn't (yet?) vectorized, so before using it with an array argument you must add @vectorize_1arg Number Base.cis On Friday, July 18, 2014 12:46:54 PM UTC-7, Steven G. Johnson wrote: On Friday, July 18, 2014 3:08:34 PM UTC-4, Hans W Borchers wrote: julia x = linspace(0, 2*pi); julia points = exp(x*1im); julia quadgk(z - 1 ./ z, points...) (-2.1630337481358435e-17 + 6.283185307179583im,5.957695483631124e-16) Note that linspace defaults to interpolating 100 points; this is pretty inefficient here, because you want quadgk to do the subdivision more efficiently for you. I would just integrate over a triangle (4 points, since the beginning point is included twice): julia quadgk(z - 1 / z, exp(linspace(0, 2pi, 4) * im)...)[1] - 2pi*im -1.1102230246251565e-16 - 8.881784197001252e-16im Note also that it is fine to use 1/z rather than 1./z (although the latter is harmless). Unlike Matlab, quadgk evaluates your function for one point at a time, so you don't need to vectorize it.(In Matlab, evaluating your integrand for a vector of arguments is essential for performance, but it is an unnecessary annoyance in Julia.)
[julia-users] Re: finding all roots of complex function
On Friday, July 18, 2014 4:11:34 PM UTC-4, Peter Simon wrote: Also, there is the cis(x) function which does the same thing as exp(x*im). But I notice that it isn't (yet?) vectorized, Should be fixed now in git master.
[julia-users] Re: finding all roots of complex function
Note that linspace defaults to interpolating 100 points; this is pretty inefficient here, ... Also, there is the cis(x) function which does the same thing as exp(x*im). I know. Please not that this was meant as a test case only. And I wanted to the integration along the unit circle, not a triangle or square. Much better for testing... By the way: Sometimes I feel there are too many functions in the base system. Wouldn't it be more appropriate to put more functions of special areas into their own packages, as for instance Python does with the SciPy module (and others). On Friday, July 18, 2014 10:36:07 PM UTC+2, Steven G. Johnson wrote: On Friday, July 18, 2014 4:11:34 PM UTC-4, Peter Simon wrote: Also, there is the cis(x) function which does the same thing as exp(x*im). But I notice that it isn't (yet?) vectorized, Should be fixed now in git master.
[julia-users] Re: finding all roots of complex function
On Friday, July 18, 2014 6:01:07 PM UTC-4, Hans W Borchers wrote: Wouldn't it be more appropriate to put more functions of special areas into their own packages, as for instance Python does with the SciPy module (and others). I think that the thinking goes: 1) breaking Base into submodules doesn't seem to actually solve many real problems, so far (since multiple-dispatch and the module system for external stuff avoids most problems) 2) the Matlab style of having a lot of functionality in one flat namespace may be ugly from a purist CS viewpoint, but it is really convenient for users: basic things just work, and tab completion finds them (if the names are guessable/standard enough), without having to hunt for the module that has inv or besselj or replace or readdlm or rand etc.
[julia-users] Re: finding all roots of complex function
I should perhaps mention that this is part of a bigger scheme, to first find all the poles of G(x,y)/F(x,y) and then use the residue theorem for solving a complex integral of the type integral( G(x,y)/F(x,y), (x,y)) On Thursday, July 17, 2014 3:15:45 PM UTC+2, Andrei Berceanu wrote: Hi guys, I would like to know if Julia (itself or though some existing package) can be used to find all the roots of a *complex*, 2 variable function F(x,y). Here x and y are real, so F:R - C. Thanks, Andrei
[julia-users] Re: finding all roots of complex function
On Thursday, July 17, 2014 9:19:04 AM UTC-4, Andrei Berceanu wrote: I should perhaps mention that this is part of a bigger scheme, to first find all the poles of G(x,y)/F(x,y) and then use the residue theorem for solving a complex integral of the type integral( G(x,y)/F(x,y), (x,y)) Unless F(x,y) is very special (e.g. a polynomial), I suspect that it would be much faster to just do the integral. Since you have analytic functions, 1d/contour integration is very efficient (with an appropriate algorithm, e.g. the built-in quadgk function) unless you have poles lying very close to your integration contour (and even then it is not too bad with an adaptive scheme like quadgk.) In contrast, finding *all* the zeros of an arbitrary analytic function is hard, usually harder than integrating it unless you have a good idea of where the zeros are. In general, it's not practical to guarantee that you have found all the zeros unless you can restrict your search to some finite portion of the complex plane. For finding the roots of analytic functions inside a given contour, some of the best algorithms actually involve doing a sequence of integrals (http://www.chebfun.org/examples/roots/ComplexRoots.html) that are just as hard as your original integral above. So, you might as well just do the integral to begin with.
