## RootsAndPoles.jl

Julia implementation of the global complex root and pole finding (GRPF) algorithm.
Author fgasdia
Popularity
18 Stars
Updated Last
1 Year Ago
Started In
October 2018

# RootsAndPoles.jl: Global complex Roots and Poles Finding in Julia

A Julia implementation of GRPF by Piotr Kowalczyk.

## Description

`RootsAndPoles.jl` attempts to find all the zeros and poles of a complex valued function with complex arguments in a fixed region. These types of problems are frequently encountered in electromagnetics, but the algorithm can also be used for similar problems in e.g. optics, acoustics, etc.

The GRPF algorithm first samples the function on a triangular mesh through Delaunay triangulation. Candidate regions to search for roots and poles are determined and the discretized Cauchy's argument principle is applied without needing the derivative of the function or integration over the contour. To improve the accuracy of the results, a self-adaptive mesh refinement occurs inside the identified candidate regions.

## Usage

### Installation

`]add RootsAndPoles`

### Example Problem

Consider a simple rational (complex) argument function `simplefcn(z)` for which we seek roots and poles.

```function simplefcn(z)
w = (z - 1)*(z - im)^2*(z + 1)^3/(z + im)
end```

Next, we define parameters for the initial complex domain over which we would like to search.

```xb = -2  # real part begin
xe = 2  # real part end
yb = -2  # imag part begin
ye = 2  # imag part end
r = 0.1  # initial mesh step
tolerance = 1e-9```

This package includes convenience functions for rectangular and disk shaped domains, but any "shape" can be used. `origcoords` below is simply a vector of complex numbers containing the original mesh coordinates which will be Delaunay triangulated. For maximum efficiency, the original mesh nodes should form equilateral triangles.

```using RootsAndPoles

origcoords = rectangulardomain(complex(xb, yb), complex(xe, ye), r)```

Roots and poles can be obtained with the `grpf` function. We only need to pass the handle to our `simplefcn` and the `origcoords`.

`zroots, zpoles = grpf(simplefcn, origcoords)`

Additional parameters can be provided to the tessellation and GRPF algorithms by explicitly passing a `GRPFParams` struct.

The two most useful parameters are `tess_sizehint` for the final total number of nodes in the internal `DelaunayTessellation2D` object and the root finder `tolerance` at which the mesh refinement stops. Just like `sizehint!` for other collections, setting `tess_sizehint` to a value approximately equal to the final number of nodes in the tessellation can improve performance. `tolerance` is the largest triangle edge length of the candidate edges defined in the `origcoords` domain. In practice, the root and pole accuracy is a larger value than `tolerance`, so `tolerance` needs to be set smaller than the desired tolerance on the roots and poles.

By default, the value of `tess_sizehint` is 5000 and the `tolerance` is 1e-9, but they can be specified by providing the `GRPFParams` argument

`zroots, zpoles = grpf(simplefcn, origcoords, GRPFParams(5000, 1e-9))`

Beginning version `v1.1.0`, calls to the provided function, e.g. `simplefcn`, can be multithreaded using Julia's `@threads` capability. The function is called at every node of the triangulation and the results should be independent of one another. For fast-running functions like `simplefcn`, the overall runtime of `grpf` is dominated by the Delaunay Triangulation itself, but for more complicated functions, threading can provide a significant advantage. To enable multithreading of the function calls, specify so as a `GRPFParams` argument

`zroots, zpoles = grpf(simplefcn, origcoords, GRPFParams(5000, 1e-9, true))`

By default, `multithreading = false`.

Additional parameters which can be controlled are `maxiterations`, `maxnodes`, and `skinnytriangle`. `maxiterations` sets the maximum number of mesh refinement iterations and `maxnodes` sets the maximum number of nodes allowed in the `DelaunayTessellation2D` before returning. `skinnytriangle` is the maximum allowed ratio of the longest to shortest side length in a tessellation triangle before the triangle is automatically subdivided in the mesh refinement step. Default values are

• `maxiterations`: 100
• `maxnodes`: 500000
• `skinnytriangle`: 3

These can be specified along with the `tess_sizehint`, `tolerance` and `multithreading` as, e.g.

`zroots, zpoles = grpf(simplefcn, origcoords, GRPFParams(100, 500000, 3, 5000, 1e-9, true))`

### Plot data

If mesh node `quadrants` and `phasediffs` are wanted for plotting, simply pass a `PlotData()` instance.

`zroots, zpoles, quadrants, phasediffs, tess, g2f = grpf(simplefcn, origcoords, PlotData())`

A `GRPFParams` can also be passed.

`zroots, zpoles, quadrants, phasediffs, tess, g2f = grpf(simplefcn, origcoords, PlotData(), GRPFParams(5000, 1e-9))`

In `v1.4.0`, a function `getplotdata` makes plotting the tessellation results more convenient. The code below was used to create the figure shown at the top of the page.

```zroots, zpoles, quadrants, phasediffs, tess, g2f = grpf(simplefcn, origcoords, PlotData());
z, edgecolors = getplotdata(tess, quadrants, phasediffs, g2f);

using Plots

pal = ["yellow", "purple", "green", "orange", "black", "cyan"]
plot(real(z), imag(z), group=edgecolors, palette=pal,
xlabel="Re(z)", ylabel="Im(z)",
xlims=(-2, 2), ylims=(-2, 2),
legend=:outerright, legendtitle="f(z)", title="Simple rational function",
label=["Re > 0, Im ≥ 0" "Re ≤ 0, Im > 0" "Re < 0, Im ≤ 0" "Re ≥ 0, Im < 0" "" ""])```

## Limitations

This package uses VoronoiDelaunay.jl to perform the Delaunay tessellation. `VoronoiDelaunay` is numerically limited to the range of `1.0+eps(Float64)` to `2.0-2eps(Float64)` for its point coordinates. `RootsAndPoles.jl` will accept functions and `origcoords` that aren't limited to `Complex{Float64}`, for example `Complex{BigFloat}`, but the internal tolerance of the root finding is limited to `Float64` precision.

## Citing

Please consider citing Piotr's publications if this code is used in scientific work:

1. P. Kowalczyk, “Complex Root Finding Algorithm Based on Delaunay Triangulation”, ACM Transactions on Mathematical Software, vol. 41, no. 3, art. 19, pp. 1-13, June 2015. https://dl.acm.org/citation.cfm?id=2699457

2. P. Kowalczyk, "Global Complex Roots and Poles Finding Algorithm Based on Phase Analysis for Propagation and Radiation Problems," IEEE Transactions on Antennas and Propagation, vol. 66, no. 12, pp. 7198-7205, Dec. 2018. https://ieeexplore.ieee.org/document/8457320

We also encourage you to cite this package if used in scientific work. Refer to the Zenodo DOI at the top of the page or CITATION.bib.