Robustness in Numerical Computation I Root Finding

1. Robustness in Numerical Computation IRoot Finding KwangheeKo School of Mechatronics Gwnagju Institute of Science and Technology

2. Motivation • Why do we need a nonlinear polynomial solver? • Many geometric problems are formulated as systems of nonlinear polynomial equations • Distance function • Intersections • Curvature extrema, etc.

3. Solution Methods • Local Solution Methods • Newton-type method • good initial approximation • No assurance that all roots have been found • Global Solution Methods • Algebraic Type Methods • Homotopy (Continuation) Methods • Subdivision Methods

4. Subdivision Methods • Advantages • Efficient and stable • Easy to implement • Disadvantages • No certainty that each root has been isolated • No explicit information about root multiplicities • Example Algorithm • Projected Polyhedron Algorithm

5. Projected Polyhedron Algorithm • Transform algebraic problems (finding roots) to geometric problems (computing intersections). • Use a powerful geometric property, the convex hullproperty in the algorithm. • Input : a system of polynomial equations in Bernstein form. • Output : intervals (regions) which contain roots. Note : Conversion from power basis to Bernstein basis is unstable. So we need to use Bernstein polynomials from the beginning to formulate a problem.

6. Projected Polyhedron Algorithm(Univariate Polynomial Case) f(t) • Make an transformation such that the range t of a function f(t) is from 0 to 1. • Construct a graph (t,f(t)) 0 t1 t2 1 t

7. Projected Polyhedron Algorithm (Univariate Polynomial Case) Construct the convex hull of the Bezier curve

8. Projected Polyhedron Algorithm (Univariate Polynomial Case) Intersect the convex hull with the parameter axis

9. Projected Polyhedron Algorithm (Univariate Polynomial Case) Discard the regions which do not contain roots after applying the de Casteljau subdivision algorithm Eliminate

10. Projected Polyhedron Algorithm (Multivariate Polynomial Case) • x = (x1,x2,…,xm) independent variables • f1(x)=f2(x)=…=fn(x)=0 • Algorithm • Affine Transformation : 0  xi  1, i=1,…,m • Construct a graph for each fj, j=1,…,n. • Project the control polygon of each fj onto m different coordinate (2D) plane, i.e. (x1,xm+1)-plane, …, (xm,xm+1)-plane. • Construct n two dimensional convex hulls in each plane. • Intersect each convex hull with the horizontal axis : x1,…,xm. • Compute intersection of all the intervals obtained at step 5. • If empty, no solution. If non-empty, discard the portions that do not contain roots. • Repeat.

11. Robustness Issues • Floating Point Arithmetic • Limited precision / errors due to rounding : possibility of missing roots

12. Robustness Issues How to Guarantee Robustness? • Rational / Exact Arithmetic • Memory intensive and time consuming • Interval Arithmetic • Inexpensive and robust

13. Robustness Issues Results in Rounded Interval Arithmetic

14. Polynomials with Multiple Roots • Loss of accuracy as a result of high multiplicity m. y = (x-0.1)m

15. Computing Multiplicity using Topological Degree Method • Let p(x,y), q(x,y) be polynomials with rational coefficients without common factors, of degrees n1 and n2, and let F=(p, q). • Let A be a rectangle in the plane defined by so that no zero of F lies its boundary , and does not vanish at its vertices. • Gauss map where S1 is the unit circle. • G is continuous ( ). • and S1 carry the counterclockwise orientation • Degree d of G : an integer indicating how many times is wrapped around S1 by G.

16. Illustration of the Gauss Map q(x,y) F=(p(x,y),q(x,y)) F1 y (0,0) p(x,y) a4 z=(x0,y0) a3 Y x (0,0) a1 a2 F1 / | F1| G=F/||F|| Gauss Map X S1

17. Example • p(x) = (x-1/2)5 = 0 ->

18. Computing Multiplicity using Topological Degree Method • Two approaches for computing the topological degree • Cauchy index computation method • Winding number computation method

