A subdivision IA-based algorithm to find roots of a simple polynomial based on the Krawczyk operator (follows

*Subdivision Algorithms for Complex Root Isolation: Empirical Comparisons, N. Kamath, master thesis, 2010.*).

\[

K(A) = m - Y\cdot f(m) - (I - Y\cdot J(A))\cdot(A-m)

\]

where

- f is the polynomial function to find zeros for
- A is the input complex square to judge if it contains a root
- m is an arbitrarily chosen point inside A (I use the center)
- Y is an abritrary invertible matrix - real valued or real interval-valued (can be the inverse of the Jacobian, I use {{2,-0.25},{-0.5,1.5}}
- J(A) is the Jacobian of f in its interval extension at square A

Algebraic operations involving different types (A-m e.g.) are enhanced appropriately (m is treated as a degenerate complex interval with zero width in its components).

The algorithm starts with a list of complex squares

*A* (here the entire 8-square) and goes over it repeatedly, computing the operator

*K(A)*, then using one test:

*exclusion predicate* If \( K(A) \cap A = \emptyset \) => *A* definitely does not contain a root and is discarded.*inclusion* If \( K(A) \subseteq A \) => *A* definitely contains a root and is subdivided into 2x2 squares put into the list.- else nothing can be deduced and
*A* is subdivided as above.

Subdivision occurs for a specified number of refinement steps and the union of all nondiscarded squares after the last refinement cover all roots (provided one starts with a sufficiently large box). Usually there are more boxes left than there are roots.

The image below is a first example for

\( z^3+iz+(-2+2.5i) \) for level 50 refinement (for the output list squares, not the image) - and it resulted in only a dozen possible root regions. White regions do not contain a root, the circumferencing red rectangle indicates the size at which those regions could be discarded. Their area decrease leads visually the way to the possible root positions - single gray pixels in a white desert. The right image is a point-sampled classic Newton control, using a modified universal set of starting points

forum link.

It was nice to see that the subdivision algorithm also judges the two bottom roots as slightly not on an axis parallel line. I like those subtle features as a plausibility test.

All calculations were done error-free using my custInt library and outward rounding with double precision at interval end points. The image itself is however only a visualization and not error-controlled (pixels can therefore be off by +-1 in position).

The exclusion predicate works quite well, but so far I haven't had a positive inclusion result. I am currently exploring different

*Y* matrices and try to find out whether the outward rounding (computing K(A) involves a lot of operations) always enlarges K(A) to such an extent that it hardly can get smaller than A.

Has anyone experience with subdivision/IA and especially inclusion tests? For my purposes, speed is not a concern, neither is finding all roots, but my subsequent analysis relies on the information that a region in the complex plane definitely contains a root.

Linkback: https://fractalforums.org/fractal-mathematics-and-new-theories/28/root-finding-via-subdivisionia/3658/