Computational Aspects of Quaternionic Polynomials

Part II: Root-Finding Methods

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

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.

Introduction

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

P(x)=(Subscript[a, n])(x^n)+(Subscript[a, n-1])(x^(n-1))+…+(Subscript[a, 1])x+(Subscript[a, 0,]) with (Subscript[a, n])≠0 and (Subscript[a, k])∈ℍ (the quaternions),

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 ([]) 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. [–]).

This article discusses two numerical methods to approximate the zeros (or roots) of polynomials of the form (). 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 [–].

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

Newton-Like Methods

Theoretical Framework

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

The real vector space (ℝ^4) can be identified with ℍ by means of

x:=((Subscript[x, 0]),(Subscript[x, 1]),(Subscript[x, 2]),(Subscript[x, 3]))∈(ℝ^4)⟷ x=(Subscript[x, 0])+i (Subscript[x, 1])+j (Subscript[x, 2])+k (Subscript[x, 3])∈ℍ,

where i, j and k are Hamilton’s imaginary units. Thus, throughout the article, we do not distinguish an element in (ℝ^4) from the corresponding quaternion in ℍ, unless we need to stress the context.

Using the simplified notation (UnderBar[x]):=i (Subscript[x, 1])+j (Subscript[x, 2])+k (Subscript[x, 3]) for the vector part of x, any arbitrary nonreal quaternion x can be written as

x=(Subscript[x, 0])+(UnderBar[x])=(Subscript[x, 0])+ω((UnderBar[x]))|(UnderBar[x])|,

where |(UnderBar[x])| is the norm of (UnderBar[x]) and ω((UnderBar[x])) is the quaternion

ω((UnderBar[x]))=((UnderBar[x])/(|(UnderBar[x])|)),

also referred to as the sign of (UnderBar[x]). In addition, since  ((ω((UnderBar[x])))^2)=-1 and |ω((UnderBar[x]))|=1, one can say that ω((UnderBar[x])) behaves like the complex imaginary unit, and for this reason we call () the complex-like form of the quaternion x.

In what follows, we consider domains Ω⊂(ℝ^4)≅ℍ and functions f:Ω→ℍ that can be written in the form

f(x)=f((Subscript[x, 0])+ω r)=u((Subscript[x, 0]),r)+ω v((Subscript[x, 0]),r),

where r:=|(UnderBar[x])|, ω:=ω((UnderBar[x])) and u and v are real-valued functions. Continuity and differentiability are defined coordinate-wise.

We define on the set (𝒞^1)(Ω,ℍ) the so-called radial operators

(∂_rad):=(1/2)((∂_0)-ω(∂_r)) and (Subscript[(OverBar[∂]), rad]):=(1/2)((∂_0)+ω(∂_r)),

where (∂_0):=(∂/(∂(Subscript[x, 0]))) and (∂_r):=(∂/(∂(Subscript[x, r]))).

We introduce the following concept.

Definition 1

Let f be a function of the form (), x∈Ω and h=(Subscript[h, 0])+ω((UnderBar[x]))(Subscript[h, r]), with (Subscript[h, 0]),(Subscript[h, r])∈ℝ. Such a function f is called radially holomorphic (or radially regular) in Ω if

(lim_(h→0))(f(x+h)-f(x))(h^(-1))

Theorem 1

A function f of the form () is radially holomorphic iff (Subscript[(OverBar[∂]), rad])f=0. In that case, we have f'=(∂_rad)f=(∂_0)f=(∂_0)u((Subscript[x, 0]),r)+ω (∂_0)v((Subscript[x, 0]),r).

It follows at once that any quaternionic polynomial of the form () but with (Subscript[a, k])∈ℝ is radially holomorphic and its radial derivative is

P'(x)=n (Subscript[a, n])(x^(n-1))+(n-1) (Subscript[a, n-1])(x^(n-2))+…+(Subscript[a, 1]).

Quaternionic Newton Method

For holomorphic complex functions f of one complex variable, the well-known Newton method for finding a zero (SuperStar[x]) consists of approximating (SuperStar[x])∈ℂ by means of the iterative process

