Weierstrass method for quaternionic polynomial root-finding

TL;DR

This paper introduces a quaternionic Weierstrass-like method for efficiently finding all zeros of unilateral quaternionic polynomials, supported by convergence analysis and numerical experiments.

Contribution

It presents the first quaternionic Weierstrass method with convergence analysis and numerical validation for polynomial root-finding.

Findings

The method converges reliably for quaternionic polynomials.

Numerical examples demonstrate the method's effectiveness.

The approach extends complex root-finding techniques to quaternions.

Abstract

Quaternions, introduced by Hamilton in 1843 as a generalization of complex numbers, have found, in more recent years, a wealth of applications in a number of different areas which motivated the design of efficient methods for numerically approximating the zeros of quaternionic polynomials. In fact, one can find in the literature recent contributions to this subject based on the use of complex techniques, but numerical methods relying on quaternion arithmetic remain scarce. In this paper we propose a Weierstrass-like method for finding simultaneously {\sl all} the zeros of unilateral quaternionic polynomials. The convergence analysis and several numerical examples illustrating the performance of the method are also presented.

Tables

Table 1. Table 1 : Quaternionic Weierstrass method for Example  1

$k$	${ϵ}_1^{\left(k\right)}$	${ϵ}_2^{\left(k\right)}$	${ϵ}_3^{\left(k\right)}$	${ϵ}_4^{\left(k\right)}$	${ϵ}_5^{\left(k\right)}$	${ϵ}_6^{\left(k\right)}$	${\rho }^{\left(k\right)}$
$0$	$2.7\text{e}-1$	$6.0\text{e}-2$	$4.5\text{e}-1$	$2.0\text{e}-2$	$8.0\text{e}-2$	$3.3\text{e}-1$	–
$1$	$9.3\text{e}-2$	$1.8\text{e}-2$	$7.7\text{e}-2$	$7.2\text{e}-3$	$6.0\text{e}-2$	$3.1\text{e}-2$	1.36
$2$	$7.9\text{e}-3$	$1.9\text{e}-3$	$5.7\text{e}-3$	$5.6\text{e}-4$	$5.3\text{e}-3$	$1.1\text{e}-3$	2.03
$3$	$6.0\text{e}-5$	$2.0\text{e}-5$	$4.0\text{e}-5$	$4.5\text{e}-6$	$9.2\text{e}-6$	$3.1\text{e}-7$	2.17
$4$	$2.4\text{e}-9$	$1.5\text{e}-9$	$3.3\text{e}-9$	$3.7\text{e}-10$	$2.0\text{e}-9$	$2.2\text{e}-13$	2.04
$5$	$1.5\text{e}-17$	$5.3\text{e}-18$	$1.6\text{e}-17$	$6.1\text{e}-19$	$7.5\text{e}-18$	$8.1\text{e}-26$	2.06

Table 2. Table 2 : Parallel version of Weierstrass method for Example  1

$k$	${ϵ}_1^{\left(k\right)}$	${ϵ}_2^{\left(k\right)}$	${ϵ}_3^{\left(k\right)}$	${ϵ}_4^{\left(k\right)}$	${ϵ}_5^{\left(k\right)}$	${ϵ}_6^{\left(k\right)}$	${\rho }^{\left(k\right)}$
$0$	$2.7\text{e}-1$	$6.0\text{e}-2$	$4.5\text{e}-1$	$2.0\text{e}-2$	$8.0\text{e}-2$	$3.3\text{e}-1$	–
$1$	$9.3\text{e}-2$	$2.3\text{e}-2$	$8.6\text{e}-2$	$5.9\text{e}-3$	$5.3\text{e}-3$	$9.8\text{e}-2$	3.26
$2$	$1.7\text{e}-2$	$4.4\text{e}-3$	$2.8\text{e}-2$	$2.0\text{e}-3$	$2.0\text{e}-2$	$3.3\text{e}-2$	1.47
$3$	$1.0\text{e}-3$	$4.3\text{e}-4$	$3.1\text{e}-3$	$3.1\text{e}-4$	$2.9\text{e}-3$	$4.2\text{e}-3$	1.60
$4$	$5.0\text{e}-5$	$1.9\text{e}-5$	$1.5\text{e}-4$	$1.2\text{e}-5$	$2.4\text{e}-4$	$3.0\text{e}-4$	1.48
$5$	$2.2\text{e}-7$	$1.1\text{e}-7$	$2.1\text{e}-6$	$1.0\text{e}-7$	$7.0\text{e}-6$	$7.9\text{e}-6$	1.44
$6$	$4.5\text{e}-10$	$2.2\text{e}-10$	$5.6\text{e}-8$	$2.5\text{e}-10$	$6.0\text{e}-7$	$1.9\text{e}-6$	1.12
$7$	$2.5\text{e}-13$	$1.2\text{e}-13$	$9.4\text{e}-11$	$2.0\text{e}-13$	$1.2\text{e}-8$	$8.6\text{e}-9$	1.39
$8$	$1.6\text{e}-18$	$1.6\text{e}-18$	$1.6\text{e}-15$	$3.7\text{e}-18$	$2.4\text{e}-12$	$3.0\text{e}-12$	1.45
$9$	$2.0\text{e}-25$	$1.7\text{e}-25$	$7.9\text{e}-21$	$5.5\text{e}-25$	$1.5\text{e}-17$	$2.3\text{e}-17$	1.44

Table 3. Table 3 : Weierstrass method for spherical roots - Example  2

$k$	${ϵ}_1^{\left(k\right)}$	${ϵ}_2^{\left(k\right)}$	${ϵ}_3^{\left(k\right)}$	${ϵ}_4^{\left(k\right)}$	${\rho }^{\left(k\right)}$
$0$	$1.3\text{e}-2$	$7.1\text{e}-2$	$7.6\text{e}-2$	$1.3\text{e}-1$	–
$1$	$6.3\text{e}-3$	$4.5\text{e}-2$	$6.1\text{e}-3$	$1.1\text{e}-2$	1.52
$2$	$1.2\text{e}-4$	$9.8\text{e}-4$	$1.4\text{e}-3$	$9.7\text{e}-5$	2.11
$3$	$9.6\text{e}-8$	$1.1\text{e}-6$	$2.6\text{e}-6$	$1.8\text{e}-8$	1.96
$4$	$9.2\text{e}-12$	$6.1\text{e}-11$	$1.6\text{e}-11$	$1.0\text{e}-15$	1.83
$5$	$1.4\text{e}-22$	$9.7\text{e}-22$	$4.9\text{e}-21$	$8.0\text{e}-31$	1.99

Table 4. Table 4 : Weierstrass method for double roots - Example  3

