M. Irene Falcão, Fernando Miranda, Ricardo Severino, M. Joana Soares

Part II: Root-Finding Methods


This article explores the numerical mathematics and visualization capabilities of Mathematica in the framework of quaternion algebra. In this context, we discuss computational aspects of the recently introduced Newton and Weierstrass methods for finding the roots of a quaternionic polynomial.


Since Niven proved in his pioneering work [1] that every nonconstant polynomial of the form


has at least one zero in H, thereby extending the fundamental theorem of algebra to quaternionic polynomials, the use of such polynomials has been considered by different authors and in different contexts. Quaternionic polynomials ([2]) have found a wealth of applications in a number of different areas and have motivated the design of efficient methods for numerically approximating their zeros (see e.g. [38]).

This article discusses two numerical methods to approximate the zeros (or roots) of polynomials of the form (1). They can be seen as the quaternionic versions of the well-known Newton and Weierstrass iterative root-finding methods and they both rely on quaternion arithmetic. Here we explain in detail how we have used Mathematica to produce the numerical results recently presented in [911].

All the computations in this article require the package , available for download at w3.math.uminho.pt/QuaternionAnalysis (see [12] and [13]).

Newton-Like Methods

Theoretical Framework

We introduce the basic definitions and results needed; we refer to Part 1 of this article [2] for recalling the main aspects of the quaternion algebra and to [14] for details on quaternionic calculus.

The real vector space can be identified with by means of

where , and are Hamiltons imaginary units. Thus, throughout the article, we do not distinguish an element in from the corresponding quaternion in , unless we need to stress the context.

Using the simplified notation for the vector part of , any arbitrary nonreal quaternion can be written as


where is the norm of and is the quaternion


also referred to as the sign of . In addition, since and , one can say that behaves like the complex imaginary unit, and for this reason we call (2) the complex-like form of the quaternion .

In what follows, we consider domains Omega subset R^4=H and functions f:Omega→H that can be written in the form


where , and and are real-valued functions. Continuity and differentiability are defined coordinate-wise.

We define on the set the so-called radial operators

where and .

We introduce the following concept.

Definition 1

Let be a function of the form (4), and , with . Such a function is called radially holomorphic (or radially regular) in if

exists. In that case, this limit is called the radial derivative of at and is denoted by .

Theorem 1

A function of the form (4) is radially holomorphic iff . In that case, we have .

It follows at once that any quaternionic polynomial of the form (1) but with is radially holomorphic and its radial derivative is


Quaternionic Newton Method

For holomorphic complex functions of one complex variable, the well-known Newton method for finding a zero consists of approximating by means of the iterative process


with sufficiently close to and . Identifying a real quaternion with a vector in , the problem of solving any quaternionic equation can always be transformed into the problem of solving a system of four nonlinear equations, whose solutions, in turn, can be obtained by using the multivariate version of (6):


with sufficiently close to and a nonsingular Jacobian matrix . Not surprisingly, recent experiments performed by some of the authors of this article ([9], [10]) have shown the substantial gain in computational effort that can be achieved when using a direct quaternionic approach to this problem.

Newton methods in the quaternion context were formally adapted for the first time by Janovská and Opfer in [7], where the authors solved equations of the form . Later, Kalantari in [15], using algebraic-combinatorial arguments, proposed a Newton method for finding roots of special quaternionic polynomials. In [9], the equivalence between the classical multivariate Newton method (7) and quaternionic versions of Newton methods for a class of functions was established.

Due to the noncommutativity of multiplication for quaternions, the quotient of two quaternions and may be interpreted in two different ways: either as (the right quotient) or (the left quotient). This leads naturally to considering two versions of Newton iteration in the quaternionic setting:


The derivative in equations (8) and (9) has been considered in [9] and [10] as the radial derivative of a radially holomorphic function. In fact, in Corollary 2 of [9] it was proved that for such functions, equations (7), (8) and (9) produce, for each , the same sequence, provided that is nonsingular. Here is a more general result.

Theorem 2 ([9], Theorem 4)

