数理解析研究所講究録
1091 巻 1999 年 36-44

36

A new method to compute zeros of polynomials
武藤秀夫 (HIDEO MUTO) ’

鈴木智博 (TOMOHIRO SUZUKI)\dagger

鈴木俊夫 (TOSHIO SUZUKI)\ddagger

1

IntroductionLet$f(z)=z^{\mathrm{p}}+a_{1}z^{\mathrm{p}-1}+\cdots+a_{p},$ $a_{p}\neq 0$

be a polynomial with complex coefficients of order

distinct zeros, that is, $f(z)= \prod_{k=1}^{s}(Z-z_{k})^{n}k$ . The root finding problem for
$s$

$f$

$p$

which has

has along history and manykinds of methods are known. Among them the Newton’s method (or its stable version) is most popularly
used because of its rapidness and reliability. We can see thousands of papers in$\mathrm{M}\mathrm{c}\mathrm{N}\mathrm{a}\mathrm{m}\mathrm{e}\mathrm{e}’ \mathrm{S}$on-linebibliography available at http://www.elsevier.nl/local.cam, in which the general subject of polynomialroot finding is divided into 29 categories.

Putting$q(z)= \frac{f’(z)}{f(z)}$, the r-th derivative$q^{(r)}(z)= \sum_{k=1}^{s}\frac{1}{r!}\frac{(-1)’n_{k}}{(z-z_{\mathrm{k}})t+1}$

a formula to give an approximation of the nearest zero to
the Fourier coefficient$z$

and radius$|\tau|,$$T(z, \mathcal{T}, m;f)$$\mathrm{s}\mathrm{a}\mathrm{t}\mathrm{i}_{\mathrm{S}}q\mathrm{i}\mathrm{n}\mathrm{g}$

$z$

as

$rarrow\infty$

is analysed in Henrici ([3]) as. On the other hand we considerof $q(z)$ using $m$ points on a circle in the complex plane with center

the following representation formula:$T(z, \tau, m;f)=\sum_{k=1}^{l}\frac{n_{k}}{1-(^{z}\lrcorner\frac{-l}{\tau})^{m}}$, whichwill be proved in section 2. It is easy to see that $T(z, \tau, m;f)$ as well as $q^{(r)}(z)$ has information on thelocation of the nearest zero to$z$.

Our new method is considered to be derived from the error terms of the numerical integral of thefollowing integral:$\frac{1}{2\pi i}\oint_{\mathrm{r}^{\frac{f’(_{Z)}}{f(z)}d}}z$

where$\Gamma$$f$$\mathrm{m}$(1.1)

is a circle on the complex plane. It is well known that the integral (1.1) gives the number ofzeros of enclosed in . Taking
If,$\Gamma$$\mathrm{m}$points on , we use a modified form of the numerical integral of (1.1).
$\Gamma$

is small then it has a large error as a numerical integral. But we can show in this paper that it has acritical information on zeros of . The main feature of our method is in$f$$\mathrm{g}\dot{\mathrm{l}}\mathrm{o}\mathrm{b}\mathrm{a}\mathrm{l}$convergence which cannot

be realized in the simple Newton algorithm.

cDept. Mathematics and Physics,$\mathrm{Y}\mathrm{m}\mathrm{a}\mathrm{n}\mathrm{a}S\dot{\mathrm{h}}\mathrm{i}$Univ.

\dagger Dept. Computer Science, Yamanashi Univ.
$t_{\mathrm{D}\mathrm{e}\mathrm{p}}\iota$. Mathematics and Physics, Yamanashi Univ. -mail address: suzuki@grape.kkb.yamanashi.ac.jp

37

Because of the usage of the modified numerical integral, our theory m.ay look like similar to the papers($\mathrm{e}\mathrm{g}$. Delves&Lynes [2]) in the categor.y of ‘5. Integral method,

values of it are used in quite different ways; Their$\mathrm{m}\mathrm{e}\mathrm{t}\mathrm{h}\mathrm{o}\mathrm{d}_{8}$$\mathrm{e}\mathrm{s}\mathrm{p}.$Le.hmer’s’ at a glance. But theuse anyway the value of the numerical integral

while our method essentially needs the error of it. So, our new method dose notas accurate as
belong in any one of the 29 categories cited above. Among the 29 categories, only ‘2. Newton’s method’$\cdot$$\mathrm{p}\mathrm{o}8\mathrm{S}\mathrm{i}\mathrm{b}\mathrm{l}\mathrm{e}$has a relation with ours through the fact that, if the numerical integral in our method is carried out with$\mathrm{m}=1$, our method is essentially reduced to the Newton’s method.