19. The Cauchy Index • Preliminaries • R(x) : a rational function q(x)/p(x), where p, q are polynomials. • [a,b] : a closed interval, a < b. R does not become infinite at the end points. • Definition of the Cauchy index • By the Cauchy index, of R over [a,b], we mean where denotes the number of points in (a,b) at which R(x) jumps from , respectively, as x is moving from a to b. Notice that from the definition.

20. Cauchy Index (continued) • Preliminaries • A : a rectangle defined by [ a1,a2] x [ a3,a4] which encloses a zero. • F= (p,q) does not vanish on the boundary of A, • is not zero at each vertex of A. • Let Then, we set (for counterclockwise traversal of ) • Proposition* IAF is an even integer and the multiplicity * T. Sakkalis, “The Euclidean Algorithm and the Degree of the Gauss Map”, SIAM J. Computing. Vol. 19, No. 3, 1990.

21. Illustrative Example for Multiplicity Computation No. Roots of f(x,a3) = 0 in [0,1] (from the IPP) 1 [0.46922316412099, 0.46922316512099] 2 [0.49273457408967, 0.492734576204823] 3 [0.499999997363532, 0.500000001889623] 4 [0.507265424645288, 0.507265426808589] 5 [0.530776834861365, 0.530776835861365] • p(x) = (x-1/2)5 = 0 • A root of p(x), [ a ] = [0.49,0.51]. • P(z); (z = x+iy) • Create • Calculate the Cauchy index • Roots of f(x,a3) = 0 • Calculation of • Roots No. 2, 3, and 4 are selected since they lie within the interval [ a ].

22. Illustrative Example (continued) • Similarly, • Calculate • The multiplicity m of the root is

23. Winding number Computation Method

24. Polynomials with Multiple Roots Intersection between f(x) and y=0 Interval Projected Polyhedron Algorithm After verifying the existence of roots using the degree of Gauss map

25. Univariate Polynomial Equation • Assume that we have a univariate polynomial equation f(t) = 0. • Substitute t with z = x + iy. Then we have f(z) = p(x,y) + iq(x,y). • Solve f(z)=0 in the complex domain based on the topological degree algorithm. • F(x,y) = (p(x,y),q(x,y)). • We can find not only real but also complex roots of f(t)=0 at the same time using the bisection algorithm.

26. Topological Degree Bisection Algorithm y x

27. Topological Degree Bisection Algorithm y x

28. Topological Degree Bisection Algorithm y x

29. Topological Degree Bisection Algorithm y x

30. Example

31. Bivariate Polynomial Equations • Topological Degree Method • Useful for a polynomial equation of one variable • Procedure for bivariate polynomial equations • Compute the resultants of the system for each variable. • Solve the resultant equations using the topological degree method. • Select tuples which satisfy the system of equation.

32. Bivariate Polynomial Equations • Resultant of a system of polynomial equations • Let f and g be two polynomials in one variable. • When f and g have a nontrivial common factor h, say f = f1h and g=g1h, and the equation uf+vg = 0 has the solution u=g1 and v=-f1. • Conversely, the existence of nonzero solutions to this equation implies the fact that the polynomials f and g admit a common non-trivial divisor h.

33. Bivariate Polynomial Equations • Resultant of a system of polynomial equations • Given = 0 Then there exists a common non-trivial divisor h of f and g, which is a common solution of the equations f=0 and g=0.

34. Bivariate Polynomial Equations • Resultant of a system of bivariate polynomial equations • Given f(x,y) = 0 and g(x,y) = 0. • h(x) = Resy(f,g), l(y) = Resx(f,g). • Find the solutions of h(x) = 0 and l(y) = 0 to make the resultants be zero. • It implies that both equations have the common factors, namely, common roots. • Find pairs of (x,y) to satisfy f(x,y) = g(x,y) = 0 at the same time.

35. Bivariate Polynomial Equations • Example