Finding poles, branch cuts and other complex singularities
- To: mathgroup at smc.vnet.net
- Subject: [mg66808] Finding poles, branch cuts and other complex singularities
- From: "Goyder Dr HGD" <h.g.d.goyder at cranfield.ac.uk>
- Date: Wed, 31 May 2006 06:30:58 -0400 (EDT)
- Sender: owner-wri-mathgroup at wolfram.com
Dear Mathgroup,
I am attempting to develop a numerical method that finds the location of
the singularities of a function on the complex plane. The preliminary
version of the code I give below works, but I have some issues that need
to be resolved. The code integrates the function around a square which
is recursively subdivide if the integral is not zero. This approach
closes-in on singular regions reasonably efficiently. I use two stopping
criteria; one stops the subdivisions if the square side length falls
below eps. The other (more difficult) declares the integral to be zero
if its absolute value is less than a tolerance tol.
Although the code works it gives a flood of warning messages. Can these
be eliminated? I can't just switch them off because I need the message
that the contour has gone through a singularity. (This is a fair warning
because I am looking for singularities and the workaround is easy -just
perturb the origin of my starting square.)
I was expecting the warning that the integrand is zero because this will
occur when the path of integration surrounds an analytic region.
However, breaking down the integral into a sum of two or more sides of
the square does not cure this warning. See below.
I also expect a warning when the path of the integral crosses a branch
cut (integrand has a step change in value). Can this be accommodated by
any options in NIntegrate?
Can I improve my criterion for the integral around a square being zero?
A good approach would be to scale the tolerance by using the maximum and
minimum absolute values of the integrand. (I think the maximum and
minimum values for the whole square must occur on the boundary.) This
would require the method to store the function evaluations made by
NIntegrate. How could this be done?
Now for some maths questions which stretch my understanding of complex
functions my experience being more applied than mathematical. In the
plots I give below I can distinguish poles from essentially
singularities because the value of the integral around a pole is
independent of the size of my square but the integral around an
essential singularity just gets bigger (at least always changes) if I
make my squares smaller. Am I correct in thinking that this is a
definitive way of distinguishing poles from essential singularities?
What can be said about branch cuts? Other than just making my squares
smaller how can I distinguish an isolated branch-cut, which is very
short, from a pole? Are there essential branch-cuts where the jump
across a branch-cut depends on the resolution of my squares? I think
there is a theorem that states that at least some types of branch-cut
can be considered as a limit of a distribution of poles. Does this only
apply to a limited type of branch cut?
The following code sets up the preliminary scheme and plots results.
Expect warning messages.
All comments and solutions thankfully received.
Hugh Goyder
Cranfield University
Tel: +44 (0) 1793 785122
ClearAll[ContourIntegral];
ContourIntegral::usage = " ContourIntegral[f, {z, z0, ds}]Integrates the
function f around the square with side length ds and bottom left
coordinates z0";
ContourIntegral[f_, {z_, z0_, ds_}] := NIntegrate[f, {z, z0, z0 + ds, z0
+ ds + I*ds, z0 + I*ds, z0}]
ClearAll[RecursiveSubDivide];
RecursiveSubDivide::usage = "RecursiveSubDivide[f,{z,z0,ds},{esp,tol}]
recursively subdivides a square into smaller squares. Squares which
contain no singularities \
are not subdivided. Subdivision stops when the side of a square is less
than eps. A square is deemed to contain no singularities if the absolute
value of the \
integral around the square is less than tol. Information on each square
which is not subdivided is output in the form {bottom left coordinate,
side length, value \
of integral around square}";
RecursiveSubDivide[f_, {z_, z0_, ds_}, {eps_, tol_}] := Module[{int,
res, ds2, z1, z2, z3, z4, r1, r2, r3, r4},
ds2 = 0.5*ds; {z1, z2, z3, z4} = {z0, z0 + ds2, z0 + ds2 + I*ds2, z0
+ I*ds2}; If[i1 = ContourIntegral[f, {z, z1, ds2}]; Abs[i1] <= tol ||
Abs[ds2] <= eps,
r1 = {z1, ds2, i1}, r1 = RecursiveSubDivide[f, {z, z1, ds2}, {eps,
tol}]]; If[i2 = ContourIntegral[f, {z, z2, ds2}]; Abs[i2] <= tol ||
Abs[ds2] <= eps,
r2 = {z2, ds2, i2}, r2 = RecursiveSubDivide[f, {z, z2, ds2}, {eps,
tol}]]; If[i3 = ContourIntegral[f, {z, z3, ds2}]; Abs[i3] <= tol ||
Abs[ds2] <= eps,
r3 = {z3, ds2, i3}, r3 = RecursiveSubDivide[f, {z, z3, ds2}, {eps,
tol}]]; If[i4 = ContourIntegral[f, {z, z4, ds2}]; Abs[i4] <= tol ||
Abs[ds2] <= eps,
r4 = {z4, ds2, i4}, r4 = RecursiveSubDivide[f, {z, z4, ds2}, {eps,
tol}]]; {r1, r2, r3, r4}]
f1 = Log[z] + Sqrt[I*z] + 0.1/(z - 0.5 - 0.5*I) + (z - 1.01 - 1.5*I)/(z
- 1. - 1.5*I) + Log[(z - 1.5 + I)/(z - 1 + I)] + E^(1/(30*(z + 1 - I)));
Timing[e1 = RecursiveSubDivide[f1, {z, -2.001 - 2.001*I, 4.01}, {0.01,
0.0001}]; ]
e2 = Partition[Flatten[e1], 3];
ClearAll[z, dz, val];
e3 = e2 /. {z_, dz_, val_} :> Line[{{Re[z], Im[z]}, {Re[z] + dz, Im[z]},
{Re[z] + dz, Im[z] + dz}, {Re[z], Im[z] + dz}, {Re[z], Im[z]}}];
Show[Graphics[{e3}, Frame -> True, PlotRange -> All, AspectRatio ->
Automatic, TextStyle -> {FontFamily -> "Times", FontSize -> 12},
FrameLabel -> {"Real", "Imaginary"}]];
ClearAll[z, dz, val];
e3 = e2 /. {{z_, dz_ /; Abs[dz] <= 0.01, val_} :> Line[{{Re[z], Im[z]},
{Re[z] + dz, Im[z]}, {Re[z] + dz, Im[z] + dz}, {Re[z], Im[z] + dz},
{Re[z], Im[z]}}],
{z_, dz_ /; Abs[dz] > 0.0001, val_} :> {}};
Show[Graphics[{e3}, Frame -> True, AspectRatio -> Automatic, PlotRange
-> {{-2, 2}, {-2, 2}}, FrameLabel -> {"Real", "Imaginary"},
TextStyle -> {FontFamily -> "Times", FontSize -> 12}]];
e2 = Partition[Flatten[e1], 3];
e3 = e2 /. {z_, dz_, val_} :> Cuboid[{Re[z], Im[z], 0}, {Re[z] + dz,
Im[z] + dz, Abs[val]}];
Show[Graphics3D[{e3}, BoxRatios -> {1, 1, 1}, Axes -> True, PlotRange ->
All, AxesLabel -> {"Real", "Imaginary", "Abs[Integral]"},
TextStyle -> {FontFamily -> "Times", FontSize -> 12}]];
(* This version fails to remove zero integral warnings
ClearAll[ContourIntegral];
ContourIntegral[f_, {z_, z0_, ds_}] := NIntegrate[f, {z, z0, z0 +
ds}] + NIntegrate[f, {z, z0 + ds, z0 + ds + I ds}] + NIntegrate[f, {
z, z0 + ds + I ds, z0 + I ds}] + NIntegrate[f, {z, z0 + I ds, z0}]
*)
--
This message has been scanned for viruses and
dangerous content by the Cranfield MailScanner, and is
believed to be clean.