MATLAB Examples

# Roots of a complex function via Cauchy integrals

Nick Trefethen, September 2011

(Chebfun example roots/ComplexRoots.m)

```function ComplexRoots ```

Poles and zeros of complex functions can be located by the evaluation of contour integrals, as mentioned in Chapter 5 of the Chebfun Guide. For example, suppose we have a function like this one with a single root s1 in the unit disk:

```ff = @(z) (z-0.5i).*exp(z); ```

We can find the root as the value of a contour integral around the unit circle:

```s = (1/2i*pi) INT z (f'(z)/f(z)) dz
```

Since Chebfun works with real independent variables, we parametrize the unit circle by a real variable t on [-1,1]:

```z = chebfun('exp(1i*pi*t)'); ```

which gives us

```s1 = (1/2i*pi) INT z ( (df/dt)(dt/dz) / f ) (dz/dt) dt
```
`    = (1/2i*pi) INT z ( (df/dt) / f ) dt`

So here is the Chebfun evaluation:

```f = ff(z); s1 = sum(z.*diff(f)./f)/(2i*pi) ```
```s1 = 0.0000 + 0.5000i ```

There is nothing in this computation that depends on the use of the unit disk. Other contours are equally tracatable in Chebfun, as illustrated in the Example complex/KeyholeContour and in Chapter 5 of the Chebfun Guide.

This method of finding a single root goes back at least to McCune in 1966 [3]. In practice we would often want to be able to find multiple roots, and a generalized algorithm for this case was published by Delves and Lyness in 1967 [1], with mathematical origins as far back as Jackson in 1917 [2]. The idea here is that if f has more than one root in the unit disk, then the value s above comes out as the sum of all these roots. Similarly

```s2 = (1/2i*pi) INT z^2 (f'(z)/f(z)) dz
```

is the sum of the squares of the roots, the analogous formula for s3 with a factor z^3 gives the sum of the cubes, and so on. And a count of the number of roots is given by

```s0 = (1/2i*pi) INT (f'(z)/f(z)) dz
```

(this is basically the argument principle). So for example we can count the number of roots of cosh(pi*z) in the unit disk like this:

```ff = @(z) cosh(pi*z); f = ff(z); s0 = sum(diff(f)./f)/(2i*pi) ```
```s0 = 2.0000 - 0.0000i ```

Here are the sum of the roots and the sum of their squares:

```s1 = sum(z.*diff(f)./f)/(2i*pi) s2 = sum(z.^2.*diff(f)./f)/(2i*pi) ```
```s1 = 1.6179e-15 + 9.5182e-17i s2 = -0.5000 - 0.0000i ```

corresponding to roots at +- 0.5i. We can find these numbers systematically by noting that the monic polynomial p(z) with these roots has coefficients c0=(s1^2-s2)/2, c1=-s1, c2=1. So here is a calculation of the two roots in the unit disk of cosh(pi*z):

```p = [1 -s1 (s1^2-s2)/2]; roots(p) ```
```ans = 0.0000 + 0.5000i 0.0000 - 0.5000i ```

Generalization to higher numbers of roots can be done via Newton's identities. We don't pursue the general case here but instead write a code that finds three roots of an analytic function in the unit disk:

```function r = roots3(ff) % find 3 roots of ff in unit disk z = chebfun('exp(1i*pi*t)'); f = ff(z); s0 = sum(diff(f)./f)/(2i*pi); s1 = sum(z.*diff(f)./f)/(2i*pi); s2 = sum(z.^2.*diff(f)./f)/(2i*pi); s3 = sum(z.^3.*diff(f)./f)/(2i*pi); p = [1 -s1 (s1^2-s2)/2 -(s1^3-3*s1*s2+2*s3)/6]; r = roots(p); end ```

Here is an example:

```ff = @(z) cosh(exp(z)).*(z-.3).*(1+4*z.^2); roots3(ff) ```
```ans = 0.0000 + 0.5000i 0.0000 - 0.5000i 0.3000 - 0.0000i ```

Here is another:

```ff = @(z) (z.^3-1/8).*exp((-1-2i)*z); roots3(ff) ```
```ans = -0.2500 - 0.4330i -0.2500 + 0.4330i 0.5000 - 0.0000i ```
```end ```

References

[1] L. M. Delves and J. N. Lyness, A numerical method for lcoating the zeros of an analytic function, Mathematics of Computation 21 (1967), 543-560.

[2] D. Jackson, Roots and singular points of analytic functions, Annals of Matheamtics 19 (1917), 142-151.

[3] J. E. McCune, Exact inversion of dispersion relations, Physics of Fluids 9 (1966), 2082-2084.