U.U.D.M. Project Report 2016:42
Examensarbete i matematik, 15 hp Handledare: Warwick Tucker Examinator: Jörgen Östensson Augusti 2016
Enclosing Zeros for Systems of Two Analytic Functions
Joel Dahne
Department of Mathematics
Enclosing Zeros for Systems of Two Analytic Functions
Joel Dahne August 24, 2016
Abstract
We present a method for finding and validating simple zeros to systems of two analytic functions. It can also prove that all zeros have been found by calculating the total number of zeros in the domain. The last part turns out to be very time consuming and not always applicable in practical examples. To illustrate the method we apply it to a non trivial academic example and also on a physical problem for which no effective methods have been found before.
1 Introduction
Finding zeros for analytic functions in one variable is a well studied, and im- portant, problem. For polynomials there exists very fast and accurate methods, see for example [4, 9]. In the case of general analytic functions there are also methods developed. One such method, that also produces validated results, is presented in [5].
For systems of analytic functions the problem of finding approximations of the zeros is not as well studied. In [6] they present a method that can generate approximations of the zeros. They use a multidimensional version of the logarithmic integral together with other methods to do this. However they conclude that proper evaluation of the method is hard, especially evaluation of the integral.
In this paper we will present another method for finding approximations of the zeros. It will also make use of the multidimensional logarithmic integral but combined with a random search for the zeros.
In addition we will make use of validated numerics to achieve rigorous results.
The method should generate a number of small boxes which are proved to contain a zero of the function. We also wish to be able to prove that all zeros are found. A short introduction to validated numerics and how this can be achieved will be given.
2 The General Strategy
We will divide the task of finding all the zeros into two parts.
1. Calculate the number of zeros in the domain
2. Perform an exhaustive random search for these zeros
The calculation for the number of zeros will be done using the logarithmic integral and will give us a stop criterion for the random search. For the random search we will make use of Newton’s method and also a validated version of Newton’s method presented in section 3.4.2. This will give us an increasing number of verified zeros and when the number of zeros is equal to the calculated number of zeros in the domain we can stop.
2.1 Calculating the Number of Zeros
To be able to calculate the number of zeros of the function we will use that it is analytic.
In one complex variable the logarithmic integral gives us a method to find the number of zeros for a function in a domain.
Theorem 1. Let f be analytic on the closure of a bounded domain D ⊂ C with piece wise smooth boundary that does no contain any zeros of f . Then the number of zeros to f , counting multiplicity, is given by
1 2πi
Z
∂D
f
0(z)
f (z) dz. (1)
This can also be extended to meromorphic functions, but we will make no use of that here.
Since we will work with systems of analytic functions we will need a gener- alization of this formula. Such a generalization exists in the literature and is given in the following theorem.
Theorem 2. Let f = f (z
1, ..., z
n) be analytic on the closure of a bounded domain D ⊂ C
nwith piece wise smooth boundary that does not contain any zeros of f . Then f has only isolated zeros in D, and their number, counting multiplicity, is expressed by the formula
N = Z
∂D
ω(f (z), f (z))), (2)
where
ω(f, f ) = (n − 1)!
(2πi)
n1
|f |
2nn
X
j=1
(−1)
j−1f
jdf
[j]∧ df. (3)
For a proof of this we refer to [1].
The integral is given in differential forms and the notation used is df = df
1∧ ... ∧ df
n,
df
[k]= df
1∧ ... ∧ df
k−1∧ df
k+1∧ ... ∧ df
n.
We want to be able to compute the integral numerically and to prepare it for
this we will rewrite it as a Riemann integral which is more suitable for numerical
integration. To accomplish this we will follow the procedure in [6], where more
details can be found.
With the boundary of D given by
∂D = ∂D
[1]∪ ... ∪ ∂D
[n](4)
where
∂D
[k]= D
1× · · · × D
k−1× ∂D
k× D
k+1× · · · × D
n. (5) We can rewrite the formula as
n
X
k=1
Z
∂D[k]
ω
k(f, f ), (6)
where
ω
k= (−1)
k−1(n − 1)!
(2πi)
ndet J
fdet J
[k]|f |
2ndz
[k]∧ dz. (7) Here J
fdenotes the Jacobian with respect to f and J
[k]is the Jacobian matrix with the k th column replaced with [f
1, ..., f
n]
T.
Since we will only be working with systems of two analytic functions we can let n = 2, which gives us
1 4π
2Z
∂D[2]
det J
fdet J
[2]|f |
4dz
1dz
1dz
2− Z
∂D[1]
det J
fdet J
[1]|f |
4dz
2dz
1dz
2. (8) When doing the computations it will be beneficial to consider it as a function of 4 real variables instead of 2 complex. By changing variables from z
1and z
2to x
1+ iy
1and x
2+ iy
2we get
1 4π
2Z
∂D[2]
det J
fdet J
[2]|f |
4(dx
1− idy
1) ∧ (dx
1+ idy
1) ∧ (dx
2+ idy
2)−
− Z
∂D[1]
det J
fdet J
[1]|f |
4(dx
2− idy
2) ∧ (dx
1+ dy
1) ∧ (dx
2+ dy
2)
. (9) Expanding the differentials we get
(dx
1− idy
1) ∧ (dx
1+ idy
1) ∧ (dx
2+ idy
2) = (dx
1∧ dx
1+ i(dx
1∧ dy
1) − i(dy
1∧ dx
1) + dy
1∧ dy
1) ∧ (dx
2+ idy
2) = 2i(dx
1∧ dy
1) ∧ (dx
2+ idy
2) = 2i(dx
1∧ dy
1∧ dx
2+ i(dx
1∧ dy
1∧ dy
2)) (10) and
(dx
2− idy
2) ∧ (dx
1+ idy
1) ∧ (dx
2+ idy
2) =
−(dx
2− idy
2) ∧ (dx
2+ idy
2) ∧ (dx
1+ idy
1) =
−(dx
2∧ dx
2+ i(dx
2∧ dy
2) − i(dy
2∧ dx
2) + dy
2∧ dy
2) ∧ (dx
1+ idy
1) =
−2i(dx
2∧ dy
2) ∧ (dx
1+ idy
1) =
−2i(dx
2∧ dy
2∧ dx
1+ i(dx
2∧ dy
2∧ dy
1)). (11) This gives us the integral formula
i 2π
2Z
∂D[2]
det J
fdet J
[2]|f |
4dx
1dy
1dx
2+ i Z
∂D[2]
det J
fdet J
[2]|f |
4dx
1dy
1dy
2−
− Z
∂D[1]
det J
fdet J
[1]|f |
4dx
2dy
2dx
1− i Z
∂D[1]
det J
fdet J
[1]|f |
4dx
2dy
2dy
1. (12)
2.2 Random Search for Zeros
When the number of zeros is known what is left is to find them. As mentioned we will make use of Newton’s method and a validated version of it to accomplish this. Newton’s method will be used to generate possible zeros, points that are likely to be close to a zero. Then the validated version is applied to check if the points indeed is close to a zero.
For Newton’s method we will start by choosing a random point, z
0, in the domain and then generate the sequence
z
n+1= z
n− J
f(z
n)
−1f (z
n) n = 0, 1, 2, .... (13) which normally converges to a zero of f within a few iterations. This conver- gence is however not guaranteed and we have no information of how good of an approximation to a zero the new point actually is.
Because of this we will need to check and validate the generated point. Using the validated version of Newton’s method we can determine if the generated point indeed is a good approximation of a zeros, more specifically we are able to prove that a small box around the point contains a zero of the function.
If the point is determined to be a good approximation of a zero we add it to the list of zeros, if it is not then it is discarded. We can then continue to generate points and try to validate them until the number of zeros in the list is equal to the calculated number of zeros in the domain.
3 Validated Numerics
Computers are fast, they can perform billions of operations per second and with them we can accomplish tasks otherwise not possible. But when it comes to mathematics, computers are not always reliable. A computer does not know what a real number is and cannot calculate with real numbers. All the computer gives us is an approximation, often times it is good enough though. But to be able to generate mathematical proofs involving the continuum with the help of a computer we need other tools. For this we will use validated numerics and we will here give a short introduction to it, for a more complete one see for example [10].
For simplicity we will in this introduction consider validated numerics only on real numbers. The step from R to C
2is not very big and a short discussion of it will be given in the end.
3.1 Floating Points
To start with we will give an introduction to floating point numbers, the comput- ers way of handling real numbers. We will not go into how they are constructed or how the computations are done, but will instead talk about some of their properties and the problems they cause.
The floating point numbers are a finite subset of the real numbers and de-
noted by F. Around the origin they are rather dense and then become more
and more sparse the further from the origin, see figure 1. The exact number of
floating points depends on which format is used.
Figure 1: The distribution of floating point numbers.
The problem with computing with floating points is that when the exact result of a computation cannot be represented in the format we have to round off to the closest number that we can represent. Even though the round off for one computation may be small we get problems when doing many computations together, each one introducing a small error.
Example 1. Consider the first one million terms in the harmonic series. For a mathematician it is clear that
1 1 + 1
2 + 1
3 + · · · + 1 10
6= 1
10
6+ · · · + 1 2 + 1
1 , (14)
because addition is associative and commutative. But computing this on a com- puter using a single precision format we get
1 1 + 1
2 + 1
3 + · · · + 1
10
6=14.357357, 1
10
6+ · · · + 1 3 + 1
2 + 1
1 =14.392651.
Clearly they are not equal and it is hard to tell how far from the real answer they are. It is obvious that this is not a property that is very good when trying to prove something about the result.
In the following sections we will introduce a method that will help us keep track of the round off error and it will allow us to actually prove things about the result.
3.2 Interval Arithmetic
At the heart of validated numerics lies the concept of intervals and interval arithmetic, essentially the use of intervals will let us keep track of the errors introduced in the computations. We will begin by giving an elementary intro- duction to interval arithmetic on R.
Let IR denote the set of all closed intervals on R, i.e. IR = {[a, b] : a ≤ b, a, b ∈ R}. We will use the notation x = [x, x] to denote an interval and its endpoints x, x respectively.
Let ? denote one of the standard arithmetic operations, i.e. ? ∈ {+, −, ×, ÷}.
We can then define x ? y, for x, y ∈ IR, by x ? y = {x ? y : x ∈ x, y ∈ y} where
we demand that 0 6∈ y for division. What might not be immediately clear is
that this will always return a new interval. The endpoints for this new interval can be characterized by
x + y =[x + y, x + y], x − y =[x − y, x − y],
x × y =[min(xy, xy, xy, xy), max(xy, xy, xy, xy), x ÷ y =x × [1/y, 1/y].
This gives us an extension of the standard arithmetic operations on R to IR.
Next we want to be able to extend functions on R to functions on IR. We want to be able to enclose the range of the function by an interval.
Definition 1. Let D ⊂ R, and consider a function f : D → R. We define the range of f over D to be the set
R(f, D) = {f (x) : x ∈ D}.
However, except in very simple cases, finding the exact range of a function is hard. There is an entire branch in mathematics – optimization theory – devoted to finding only the smallest element in that set.
For rational functions, functions on the form f (x) = p(x)/q(x) where p and q are polynomials, we can create an extension by simply substituting all instances of the real variable x with the interval variable x. This produces a rational interval function, F (x) called the natural interval extension of f .
We also want find interval extensions of more general functions. For mono- tonic functions this is easy, for example we may define
e
x= exp x = [e
x, e
x].
Other functions can also be extended. Especially the set of standard functions S ={a
x, log
ax, x, x
p/q, abs x, sin x, cos x, tan x, ...
..., sinh x, cosh x, tanh x, arcsin x, arccos x, arctan x}
can be extended to interval functions. By doing this we can find interval ex- tensions to all functions expressed as a finite number of standard functions combined with arithmetic operations, constants and composition. These kind of functions are often called elementary functions and the set of all elementary functions is denoted by E.
Theorem 3 (The Fundamental Theorem of Interval Analysis). Given an elementary function f and a natural interval extension F such that F (x) is well-defined for some x ∈ IR, we have
(1)z ⊆ z
0⊆ x =⇒ F (z) ⊆ F (z
0)(inclusion isotonicity) (2)R(f ; x) ⊆ F (x)(range enclosure)
While finding the exact range of a function is hard we can, by evaluating
the interval extension of function, find an enclosure of the range. With the help
of Theorem 2 we can then say something about what is not in the range of the
function.
Example 2. Let f (x) = (sin x − x
2+ 1) cos x. Prove that f has no zeros in [0,
12].
By substituting x with x = [0,
12] and f with its natural interval extension F we end up with
F ([0, 1
2 ]) = (sin[0, 1
2 ] − [0, 1
2 ]
2+ 1) cos[0, 1 2 ]
= ([0, sin 1
2 ] − [0, 1
4 ] + [1, 1])[cos 1 2 , 1]
= [ 3
4 , 1][cos 1
2 , 1] = [ 3 4 cos 1
2 , 1 + sin 1 2 ]
⊆ [0.65818, 1.4795].
Since 0 6∈ F ([0,
12]), it follows that 0 6∈ R(f, [0,
12]) and f has no roots in the domain.
3.3 Directed Rounding
We have now discussed the basics of interval arithmetic using real numbers.
However, since we are doing computations on a computer, we will in practice be limited to floating point numbers. Therefor we will discuss how we can limit the domain to floating points and still keep the wanted properties such as inclusion isotonicity and range enclosure.
Instead of working with real intervals, IR, we will now work with floating point intervals, IF = {[x, x] : x ≤ x, x, x ∈ F}, recall that F is the set of floating point numbers. To do this we will use two modes of rounding, round towards minus infinity and round towards infinity denoted by 5(x) and 4(x). In general the left point of the interval will be rounded towards minus infinity and the right points towards infinity. For example the addition will be implemented as
x + y = [5(x + y), 4(x + y)]. (15)
This ensures us that the floating point interval contains the corresponding real interval.
Example 3. Lets, as in example 1, again consider the first one million terms in the harmonic series. This time we add them together using directed rounding instead and get an interval proved to contain the exact result. For the lower part of the interval all roundings were done towards minus infinity, giving us a lower bound for the result, and for the upper part towards infinity, giving us an upper bound.
1 1 + 1
2 + 1
3 + · · · + 1
10
6∈[13.971489, 14.924348], 1
10
6+ · · · + 1 3 + 1
2 + 1
1 ∈[14.350338, 14.436089].
Both of them give us an interval containing the exact result, but the later is to prefer since it gives a narrower interval.
With this directed rounding together with interval arithmetic we are able
have the computer generate a guaranteed enclosure of a function even though
we are limited to floating point numbers.
3.4 Validated Methods
Using interval arithmetic we have a method of enclosing the exact value of a computation, with the help of this we can get rigorous results. Most numerical methods however does not only suffer from floating point errors but also dis- cretization errors from the method itself. To truly get rigorous results we need adapt the methods we use to account for the discretization error as well. We will describe two validated methods, one for computing definite integrals and one for finding zeros, that will be used in the coming algorithm.
3.4.1 Computing Integrals
There are several methods developed for computing definite integrals. From the very simple midpoint method and trapezoid method to the slightly more complex Simpson’s method. Even tough all of them many times give a good approximation of the integral they are not rigorous. We present a method that can give us a guaranteed enclosure of the integral.
Let f be a real valued function which admits a well-defined interval extension F over the integration domain [a, b]. We then have the enclosure of the integral I
I ∈ w([a, b])F ([a, b]), (16)
where w([a, b]) = b − a denotes the width of the interval. In most cases this will however produce a terrible estimate of the integral. By partitioning the domain [a, b] into several smaller intervals we can produce a much better estimate. Given a partition a = x
0≤ x
1≤ · · · ≤ x
N= b we can decompose the domin,
[a, b] = [x
0, x
1] ∪ [x
1, x
2] ∪ · · · ∪ [x
N −1, x
N]. (17) Now we get the estimate
I ∈
N
X
i=1
w([x
i−1, x
i])F ([x
i−1, x
i]). (18)
It can be noted that this method is very similar to the definition of the Riemann- Stieltjes integral.
Even if this method can be used to generate a tight enclosure of the integral it converges very slowly. It can be shown to have linear convergence rate which is very poor compared to, for example, Simpson’s method which has a quartic convergence rate. Even though there are validated methods for computing in- tegrals that are much faster we will actually stick with this one in the end. The main reason for this is that it is easily adapted to work for integrals in more than one dimension and does not require that the function is differentiable.
3.4.2 A Validated Newton’s method
We have already given a short presentation of Newton’s method in section 2.2,
it gives us a way of finding approximations of the zeros of a function. Now we
will introduce a variant of Newton’s method that uses intervals instead of real
numbers. It will give us a way of finding guaranteed enclosures of zeros and
even proving that a interval contains exactly one zero of the function.
Definition 2 (The Interval Newton Operator). For a continuous and dif- ferentiable function f : x → R with the interval extension F we define the interval Newton method
N (x) = m − f (m)
F
0(x) . (19)
Where m denotes the midpoint of the interval x.
Much like for the regular Newton method we can define a sequence, only now it will be a sequence of intervals. Starting with an interval x
0for which N (x
0) is well-defined, i.e. 0 6∈ F
0(x
0), we define the sequence
x
k+1= N (x
k) ∩ x
k, k = 0, 1, 2, .... (20) This sequence has several nice properties that will allow us to find zeros to the function.
Theorem 4. Assume that N (x
0) is well-defined. If x
0contains a zero, x
∗, of f then so does x
kfor all k ∈ N. The sequence {x
k}
∞k=1forms a nested sequence converging to x
∗.
For a proof of this we refer to [10]. So like Newton’s method this sequence allows us to get better and better approximations of the zero, only this time we get a guaranteed enclosure of the zero. This is however not the only usable property of the sequence. What we will use it for the most is its ability to prove existence of a unique zero in an interval.
Theorem 5. Assume that f is twice continuous differentiable, and that N (x) is well-defined. Then the following statements hold
1. if N (x) ∩ x = ∅, then x contains no zeros of f ; 2. if N (x) ⊆ x, then x contains exactly one zero of f .
The first statement clearly follows from Theorem 3. For a proof of the second statement we again refer to [10].
One thing to note here is that the method can never prove the existence of a root of higher multiplicity. This is simply because by definition the derivative at these roots is equal to zero and the method requires that 0 6∈ F
0(x). Thus we will only be able to locate simple zeros to the function.
We at least have a method for proving the existence of a unique simple zero in an interval and in the algorithm we will use an higher order version of this for validating the zeros.
3.5 The Step to C
2We will here give a short description of the step from handling real intervals to working with intervals on C
2. For the extension to C we will represent a complex interval with two real intervals as
z = x + iy = {x + iy|x ∈ x, y ∈ y}, (21)
so an interval in C will be a rectangle. When taking the step from C to C
2we no longer talk about intervals but instead about boxes, a box will be represented by a product of intervals. So the boxes in C
2will be given by
z
1×z
2= x
1+iy
1×x
2+iy
2= {(x
1+iy
1, x
2+iy
2)|x
1∈ x
1, y
1∈ y
1, x
2∈ x
2, y
2∈ y
2}.
(22) Both the methods described above extends naturally to these kind of boxes.
For calculating the integral all we have to change is that instead of using the width in equation 16 we use the volume of the box. In the case of the validated Newton’s method we substitute the derivative with the Jacobian in equation 19 and get
N (z) = m − J
F(z)
−1f (m), (23) here z is a box in C
2and m denotes the midpoint of it.
4 The Algorithm
We are now ready to go through the algorithm for finding zeros to a system of two analytic functions.
For a domain D ⊂ C
2let f = (f
1, f
2) : D → C
2be an analytic function with no zeros on the boundary ∂D. The methods will require that we can find an enclosure to the range of both f and its first order partial derivatives
∂f1
∂z1
,
∂f∂z12
,
∂f∂z21
,
∂f∂z22
. Because of this we will also require that we can find an interval extension to all of these. In practice this will be automated by the use of different libraries.
Since we will make use of intervals we will also require that the domain D is a box in C
2, i.e. given by a product of intervals. So it will be on the form
D = z
1× z
2, (24)
where z
1= x
1+ iy
1and z
2= x
2+ iy
2with x
1, y
1, x
2, y
2∈ IR.
4.1 Integration
We want to compute the integral given in equation 12, i
2π
2Z
∂D[2]
det J
fdet J
[2]|f |
4dx
1dy
1dx
2+ i Z
∂D[2]
det J
fdet J
[2]|f |
4dx
1dy
1dy
2−
− Z
∂D[1]
det J
fdet J
[1]|f |
4dx
2dy
2dx
1+ i Z
∂D[1]
det J
fdet J
[1]|f |
4dx
2dy
2dy
1.
With D being a box in C
2the boundary will be given by 8 cuboids, D
[1]and D
[2]will represent 4 cuboids each. These will be given by
D
1= x
1+ iy
1× x
2+ iy
2, D
2= x
1+ iy
1× x
2+ iy
2, D
3= x
1+ iy
1
× x
2+ iy
2, D
4= x
1+ iy
1× x
2+ iy
2, D
5= x
1+ iy
1× x
2+ iy
2, D
6= x
1+ iy
1× x
2+ iy
2, D
7= x
1+ iy
1× x
2+ iy
2