Let be a function defined on the set such that the , , are radially holomorphic functions in and the are quaternions not all zero. If is a root of such that is nonsingular and is Lipschitz continuous on a neighborhood of , then for all sufficiently close to such that commutes with all , the Newton processes
both produce the same sequence as (7), which converges quadratically to .

Each step of the iterative schemes (10) and (11) is implemented in the function , which has as arguments the quaternion and the indication of the version: for (10) or for (11). At each step, a test of the value of is also performed. We recall again that all the functions presented here require the package .

The -Newton methods consist of the successive application of the iterative schemes (9) or (10) through the function , using a stopping criteria based on the incremental size and on the maximum number of iterations .

Example 1

Consider the radially holomorphic polynomial , whose only roots in are the real isolated roots , and . For the concepts of isolated and spherical roots, we refer the reader to [2], Definition 4.

The use of the initial guess requires nine iterations to get an approximation to the root 0 with precision The fact that both methods produce the same sequence is also confirmed.

The use of the initial guesses and requires 14 iterations to get an approximation to the roots and , respectively.

Example 2

The polynomial has a real root 0 and the sphere of zeros . Since the polynomial is radially holomorphic, both methods produce the same sequence. Here we would like to call attention to the convergence to the spherical root.

As pointed out in Example 3 of [10], the behavior of the Newton methods in case of convergence to values generating a spherical root is clear: if is the initial guess, then the Newton sequence converges to the root such that . This phenomenon can be easily seen from the preceding results or by computing the sign (3) of the vector part of the iterations.

When the initial guess is chosen in , then all the subsequent iterations belong to .

Example 3

Now consider the polynomial with the three isolated roots , and (cf. [9], Example 3). This polynomial is not radially holomorphic, which means that we cannot anticipate the behavior of Newton methods unless we choose initial guesses such that Theorem 2 applies, that is, such that commutes with . In other words, must be of the form .

What happens if the assumptions of Theorem 2 are not valid? In fact, as we next illustrate, although the left and right Newton methods do not give the same sequence, we can observe convergence in both cases.

The choice leads, in both versions, to convergence to the root .

With the choice , the right version of the Newton method converges to the root , while the left version converges to .

It is interesting that the 4D Newton method (7) gives convergence to the other root , as observed in [9].

Following [9] and [10], consider a function that gives the number of iterations required for each process to converge, within a certain precision, to one of the solutions of the problem under consideration, using as the initial guess.

We now consider different initial guesses by choosing points in special regions and we show density plots of . The white regions that may appear correspond to a choice of for which the method under consideration does not reach the level of precision with iterations. The default choices of and usually lead to realistic plots that require some minutes to be produced. A smoother density can be obtained by increasing the option .

Example 4

We consider again the polynomial of Example 3, whose roots are the isolated roots , and . The following code produces the plots corresponding to the choice of in one of the following regions:

As was already pointed out, Theorem 2 can be applied only in ; this is why both methods produce the same plots in this case.

Here is the behavior of the -Newton methods in .

Here is the behavior of the -Newton methods in .

Basin of Attraction of a Root

The plots produced by give information on the number of iterations required by each of the quaternionic Newton methods to converge within a certain precision to any of the roots of the polynomial under consideration. However, those plots do not give any information about the root and how the convergence occurs. This issue can be easily overcome by plotting the basins of attraction of the roots with respect to the iterative function. More precisely, we introduce a new input parameter in the function with the information of the root for which we want to compute the basin of attraction. A new function takes into account the existence of spheres of zeros. The functions and give the number of iterations needed to observe convergence to an isolated root or a spherical one, respectively. These functions return when the corresponding convergence test fails.

The functions that plot the basin of attraction of an isolated root or a spherical root have an input parameter associated to that root. The color coding used is the following: if the initial guess , chosen in a domain , causes the process to converge to a certain isolated root to which the color was associated, then the point is plotted with the color . For a sphere of zeros , all the points that converge to a point in have the color assigned to . Dark shades of a color mean fast convergence, while lighter-colored points lead to slower convergence. As before, white regions mean that the method does not converge.

Example 5