$P$		$Q$
${ϵ}^{\left(k\right)}$	${\rho }^{\left(k\right)}$	${ϵ}^{\left(k\right)}$	${\rho }^{\left(k\right)}$
$6.8\text{e}-10$	1.05	$1.9\text{e}-13$	1.03
$2.6\text{e}-10$	1.04	$7.3\text{e}-14$	1.03
$1.0\text{e}-10$	1.04	$2.8\text{e}-14$	1.03
$3.9\text{e}-11$	1.04	$1.1\text{e}-14$	1.03
$1.5\text{e}-11$	1.04	$4.1\text{e}-15$	1.03
$5.9\text{e}-12$	1.04	$1.6\text{e}-15$	1.03
$2.3\text{e}-12$	1.04	$6.1\text{e}-16$	1.03
$8.8\text{e}-13$	1.04	$2.4\text{e}-16$	1.03
$3.4\text{e}-13$	1.03	$9.1\text{e}-17$	1.03

Equations

$${\mathbf{i}}^2={\mathbf{j}}^2={\mathbf{k}}^2=-1,\mathbf{ij}=-\mathbf{ji}=\mathbf{k}.$$

$$[q]=\left\{q^{\prime }\in \mathbb{H}:\mathrm{Re}q=\mathrm{Re}{q}^{\prime }\text{ and }|q|=|q^{\prime }|\right\}.$$

$$P(x)=a_n{x}^n+a_{n-1}{x}^{n-1}+\cdots +a_{1}x+a_0,a_i\in \mathbb{H};i=0,\dots ,n.$$

$$P(x)+Q(x):=$$

$$P(x)\ast Q(x):=$$

$$P(q)=a_n{q}^n+a_{n-1}{q}^{n-1}+\cdots +a_{1}q+a_0,$$

$$P(q)=L(\tilde{q})R(q)\text{with}\tilde{q}=hq{h}^{-1}.$$

$$P(q)=R(q)L(q).$$

$$P(x)=Q(x)\ast (x-q).$$

$$P(x)=(x-x_n)\ast (x-x_{n-1})\ast \cdots \ast (x-x_1).$$

$$P(x)=(x-y_n)\ast (x-y_{n-1})\ast \cdots \ast (x-y_1)$$

$$y_{\pi (i)}=h_i{x}_i{h}_i^{-1}.$$

$$\overline{P}(x)=\sum _{k=0}^n\overline{a_k}{x}^k.$$

$$\overline{P\ast Q}$$

$$P\ast \overline{P}$$

$${\mathcal{Q}}_q(x):=(x-q)\ast (x-\overline{q})=x^2-2\mathrm{Re}qx+|q{|}^2,$$

$$z_i^{(k+1)}=z_i^{(k)}-\frac{P(z_i^{(k)})}{\prod _{\genfrac{}{}{0ex}{}{j=1}{j\ne i}}^n(z_i^{(k)}-z_j^{(k)})};i=1,\dots ,n;k=0,1,2,\dots$$

$$z_i^{(k+1)}=z_i^{(k)}-\frac{P(z_i^{(k)})}{\prod _{j=1}^{i-1}(z_i^{(k)}-z_j^{(k+1)})\prod _{j=i+1}^n(z_i^{(k)}-z_j^{(k)})};i=1,\dots ,n;k=0,1,2,\dots$$

$$P(x)=(x-x_n)\ast (x-x_{n-1})\ast \cdots \ast (x-x_1),$$

$$\mathop{\mathrel{\mathchoice{\ooalign{\ooalign{$\displaystyle\prod$\cr\hfil$\displaystyle\ast$\hfil\cr}\cr\hfil$\displaystyle\phantom{\prod}$\hfil\cr}{}}{\ooalign{\ooalign{$\textstyle\prod$\cr\hfil$\textstyle\ast$\hfil\cr}\cr\hfil$\textstyle\phantom{\prod}$\hfil\cr}{}}{\ooalign{\ooalign{$\scriptstyle\prod$\cr\hfil$\scriptstyle\ast$\hfil\cr}\cr\hfil$\scriptstyle\phantom{\prod}$\hfil\cr}{}}{\ooalign{\ooalign{$\scriptscriptstyle\prod$\cr\hfil$\scriptscriptstyle\ast$\hfil\cr}\cr\hfil$\scriptscriptstyle\phantom{\prod}$\hfil\cr}{}}}}_{i=k}^{m}(x-\alpha_{i}):=(x-\alpha_{m})\ast(x-\alpha_{m-1})\ast\cdots\ast(x-\alpha_{k}).$$

$$\mathop{\mathrel{\mathchoice{\ooalign{\ooalign{$\displaystyle\prod$\cr\hfil$\displaystyle\ast$\hfil\cr}\cr\hfil$\displaystyle\phantom{\prod}$\hfil\cr}{}}{\ooalign{\ooalign{$\textstyle\prod$\cr\hfil$\textstyle\ast$\hfil\cr}\cr\hfil$\textstyle\phantom{\prod}$\hfil\cr}{}}{\ooalign{\ooalign{$\scriptstyle\prod$\cr\hfil$\scriptstyle\ast$\hfil\cr}\cr\hfil$\scriptstyle\phantom{\prod}$\hfil\cr}{}}{\ooalign{\ooalign{$\scriptscriptstyle\prod$\cr\hfil$\scriptscriptstyle\ast$\hfil\cr}\cr\hfil$\scriptscriptstyle\phantom{\prod}$\hfil\cr}{}}}}_{i=k}^{m}(x-\alpha_{i}):=1,\quad{\text{whenever\ }}\ k>m.$$

$$\zeta_{i}=\overline{R_{i}}(x_{i})\,x_{i}\,\bigl{(}\overline{R_{i}}(x_{i})\bigr{)}^{-1};\ i=1,2,\ldots,n,$$

$$R_{i}:=\mathop{\mathrel{\mathchoice{\ooalign{\ooalign{$\displaystyle\prod$\cr\hfil$\displaystyle\ast$\hfil\cr}\cr\hfil$\displaystyle\phantom{\prod}$\hfil\cr}{}}{\ooalign{\ooalign{$\textstyle\prod$\cr\hfil$\textstyle\ast$\hfil\cr}\cr\hfil$\textstyle\phantom{\prod}$\hfil\cr}{}}{\ooalign{\ooalign{$\scriptstyle\prod$\cr\hfil$\scriptstyle\ast$\hfil\cr}\cr\hfil$\scriptstyle\phantom{\prod}$\hfil\cr}{}}{\ooalign{\ooalign{$\scriptscriptstyle\prod$\cr\hfil$\scriptscriptstyle\ast$\hfil\cr}\cr\hfil$\scriptscriptstyle\phantom{\prod}$\hfil\cr}{}}}}_{j=1}^{i-1}(x-x_{j}).$$

$$z_i^{(k+1)}=z_i^{(k)}-\left(\overline{{\mathcal{L}}_i^{(k)}}\ast P\ast \overline{{\mathcal{R}}_i^{(k)}}\right)(z_i^{(k)}){\left({\mathcal{Q}}_i^{(k)}(z_i^{(k)})\right)}^{-1},$$