In$8\dot{\mathrm{e}}\mathrm{c}\mathrm{t}\mathrm{i}\mathrm{o}\mathrm{n}2$, we study the basic theory, from which algorithms for finding zeros of$f$

are derived. We

can prove that our fundamental algorithm is generally convergent. Our algorithm for computing zeros
are stated in two ways in section 3 with some remarks. The results of numerical experiments are shownin section 4 with the fact that our method is comparable with the Newton’s method in its effectivity.
Concluding remarks are stated in section 5.

Basic theory

2
Let

$f$

be a polynomial of degree

a positive integer$\mathrm{m}$

, we

$p$

with zeros$8 \mathrm{e}\mathrm{t}\theta=\theta_{m}=\frac{2\pi}{m}$$z_{k},$and define$T( \lambda, \mathcal{T}, m;f)=\frac{\tau}{m}\sum_{j=0}^{m-1}\frac{f’(\lambda+\tau e^{1\rho}j)}{f(\lambda+re^{i\theta j})}.et\theta j$

Denoting$\nu_{k}=z_{k}-\lambda$

Lemma 2.1 If$k=1,$,$\cdots,$$s$, of multiplicity$T(\lambda, \tau, m;f)$$n_{k}$. For$\lambda\in \mathrm{C},$$\tau\in \mathrm{C}$andand $S(\lambda, \tau, m;f)$ as follows:.

$S( \lambda, \mathcal{T}, m;f)=\frac{\tau^{2}}{m}\sum_{i=0}^{m-1}\frac{f’(\lambda+\tau ej)i\theta}{f(\lambda+\tau ej)i\theta}e2\iota\theta j$, we have the following lemma.

$\tau^{m}\neq\nu_{k}^{m}$$\sum_{k=1}^{l}\frac{n_{k}}{1-(_{\tau}^{\nu}\alpha)^{m}}(z_{k}-\lambda)$

for all , then we have$k$$T( \lambda, \tau, m;f)=\sum_{k=1}^{l}\frac{n_{k}}{1-(\begin{array}{l}\perp\nu\tau\end{array})}$

.

Proof
$T( \lambda, \tau, m;f)=\frac{\tau}{m}m-1j=0\sum k=1\sum^{s}\frac{n_{k}e^{*\theta j}}{\lambda+\tau e^{|\theta}j-z_{k}}$$= \frac{r}{m}\sum_{k=1\mathrm{j}}\sum^{-1}\frac{n_{k}e^{\epsilon\theta j}}{-\nu_{k}+\tau e\iota\theta j}sm=0$

and$\frac{\tau}{m}\sum_{j=0}^{1}\frac{n_{k}e^{\iota\theta \mathrm{j}}}{-\nu_{k}+\tau e^{l\theta j}}=m-\frac{n_{k}}{m}\sum_{\mathrm{j}=0}^{m-1}\frac{\sum_{l_{-}0}^{m_{\vee}-}1(\begin{array}{l}-\nu\perp\tau\end{array})e^{-}i\theta jl}{(1-^{\underline{\nu}}\tau\iota_{e}-l\theta j)\sum^{m-1}l=0(\begin{array}{l}\lrcorner^{\nu}\tau\end{array})le^{-i}\theta jl}$$= \frac{n_{k}}{m}\sum_{0j=}^{1}\frac{\sum_{l_{-^{0}}^{-}}^{m-1}(_{\tau}\underline{\nu}_{\mathrm{A}})\iota_{6}-i\theta j/}{1-(_{\tau}^{\nu}\lrcorner)^{m}}m-$$= \frac{n_{k}}{1-(\begin{array}{l}\underline{\nu}_{\mathrm{A}}\tau\end{array})}$and $S(\lambda, \tau, m, f)=$

38

We can derive the second equality by the similar way. 1
Now we consider zer$T= \frac{n}{1-(\ _{r})m},+\delta$$o\mathrm{s}$

of

$f$

. Put

, we can solve$T=T(\lambda, \tau, m;f)$. Let$z_{1}= \lambda+\tau(\frac{\mathcal{T}-n}{T-}\frac{-\delta}{\delta})^{\perp}m$$z_{1}$

be the nearest zero to

So if the parametersproperly so that 6 is small enough, then we have an approximation formulaSince$\lambda,$$\tau$$\lambda$. Expressing that

and$\lambda+\tau(^{\underline{T}_{\frac{-n}{T}}})^{\perp}n$can be setfor$z_{1}$.