2D-NM: (Subscript[z, k+1])=(Subscript[z, k])-((f((Subscript[z, k])))/(f'((Subscript[z, k])))),k=0,1,2,…,

with (Subscript[z, k])∈ℂ sufficiently close to (SuperStar[x])∈(ℝ^2) and f'((Subscript[z, k]))≠0. Identifying a real quaternion with a vector in (ℝ^4), 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 ():

4D-NM: (Subscript[z, k+1])=(Subscript[z, k])-(((J f((Subscript[z, k]))))^(-1))f((Subscript[z, k])), k=0,1,2,…,

with (Subscript[z, k])∈(ℝ^4) sufficiently close to (SuperStar[x])∈(ℝ^4) and a nonsingular Jacobian matrix J f((Subscript[z, k])). Not surprisingly, recent experiments performed by some of the authors of this article ([], []) 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 [], where the authors solved equations of the form (x^n)=a,a∈ ℍ. Later, Kalantari in [], using algebraic-combinatorial arguments, proposed a Newton method for finding roots of special quaternionic polynomials. In [], the equivalence between the classical multivariate Newton method () 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 p and q may be interpreted in two different ways: either as p (q^(-1)) (the right quotient) or (q^(-1))p (the left quotient). This leads naturally to considering two versions of Newton iteration in the quaternionic setting:

ℍ-(Subscript[NM, right]): (Subscript[z, k+1])=(Subscript[z, k])-f((Subscript[z, k])) (((f'((Subscript[z, k]))))^(-1)), k=0,1,2,…,

ℍ-(Subscript[NM, left]): (Subscript[z, k+1])=(Subscript[z, k])-(((f'((Subscript[z, k]))))^(-1))f((Subscript[z, k])) , k=0,1,2,….

The derivative in equations () and () has been considered in [] and [] as the radial derivative of a radially holomorphic function. In fact, in Corollary 2 of [] it was proved that for such functions, equations (), () and () produce, for each (Subscript[z, 0]), the same sequence, provided that J f((Subscript[z, k])), k=0,1,… is nonsingular. Here is a more general result.

Theorem 2 ([], Theorem 4)

Let f(x)=(Underoverscript[∑ , (i=0), s])(Subscript[α, i])(Subscript[f, i])(x) be a function defined on the set (𝒞^1)(Ω,ℍ) such that the (Subscript[f, i]),  i=0,…,s, are radially holomorphic functions in Ω and the (Subscript[α, i]) are quaternions not all zero. If (SuperStar[z]) is a root of f such that J f((SuperStar[z])) is nonsingular and J f is Lipschitz continuous on a neighborhood of (SuperStar[z]), then for all (Subscript[z, 0])∈ℍ sufficiently close to (SuperStar[z]) such that ω((UnderBar[(Subscript[z, 0])])) commutes with all ω((UnderBar[(Subscript[α, i])])), the Newton processes

ℍ-(Subscript[NM, right]): (Subscript[z, k+1])=(Subscript[z, k])-f((Subscript[z, k]))((((∑_(i=0)^s)(Subscript[α, i])(Subscript[f, i])'((Subscript[z, k]))))^(-1)),k=0,1,2,…

ℍ-(Subscript[NM, left]): (Subscript[z, k+1])=(Subscript[z, k])-((((∑_(i=0)^s)(Subscript[α, i])(Subscript[f, i])'((Subscript[z, k]))))^(-1))f((Subscript[z, k])),k=0,1,2,…

Each step k of the iterative schemes () and () is implemented in the function NewtonIterativeFunction, which has as arguments the quaternion (Subscript[z, k]) and the indication of the version: “right” for () or “left” for (). At each step, a test of the value of f'((Subscript[z, k])) is also performed. We recall again that all the functions presented here require the package QuaternionAnalysis.

Needs["QuaternionAnalysis`"]

NewtonIterativeFunction[f_,Df_,q_,"right",t_:10^-16]:=If[Abs@(Df@@q)>t,q- f@@q**(1/Df@@q),Null]

NewtonIterativeFunction[f_,Df_,q_,"left",t_:10^-16]:=If[Abs@(Df@@q)>t,q- (1/Df@@q)**f@@q,Null]

The ℍ-Newton methods consist of the successive application of the iterative schemes () or () through the function NewtonIterations, using a stopping criteria based on the incremental size |(Subscript[z, k+1])-(Subscript[z, k])| and on the maximum number of iterations ItMax.

TestConvergence[z1_,z2_,eps_]:=Norm[z1-z2] >eps

NewtonIterations[f_,Df_,z0_,method_,eps_:10^-12,ItMax_:20]:=NestWhileList[NewtonIterativeFunction[f,Df,#,method]&,z0,TestConvergence[##,eps]&,2,ItMax]

Example 1

Consider the radially holomorphic polynomial p(x)=(x^3)-x, whose only roots in ℍ are the real isolated roots -1, 1 and 0. For the concepts of isolated and spherical roots, we refer the reader to [], Definition 4.

f1[X0_,X1_,X2_,X3_]=QPower[Quaternion[X0,X1,X2,X3],3]-Quaternion[X0,X1,X2,X3];
Df1=Function[{X0,X1,X2,X3},Evaluate[Map[D[#,X0]&,f1[X0,X1,X2,X3]]]];

NewtonIterations[f1,Df1,Quaternion[1.,-1.,1.,-1.],"right"]//Column

NewtonIterations[f1,Df1,Quaternion[1.,-1.,1.,-1.],"left"]//Column

The use of the initial guesses (Subscript[z, 0])=1+j and (Subscript[z, 0])=1-j requires 14 iterations to get an approximation to the roots 1 and -1, respectively.

NewtonIterations[f1,Df1,Quaternion[1.,0.,1.,0.],"right"]//Column

NewtonIterations[f1,Df1,Quaternion[1.,0.,1.,0.],"left"]//Column

NewtonIterations[f1,Df1,Quaternion[-1.,0.,1.,0.],"right"]//Column

NewtonIterations[f1,Df1,Quaternion[-1.,0.,1.,0.],"left"]//Column

Example 2

The polynomial p(x)=(x^3)+x has a real root 0 and the sphere of zeros [i]. 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.

f2[X0_,X1_,X2_,X3_]=QPower[Quaternion[X0,X1,X2,X3],3]+Quaternion[X0,X1,X2,X3];
Df2=Function[{X0,X1,X2,X3},Evaluate[Map[D[#,X0]&,f2[X0,X1,X2,X3]]]];

(s1=NewtonIterations[f2,Df2,Quaternion[1.,-1.,1.,-1.],"right"])//Last

NewtonIterations[f2,Df2,Quaternion[-1.,1.,0.,1.],"right"]//Last

NewtonIterations[f2,Df2,Quaternion[-1.,2.,3.,4.],"right"]//Last

W/@s1//Column

Quaternion[0,𝓍,𝓎,0]- f2[0,𝓍,𝓎,0]**(1/Df2[0,𝓍,𝓎,0])

InitialGuess=
{
Quaternion[0.,1.,2.,0.],
Quaternion[0.,-1.,-1.5,0.],
Quaternion[0.,-.3,.3,0.],
Quaternion[0.,.02,.7,0.],
Quaternion[0.,1.2,.8,0.],
Quaternion[0.,.4,.5,0.],
Quaternion[0.,-.4,.5,0.]
};
seqs4D=Apply[List,NewtonIterations[f2,Df2,#,"right"]&/@InitialGuess,{2}];
seqs2D=Map[Take[#,{2,3}]&,seqs4D,{2}];
seqs2D=Map[(Transpose@{#1,Range[Length@#1]})&,seqs2D];

Show[
Plot[{Sqrt[1-x^2],-Sqrt[1-x^2]},{x,-2,2},PlotRange->All,PlotStyle->{Red}],
Graphics[{Circle[#,.03]&/@(List@@@(Rest@Most@#&/@InitialGuess)),
Apply[{AbsolutePointSize[8/#2],Point[#1]}&,seqs2D,{2}],
Map[Line,Map[First,seqs2D,{2}]]}],
Frame->True,FrameLabel-> {x,y},AspectRatio->Automatic,Ticks->None,Axes->False
]

Example 3

Now consider the polynomial p(x)=(x^3)-j (x^2)- x+j  with the three isolated roots (Subscript[r, 1])=1, (Subscript[r, 2])=-1 and (Subscript[r, 3])=j (cf. [], Example 3). This polynomial is not radially holomorphic, which means that we cannot anticipate the behavior of Newton methods unless we choose initial guesses (Subscript[z, 0]) such that Theorem 2 applies, that is, such that ω((UnderBar[(Subscript[z, 0])])) commutes with ω(j)=j. In other words, (Subscript[z, 0]) must be of the form a+j b.

f3[X0_,X1_,X2_,X3_]=QPower[Quaternion[X0,X1,X2,X3],3]-Quaternion[0,0,1,0]**QPower[Quaternion[X0,X1,X2,X3],2]-Quaternion[X0,X1,X2,X3]+Quaternion[0,0,1,0];
Df3=Function[{X0,X1,X2,X3},Evaluate[Map[D[#,X0]&,f3[X0,X1,X2,X3]]]];

NewtonIterations[f3,Df3,Quaternion[1.,0.,2.,0.],"right"]//Last

NewtonIterations[f3,Df3,Quaternion[1.,1.,2.,0.],"right"]//Column

NewtonIterations[f3,Df3,Quaternion[1.,1.,2.,0.],"left"]//Column

NewtonIterations[f3,Df3,Quaternion[1.31,2.,0.,0.],"right"]//Column

NewtonIterations[f3,Df3,Quaternion[1.31,2.,0.,0.],"left"]//Column

Following [] and [], consider a function Iterations 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 (Subscript[z, 0]) as the initial guess.

Iterations[f_,Df_,z0_,method_,eps_,ItMax_]:=Module[
{seq=NewtonIterations[f,Df,z0,method,eps,ItMax],nseq},
nseq=Length@seq;
If[Last@seq===Null||nseq>ItMax,Null,nseq]
]

We now consider different initial guesses (Subscript[z, 0]) by choosing points in special regions Ω=Ω(x,y)⊂(ℝ^4) and we show density plots of Iterations. The white regions that may appear correspond to a choice of (Subscript[z, 0])∈Ω for which the method under consideration does not reach the level of precision eps with ItMax iterations. The default choices of eps=10^-2 and ItMax=20 usually lead to realistic plots that require some minutes to be produced. A smoother density can be obtained by increasing the option PlotPoints.

ViewSolution[f_,Df_,domain_,method_,eps_:10^-2,ItMax_:20]:=
DensityPlot[Iterations[f,Df,N@First@domain,method,eps,ItMax],Evaluate[Sequence@@(Rest@domain)],
PlotRange->{0,ItMax},PlotPoints ->50,AspectRatio -> Automatic,MaxRecursion->1,FrameTicksStyle->Directive[7],ColorFunctionScaling-> False,ColorFunction->(ColorData["Rainbow"][#/ItMax]&),PlotLegends->Automatic]

Example 4

We consider again the polynomial p(x)=(x^3)-j (x^2)- x+j  of Example 3, whose roots are the isolated roots -1, 1 and j. The following code produces the plots corresponding to the choice of (Subscript[z, 0]) in one of the following regions:

(Subscript[Ω, 1])={(x,0,y,0):x,y∈[-2,2]}},

(Subscript[Ω, 2])={(x,y,0,0):x,y∈[-2,2]}},

(Subscript[Ω, 3])={(0,x,y,0):x,y∈[-2,2]}}.

As was already pointed out, Theorem 2 can be applied only in (Subscript[Ω, 1]); this is why both methods produce the same plots in this case.

Show[ViewSolution[f3,Df3,{Quaternion[x,0,y,0],{x,-2,2},{y,-2,2}},"left"]]
Show[ViewSolution[f3,Df3,{Quaternion[x,0,y,0],{x,-2,2},{y,-2,2}},"right"]]

Here is the behavior of the ℍ-Newton methods in (Subscript[Ω, 2]).

Show[ViewSolution[f3,Df3,{Quaternion[x,y,0,0],{x,-2,2},{y,-2,2}},"left"]]
Show[ViewSolution[f3,Df3,{Quaternion[x,y,0,0],{x,-2,2},{y,-2,2}},"right"]]

Here is the behavior of the ℍ-Newton methods in (Subscript[Ω, 3]).

Show[ViewSolution[f3,Df3,{Quaternion[0,x,y,0],{x,-2,2},{y,-2,2}},"left"]]
Show[ViewSolution[f3,Df3,{Quaternion[0,x,y,0],{x,-2,2},{y,-2,2}},"right"]]

Basin of Attraction of a Root

The plots produced by ViewSolution 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 NewtonIterations with the information of the root for which we want to compute the basin of attraction. A new function NewtonIterationsSphere takes into account the existence of spheres of zeros. The functions Iterations and IterationsSphere give the number of iterations needed to observe convergence to an isolated root or a spherical one, respectively. These functions return Null when the corresponding convergence test fails.

TestConvergenceSphere[z1_,z2_,eps_]:=Abs[Re[z1]-Re[z2]]>eps||Abs[AbsVec[z1]-AbsVec[z2]]>eps

NewtonIterations[f_,Df_,z0_,root_,method_,eps_,ItMax_]:=NestWhileList[NewtonIterativeFunction[f,Df,#,method]&,z0,TestConvergence[#,root,eps]&,1,ItMax]

NewtonIterationsSphere[f_,Df_,z0_,root_,method_,eps_,ItMax_]:=NestWhileList[NewtonIterativeFunction[f,Df,#,method]&,N[z0],TestConvergenceSphere[#,root,eps]&,1,ItMax]

Iterations[f_,Df_,z0_,root_,method_,eps_,ItMax_]:=Module[{seq=NewtonIterations[f,Df,z0,root,method,eps,ItMax],nseq},nseq=Length@seq;
If[Last@seq===Null||nseq>ItMax,Null,nseq]]

IterationsSphere[f_,Df_,z0_,root_,method_,eps_,ItMax_]:=Module[{seq=NewtonIterationsSphere[f,Df,z0,root,method,eps,ItMax],nseq},nseq=Length@seq;
If[Last@seq===Null||nseq>ItMax,Null,nseq]]

The functions that plot the basin of attraction of an isolated root Basin or a spherical root BasinSphere have an input parameter color associated to that root. The color coding used is the following: if the initial guess (Subscript[z, 0]), chosen in a domain Ω=Ω(x,y)⊂(ℝ^4), causes the process to converge to a certain isolated root r to which the color (Subscript[c, r]) was associated, then the point (Subscript[z, 0]) is plotted with the color (Subscript[c, r]). For a sphere of zeros [s], all the points (Subscript[z, 0]) that converge to a point in [s] have the color assigned to [s]. 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.

Basin[f_,Df_,domain_,root_,method_,color_,eps_:10^-2,ItMax_:20]:=
DensityPlot[Iterations[f,Df,N@First@domain,root,method,eps,ItMax],Evaluate[Sequence@@(Rest@domain)],PlotRange->{0,ItMax},PlotPoints ->50,AspectRatio -> Automatic,MaxRecursion->1,FrameTicksStyle->Directive[7],ColorFunctionScaling-> False,ColorFunction->(Lighter[color,#/ItMax]&)]

BasinSphere[f_,Df_,domain_,root_,method_,color_,eps_:10^-2,ItMax_:20]:=
DensityPlot[IterationsSphere[f,Df,N@First@domain,root,method,eps,ItMax],Evaluate[Sequence@@(Rest@domain)],PlotRange->{0,ItMax},PlotPoints ->50,AspectRatio -> Automatic,MaxRecursion->1,FrameTicksStyle->Directive[7],ColorFunctionScaling-> False,ColorFunction->(Lighter[color,#/ItMax]&)]

Example 5

We consider once more the polynomial p(x)=(x^3)-j (x^2)- x+j of Example 4, now from the perspective of the basins of attraction of each of the roots -1, 1 and j. We associate with these roots the colors red, blue and green, respectively, and consider the domains (Subscript[Ω, 1]), (Subscript[Ω, 2]) and (Subscript[Ω, 3]), 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 (Subscript[Ω, 1]) (left).

Show[Basin[f3,Df3,{Quaternion[x,0,y,0],{x,-2,2},{y,-2,2}},1,"left",Darker@Red],Basin[f3,Df3,{Quaternion[x,0,y,0],{x,-2,2},{y,-2,2}},-1,"left",Darker@Blue],Basin[f3,Df3,{Quaternion[x,0,y,0],{x,-2,2},{y,-2,2}},Quaternion[0,0,1,0],"left",Darker@Green],PlotRange->All]

Here are the basins of attraction in (Subscript[Ω, 2]) (left and right).

Quiet@Show[Basin[f3,Df3,{Quaternion[x,y,0,0],{x,-2,2},{y,-2,2}},1,"left",Darker@Red],Basin[f3,Df3,{Quaternion[x,y,0,0],{x,-2,2},{y,-2,2}},-1,"left",Darker@Blue],Basin[f3,Df3,{Quaternion[x,y,0,0],{x,-2,2},{y,-2,2}},Quaternion[0,0,1,0],"left",Darker@Green],PlotRange->All]

Show[Basin[f3,Df3,{Quaternion[x,y,0,0],{x,-2,2},{y,-2,2}},1,"right",Darker@Red],Basin[f3,Df3,{Quaternion[x,y,0,0],{x,-2,2},{y,-2,2}},-1,"right",Darker@Blue],Basin[f3,Df3,{Quaternion[x,y,0,0],{x,-2,2},{y,-2,2}},Quaternion[0,0,1,0],"right",Darker@Green],PlotRange->All]

Here are the basins of attraction in (Subscript[Ω, 3]) (left and right).

Show[Basin[f3,Df3,{Quaternion[0,x,y,0],{x,-2,2},{y,-2,2}},1,"left",Darker@Red],Basin[f3,Df3,{Quaternion[0,x,y,0],{x,-2,2},{y,-2,2}},-1,"left",Darker@Blue],Basin[f3,Df3,{Quaternion[0,x,y,0],{x,-2,2},{y,-2,2}},Quaternion[0,0,1,0],"left",Darker@Green],PlotRange->All]

Show[Basin[f3,Df3,{Quaternion[0,x,y,0],{x,-2,2},{y,-2,2}},1,"right",Darker@Red],Basin[f3,Df3,{Quaternion[0,x,y,0],{x,-2,2},{y,-2,2}},-1,"right",Darker@Blue],Basin[f3,Df3,{Quaternion[0,x,y,0],{x,-2,2},{y,-2,2}},Quaternion[0,0,1,0],"right",Darker@Green],PlotRange->All]

Example 6

This example concerns the polynomial p(x)=(x^3)+x studied in Example 2, which has an isolated root 0 (red) and a sphere of zeros [i] (blue). The corresponding plots can be obtained as follows.

Here are the basins of attraction in (Subscript[Ω, 1]) (left).

Show[Basin[f2,Df2,{Quaternion[x,0,y,0],{x,-2,2},{y,-2,2}},0,
"left",Darker@Red],BasinSphere[f2,Df2,{Quaternion[x,0,y,0],{x,-2,2},{y,-2,2}},Quaternion[0,1,0,0],"left",Darker@Blue],PlotRange->All]

Here are the basins of attraction in (Subscript[Ω, 2]) (left); as expected, the behavior is similar to that in (Subscript[Ω, 1]), since [j]=[i].

Show[Basin[f2,Df2,{Quaternion[x,y,0,0],{x,-2,2},{y,-2,2}},0,"left",Darker@Red],
BasinSphere[f2,Df2,{Quaternion[x,y,0,0],{x,-2,2},{y,-2,2}},Quaternion[0,1,0,0],"left",Darker@Blue],PlotRange->All]

Here are the basins of attraction in (Subscript[Ω, 3]) (left).

Show[Basin[f2,Df2,{Quaternion[0,x,y,0],{x,-2,2},{y,-2,2}},0,"left",Darker@Red],BasinSphere[f2,Df2,{Quaternion[0,x,y,0],{x,-2,2},{y,-2,2}},Quaternion[0,1,0,0],"left",Darker@Blue],PlotRange->All]

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 [] in 1891 and later rediscovered and derived in different ways by Durand [] in 1960, Dočev [] in 1962 and Kerner [] and Prešić [] in 1966.

Let P be a complex monic polynomial of degree n with roots (Subscript[ζ, 1]),…, (Subscript[ζ, n]) and let ((Subscript[z, 1])^((0))),…,((Subscript[z, n])^((0))) be n distinct numbers. The classical Weierstrass method for approximating the roots (Subscript[ζ, i]) is defined by the iterative scheme:

((Subscript[z, i])^((k+1)))=((Subscript[z, i])^((k)))-((P(((Subscript[z, i])^((k)))))/((∏_(j=1,j≠i)^n)((Subscript[z, i])^((k)))-((Subscript[z, j])^((k))))); i=1,… ,n; k=0,1,….

If the roots (Subscript[ζ, 1]),… ,(Subscript[ζ, n]) are distinct and ((Subscript[z, 1])^((0))),… ,((Subscript[z, n])^((0))) are sufficiently good initial approximations to these roots, then the method converges at a quadratic rate, as was first proved by Dočev []. The iteration procedure () 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:

((Subscript[z, i])^((k+1)))=((Subscript[z, i])^((k)))-((P(((Subscript[z, i])^((k)))))/((∏_(j=1)^(i-1))(((Subscript[z, i])^((k)))-((Subscript[z, j])^((k+1))))(∏_(j=i+1)^n)(((Subscript[z, i])^((k)))-((Subscript[z, j])^((k)))))); i=1,…,n; k=0,1,….

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

Theorem 3—Factorization into linear terms

Any monic polynomial P of degree n≥1 in ℍ[x] admits a factorization into linear factors; that is, there exist (Subscript[x, 1]),… ,(Subscript[x, n]) such that

P(x)=(x-(Subscript[x, n]))… (x-(Subscript[x, 1])).

Theorem 4—Zeros from factors

Consider a polynomial P whose factor terms are (Subscript[x, 1]),… ,(Subscript[x, n]); that is, P admits a factorization of the form (). If the similarity classes [(Subscript[x, i])], i=1,… ,n, are distinct, then P has exactly n zeros (Subscript[ζ, 1]),… ,(Subscript[ζ, n]), which are given by:

(((Subscript[ζ, i])=(OverBar[ℛ])((Subscript[x, i]))(Subscript[x, i])ℛ((Subscript[x, i]))))^(-1)), i=1,…,n

(Subscript[ℛ, i])(x)=({{{, ({{1,, if i=1}, {(x-(Subscript[x, i-1]))…(x-(Subscript[x, 1])),, otherwise}})}})

Weierstrass Algorithm

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

Theorem 5 ([])

Let P be a polynomial in ℍ[x] of degree n with simple roots and, for i=1,…, n; k=0,1,2,… , let

((Subscript[z, i])^((k+1)))=((Subscript[z, i])^((k)))-((OverBar[((Subscript[ℒ, i])^((k)))])P (OverBar[((Subscript[ℛ, i])^((k)))]))(((Subscript[z, i])^((k))))(((((Subscript[𝒞, i])^((k)))((Subscript[z, i])^((k)))))^(-1)),

((Subscript[ℒ, i])^((k)))(x):=(x-((Subscript[z, n])^((k))))(x-((Subscript[z, n-1])^((k))))…(x-((Subscript[z, i+1])^((k)))), i=1,… ,n-1 and ((Subscript[ℒ, n])^((k)))(x):=1,

((Subscript[ℛ, i])^((k)))(x):=(x-((Subscript[z, i-1])^((k+1))))…(x-((Subscript[z, 2])^((k+1))))…(x-((Subscript[z, 1])^((k+1)))), i=2,… ,n and ((Subscript[ℛ, 1])^((k)))(x):=1

((Subscript[𝒞, i])^((k)))(x):=(∏_(j=1)^(i-1))(Subscript[Ψ, ((Subscript[z, j])^((k+1)))])(x)(∏_(j=i+1)^n)(Subscript[Ψ, ((Subscript[z, j])^((k)))])(x),

((Subscript[ζ, i])^((k+1)))=(OverBar[((Subscript[ℛ, i])^((k)))])(((Subscript[z, i])^((k+1))))((((Subscript[z, i])^((k+1)))(((Subscript[ℛ, i])^((k)))(((Subscript[z, i])^((k+1))))))^(-1)), k=1,…,n

The functions (OverBar[((Subscript[ℒ, i])^((k)))]), (OverBar[((Subscript[ℛ, i])^((k)))]) and ((Subscript[𝒞, i])^((k))) are implemented as the functions Polynomialℒ, Polynomialℛ and Polynomial𝒞, respectively. The support file QPolynomial associated with [] needs to be loaded.

<<QPolynomial`

Polynomialℒ[fact_,i_]:=Module[
{n=Length@fact},
Which[
i<n-1,(NonCommutativeMultiply@@(Polynomial[1,-Conjugate@#]&/@(Drop[fact,i]))),
i==n-1,Polynomial[1,-Conjugate@#]&@@(Drop[fact,i]),
i==n,Polynomial[1]
]
]

Polynomialℛ[fact_,i_]:=Module[
{n=Length@fact},
Which[
i==1,Polynomial[1],
i==2,Polynomial[1,-Conjugate@#]&@@(Drop[fact,i-n-1]),
i≤n,(NonCommutativeMultiply@@(Polynomial[1,-Conjugate@#]&/@(Drop[fact,i-n-1])))
]
]

Polynomial𝒞[factors:{_,_},i_]:=Eval[CharacteristicPolynomial@@Drop[factors,{i}]]@factors[[i]]

Polynomial𝒞[factors_,i_]:=Module[
{factori=factors[[i]]},NonCommutativeMultiply@@Map[Eval[#]@factori&,CharacteristicPolynomial/@Drop[factors,{i}]]
]

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

WeierstrassIteration[pol_,fact_,___]:=Module[
{n=Length@fact,factors=fact,roots={},RPoli,RPoliz},
For[i=1,i≤n,i++,
RPoli=Polynomialℛ[factors,i];
factors[[i]]=factors[[i]]-Eval[Polynomialℒ[factors,i]**pol**RPoli][factors[[i]]]**(1/Polynomial𝒞[factors,i]);
RPoliz=Eval[RPoli][factors[[i]]];
AppendTo[roots,RPoliz**factors[[i]]**(1/RPoliz)];
];
{pol,factors,roots}
]

The quaternionic Weierstrass iterative method is implemented in the function WeierstrassMethod.

TestConvergenceW[pol_,x1_,x2_,eps_]:=Max@(Abs[Norm/@(Last@x1)-Norm/@(Last@x2)])>eps||Max@(Abs[Re/@(Last@x1)-Re/@(Last@x2)])>eps || Max@(Norm/@Eval[pol,Last@x2])>eps;

WeierstrassMethod[pol_Polynomial,fact_,eps_:10.^-12,max_Integer:50]:=Module[{out,last,errorroots,errorfact},
If[Length@pol==Length@fact+1,
out=NestWhileList[WeierstrassIteration@@#&,{pol,fact},TestConvergenceW[pol,#1,#2,eps]&,2,max];
last=Take[Reverse@out,{1,2}];
errorroots=Max[Max@Abs[Re@(last[[1,3]])-Re/@(last[[2,3]])],Max@(Abs[Norm/@(last[[1,3]])-Norm/@(last[[2,3]])])];
{If[errorroots<eps,"Convergence","Divergence"],Length@out-1,N@errorroots,N@last[[1,3]],N@last[[1,2]]},
Message[WeierstrassMethod::Dimensions]
]
];

The usual convergence test (max_i)|((Subscript[ζ, i])^((k+1)))-((Subscript[ζ, i])^((k)))|<(Subscript[ε, 1]) has been replaced in TestConvergenceW by

(max_i)|Re(((Subscript[ζ, i])^((k+1))))-Re(((Subscript[ζ, i])^((k))))|<(Subscript[ε, 1]) and (max_i)| |((Subscript[ζ, i])^((k+1)))|-|((Subscript[ζ, i])^((k)))| |<(Subscript[ε, 1])

in order to let the function WeierstrassMethod recognize a sphere of zeros. Since we also include a test on the value of (max_i)|P(((Subscript[ζ, i])^((k+1))))|, 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 p(x)=(x^3)-j (x^2)- x+j  of Example 3, which we recall are (Subscript[r, 1])=1, (Subscript[r, 2])=-1 and (Subscript[r, 3])=j. All of the initial approximations ((Subscript[z, 1])^((0))), ((Subscript[z, 2])^((0))) and ((Subscript[z, 3])^((0))) have to lie in distinct congruence classes.

WeierstrassMethod[Polynomial[1,Quaternion[0,0,-1,0],-1,Quaternion[0,0,1,0]],{1.,2.,Quaternion[1.,0.,1.,0.]}]

Example 8

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

p(x)=(x+2i)(x+1+k)(x-2)(x-1)(x-2+j)(x-1+i)

(Subscript[ζ, 1])=1-i, (Subscript[ζ, 2])=2-(2/3)i-(1/3)j+(2/3)k, (Subscript[ζ, 3])=1,(Subscript[ζ, 4])=2, (Subscript[ζ, 5])=-1-(29/39)i+(14/39)j-(22/39)k, (Subscript[ζ, 6])=-(224/113)i-(30/113)k .

factors={Quaternion[1,-1,0,0],Quaternion[2,0,-1,0],1,2,Quaternion[-1,0,0,-1],Quaternion[0,-2,0,0]};
zeros=ZerosFromChain[factors]

WeierstrassMethod[PolynomialFromChain[factors],{Quaternion[0.5,0.,0.,0.],Quaternion[1.5,0.,-1.,0.],Quaternion[1.5,1.,-1.,1.],Quaternion[1.5,1.,-1.,0.],Quaternion[-0.5,0.,0.,0.],Quaternion[-1.,-2.,0.,0.]}]

N[zeros]

Example 9

The polynomial p(x)=(x^3)+(1+j)(x^2)+x+1+j has an isolated root 1+j and a sphere of zeros [i]. 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 ((Subscript[𝒬, i])^((k)))(((Subscript[z, i])^((k)))) in () become close to zero and some care is required.

WeierstrassMethod[Polynomial[1,Quaternion[1,0,1,0],1,Quaternion[1,0,1,0]],{Quaternion[-1.,1.,0,0],Quaternion[2,0,0,0],Quaternion[1,0,0,0]}]

out=WeierstrassMethod[Polynomial[1,Quaternion[1,0,1,0],1,Quaternion[1,0,1,0]],N[{Quaternion[-1,1,0,0],Quaternion[2,0,0,0],Quaternion[1,0,0,0]},200]]

The spherical root can be identified at once by observing that, up to the required precision, we have [(Subscript[ζ, 1])]=[(Subscript[ζ, 3])].

ζ=out[[5]];

Abs[Re[ζ[[1]]]-Re[ζ[[3]]]]

Abs[Norm[ζ[[1]]]-Norm[ζ[[3]]]]

Conclusion

This is the second article on several computational aspects of polynomials in the ring ℍ[x]. 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

I. Niven, “Equations in Quaternions,” The American Mathematical Monthly, 48(10), 1941 pp. 654–661. www.jstor.org/stable/2303304.

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.

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.

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.

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.

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.

D. Janovská and G. Opfer, “Computing Quaternionic Roots by Newton's Method,” Electronic Transactions on Numerical Analysis, 26, 2007 pp. 82–102.

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.

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.

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.

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.

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.

F. Miranda and M. I. Falcão. “QuaternionAnalysis Mathematica Package.” w3.math.uminho.pt/QuaternionAnalysis.

K. Gürlebeck, K. Habetha and W. Sprössig, Holomorphic Functions in the Plane and n-Dimensional Space, Basel: Birkhäuser, 2008. doi.org/10.1007/978-3-7643-8272-8.

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.

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.

E. Durand, Solutions numériques des équations algébriques. Tome I: Equations du type F(x); racines d’un polynôme, Paris: Masson,1960.

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.

I. O. Kerner, “Ein Gesamtschrittverfahren zur Berechnung der Nullstellen von Polynomen,” Numerische Mathematik, 8(3), 1966 pp. 290–294. doi.org/10.1007/BF02162564.

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.

The package QuaternionAnalysis.

Available at: w3.math.uminho.pt/QuaternionAnalysis

The file QPolynomial.m.

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