NBThis 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.
Introduction
Since Niven proved in his pioneering work [1] that every nonconstant polynomial of the form
 |
(1) |
has at least one zero in 
, 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. [3–8]).
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 [9–11].
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 Hamilton’s 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
 |
(2) |
where 
is the norm of 
and 
is the quaternion
 |
(3) |
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 
and functions 
that can be written in the form
 |
(4) |
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
be a function of the form (4), 
and 
, with 
. Such a function 
is called radially holomorphic (or radially regular) in 
if
at 
and is denoted by 
.Theorem 1
of the form (
. 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
 |
(5) |
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
 |
(6) |
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):
 |
(7) |
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:
 |
(8) |
 |
(9) |
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)
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 |
(10) |
 |
(11) |
.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.
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.
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.
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 
.
leads, in both versions, to convergence to the root 
.
, the right version of the Newton method converges to the root 
, while the left version converges to 
.
, 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:
 |
(12) |
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 variant—the so-called single-step, serial or sequential mode—that makes use of the most recent updated approximations to the roots as soon as they are available:
 |
(13) |
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 3—Factorization into linear terms
of degree 
in 
admits a factorization into linear factors; that is, there exist 
such that |
(14) |
Theorem 4—Zeros from factors
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: |
(15) |
 |
(16) |
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])
be a polynomial in 
of degree 
with simple roots and, for 
, let |
(17) |
 |
(18) |
 |
(19) |
 |
(20) |
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 |
(21) |
of 
.The functions 
, 
and 
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.
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
 |
(22) |
 |
(23) |
included in 
can be used.
produces the following results.
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.
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 
.
Conclusion
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.
Acknowledgments
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/MEC’s (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.
References
| [1] | I. Niven, “Equations in Quaternions,” The American Mathematical Monthly, 48(10), 1941 pp. 654–661. 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. 1–24. doi.org/10.1007/s10444-016-9473-0. |
| [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. 121–140. 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. 1229–1237. 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. 297–313. doi.org/10.13001/1081-3810.1240. |
| [7] | D. Janovská and G. Opfer, “Computing Quaternionic Roots by Newton’s Method,” Electronic Transactions on Numerical Analysis, 26, 2007 pp. 82–102. |
| [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. 244–256. 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. 458–470. 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. 146–161. 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. 1–15. 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. 200–214. 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(3–4) 2013 pp. 302–322. 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 d’un 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. 136–139. |
| [19] | I. O. Kerner, “Ein Gesamtschrittverfahren zur Berechnung der Nullstellen von Polynomen,” Numerische Mathematik, 8(3), 1966 pp. 290–294. doi.org/10.1007/BF02162564. |
| [20] | S. B. Prešić, “Un procédé itératif pour la factorisation des polynômes,” Comptes Rendus de l’Académie des Sciences Paris Série A, 262, 1966 pp. 862–863. |
| 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
.Available at: w3.math.uminho.pt/QuaternionAnalysis
.Available at: content.wolfram.com/sites/19/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
Portugal
mif@math.uminho.pt
Fernando Miranda
CMAT – Centre of Mathematics
DMA – Department of Mathematics and Applications
University of Minho
Campus de Gualtar, 4710-057 Braga
Portugal
fmiranda@math.uminho.pt
Ricardo Severino
DMA – Department of Mathematics and Applications
University of Minho
Campus de Gualtar, 4710-057 Braga
Portugal
ricardo@math.uminho.pt
M. Joana Soares
NIPE – Economics Politics Research Unit
DMA – Department of Mathematics and Applications
University of Minho
Campus de Gualtar, 4710-057 Braga
Portugal
jsoares@math.uminho.pt