We consider once more the polynomial of Example 4, now from the perspective of the basins of attraction of each of the roots , and . We associate with these roots the colors red, blue and green, respectively, and consider the domains , and , described in Example 4. The corresponding plots can be obtained as follows (it can take some time to produce the figures).

Here are the basins of attraction in (left).

Here are the basins of attraction in (left and right).

Here are the basins of attraction in (left and right).

Example 6

This example concerns the polynomial studied in Example 2, which has an isolated root 0 (red) and a sphere of zeros (blue). The corresponding plots can be obtained as follows.

Here are the basins of attraction in (left).

Here are the basins of attraction in (left); as expected, the behavior is similar to that in , since .

Here are the basins of attraction in (left).

Weierstrass Method

Theoretical Framework

The Weierstrass method is one of the most popular iterative methods for obtaining simultaneously approximations to all the roots of a polynomial with complex coefficients. The method was first proposed by Weierstrass [16] in 1891 and later rediscovered and derived in different ways by Durand [17] in 1960, Dočev [18] in 1962 and Kerner [19] and Prešić [20] in 1966.

Let be a complex monic polynomial of degree with roots and let be distinct numbers. The classical Weierstrass method for approximating the roots is defined by the iterative scheme:


If the roots are distinct and are sufficiently good initial approximations to these roots, then the method converges at a quadratic rate, as was first proved by Dočev [18]. The iteration procedure (12) computes one approximation at a time based on the already computed approximations. For this reason, it is usually referred to as the total-step or parallel mode. The convergence of the method can be accelerated by using a variantthe so-called single-step, serial or sequential modethat makes use of the most recent updated approximations to the roots as soon as they are available:


In a recent article [11], we adapted the Weierstrass method to the quaternion algebra setting. We refer to [2] and references therein to recall the main concepts and properties of the ring of unilateral quaternionic polynomials. In particular, we recall the factorization of polynomials in into linear terms and the relation between zeros and factors of .

Theorem 3Factorization into linear terms

Any monic polynomial of degree in admits a factorization into linear factors; that is, there exist such that

Theorem 4Zeros from factors

Consider a polynomial whose factor terms are ; that is, admits a factorization of the form (14). If the similarity classes , , are distinct, then has exactly zeros , which are given by:

Weierstrass Algorithm

Following the idea of the Weierstrass method in its sequential version (13), the next results show how to obtain sequences converging, at a quadratic rate, to the factor terms in (14) of a given polynomial . Moreover, by making use of Theorem 4, it is possible to construct sequences converging quadratically to the roots of .

Theorem 5 ([11])

Let be a polynomial in of degree with simple roots and, for , let
with denoting the characteristic polynomial of , that is, . If the initial approximations are sufficiently close to the factor terms in a factorization of in the form (14), then the sequences converge quadratically to . Moreover, the sequences defined by
converge quadratically to the roots of .

The functions , and C_i^((k)) are implemented as the functions , and , respectively. The support file associated with [2] needs to be loaded.

The iterative functions associated with (17) and (21) are built into the function.

The quaternionic Weierstrass iterative method is implemented in the function .

The usual convergence test has been replaced in by

in order to let the function recognize a sphere of zeros. Since we also include a test on the value of , there is no risk of misidentifying an isolated root.

Example 7

We consider now the application of the Weierstrass method to the computation of the roots of the polynomial of Example 3, which we recall are , and . All of the initial approximations , and have to lie in distinct congruence classes.

Some explanation of the output is needed. The first entry indicates the convergence or divergence of the method. The second entry is the error in the approximations to the zeros. The last two entries contain approximations to the roots and factors terms. Since there are two real roots and just one nonreal root, the roots and factor terms coincide.

Example 8

Our next test example is a polynomial that also fulfills the assumptions of Theorem 5 and has simple zeros (see [11], Example 1). First, we check that the polynomial

has the roots
Recall that the function included in can be used.

The Weierstrass method applied to the extended form of produces the following results.

The convergence to the roots is in a order different from the one given in (22) because the convergence to the factor terms also occurs in a sequence different from the one given in (23).

Example 9