$${\mathcal{L}}_i^{(k)}(x):=$$

$${\mathcal{R}}_i^{(k)}(x):=$$

$${\mathcal{Q}}_i^{(k)}(x):=$$

$${\epsilon }_i^{(k)}:=x_i-z_i^{(k)},i=1,\dots ,n,$$

$${\epsilon }^{(k)}:=\underset{i}{\max }|{\epsilon }_i^{(k)}|.$$

$$\varepsilon_{i}^{(k+1)}={\mathcal{O}}\bigl{(}(\varepsilon^{(k)})^{2}\bigr{)}.$$

$$P(x)$$

$$\displaystyle=\Bigl{(}\mathop{\mathrel{\mathchoice{\ooalign{\ooalign{$\displaystyle\prod$\cr\hfil$\displaystyle\ast$\hfil\cr}\cr\hfil$\displaystyle\phantom{\prod}$\hfil\cr}{}}{\ooalign{\ooalign{$\textstyle\prod$\cr\hfil$\textstyle\ast$\hfil\cr}\cr\hfil$\textstyle\phantom{\prod}$\hfil\cr}{}}{\ooalign{\ooalign{$\scriptstyle\prod$\cr\hfil$\scriptstyle\ast$\hfil\cr}\cr\hfil$\scriptstyle\phantom{\prod}$\hfil\cr}{}}{\ooalign{\ooalign{$\scriptscriptstyle\prod$\cr\hfil$\scriptscriptstyle\ast$\hfil\cr}\cr\hfil$\scriptscriptstyle\phantom{\prod}$\hfil\cr}{}}}}_{j=2}^{n}(x-z_{j})+\mathscr{E}_{1}(x)\Big{)}\ast(x-z_{1}-\varepsilon_{1}),$$

$$-(x-z_n)\ast \cdots \ast (x-z_{j-1})\ast {\epsilon }_j\ast (x-z_{j+1})\ast \cdots \ast (x-z_2),$$

Full text

Keywords: Quaternionic polynomials $\cdot$ Root-finding methods $\cdot$ Weierstrass algorithm

1 Introduction

The increasing interest in using quaternions and their applications in areas as diverse as number theory, robotics, virtual reality or image processing (see e.g. [4, 10, 24, 25, 26]), motivated several authors to consider extending well-known (complex) numerical methods, in particular root-finding methods, to the quaternion algebra framework. However, the problem of finding the zeros of quaternionic polynomials turns out to be much more demanding than the analogous problem over the real and complex fields. Niven, in his pioneering work, [23], gave a first extension of the Fundamental Theorem of Algebra for the quaternion context, proving that any quaternionic polynomial of positive degree whose coefficients are located only on one side of the powers must have at least one quaternionic root. In the aforementioned paper, Niven also proposed a method for computing the roots of such polynomials. This algorithm is, however, as stated by Niven a “not very practical” one, due to the need of solving two coupled nonlinear equations for the determination of pairs of real constants. Later, in [30], the authors, by making use of (a complexified version of) the companion matrix of the polynomial, turned the ideas of Niven into what can be considered as the first really usable numerical algorithm.

Nowadays, other quaternionic root-finding algorithms are available which essentially replace the problem of computing the roots of a quaternionic polynomial of degree $n$, by the problem of determining the roots of a real or complex polynomial of degree $2n$ (usually with multiple roots), relying in this way on algorithms for complex polynomial root-finding (see [5, 31] and the references therein). Several experiments performed by two of the authors of this paper ([8, 22]) have shown the substantial gain in computational effort that can be achieved when using a direct quaternionic approach to this problem.

The Weierstrass method, also known in the literature as the Durand-Kerner method or Dochev method, is one of the most popular iterative methods for obtaining simultaneously approximations to all the roots of a given polynomial with complex coefficients (for a survey on most of the traditional methods for root-finding we refer to [20]). The formula involved in the method was first proposed by Weierstrass [33], in connection with a constructive proof of the Fundamental Theorem of Algebra, and later rediscovered and derived in different ways by Durand [7], Dochev [6] and Kerner [18], among others.

The main purpose of this paper is to present an adaptation of the Weierstrass method to the case of quaternionic polynomials. By making use of the so-called Factor Theorem for quaternions we derive an iterative method which shows fast convergence and robustness with respect to the initial approximations.

The paper is organized as follows: in Section 2 we review some basic results on the algebra of real quaternions and on quaternionic polynomials; Section 3 contains the main results of the paper; after revisiting the classical (complex) Weierstrass method we derive a generalization to the quaternionic case and prove, under some natural assumptions, its quadratic order of convergence; in Section 4 we present several numerical experiments illustrating the results obtained in Section 3; finally, in Section 5 we draw some conclusions and indicate some future work.

2 Basic results on quaternions

In this section we present a brief summary on the main results on the algebra of real quaternions and on the ring of polynomials over the quaternions needed in the sequel.

2.1 The algebra of real quaternions

Let $\{1,\mathbf{i},\mathbf{j},\mathbf{k}\}$ be an orthonormal basis of the Euclidean vector space $\mathbb{R}^{4}$ with a product given according to the multiplication rules

[TABLE]

This non-commutative product generates the well known algebra of real quaternions $\mathbb{H}$.

Given a quaternion $q=q_{0}+q_{1}\mathbf{i}+q_{2}\mathbf{j}+q_{3}\mathbf{k}\in\mathbb{H}$, its conjugate $\overline{q}$ is defined as $\overline{q}=q_{0}-q_{1}\mathbf{i}-q_{2}\mathbf{j}-q_{3}\mathbf{k}$; the number $q_{0}$ is called the real part of $q$ and denoted by $\operatorname{Re}q$ and the vector part of $q$, denoted by $\operatorname{Vec}q$, is given by $\operatorname{Vec}q=q_{1}\mathbf{i}+q_{2}\mathbf{j}+q_{3}\mathbf{k}$; the norm of $q$, $|q|$, is given by $|q|=\sqrt{q\overline{q}}=\sqrt{q_{0}^{2}+q_{1}^{2}+q_{2}^{2}+q_{3}^{2}};$ the inverse of $q$ (if $q\neq 0$), denoted by $q^{-1}$ is the (unique) quaternion such that $qq^{-1}=q^{-1}q=1$ and is given by $q^{-1}=\dfrac{\overline{q}}{|q|^{2}}$.

We say that a quaternion $q$ is congruent to a quaternion $q^{\prime}$, and write $q\sim q^{\prime}$, if there exists a non-zero quaternion $h$ such that $q^{\prime}=hqh^{-1}$. This is an equivalence relation in $\mathbb{H}$, partitioning $\mathbb{H}$ in the so-called congruence classes. We denote by $[q]$ the congruence class containing a given quaternion $q$. It can be shown (see, e.g. [35]) that