$T= \sum_{k=1}^{\delta}\frac{n_{k}}{1-(\begin{array}{l}\sim\nu_{\mathrm{A}}\tau\end{array})}=\frac{n_{1}}{1-(_{\tau}^{\underline{\nu}_{1}}\rangle m}\{1+\sum_{k=2}^{s}\frac{n_{k}}{n_{1}}\frac{1-(\begin{array}{l}\lrcorner\nu\tau\end{array})}{1-(\begin{array}{l}\Delta\nu\tau\end{array})}\}$$= \frac{n_{1}}{1-(\begin{array}{l}\lrcorner\nu\tau\end{array})}\{1+\sum_{k=2}^{\theta}\frac{n_{k}}{n_{1}}(\frac{\nu_{1}}{\nu_{k}})^{m}\frac{1-(\frac{\tau}{\nu_{1}})^{m}}{1-(\frac{\tau}{\nu_{k}})^{m}}\}$

putting$\mathrm{m}$$\Delta=\sum_{k=2}’\frac{n_{k}}{n_{1}}(\frac{\nu_{1}}{\nu_{k}})^{m}\frac{1-(\frac{\tau}{\nu_{1}}\rangle^{m}}{1-(\frac{\tau}{\nu_{k}})^{m}}$,, we have.

$\frac{T-n_{1}}{T}=(\frac{\nu_{1}}{\tau})^{m}\{$ $1+4( \frac{\tau}{\nu_{1}})^{m}\Delta\}\{1+\Delta\}-1$

Choosing a proper branch in the complex plane, the$\mathrm{f}\mathrm{o}\mathrm{U}\mathrm{o}\mathrm{W}\mathrm{i}\mathrm{n}\mathrm{g}$proposition is derived$\mathrm{f}\mathrm{o}\mathrm{r}\mathrm{m}\mathrm{a}\mathrm{u}_{\mathrm{y}}$from thisequality.
Proposition 2.1$\lambda+\tau(\frac{T-n_{1}}{T})^{m}-z_{1}\perp=\nu_{1}[(1+(\frac{\tau}{\nu_{1}})^{m}\Delta)^{m}\perp(\frac{1}{1+\Delta})^{\perp}n-1]$

holds.

Throughout this paper, let each zero of be numbered according to the distance from in the complex$f$plane, that is,$|\nu_{1}|\leq|\nu_{2}|\leq\cdots\leq|\nu_{s}|$$\lambda$. Here we have the following theorem.

Theorem 2.1 $If|\tau|<|\nu_{1}|<|\nu_{2}|$ then there exists a positive integer
inequality

If$| \lambda+\tau(\frac{T-n}{T}\iota)^{\frac{1}{n}}-z_{1}|\leq\sim_{m}164_{\frac{|\nu_{1}|}{1-1_{\nu}^{\mathrm{n}_{2}}1^{m}}}$$|\nu_{1}|<|\tau|<|\nu_{2}|$Proof It is easy to derive the inequality
such that$| \Delta|<\frac{1}{2}$that, for all $m>N$ , theholds.

then there exists a positive integer$| \lambda+\tau(\frac{T-n}{T})^{\frac{1}{m}}-z_{1}|\leq\underline{1}6m4_{\frac{|\nu_{1}|}{1-|\frac{r}{\nu_{2}}|^{m}}}|\frac{\tau^{2}}{\nu_{1}\nu_{2}}|^{m}$$Nsu\mathrm{c}h$$N$such that, for all $m>N$ , the inequality

holds, provided that$|\Delta|\leq p|_{\nu_{2}}\nu\lrcorner\downarrow m_{\frac{2}{1-|\frac{r}{\nu_{2}}|^{m}}}$)$| \frac{\tau^{2}}{\nu_{1}\nu_{\mathit{2}}}|<1$.

so we can choose an positive integer $N$

for all $m>N$ . First, we consider the case $|\tau[<|\nu_{1}|<|\nu_{2}|$ . Since$( \frac{1+(\frac{\tau}{\nu_{1}})^{m}\Delta}{1+\Delta})^{\frac{1}{m}}-1=\frac{1}{(1+\Delta)^{\perp}n}\{(1+(\frac{\tau}{\nu_{1}})m_{\Delta})^{\perp}m-(1+\Delta)m\}\perp$$= \frac{1}{(1+\Delta)^{\frac{1}{n}}}\{[1+\frac{1}{m}(\frac{\tau}{\nu_{1}})m\Delta+\frac{1}{2}\frac{1}{m}(\frac{1}{m}-1)(\frac{\tau}{\nu_{1}})^{2}m\Delta 2+\cdots]$$-[1+ \frac{1}{m}\Delta+\frac{1}{2}\frac{1}{m}(\frac{1}{m}-1)\Delta^{2}+\cdots]\}$$= \frac{-1}{(1+\Delta)^{\frac{1}{m}}}\{\frac{1}{m}(1-(\frac{\tau}{\nu_{1}})^{m})\Delta+\sim\frac{1}{m}(\frac{1}{m}-1)21(1-(\frac{\tau}{\nu_{1}})^{2}m)\Delta^{\mathrm{z}}+\cdots\}$

we have, for $m>N$ ,$| \lambda+\tau(\frac{T-n_{1}}{T})n-z_{1}|\perp=|\nu_{1}||(1+(\frac{\tau}{\nu_{1}})^{m_{\Delta}})^{n}\perp(\frac{1}{1+\Delta})^{n}-1|\perp$,

39$\leq|\nu_{1}|\frac{1}{(1-|\Delta|)^{\perp}m}\{\frac{1}{m}2|\Delta|+\frac{1}{m}2|\Delta|2+\cdots\}$$\leq|\nu_{1}|\frac{1}{(1-|\Delta|)^{1+\frac{1}{m}}}\frac{1}{m}2|\Delta|$$\leq\frac{16p}{m}\frac{|\nu_{1}|}{1-|\frac{\tau}{\nu_{2}}|^{m}}|\frac{\nu_{1}}{\nu_{2}}|^{m}$

As for the case$2| \frac{\tau}{\nu_{1}}|^{m}|\Delta|$$|\nu_{1}|<|\tau|<|\nu_{2}|$, considering that

and$|(1-( \frac{\tau}{\nu_{1}})^{m}\Delta|<$, we can prove the second inequality by the similar way for the first case. 1

Remark 2.1 Since the measure of the set$\lambda+\tau_{m}(\frac{T(\lambda\tau mmi\int 1)-n1}{T(\lambda\tau_{m},m,f1)},.)^{\frac{1}{n}}$Remark 2.2 Let

$C$

with

$|\tau_{m}|<|\nu_{1}|$$\{\lambda\in \mathrm{C}||\nu_{1}|=|\nu_{2}|\}$

for each

be a circle with center$\lambda$$\lim_{marrow\infty}\tau(\lambda,\tau,m;f)=\frac{1}{2\pi i}\oint_{C}\frac{f’(_{Z)}}{f(z)}d_{Z}$provided that no zero ofto the cases:
Put$\Delta=\sum_{k=2}’\frac{n_{k}}{n_{1}}(\frac{\tau}{\nu_{k}})^{m}\frac{(\frac{\nu_{1}}{\tau})^{m}-1}{1-(\frac{\tau}{\nu_{k}})^{m}}$$|\tau|<|\nu_{1}|$

$f$

or

is situated on$C$$|\nu_{1}|<|\tau|<|\nu_{2}|$$A=- \sum_{k\neq 1}\frac{n_{k}}{1-(\begin{array}{l}-\nu_{\mathrm{A}}\tau\end{array})}$,$m$} converges toand radiusand,is equal to zero, the$|\tau|$$z_{1}$

for almost all$8\mathrm{e}\mathrm{q}\mathrm{u}\mathrm{e}\mathrm{n}\mathrm{c}\mathrm{e}\{\lambda_{m}=$$\lambda\in \mathrm{C}$.

, then we have$\varliminf_{m\infty}S\mathrm{t}\lambda,$$\tau,m;f$)$= \frac{1}{2\pi i}\oint_{C}\frac{f’(Z)(Z-\lambda)}{f(z)}d_{Z}$. We should note that the first one is zero or$n_{1}$,

corresponding.

$B=- \sum_{k\neq 1}\frac{n_{k}\nu_{k}}{1-(\begin{array}{l}\lrcorner\nu\tau\end{array})}$and$C= \frac{1-(\begin{array}{l}\lrcorner\nu\tau\end{array})}{n_{1}}$. Then we have th. following$\mathrm{e}$proposition directly from Lemma 2.1.
Proposition 2.2$z_{1}-( \lambda+\frac{S}{T})=C\frac{\nu_{1}A-B}{1-CA}$holds.

Theorem 2.2 Assume that a polynomial $f(z)$ satisfies the following conditions (1)$\lambda|,$$| \tau|<r_{J}(\mathit{2})\min_{\mathrm{k}\neq 1}|z_{k}-\lambda|\geq R_{2}>r,$$\kappa=\max\{|\tau|,$ $|z_{1}-\lambda|1$(3) $R_{3} \geq\max_{k}|z_{k}-\lambda|,$ (4)$( \frac{f}{R_{2}})^{m}<\frac{1}{(2p-1)}$, we have$|z_{1}-( \lambda+\frac{S}{T})|\leq\frac{2p(R_{3}+\kappa)}{1-(2p-1)(\frac{\kappa}{R_{2}}\rangle^{m}}(\frac{\kappa}{R_{2}})^{m}$

Proof From the definitions of $A,$

$B$

and

$C$

we have$|A| \leq\sum_{k\neq 1}\frac{n_{k}}{1-|\frac{\tau}{z_{k}-\lambda}|^{m}}|\frac{\tau}{z_{k}-\lambda}|^{m}$$\leq\sum_{k\neq 1}\frac{n_{k}}{1-(_{R\mathrm{a}}^{\tau}\mathrm{u})m}(\frac{|\tau|}{R_{2}})m$$=(p-n1) \frac{1}{1-(_{R_{2}}r\mathrm{u})m}1\frac{|\tau|}{R_{2}})^{m}$.

,$|B| \leq\sum k\neq 1\frac{n_{k}|Z_{k}-\lambda|}{1-|\frac{\tau}{z_{k}-\lambda}|^{m}}|\frac{r}{z_{k}-\lambda}|m$$\leq\sum_{k\neq 1}\frac{n_{k}R_{3}}{1-(_{R_{2}}\mathrm{u}\tau)^{m}}(\frac{|\tau|}{R_{2}})m$$=(p-n_{1}) \frac{R_{3}}{1-(_{R_{2}}^{\tau}\mathrm{u}_{)}m}(\frac{|\tau|}{R_{2}})^{m}$.

$|C| \leq|\frac{\tau^{m}-\nu_{1}^{m}}{n_{1}\tau^{m}}|\leq 2\frac{\kappa^{m}}{|\tau|^{m}}$,$-(\mathit{4}):(\mathit{1})0<|z_{1}-$. Then, putiing

40

And
$|z_{1}-( \lambda+\frac{S}{T})|\leq \mathrm{I}C|\frac{|A\nu_{1}-B|}{1-|C||A|}$$\leq\frac{|\nu_{1}|+R_{3}}{1-|c|(p-n1)\frac{1}{1-(\mathrm{P}^{r_{2}})m}(R\mathrm{u}\tau)^{m}2}|C|(p-n1)\frac{1}{1-(_{R}^{\tau_{2}}\mathrm{u}_{)}m}(\frac{|\tau|}{R_{2}})m$$\leq\frac{2p(|\nu 1|+R3)}{1-(2p-1)(\frac{\kappa}{R_{2}})^{m}}(\frac{\kappa}{R_{2}})^{m}1$

Remark 2.3 The second case of Theorem 2.1 is not easy to check, so$S$$\lambda+\overline{T}$can be used as a comple-

mentary value for that cas .
$e$Algorithms for finding zeros of3In$\mathrm{t}\mathrm{h}\mathrm{i}_{8}$section we$f$algorithms to compute zeros of , based on the theory developed in the$f$$\mathrm{P}^{\mathrm{r}\mathrm{o}}\mathrm{P}^{\mathrm{o}8\mathrm{e}}$preceding section. We assume hereafter that the numerical computation in our algorithm is carried out

under double precision arithmetic of 16 decimal places. We set a value

considered to be a zero of

For almost all fixed$f$such that if $|f(z)|<\epsilon$ then

$z$

is

numerically.

$\lambda\in \mathrm{C}$, we can get$z_{1}$

by the following algorithm as is stated in Remark 2.1.

Algorithm 1: Fundamental algorithm for computing(step 1) Check that$\epsilon$$|f(\lambda)|>\epsilon$. Put$z_{1}$.

$R_{1}= \min(p|,\frac{f(\lambda)}{f(\lambda)}|, |f(\lambda)|^{1/p})$.

(step 2) Set initial parameters: $m=2,$ $\tau_{\max}=R_{1},$ $\tau_{\min}=0,$ $\tau=R_{1}/p$(s.tep 3) Compute $T=T(\lambda, \tau,m;f)$ .
(step 4) If $|T|\leq 10^{-5}$ then$\tau_{\min}=\tau,$$\tau=\ovalbox{\tt\small REJECT}_{2}\tau\tau$and go to (step 3).

If $|T|\geq.99$ then $\tau_{\max}=r,$ $\tau=\frac{\tau+\tau}{2}$ and go to (step 3).
(step 5) Put $n_{1}=1$ . Get$m\mathrm{v}\mathrm{a}\mathrm{l}\mathrm{U}\mathrm{e}8$$|f( \lambda+\tau x)|=\min_{k}|f(\lambda+\tau x_{k})|$$|f(Z)|\geq\epsilon$(step 7)

is a zero of . End.

Remark 3.1

$0$

as

$\{x_{1}, \cdots, xm\}$. Choose$x$from $\{x_{1}, \cdots, xm\}$ such that

. Put $z=\lambda+\tau x$ .

$f$is set as a value such that $|z_{1}-\lambda|<R_{1}$ (cf. [3], p454 Theorem 6. ).
$4\mathrm{e}$$R_{1}$&mark 3.2 For fixedffom$( \frac{T-n}{T})^{\frac{1}{m}}$then $m=m\cross 2$ and go to (step 3).

(step 6) If$z$ofto over$|\nu_{1}|,$$\lambda,$$|T(\tau)|$$m$

and $f,$ $T(\tau)=T(\lambda, \tau, m;f)$ is a continuous function of . Whendose from$\tau$$0$

to over$n_{1}$$|\tau|$variescontinuously as is seen fiom Remark 2.2. So (step 4)

can be passed after at most several times returns to (step 3). Practically$|\tau|$can be any value as far as

it can pass the criterion in (step 4).
Remark 3.3 If
be better thatvalues of$z_{1}$$n_{1}$

may be of multiplicity greater than 1, for example, of multiplicity 3, then it wouldof (step 5) is set to range over three cases: $n_{1}=1,2,3$ and that$( \underline{\tau}_{\overline{T}^{n_{1}}})\frac{1}{m}$, $n_{1}=1,2,3$ .

$x$is chosen from $3m$

41

Remark 3.4 Since the lessis the more effective this algorithm is, we would be able to get the$|_{\nu_{2}}^{\nu}\lrcorner|$more practical algorithm by replacing$\lambda$with a better approximation. So we propose the following morepractical algorithm, though the convergence of which is not proved yet.

Algorithm 2: Applied algorithm which converge to a zero of(step 1) Take A$\tau\max\sim-R_{1},$C. Check that$\in$$\mathcal{T}_{\min}=0$$|f(\lambda)|>\epsilon$

. Let

$f$

.

$R_{1}= \min(p|,\frac{f(\lambda)}{f(\lambda)}|, |f(\lambda)|^{1/\mathrm{p}}))m=5,$$L=1$ ,

and $\tau=R_{1}/p$ .

tsteo 2) Comoute $T=T(\lambda.\tau.m:f)$ .
it$|\mathit{1}|\geq.9\mathrm{W}$tnen $\tau_{\max}=\tau,$$\tau=\frac{\mathrm{T}’}{2}\mathrm{a}\mathrm{n}a$(step 4) Get the $3m$ values ofsuch that$|f( \lambda+\tau x)|=\min|f(\lambda+\tau x_{k}^{n_{1}})|$

If $|f(\lambda+\tau x)|>|f(\lambda)|$ then let$\mathrm{I}\mathrm{f}|f(\lambda)|<\epsilon$

else letas$(^{\underline{T}-}\tau^{n}arrow)^{\frac{1}{m}}$go to$(\mathrm{s}\mathrm{t}\mathrm{e}_{\mathrm{P}^{l}})$$\{x_{1}^{n_{1}}, \cdots , x_{m}^{n_{1}}\},$.

$n_{1}=1,2,3$ .

Choose$x$

from$\{X_{1}, \cdots, X_{m^{1}}\}n_{1n}$.
go , (step 2) else let$\mathrm{t}\mathrm{o}$$m_{\overline{\sim}}2m\mathrm{a}\mathrm{n}..\mathrm{d}$$\lambda=\lambda+\tau x$.

then stop$R_{1}= \min(p.

|,\frac{f(\lambda)}{f(\lambda)}| , |f(\lambda)|^{1}/p),$$\mathcal{T}maG=R_{1},$ $\tau_{\min}=0$

if

Let $m=5+M$ ,or let $m=3+M$ , ifor let$m=1$ ,

ifand $\tau=R_{1}/p$ .

$\tau>10^{-2}$$10^{-2}\geq\tau>10^{-9}$$10^{-9}\geq\tau$.

(step 5) Go to (step 2).

Remark 3.5 We$\mathrm{s}\mathrm{h}_{0}\mathrm{u}\mathrm{l}\mathrm{d}\backslash$remark that$\tau$

may be decided easily, satisfying

small. In fact for polynomials of lower order, we can set
to determine the parameters $m$ and$\tau$$\tau$$|\tau|<|\nu_{1}|$more roughly, say$\mathrm{m}$$\mathrm{g}\mathrm{r}e$. If the multiplicity or the number of close zeros is guessed to beeffective to set $n_{1}=1,2,$$\cdots$, . If$t_{0}$$x$. The criterion

larger than this case for polynomials of higher order.

Remark 3.6 We assume here that the multiplicities of zeros are rarely$n_{1}\geq 2$$|\tau|=|f(\lambda)|$in (step 3) is suitable for the polynomials of order less than 50. As

is seen from Theorem 2.1, we had better to take$n_{1}=1,2,3$, as far as $|T|$ is not too

is chosen from$\{x_{1}^{1}, \cdots , x_{m}^{1}\}$ater than 3. So we set$t_{0}$

then it will be morethen it does not need to compute forin the following iterations.

Remark 3.7 The case of $m=1$ is essentially coincide with the Newton’s method.

Remark 3.8 When and$\tau$$m$are fixed, this algorithm is not generally convergent as is proved in[4] that there is no generally convergent purely iterative algorithm, rational overpolynomials of degree$\geq 4$$\mathrm{C}$$\mathrm{M}\mathrm{c}\mathrm{M}\mathrm{u}\mathrm{l}\mathrm{l}\mathrm{e}\mathrm{n}$, for finding roots of. There are examples which does not converge to a zero of . In fact Professor$f$each of which guarantees the existence of

Sugiura of Nagoya Univ. gave examples of triples$(f(z), m, \mathcal{T})$a regionthen the iteration of this aJgorithm- for fixed$D$such that if the$\mathrm{i}\mathrm{n}_{l}\mathrm{i}\mathrm{t}\mathrm{i}\mathrm{a}1\lambda$

is taken in$D$gives a sequence that approaches to two different values periodically.

$m$

and

$\tau$

42

Remark 3.9 When a cluster of very close zeros apart from other zeros exists, we can easily see from
Lemma 2.1 that$\lambda+\frac{S}{T}$circle with center
value of$T$$\lambda$

yields a certain weighted mean of zeros contained in the cluster by taking aand radius

which enclose the cluster and exclude all the other zeros. Then the$|\tau|$gives the number of zeros contained in it. Moreover considering the weighted mean like amultiple zero

$z_{1}$, the estimate in Theor$e\mathrm{m}2.2$

can be applied to this case. That is, a set of very close

zeros is treated as a multiple zero in the formula of

4$8\mathrm{m}\mathrm{a}\mathrm{U}$$\lambda+\frac{S}{T}$.

Results of numerical experiments

Numerical experiments were carried out for many examples. Algorithm 1 was testified to work$\mathrm{w}\mathrm{e}\mathrm{l}\mathrm{l}|$Thatis, Algorithm 1 is generally convergent and the errors of test computations are theoretically reasonablevalues, except machine$\mathrm{e}\mathrm{r}\mathrm{r}\mathrm{o}\mathrm{r}8$. So we show only a$\mathrm{f},\mathrm{e}\mathrm{w}$results which characterize how Algorithm 1 works.

The polynomial for the experiment to test Algorithm 1 is $f(z)=z^{3}-3z+3$ whose zeros are $z_{1}=$$1.051\cdots+0.565\cdots i,$ $z_{2}=1.0’ 51\cdots-0.565\cdots i$bound stated in Theorem 2.1 by$\lambda$

We can seee in Table 1-2 that if$\mathrm{v}e$$E= \frac{16}{m}\frac{|\nu_{1}|}{1-1_{\dot{\nu}_{2}}^{L}|^{n}}$. Let$z= \lambda+\tau(\frac{T-1}{T})^{\perp}m$.

$f[_{Z})=\mathrm{z}^{3}-.\mathrm{q}_{z}\neq.\mathrm{q}$is a good approximation of zero then Algorithm 1 converges ven for small$r\mathrm{p}_{B*}\mathrm{h}\iota‘*1_{-}?arrow$converges, i.e.

hv Alnnrit.hm 1

Tahlr 1-1

Table 1-1 shows that if$E$

and $z_{3}=-2.103\cdots+0i$ . We denote the theoretical error

$e$

hv

$\mathrm{A}1\sigma\cap 1\mathrm{i}\mathrm{t}.\mathrm{h}\mathfrak{m}1$$\frac{|z_{1}-\lambda 1}{\mathrm{I}^{z_{2^{-}}\lambda}1}\sim 1$

then, though$f[,)=r^{3}.-.\mathrm{q}_{r}$.

$\mathrm{i}\mathrm{t}$$\mathrm{m}$.

$\lrcorner_{-}.\mathrm{q}$needs a large amount of computation, it

ry slowly.

As for the examples for Algorithm 2, we tried some polynomials of order higher than 100 withsatisfactory results, added to many ones of lower order. Here we show some results which illustrate
how the convergent process of Algorithm 2 is. In these experiments roughly determined$\tau$

such that43$\}\tau|=O(|f(\lambda)|)$is used. Table 2-1 shows that Algorithm 2 converges even if$|\nu_{1}|=|\nu_{2}|$. Table 2-2 is theresult for $f(z)=z^{20}+1$ which has simple zeros distributed on the unit circle. Table 2-3 is the result fora polynomial which has double zeros. We set the range of $n_{1}\{1,2,3\}$ in these experiments. It is easily
seen from our theorems that if$|z_{1}-\dot{\lambda}|<|\tau|$then