The polynomial has an isolated root and a sphere of zeros . The assumptions of Theorem 5 do not apply to this polynomial, but we can observe convergence to the roots as we increase the precision of the computations. When a polynomial has a spherical root, two of its factor terms are in the same congruence class. Therefore, as the iteration proceeds, the values in (17) become close to zero and some care is required.

Using the usual precision, it was not possible to reach the required tolerance. However, performing the calculations with more decimal places causes a fast convergence, under the same assumptions.

The spherical root can be identified at once by observing that, up to the required precision, we have .


This is the second article on several computational aspects of polynomials in the ring . One can find in the literature methods for numerically approximating the zeros of quaternionic polynomials based on the use of complex techniques, but numerical methods relying on quaternion arithmetic remain scarce, with the exceptions of the Newton and Weierstrass methods discussed in this article. We developed several functions to implement those methods and we also added some visualization tools.


Research at the Centre of Mathematics (CMAT) was financed by Portuguese Funds through FCT – Fundação para a Ciência e a Tecnologia, within the Project UID/MAT/00013/2013. Research at the Economics Politics Research Unit (NIPE) was carried out within the funding with COMPETE reference number POCI-01-0145-FEDER-006683 (UID/ECO/03182/2013), with the FCT/MECs (Fundação para a Ciência e a Tecnologia, I.P.) financial support through national funding and by the European Regional Development Fund (ERDF) through the Operational Programme on Competitiveness and Internationalization – COMPETE 2020 under the PT2020 Partnership Agreement.


[1] I. Niven, Equations in Quaternions, The American Mathematical Monthly, 48(10), 1941
pp. 654661. www.jstor.org/stable/2303304.
[2] M. I. Falcão, F. Miranda, R. Severino, and M. J. Soares, Computational Aspects of Quaternionic Polynomials: Part 1, The Mathematica Journal, 20(4), 2018. doi.org/10.3888/tmj.20-4.
[3] R. Farouki, G. Gentili, C. Giannelli, A. Sestini and C. Stoppato, A Comprehensive Characterization of the Set of Polynomial Curves with Rational Rotation-Minimizing Frames, Advances in Computational Mathematics, 43(1), 2017 pp. 124.
[4] R. Pereira, P. Rocha and P. Vettori, Algebraic Tools for the Study of Quaternionic Behavioral Systems, Linear Algebra and Its Applications, 400, 2005 pp. 121140. doi.org/10.1016/j.laa.2005.01.008.
[5] R. Serôdio, E. Pereira and J. Vitória, Computing the Zeros of Quaternion Polynomials, Computers and Mathematics with Applications, 42(8-9) 2001 pp. 12291237. doi.org/10.1016/S0898-1221(01)00235-8.
[6] S. De Leo, G. Ducati, and V. Leonardi, Zeros of Unilateral Quaternionic Polynomials, The Electronic Journal of Linear Algebra, 15(1), 2006 pp. 297313.
[7] D. Janovská and G. Opfer, Computing Quaternionic Roots by Newton’s Method, Electronic Transactions on Numerical Analysis, 26, 2007 pp. 82102.
[8] D. Janovská and G. Opfer, A Note on the Computation of All Zeros of Simple Quaternionic Polynomials, SIAM Journal on Numerical Analysis, 48(1), 2010 pp. 244256. doi.org/10.1137/090748871.
[9] M. I. Falcão, Newton Method in the Context of Quaternion Analysis, Applied Mathematics and Computation, 236, 2014 pp. 458470. doi.org/10.1016/j.amc.2014.03.050.
[10] F. Miranda and M. I. Falcão, Modified Quaternion Newton Methods, in Computational Science and Its Applications (ICCSA 2014), Guimarães, Portugal, Lecture Notes in Computer Science, 8579 (B. Murgante et al., eds.), Berlin, Heidelberg: Springer, 2014 pp. 146161. doi.org/10.1007/978-3-319-09144-0_ 11.
[11] M. I. Falcão, F. Miranda, R. Severino and M. J. Soares, Weierstrass Method for Quaternionic Polynomial Root-Finding, Mathematical Methods in the Applied Sciences, 2017 pp. 115. doi:10.1002/mma.4623.
[12] M. I. Falcão and F. Miranda, Quaternions: A Mathematica Package for Quaternionic Analysis, in Computational Science and Its Applications (ICCSA 2011), Lecture Notes in Computer Science, 6784 (B. Murgante, O. Gervasi, A. Iglesias, D. Taniar and B. O. Apduhan, eds.), Berlin, Heidelberg: Springer, 2011 pp. 200214. doi:10.1007/978-3-642-21931-3_17.
[13] F. Miranda and M. I. Falcão. QuaternionAnalysis Mathematica Package. w3.math.uminho.pt/QuaternionAnalysis.
[14] K. Gürlebeck, K. Habetha and W. Sprössig, Holomorphic Functions in the Plane and
Dimensional Space, Basel: Birkhäuser, 2008. doi.org/10.1007/978-3-7643-8272-8.
[15] B. Kalantari, Algorithms for Quaternion Polynomial Root-Finding, Journal of Complexity, 29(34) 2013 pp. 302322. doi.org/10.1016/j.jco.2013.03.001.
[16] K. Weierstrass, Neuer Beweis des Satzes, dass jede ganze rationale Function einer Veränderlichen dargestellt werden kann als ein Product aus linearen Functionen derselben Veränderlichen, in Sitzungsberichte der Königlich Preussischen Akademie der Wissenschaften zu Berlin, 1891.
[17] E. Durand, Solutions numériques des équations algébriques. Tome I: Equations du type F(x); racines dun polynôme, Paris: Masson,1960.
[18] K. Dočev, A Variant of Newton’s Method for the Simultaneous Approximation of All Roots of an Algebraic Equation, Fiziko-Matematichesko Spisanie. Bulgarska Akademiya na Naukite, 5(38), 1962 pp. 136139.
[19] I. O. Kerner, Ein Gesamtschrittverfahren zur Berechnung der Nullstellen von Polynomen, Numerische Mathematik, 8(3), 1966 pp. 290294. doi.org/10.1007/BF02162564.
[20] S. B. Prešić, Un procédé itératif pour la factorisation des polynômes, Comptes Rendus de lAcadémie des Sciences Paris Série A, 262, 1966 pp. 862863.
M. I. Falcão, F. Miranda, R. Severino and M. J. Soares, Computational Aspects of Quaternionic Polynomials: Part 2, The Mathematica Journal, 2018. dx.doi.org/doi:10.3888/tmj.20-5.

