Finding roots of univariate polynomials is one of the fundamental tasks of numerics, and there is still a wide gap between root finders that are well understood in theory and those that perform well in practice. We investigate the root-finding method of Weierstrass, also known as the Durand–Kerner-method: this is a root finder that tries to approximate all roots of a given polynomial in parallel. This method has been introduced 130 years ago and has since then a good reputation for finding all roots in practice except in obvious cases of symmetry. Nonetheless, very little is known about its global dynamics and convergence properties. We show that the Weierstrass method, like the well-known Newton method, is not generally convergent: there are open sets of polynomials p p of every degree d ≥ 3 d \ge 3 such that the dynamics of the Weierstrass method applied to p p exhibits attracting periodic orbits. Specifically, all polynomials sufficiently close to Z 3 + Z + 175 Z^3 + Z + 175 have attracting cycles of period 4 4 . Here, period 4 4 is minimal: we show that for cubic polynomials, there are no periodic orbits of length 2 2 or 3 3 that attract open sets of starting points. We also establish another convergence problem for the Weierstrass method: for almost every polynomial of degree d ≥ 3 d\ge 3 there are orbits that are defined for all iterates but converge to ∞ \infty ; this is a problem that does not occur for Newton’s method. Our results are obtained by first interpreting the original problem coming from numerical mathematics in terms of higher-dimensional complex dynamics, then phrasing the question in algebraic terms in such a way that we could finally answer it by applying methods from computer algebra. The main innovation here is the translation into an algebraic question, which is amenable to (exact) computational methods close to the limits of current computer algebra systems.
more »
« less
High-Performance Polynomial Root Finding for Graphics
We present a computationally-efficient and numerically-robust algorithm for finding real roots of polynomials. It begins with determining the intervals where the given polynomial is monotonic. Then, it performs a robust variant of Newton iterations to find the real root within each interval, providing fast and guaranteed convergence and satisfying the given error bound, as permitted by the numerical precision used. For cubic polynomials, the algorithm is more accurate and faster than both the analytical solution and directly applying Newton iterations. It trivially extends to polynomials with arbitrary degrees, but it is limited to finding the real roots only and has quadratic worst-case complexity in terms of the polynomial's degree. We show that our method outperforms alternative polynomial solutions we tested up to degree 20. We also present an example rendering application with a known efficient numerical solution and show that our method provides faster, more accurate, and more robust solutions by solving polynomials of degree 10.
more »
« less
- Award ID(s):
- 1764071
- NSF-PAR ID:
- 10382084
- Date Published:
- Journal Name:
- Proceedings of the ACM on Computer Graphics and Interactive Techniques
- Volume:
- 5
- Issue:
- 3
- ISSN:
- 2577-6193
- Page Range / eLocation ID:
- 1 to 15
- Format(s):
- Medium: X
- Sponsoring Org:
- National Science Foundation
More Like this
-
-
Abstract This paper presents an implementation of a homotopy path tracking algorithm for polynomial numerical continuation on a graphical processing unit (GPU). The goal of this algorithm is to track homotopy curves from known roots to the unknown roots of a target polynomial system. The path tracker solves a set of ordinary differential equations to predict the next step and uses a Newton root finder to correct the prediction so the path stays on the homotopy solution curves. In order to benefit from the computational performance of a GPU, we organize the procedure so it is executed as a single instruction set, which means the path tracker has a fixed step size and the corrector has a fixed number iterations. This trade-off between accuracy and GPU computation speed is useful in numerical kinematic synthesis where a large number of solutions must be generated to find a few effective designs. In this paper, we show that our implementation of GPU-based numerical continuation yields 85 effective designs in 63 s, while an existing numerical continuation algorithm yields 455 effective designs in 2 h running on eight threads of a workstation.more » « less
-
Until recently, the only known method of finding the roots of polynomials over prime power rings, other than fields, was brute force. One reason for this is the lack of a division algorithm, obstructing the use of greatest common divisors. Fix a prime integer p and f in (Z/pnZ)[x] any nonzero polynomial of degree d whose coefficients are not all divisible by p. For the case n=2, we prove a new efficient algorithm to count the roots of f in Z/p2Z within time (d+size(f)+log p),2+o(1), based on a formula conjectured by Cheng, Gao, Rojas, and Wan.more » « less
-
more » « less
Abstract Strain localization and resulting plasticity and failure play an important role in the evolution of the lithosphere. These phenomena are commonly modeled by Stokes flows with viscoplastic rheologies. The nonlinearities of these rheologies make the numerical solution of the resulting systems challenging, and iterative methods often converge slowly or not at all. Yet accurate solutions are critical for representing the physics. Moreover, for some rheology laws, aspects of solvability are still unknown. We study a basic but representative viscoplastic rheology law. The law involves a yield stress that is independent of the dynamic pressure, referred to as von Mises yield criterion. Two commonly used variants, perfect/ideal and composite viscoplasticity, are compared. We derive both variants from energy minimization principles, and we use this perspective to argue when solutions are unique. We propose a new stress‐velocity Newton solution algorithm that treats the stress as an independent variable during the Newton linearization but requires solution only of Stokes systems that are of the usual velocity‐pressure form. To study different solution algorithms, we implement 2‐D and 3‐D finite element discretizations, and we generate Stokes problems with up to 7 orders of magnitude viscosity contrasts, in which compression or tension results in significant nonlinear localization effects. Comparing the performance of the proposed Newton method with the standard Newton method and the Picard fixed‐point method, we observe a significant reduction in the number of iterations and improved stability with respect to problem nonlinearity, mesh refinement, and the polynomial order of the discretization.
-
Suppose $F:=(f_1,\ldots,f_n)$ is a system of random $n$-variate polynomials with $f_i$ having degree $\leq\!d_i$ and the coefficient of $x^{a_1}_1\cdots x^{a_n}_n$ in $f_i$ being an independent complex Gaussian of mean $0$ and variance $\frac{d_i!}{a_1!\cdots a_n!\left(d_i-\sum^n_{j=1}a_j \right)!}$. Recent progress on Smale's 17$\thth$ Problem by Lairez --- building upon seminal work of Shub, Beltran, Pardo, B\"{u}rgisser, and Cucker --- has resulted in a deterministic algorithm that finds a single (complex) approximate root of $F$ using just $N^{O(1)}$ arithmetic operations on average, where $N\!:=\!\sum^n_{i=1}\frac{(n+d_i)!}{n!d_i!}$ ($=n(n+\max_i d_i)^{O(\min\{n,\max_i d_i)\}}$) is the maximum possible total number of monomial terms for such an $F$. However, can one go faster when the number of terms is smaller, and we restrict to real coefficient and real roots? And can one still maintain average-case polynomial-time with more general probability measures? We show the answer is yes when $F$ is instead a binomial system --- a case whose numerical solution is a key step in polyhedral homotopy algorithms for solving arbitrary polynomial systems. We give a deterministic algorithm that finds a real approximate root (or correctly decides there are none) using just $O(n^3\log^2(n\max_i d_i))$ arithmetic operations on average. Furthermore, our approach allows Gaussians with arbitrary variance. We also discuss briefly the obstructions to maintaining average-case time polynomial in $n\log \max_i d_i$ when $F$ has more terms.more » « less
- Free Publicly Accessible Full Text
-
Accepted Manuscript1.0
- Journal Article:
- https://doi.org/10.1145/3543865