[TABLE]

It follows that $[q]$ reduces to a single element if and only if $q$ is a real number. If $q=q_{0}+q_{1}\mathbf{i}+q_{2}\mathbf{j}+q_{3}\mathbf{k}$ is not real, its congruence class can be identified with the three-dimensional sphere in the hyperplane $\{(x_{0},x_{1},x_{2},x_{3})\in\mathbb{R}^{4}:x_{0}=q_{0}\}$, with center $(q_{0},0,0,0)$ and radius $\sqrt{q_{1}^{2}+q_{2}^{2}+q_{3}^{2}}$.

2.2 Ring of left quaternionic polynomials

Because of the non-commutativity of quaternion multiplication, one can consider different classes of polynomials in one quaternion variable, depending on whether the variable commutes with the polynomial coefficients or not. General polynomials in the indeterminate $x$ are defined as finite sums of non-commutative monomials of the form $a_{0}xa_{1}\dots xa_{j}$. In this work we restrict our attention to polynomials whose coefficients are located only on the left-hand side of the powers of $x$, i.e. have the special form

[TABLE]

These polynomials are usually called in the literature one-sided or unilateral (left) polynomials. As usual, if $a_{n}\neq 0$, we will say that the degree of the polynomial $P(x)$ is $n$ and refer to $a_{n}$ as the leading coefficient of the polynomial. When $a_{n}=1$, we say that $P(x)$ is monic. If the coefficients $a_{i}$ in (2) are real, then we say that $P(x)$ is a real polynomial and write $P(x)\in\mathbb{R}[x].$

The set of polynomials of the form (2) is a ring with respect to the operations of addition and multiplication defined as in the commutative case: for any two polynomials $P(x)=\sum_{i=0}^{n}a_{i}x^{i}$ and $Q(x)=\sum_{j=0}^{m}b_{j}x^{j}$,

[TABLE]

with the implicit assumption that $a_{k}=0$ for $k>n$ and $b_{k}=0$ for $k>m$. We will denote this ring of polynomials by $\mathbb{H}[x]$. Naturally, due to the non-commutativity of the quaternionic multiplication, $\mathbb{H}[x]$ is a non-commutative ring. However, if $P(x)$ is a real polynomial, then $P(x)$ commutes with any polynomial in $\mathbb{H}[x]$.

We should also observe that the evaluation map at a given quaternion $q$, defined, for the polynomial $P(x)$ given by (2), by

[TABLE]

is not a homomorphism from the ring $\mathbb{H}[x]$ into $\mathbb{H}$. In fact, $P(x)=L(x)\ast R(x)$ does not lead, in general, to $P(q)=L(q)R(q)$.

Remark 1**.**

Since all the polynomials considered will be in the indeterminate $x$, we will usually omit the reference to this variable and write simply $P$ when referring to an element $P(x)\in\mathbb{H}[x]$, the expression $P(q)$ being preferably reserved for the evaluation of $P$ at a specific value $q\in\mathbb{H}$.

We say that a quaternion $q$ is a zero of a polynomial $P$, if $P(q)=0$, and we use the notation ${\mathbf{Z}}_{P}$ to denote the zero-set of $P$, i.e. the set of all the zeros of $P$. Since this work is concerned with the computation of zeros of polynomials, there is no loss of generality in assuming that the polynomials are monic and we will do so in what follows.

We now review some basic properties of unilateral (left) quaternion polynomials needed in the sequel.

The next theorem shows a way of evaluating the product of two polynomials at a given quaternion, without explicitly performing their product. The proof of the first two results can be seen in e.g. [19] and the last result is a simple consequence of the definition of the product of polynomials and of the fact that any real number commutes with a quaternion.

Theorem 1**.**

Let $P=L\ast R$ with $L,R\in\mathbb{H}[x]$, $q\in\mathbb{H}$ and $h=R(q)$.

1. (i)

If $h=0$, then $P(q)=0$ (i.e. if $q$ is a zero of the right factor $R$, then $q$ is also a zero of the product $P$). 2. (ii)

If $h\neq 0$, then

[TABLE]

In particular, if $q$ is a zero of $P$ which is not a zero of $R$, then $\tilde{q}$ is a zero of $L$. 3. (iii)

If $L\in\mathbb{R}[x]$, then

[TABLE]

The following result, first proved by Gordon and Motzkin [15], can also be seen in [19].

Theorem 2** (Factor Theorem).**

Let $P\in\mathbb{H}[x]$ and $q\in\mathbb{H}$. Then, $q$ is a zero of $P$ if and only if there exists $Q\in\mathbb{H}[x]$ such that

[TABLE]

In 1941, Niven [23] proved the Fundamental Theorem of Algebra for unilateral quaternionic polynomials, establishing that any non-constant polynomial in $\mathbb{H}[x]$ always has a zero in $\mathbb{H}$. More general results are contained in the following theorem.

Theorem 3**.**

Let $P$ be a monic polynomial of degree $n\,(n\geq 1)$ in $\mathbb{H}[x]$. Then:

1. (i)

$P$* admits a factorization into linear factors, i.e. there exist $x_{1},\ldots,x_{n}\in\mathbb{H}$, such that*

[TABLE] 2. (ii)

For the factor terms $x_{i}$ referred in (i), we have:

1. (a)

${\mathbf{Z}}_{P}\subseteq\displaystyle{\bigcup_{i=1}^{n}[x_{i}]}$. 2. (b)

Each of the congruence classes $[x_{i}];i=1,\ldots,n,$ contains (at least) a zero of $P$. 3. (iii)

If

[TABLE]

is another factorization of $P$ into linear factors, then there exists a permutation $\pi$ of $(1,2,\ldots,n)$ and $h_{i}\in\mathbb{H};i=1,\dots,n$, such that

[TABLE]

The first result in the above theorem is an immediate consequence of the Fundamental Theorem of Algebra for quaternion polynomials and of the Factor Theorem; the proof of the other results can be found in [19] and [31].

Given a polynomial $P(x)=\sum_{k=0}^{n}a_{k}x^{k}$, its conjugate polynomial, denoted by $\overline{P}(x)$, is given by

[TABLE]

It is very simple to verify that, for all $P,Q\in\mathbb{H}[x]$:

[TABLE]

To each quaternion $q$, we will associate the following polynomial

[TABLE]

called the characteristic polynomial of $q$. Since the characteristic polynomial of $q$ only depends on the real part and norm of $q$ and recalling (1), we immediately conclude that ${\mathcal{Q}}_{q}={\mathcal{Q}}_{q^{\prime}}$ if and only if $[q]=[q^{\prime}]$. Note that ${\mathcal{Q}}_{q}$ is a quadratic polynomial with real coefficients. It can also be shown that the zero-set of ${\mathcal{Q}}_{q}$ is the congruence class of $q$, i.e. $\mathbf{Z}_{{\mathcal{Q}}_{q}}=[q]$; see, e.g. [35]. This result already shows that, in what concerns the number of zeros, polynomials in $\mathbb{H}[x]$ can behave very differently from complex polynomials: a polynomial in $\mathbb{H}[x]$ can have an infinite number of zeros. However, as Theorem 3 shows, the zeros of a polynomial of degree $n$ belong to, at most, $n$ congruence classes in $\mathbb{H}$.