Additional Material

  • The package .
  • Available at: w3.math.uminho.pt/QuaternionAnalysis

  • The file .
  • Available at: www.mathematica-journal.com/data/uploads/2018/05/QPolynomial.m

    About the Authors

    M. Irene Falcão is an associate professor in the Department of Mathematics and Applications of the University of Minho. Her research interests are numerical analysis, hypercomplex analysis and scientific software.

    Fernando Miranda is an assistant professor in the Department of Mathematics and Applications of the University of Minho. His research interests are differential equations, quaternions and related algebras and scientific software.

    Ricardo Severino is an assistant professor in the Department of Mathematics and Applications of the University of Minho. His research interests are dynamical systems, quaternions and related algebras and scientific software.

    M. Joana Soares is an associate professor in the Department of Mathematics and Applications of the University of Minho. Her research interests are numerical analysis, wavelets mainly in applications to economics, and quaternions and related algebras.

    M. Irene Falcão
    CMAT – Centre of Mathematics
    DMA – Department of Mathematics and Applications
    University of Minho
    Campus de Gualtar, 4710-057 Braga


    Fernando Miranda
    CMAT – Centre of Mathematics
    DMA – Department of Mathematics and Applications
    University of Minho
    Campus de Gualtar, 4710-057 Braga



    Ricardo Severino
    DMA – Department of Mathematics and Applications
    University of Minho
    Campus de Gualtar, 4710-057 Braga



    M. Joana Soares
    NIPE – Economics Politics Research Unit
    DMA – Department of Mathematics and Applications
    University of Minho
    Campus de Gualtar, 4710-057 Braga