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.