The zeros of an unilateral quaternionic polynomial can be of two distinct types, the so-called isolated zeros and spherical zeros, whose definitions we now recall. Let $q$ be a zero of a given polynomial $P$. We say that $q$ is an isolated zero of $P$ if the congruence class of $q$ contains no other zero of $P$. If $q$ is not an isolated zero of $P$, we call it a spherical zero of $P$. Note that, according to the definition, real zeros are always isolated zeros. The next theorem gives conditions under which a non-real zero is a spherical zero (see e.g. [29]).

Theorem 4**.**

Let $q$ be a non-real zero of a given polynomial $P\in\mathbb{H}[x]$. Then, $q$ is a spherical zero of $P$ if and only if any of the following equivalent conditions hold:

1. (i)

$q$* and $\overline{q}$ are both zeros of $P$.* 2. (ii)

$[q]\subseteq{\mathbf{Z}}_{P}$. 3. (iii)

The characteristic polynomial of $q$, ${\mathcal{Q}}_{q}$, is a divisor of $P$, i.e. there exists a polynomial $Q\in\mathbb{H}[x]$ such that $P=Q\ast{\mathcal{Q}}_{q}.$

Recalling that the congruence classes of non-real quaternions can be identified with spheres, condition (ii) justifies the choice of the term spherical to designate this type of zeros. When $q$ is a spherical zero, we also say that $q$ generates the sphere of zeros $[q]$.

3 The Weierstrass method in $\mathbb{H}[x]$

Let $P$ be a complex monic polynomial of degree $n$ with roots $\zeta_{1},\ldots,\zeta_{n}$ and let $z_{1}^{(0)},\ldots,z_{n}^{(0)}$ be $n$ given distinct numbers. The (classical) Weierstrass method for approximating the roots $\zeta_{i}$ is defined by the iterative scheme:

[TABLE]

If the roots $\zeta_{1},\ldots,\zeta_{n}$ are distinct and $z_{1}^{(0)},\ldots,z_{n}^{(0)}$ are sufficiently good initial approximations to these roots, then the method converges at a quadratic rate, as was firstly proven by Dochev [6] (see also [1, 34]). For multiple roots, the method still converges (locally) but the quadratic convergence is lost; see e.g. [11].

Formula (5) is realized in parallel mode and is often called the total-step mode. The convergence of the method can be accelerated by using a different variant that makes use of the most recent updated approximations to the roots as soon as they are available, as follows:

[TABLE]

The above variant of the Weierstrass method is usually referred to as the serial, sequential or single-step mode (see [28] and references therein).

3.1 A quaternionic Weierstrass-like scheme

Our purpose is to adapt the idea of the Weierstrass method to the computation of the zeros of quaternionic polynomials. So, let $P$ be a given monic polynomial of degree $n$ in $\mathbb{H}[x].$ Corresponding to the assumption imposed in the complex case to guarantee the quadratic convergence of the method – i.e. that the zeros of the polynomial are simple – we will now assume that the polynomial $P$ has $n$ distinct isolated roots. By analogy with the complex case, in this situation, we will still say that $P$ has only simple roots. As stated in the previous section, $P$ can be factorized in the form

[TABLE]

with the factor terms $x_{i}\in\mathbb{H}.$ For simplicity, we introduce the following convenient notation, which we borrow and adapt from [13],

[TABLE]

Remark 2**.**

Note that the order of the factors, due to the non-commutativity of the product in $\mathbb{H}[x]$, is important. We also adopt the convention that

[TABLE]

We first present a simple lemma, relating the roots of $P$ with the quaternions involved in any of its factorizations.

Lemma 1**.**

Let $P$ be a (monic) polynomial of degree $n$ in $\mathbb{H}[x]$ with simple roots and let (6) be one of its factorizations. Then:

1. (i)

The congruence classes of the elements $x_{j};j=1,\dots,n,$ in (6) are distinct. 2. (ii)

The roots $\zeta_{1},\ldots,\zeta_{n}$ of $P$ can be obtained from the quaternions $x_{1},\ldots,x_{n}$ as follows:

[TABLE]

where $R_{i}$ are the polynomials given by

[TABLE]

Proof.

The fact that the congruence classes $[x_{j}];j=1,\ldots,n$, are distinct is an immediate consequence of the results in Theorem 3 and of the assumption that $P$ has only simple roots, i.e. it has exactly $n$ isolated roots. The proof that the roots of $P$ are given by (7) is a simple adaptation of the proof of [19, Proposition 16.3]. ∎

Following the idea of the Weierstrass method in its sequential version, we will now show how to obtain sequences converging, at a quadratic rate, to the factor terms in (6) of a given polynomial $P$. Then, we will show how these sequences can be used to estimate the zeros of $P$.

Theorem 5**.**

Let $P$ be a polynomial of degree $n$ in $\mathbb{H}[x]$ with simple roots and, for $i=1,\ldots,n;$ $k=0,1,2,\ldots$, let

[TABLE]

where

[TABLE]

with ${\mathcal{Q}}_{q}$ denoting the characteristic polynomial of $q$. If the initial approximations $z_{i}^{(0)}$ are sufficiently close to the factor terms $x_{i}$ in a factorization of $P$ in the form (6), then the sequences $\{z_{i}^{(k)}\}$ converge quadratically to $x_{i}$.

Proof.

Let $z_{i}^{(k)}$ be approximations to $x_{i}$ with errors $\varepsilon_{i}^{(k)}$, i.e.

[TABLE]

and let

[TABLE]

We assume that $\varepsilon^{(k)}$ is small enough, i.e. that $z_{i}^{(k)}$ are sufficiently good approximations to $x_{i}$. We want to show that the next iterates $z_{i}^{(k+1)}$ are approximations to $x_{i}$ with errors $\varepsilon_{i}^{(k+1)}$ such that

[TABLE]

We will do this by induction on $i$. For simplicity, we will omit the iteration superscript $(k)$, writing simply $z_{i}$ for $z_{i}^{(k)}$, $\varepsilon_{i}$ for $\varepsilon_{i}^{(k)},{\mathcal{L}}_{i}$ for ${\mathcal{L}}_{i}^{(k)}$ etc. and will replace the superscript $(k+1)$ by a tilde symbol, using $\tilde{z_{i}}$ for $z_{i}^{(k+1)}$, $\tilde{\varepsilon_{i}}$ for $\varepsilon_{i}^{(k+1)}$, etc.