we added to the result of Algorithm 2 the value$\lambda_{6}$$\lambda+\frac{s}{T}$approximates$|f( \lambda+\frac{s}{T})|$

better than$\lambda+r(\frac{T-n}{T}\mathrm{L})^{\frac{1}{m}}$. So

in the following tabl es to compare two values.

in Table 2-2 is eqaul to the exact zero to the full digits. Note that$\mathrm{T}\mathrm{a}\mathrm{b}1\tilde{e}2,2$$z_{1}$$\lambda_{2}+\#$in Table 2-1 and’$\lambda_{4}+\dot{F}$incan be considered to be the numerical solutions.

Remark 4.1 Since$\lambda_{k-1}+\frac{s}{T}$can be computed cheaply by adding one process to Algorithm 2, we can

expect that such a device combining two algorithms will give another more efficient algorithm.

44

Table 2-3$‘ \mathrm{s}\mathrm{h}_{0}\mathrm{w}8$that in the case of zeros of multiplicity greater than 1, our method has the advantage of

rapidity and of accuracy over Newton’s method.

5

Concluding remarks&mark 5.1 The feature of our algorithm 1 and 2$\mathrm{i}_{8}$in the usage of the error of the numerical integralof (1.1). Our root finding method is a new one which we cannnot find in preceding papers cited in$\mathrm{M}\mathrm{c}\mathrm{N}\mathrm{a}\mathrm{m}\mathrm{e}\mathrm{e}’ \mathrm{S}$biblography.