Step 1: We first prove that the result is true for $i=1$, i.e. that we have $\tilde{\varepsilon}_{1}={\mathcal{O}}(\varepsilon^{2}).$

By making use of (61), we can rewrite the polynomial $P(x)$ as

[TABLE]

where $\mathscr{E}_{1}(x)$ designates a remainder polynomial consisting of a sum of $n-1$ terms of the form

[TABLE]

($j=2,\ldots,n$) with terms with $\ast$-products involving at least two $\varepsilon_{j}$’s. By using the definition (26) of the polynomial ${\mathcal{L}}_{1}$, we can write $P(x)$ in the following form

[TABLE]

Let $\overline{{\mathcal{L}}}_{1}$ be the conjugate of ${\mathcal{L}}_{1}$ and note that $\overline{{\mathcal{L}}}_{1}\ast{\mathcal{L}}_{1}$ is precisely the real polynomial ${\boldsymbol{\mathcal{Q}}}_{1}$ defined by (60). Hence, if we multiply $P(x)$ on the left by $\overline{{\mathcal{L}}}_{1}$ and evaluate the resulting polynomial at the point $x=z_{1}$, we obtain, recalling the results (3) and (4) in Theorem 1,

[TABLE]

where $\hat{z}_{1}=\varepsilon_{1}z_{1}\varepsilon_{1}^{-1}$. Observing that we may assume that we are working in a bounded domain ${\mathcal{D}}$ of $\mathbb{H}$ (a sufficiently large disk containing all $z_{i}$) and recalling the definition of $\mathscr{E}_{1}$, it is easily seen that we have

[TABLE]

and therefore

[TABLE]

Since we are assuming that the congruence classes $[x_{j}]$ are distinct, then, for sufficiently small $\varepsilon$, $|{\boldsymbol{\mathcal{Q}}}_{1}(z_{1})|$ is bounded away from zero and so, by multiplying both sides of the above equality on the right by \bigl{(}{\boldsymbol{\mathcal{Q}}}_{1}(z_{1})\bigr{)}^{-1}, we obtain

[TABLE]

or, in other words (cf. (61)),

[TABLE]

which means that the next approximation to $x_{1}$

[TABLE]

is such that

[TABLE]

Step $\boldsymbol{i}$: We now assume that, for $j=1,\ldots,i-1$, ${\tilde{z}_{j}}$ approximates $x_{j}$ with an error $\tilde{\varepsilon}_{j}$ such that $\tilde{\varepsilon}_{j}={\mathcal{O}}(\varepsilon^{2})$ and prove that ${\tilde{z}_{i}}$ is also an ${\mathcal{O}}(\varepsilon^{2})$ approximation to $x_{i}$.

Using the polynomials

[TABLE]

we can write

[TABLE]

where $\mathscr{E}_{i}$ and $\tilde{\mathscr{E}}_{i}$ are remainder polynomials defined in an analogous manner to $\mathscr{E}_{1}$, with the obvious modifications. Note that $\mathscr{E}_{i}$ is a sum of terms, all of which involve at least the product by a $\varepsilon_{j}$ ($j\in\{i+1,\ldots,n\}$) and $\tilde{\mathscr{E}}_{i}$ a sum of terms, all of which involve at least the product by an $\tilde{\varepsilon}_{j}$ ($j\in\{1,\ldots,i-1\}$). Therefore

[TABLE]

Hence the polynomial $P$ can be written as

[TABLE]

Multiplying both sides of the last equality on the left by $\overline{{\mathcal{L}}}_{i}$ and on the right by $\overline{{\mathcal{R}}}_{i}$ and evaluating at $x=z_{i}$, we obtain

[TABLE]

where we made use of the fact that ${\mathcal{R}}_{i}\ast\overline{{\mathcal{R}}}_{i}$ is a real polynomial and hence commutes with any other polynomial. Observing that $\overline{{\mathcal{L}}}_{i}\ast{\mathcal{L}}_{i}\ast{\mathcal{R}}_{i}\ast\overline{{\mathcal{R}}}_{i}$ is the real polynomial ${\boldsymbol{\mathcal{Q}}}_{i}$, using again the results (3) and (4) in Theorem 1 and having in mind the form of the remainder polynomials ${\mathscr{E}}_{i}$ and $\tilde{\mathscr{E}}_{i}$, we can write

[TABLE]

where $\hat{z}_{i}=\varepsilon_{i}z_{i}\varepsilon_{i}^{-1}$. Multiplying by $\left({\boldsymbol{\mathcal{Q}}}_{i}(z_{i})\right)^{-1}$ on the right and observing, once more, that $|\boldsymbol{{\mathcal{Q}}}_{i}(z_{i})|$ is bounded away from zero, we obtain

[TABLE]

or, equivalently, recalling the definition of the errors $\varepsilon_{i}$,

[TABLE]

showing that

[TABLE]

is an ${\mathcal{O}}(\varepsilon^{2})$ approximation to $x_{i}$, which is precisely the result that we wanted to establish. ∎

Remark 3**.**

We should observe that, for each $i=1,\ldots,n$, formula (9) for the computation of the approximation $z_{i}^{(k+1)}$ to $x_{i}$ involves the polynomials ${\mathcal{R}}_{i}^{(k)}$ and ${\boldsymbol{\mathcal{Q}}}^{(k)}_{i}$ which make use of the already computed $z_{1}^{(k+1)},\ldots,z_{i-1}^{(k+1)}$, i.e. the method here described can be seen as a generalization of the sequential version of the Weierstrass method. A careful analysis of the proof, namely the deduction of formula (3.1), shows that the use of the updated $z_{j}^{(k+1)};j=1,\ldots,i-1,$ when computing $z_{i}^{(k+1)}$, is essential for establishing the quadratic order of convergence of the method.

We now show how, with some additional little effort, one can use the iterative scheme (9)–(60) to produce, not only the factor terms, but also the roots of the polynomial.

Theorem 6**.**

Let $P$ be a monic polynomial of degree $n$ in $\mathbb{H}[x]$ with simple roots and let $\{z_{i}^{(k)}\}$ be the sequences defined by the Weierstrass iterative scheme (9)–(60) under the assumptions of Theorem 5. Finally, let $\{\zeta_{i}^{(k)}\}$ be the sequences defined by

[TABLE]

where ${\mathcal{R}}_{i}^{(k)}$ are the polynomials given by (43). Then, $\{\zeta_{1}^{(k)}\},\ldots,\{\zeta_{n}^{(k)}\}$ converge quadratically to the roots of $P$.

Proof.

We start by first recalling that the roots $\zeta_{i}$; $i=1,\ldots,n,$ of $P$ are related to $x_{i}$ in (6) through (7), i.e.

[TABLE]

with $R_{i}$ defined by (8).

Next, denote by $\varepsilon_{i}^{(k+1)}$ the errors in the approximations $z_{i}^{(k+1)}$ to $x_{i}$ and let $\varepsilon^{(k+1)}:=\max_{i}|\varepsilon_{i}^{(k+1)}|$. We will show that