Remark 5.2 Our fundamental algorithm is proved to be globally convergent for almost all initial values,
while our applied algorithm is not proved, yet. But we have the conjecture that Algorithm 2 is globally
convergent from the theoretical situation and the results of the numerical.

$\mathrm{e}\mathrm{x}\mathrm{p}\mathrm{e}\mathrm{r}\mathrm{i}\mathrm{m}\mathrm{e}\mathrm{n}\mathrm{t}8$Remark 5.3 Comparing with the Newton’s metho,’ Algorithm 2 is comparable in rapidity and is equal
or better in accuracy, especially for zeros of multiplicity greater than 1.
$\mathrm{d}$Remark 5.4 Connecting Algorithm 2 and the value$\lambda+\frac{S}{T}$, the more rapid and reliable algorithm will

be realized. Comparing the two kinds of estimates in twc theorems, we can expect to find the method todistinguish the multiple zeros from the close ones.

References[1] Atkinson, Kendall E. An introduction to NUMERICAL ANALYSIS, 2nd ed., John Wiley &Sons,
New York (1987).
[2] Delves, L. M. &Lyness, J. N., A Numerical Method for Locating the Zeros of an Analytic Function,

Math. Comput., 21 (1967) 543-560.
[3] Henrici, P., Applied and computational complex analysis Vol l, Wiley$\mathrm{I}\mathrm{n}\mathrm{t}\mathrm{e}\mathrm{r}\mathrm{S}\mathrm{C}\mathrm{i}\mathrm{e}^{r}\mathrm{n}$ce, New York (1974).

[4] McMullen, C., Families of rational maps and iterative root-finding algorithms, Ann. Math. 125 (1987)

467-793.
[5] Noda, M. &Sasaki, T., Approximated GCD and its application to ill-conditioned algebraic equations,J. Comp. Appl. Math. 38 (1991) 335-351.