[TABLE]

This, conjugated with the results of Theorem 5, will prove the assertion of the theorem.

Similarly to what we did in the proof of Theorem 5, we will simply write $\tilde{z}_{i}$ for $z_{i}^{(k+1)}$, $\tilde{\varepsilon}_{i}$ for $\varepsilon_{i}^{(k+1)}$, $\tilde{\varepsilon}$ for $\varepsilon^{(k+1)}$, $\tilde{\zeta}_{i}$ for $\zeta_{i}^{(k+1)}$ and ${\mathcal{R}}_{i}$ for ${\mathcal{R}}_{i}^{(k)}$. Taking into account that the polynomials $R_{i}$ in (8) are exactly the same polynomials presented in (127), we can write

[TABLE]

and therefore

[TABLE]

Expressing $\overline{{\mathcal{R}}_{i}}(x)$ in the expanded form $\sum_{j=1}^{i-1}\bar{r}_{j}x^{j}$, it follows at once that

[TABLE]

Combining the fact that both $|\overline{R_{i}}(x_{i})|$ and $|\overline{\mathcal{R}}_{i}(\tilde{z}_{i}|$ are bounded away from zero with the result (130), we can conclude that

[TABLE]

Finally, result (7) together with (130), (131) and the assumption (129) gives

[TABLE]

which is precisely the result we want to prove. ∎

3.2 Computational details

We now summarize the proposed algorithm for computing the roots of a given quaternionic unilateral polynomial $P$ of degree $n$ and make some practical comments regarding its implementation.

Quaternionic-Weierstrass algorithm

Input:

polynomial coefficients

• -

initial values $z_{i}^{(0)}$

• -

error tolerances $\varepsilon_{1},\varepsilon_{2}$

• -

maximum number of iterations $kmax$

Set $\zeta_{i}^{(0)}=z_{i}^{(0)}$ 2. 2.

For $k=1,2,\ldots$ until Stopping Criterion is true

1. (a)

Compute $z_{i}^{(k)}$, by means of (9)-(60). 2. (b)

Compute $\zeta_{i}^{(k)}$, by means of (129) and (43).

Stopping Criterion:

(\displaystyle{\max_{i}}\big{|}\zeta_{i}^{(k)}-\zeta_{i}^{(k-1)}\big{|}<\varepsilon_{1} and \displaystyle{\max_{i}}\big{|}P(\zeta_{i}^{(k)})\big{|}<\varepsilon_{2}\Big{)} or $k=kmax$.

Ouput: Factors $\tilde{x}_{i}=z_{i}^{(k)}$ and roots $\tilde{\zeta}_{i}=\zeta_{i}^{(k)}$.

Choice of initial approximations

In the classical case, Weierstrass method seems in practice to converge from nearly all starting points (see [20] and the references therein for details). The numerical experiments that we have conducted also show the robustness of the quaternionic version of the method in what concerns the choice of initial approximations. In any case, there are some aspects that should be taken into account.

First, for formula (9) to be meaningful, a first requirement one has to have in mind when choosing the initial approximations $z_{1}^{(0)},\ldots,z_{n}^{(0)}$ is that all of them belong to distinct congruence classes. This does not necessarily guarantee that, in the course of the computations, two approximations do not fall into the same congruence class, although this is very unlikely to happen. In such a case, a small perturbation of the initial guesses should be sufficient to regain convergence.

Second, it is, naturally, convenient to select the initial approximations from a region where the $x_{i}$ in any factorization of the polynomial $P$ are known to lie. Since the $x_{i}$ and the roots $\zeta_{i}$ of $P$ have the same norm, bounds on $|\zeta_{i}|$ are also valid for $|x_{i}|$. Moreover, since $P(x)\ast\overline{P}(x)$ is a real polynomial111The use of the polynomial $P(x)\ast\overline{P}(x)$ goes back to the work of Niven [23]. , whose roots $r_{i}$ also have the same norm as the roots $\zeta_{i}$ of $P$ (this is an immediate consequence of Theorem 2 in [31] and the characterization of the congruence classes given by (1)), one can use any known result on bounds on (complex) polynomial roots to obtain a region from where the initial approximations should be selected.

Non simple zeros

The proof of Theorem 5 was done under the assumption that the roots $\zeta_{1},\dots,\zeta_{n}$ of the polynomial (6) are simple, i.e. that $[\zeta_{i}]\neq[\zeta_{j}]$ for all $i\neq j$, or, equivalently, $[x_{i}]\neq[x_{j}]$ (cf. Lemma 6). When $[x_{i}]=[x_{j}]$ for some $i\neq j$, the characterization of the zero-set of the polynomial can be done taking into account the following two results.

Lemma 2**.**

If $x_{1},x_{2}\in\mathbb{H}$ and $h=\overline{x}_{2}-x_{1}$, then

[TABLE]

Proof.

The result follows by simple manipulation; see also [31] for a different, but equivalent result. ∎

Lemma 3**.**

Consider a quadratic polynomial factorized in the form

[TABLE]

where $x_{1},x_{2}\in\mathbb{H}\setminus\mathbb{R}$ and $[x_{1}]=[x_{2}]$.

1. (i)

If $x_{1}\neq\overline{x}_{2}$, then the only zero of $P$ is $x_{1}$. 2. (ii)

If $x_{1}=\overline{x}_{2}$, then $x_{1}$ generates the sphere of zeros $[x_{1}]$, i.e. $x_{1}$ is a spherical zero.

Proof.

See e.g. [13]. ∎

The problem of finding a natural definition of multiplicity for zeros of quaternionic polynomials is a rather complicated task and as a consequence one can find in the literature different (not always equivalent, see [9]) concepts of multiplicity [2, 3, 12, 14, 26, 32]. In the case (i) above, we will say that $x_{1}$ is a root with (isolated) multiplicity equal to two. For example, the polynomials $(x+1-\mathbf{i})\ast(x+1+\mathbf{k})$ and $(x+1+\mathbf{k})\ast(x+1+\mathbf{k})$ both have $\zeta=-1-\mathbf{k}$ as a root with multiplicity two.

Returning to the case of a general polynomial of degree $n$ of the form (6), we consider, for simplicity, that $[x_{i}]$ and $[x_{j}]$ are the only non-distinct congruent classes. Using Lemma 2, we can freely move the factors $(x-x_{i})$ and $(x-x_{j})$ to the right of the factorization without changing the set of congruence classes so that in the new factorization

[TABLE]

we have $[y_{1}]=[y_{2}]$. Observe that all the roots of the $n-2$ degree polynomial $Q$ are simple and, therefore, the complete characterization of the roots of $P$ can be done applying Lemma 3 to the quadratic polynomial $(x-y_{2})\ast(x-y_{1})$.

We considered the application of the quaternionic Weierstrass method to several examples of polynomials having double (isolated) or spherical roots and, in all the cases, we have observed the following: when $y_{1}$ is a double isolated root ($y_{1}\neq\overline{y}_{2}$), the behavior is analogous to the one observed in the classical case, i.e. the rate of convergence drops to one; on the other hand, if $y_{1}$ is a spherical root, the iterative scheme produces two distinct roots $\zeta_{1}$ and $\zeta_{2}$ belonging to the congruence class $[y_{1}]$ and still shows a quadratic order of convergence.

4 Numerical examples

In this section we present several examples illustrating the performance of the quaternionic Weierstrass method introduced in Section 3.

All the numerical experiments here reported were obtained by the use of the Mathematica add-on application QuaternionAnalysis [21] designed by two of the authors of this paper for symbolic manipulation of quaternion valued functions. A collection of new functions, including an implementation of the Weierstrass method described in this paper, has been recently developed in order to endow the aforementioned package with the ability to perform operations in the non-commutative ring of polynomials $\mathbb{H}[x]$.

Example 1**.**

Our first test example is a polynomial which fulfills the assumptions of Theorem 6. In fact, it is easy to see that the polynomial

[TABLE]

has only simple roots, namely

[TABLE]

Since, in this case, the polynomial roots $\zeta_{i}$ are known exactly, we replace the stopping criterion based on the incremental size of the iterations by the following one:

[TABLE]

where $\pi_{k}$ is an appropriate permutation of $\{1,\dots,6\}$. Here, we considered $\varepsilon_{1}=\varepsilon_{2}=10^{-16}$ and chose initial approximations so that $\epsilon^{(0)}\leq 0.5$.

The Weierstrass method applied to the extended form of $P$ produced, after 5 iterations, the following approximations to the factor terms (with 15 decimal places222The notation $(0)$ after the decimal point represents a sequence of 15 zeros.)

[TABLE]

corresponding to the approximate roots

[TABLE]

It is interesting to observe that the approximations $x_{i}^{(5)}$ to the factor terms lead to a factorization of $P$ different from (132), but of course in line with Theorem 3-(iii).

Table 1 contains the relevant information concerning the errors in the successive approximations $\zeta_{i}^{(k)}$ ($k=0,\dots,5$, $i=1,\dots,6$) to the roots $\zeta_{i}$ of $P$. Estimates $\rho$ for the computational local order of convergence of the method, based on the use of (see e.g. [16] for details).

[TABLE]

were also computed and are included in the last column of the table.

In order to illustrate Remark 3 we have also implemented the parallel version of Weierstrass method. In this case, using the same initial guesses, 9 iterations were required to achieve the same precision. The results presented in Table 2 clearly indicate the deterioration of the speed of convergence of this version of the method.

Our next examples concern situations where the polynomials under consideration have zeros which are not simple.

Example 2**.**

The polynomial

[TABLE]

has, apart from the isolated zeros $-\mathbf{i}+\mathbf{k}$ and $1-\mathbf{k}$, a whole sphere of zeros, $[\,\mathbf{i}\,]$. In this case, since all the spherical roots have the same real part and modulus, we replaced the stopping criterion used in the previous example by the following one:

[TABLE]

where

[TABLE]

Starting with an initial guess chosen so that $\epsilon^{(0)}\leq 0.15$, we obtained, after 5 iterations, the following approximations:

[TABLE]

The spherical root can be identified at once by observing that, up to the required precision, we have $[\zeta_{1}^{(5)}]=[\zeta_{2}^{(5)}]$, since $\operatorname{Re}\zeta_{1}^{(5)}=\operatorname{Re}\zeta_{2}^{(5)}=0$ and $|\zeta_{1}^{(5)}|=|\zeta_{2}^{(5)}|=1$.

The numerical details related to this example are displayed in Table 3. Here the numerical computations have been carried out with the precision increased to 512 significant digits.

As we can observe from Table 3, the quaternionic Weierstrass method works, produces all the roots simultaneously with machine precision and exhibits quadratic order of convergence. As expected, for the case of the spherical root, we obtain convergence to two distinct members of the sphere of zeros.

Example 3**.**

In our last example we address the problem of using Weierstrass method in cases where the polynomial under consideration has multiple (isolated) roots. The polynomials

[TABLE]

have one non-real root with multiplicity one and $-1-\mathbf{k}$ as a double root. The approximations to the roots of $P$ obtained by the use of the quaternionic Weierstrass method are

[TABLE]

while, for the roots of $Q$, we obtained

[TABLE]

As we can observe from Table 4, the behavior of the quaternionic Weierstrass method is very similar to that one observed for the classical complex case, where the rate of convergence is linear. This table shows $\epsilon^{(k)}$ for the last 9 iterations of the method together with $\rho^{(k)}$ for both polynomials.

5 Final Remarks

In this paper we proposed a generalization to the quaternionic context of the well-known Weierstrass method for approximating all zeros of a polynomial simultaneously. Due to the structure of the zero-set of a quaternionic polynomial, the claim that the method we have proposed produces all the zeros simultaneously, requires an additional explanation. Assuming the convergence of the method to the roots $\zeta_{1},\dots,\zeta_{n}$ of a polynomial $P$ of degree $n$, it is easy to identify ${\mathbf{Z}}_{P}$, once we test if each element of $\{\zeta_{1},\dots,\zeta_{n}\}$ is an isolated or a spherical zero of $P$ (cf. Theorem 4).

The quaternionic Weierstrass algorithm is entirely based on quaternionic arithmetic and shows fast convergence for simple and spherical roots. We proved the quadratic convergence of the sequential iterative scheme, under the assumptions that all the roots of the polynomial are simple, and presented numerical examples supporting this fact. In [8], it was proved that the same rate of convergence can be achieved by quaternion versions of Newton’s method for the so-called radially holomorphic functions [17, p. 234]. None of the polynomials presented in this section are radially holomorphic or are in the less restrictive conditions of [8, Theorem 4]. As far as we are aware, the method proposed in this paper is the first numerical method entirely based on quaternionic arithmetic for which we can observe theoretical and experimental results for general unilateral quaternion polynomials.

Several authors, namely Petkovic and collaborators (see e.g. [27]), have described conditions for the safe convergence of the classical method depending only on the initial approximations. This is a very interesting question that we intend to address in the near future, in the quaternionic case.

One can find in the literature several modifications to the classical Weierstrass method which improve the speed of convergence to multiple roots (see e.g. [11]). It is also in our plans of research to consider adaptations of such strategies.

Acknowledgments

Research at 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 NIPE has been carried out within the funding with COMPETE reference number POCI-01-0145-FEDER-006683, with the FCT/MEC’s financial support through national funding and by the ERDF through the Operational Programme on “Competitiveness and Internationalization – COMPETE 2020” under the PT2020 Partnership Agreement.