DETERMINATION OF THE ZEROS AND POLES OF A MINIMUM-PHASE TRANSFER FUNCTION BY THE METHOD OF MOMENT SEQUENCES

V. M. OSIPOV

(Presented by the scientific-methodological seminar of the Department of Automation and Telemechanics)

Let the transfer function of a certain linear dynamic system be given by

\[ W(p)=\frac{H(p)}{D(p)}, \tag{1} \]

where

\[ \begin{aligned} H(p)&=a_{m+1}p^m+a_m p^{m-1}+\ldots+a_2p+a_1,\\ D(p)&=p^n+\beta_n p^{n-1}+\ldots+\beta_2p+\beta_1, \end{aligned} \tag{2} \]

with \(m<n\). We shall assume that the zeros and poles of the function \(W(p)\), that is, the roots of the equations,

\[ H(p)=0 \tag{3} \]

and

\[ D(p)=0, \tag{4} \]

are simple and lie in the left half-plane. In other words, the transfer function (1) is minimum-phase. The corresponding weighting function will evidently have the form

\[ h(t)=\sum_{i=1}^{n} c_i e^{p_i t}, \tag{5} \]

where \(c_i\) is the residue of the function \(W(p)\) at the simple pole \(p=p_i\), that is,

\[ c_i=\operatorname*{Res}_{p=p_i}\frac{H(p)}{D(p)} =\frac{H(p_i)}{D'(p_i)},\qquad (i=1,2\ldots). \]

The response of the system to the input \(F(t)\) can be found by means of the convolution integral

\[ x(t)=\int_{0}^{t} h(t-\tau)F(\tau)\,d\tau =\sum_{i=1}^{n} c_i \int_{0}^{t} e^{p_i(t-\tau)}F(\tau)\,d\tau . \]

Let us consider a numerical method for determining the zeros and poles of the transfer function (1), based on constructing a sequence of real numbers having the meaning of moments of the weighting function (5). The problem reduces to determining the roots of equations (3) and (4).

§ 1. Determination of the root of smallest modulus

Number the roots of equation (4) in order of increasing moduli,

\[ |p_1| \leq |p_2| \leq |p_3| \leq \ldots \leq |p_n| \tag{1-1} \]

and consider the sequence of moments of the weight function

\[ m_k=\int_{0}^{\infty} t^k h(t)\,dt,\quad (k=0,1,2\ldots). \tag{1-2} \]

Taking (5) into account, we obtain

\[ m_k=(-1)^k \sum_{i=1}^{n} c_i\,\frac{k!}{p_i^{k+1}} . \tag{1-3} \]

Suppose that the root of smallest modulus is real. According to (1-1), this root is \(p_1\). Obviously, as \(k\) increases, the first term in (1-3) will increasingly dominate the others and, for sufficiently large \(k\), one may approximately take

\[ \frac{m_k}{k!}\approx \frac{c_1}{p_1^{k+1}}(-1)^k, \]

and also

\[ \frac{m_{k+1}}{(k+1)!}\approx \frac{c_1}{p_1^{k+2}}(-1)^{k+1}. \]

Dividing one by the other, we obtain

\[ |p_1|\approx (k+1)\frac{|m_k|}{|m_{k+1}|}. \tag{1-4} \]

The larger \(k\) is, the more accurate formula (1-4) is, since in the limit

\[ |p_1|=\lim_{k\to\infty}(k+1)\frac{|m_k|}{|m_{k+1}|}. \]

Now suppose that the roots \(p_1\) and \(p_2\) form a complex-conjugate pair with the smallest modulus. Introduce the notation

\[ \left. \begin{aligned} p_1&=-\alpha_1+j\omega_1=\sqrt{\alpha_1^2+\omega_1^2}\,e^{j\varphi_1}=|p_1|e^{j\varphi_1},\\ p_2&=-\alpha_1-j\omega_1=\sqrt{\alpha_1^2+\omega_1^2}\,e^{-j\varphi_1}=|p_1|e^{-j\varphi_1},\\ \varphi_1&=\arccos\frac{-\alpha}{|p_1|}, \end{aligned} \right\} \tag{1-5} \]

For sufficiently large \(k\), we obtain the following approximate relations \([(-1)^k\) may be omitted]:

\[ \left. \begin{aligned} \frac{m_{k-1}}{(k-1)!}&\approx \frac{2}{|p_1|^k}\left[c_1' \cos k\varphi_1+c_1''\sin k\varphi_1\right],\\ \frac{m_k}{k!}&\approx \frac{2}{|p_1|^{k+1}}\left[c_1' \cos (k+1)\varphi_1+c_1''\sin (k+1)\varphi_1\right],\\ \frac{m_{k+1}}{(k+1)!}&\approx \frac{2}{|p_1|^{k+2}}\left[c_1' \cos (k+2)\varphi_1+c_1''\sin (k+2)\varphi_1\right], \end{aligned} \right\} \tag{1-6} \]

where

\[ c_1'=\operatorname{Re}c_1=\operatorname{Re}c_2;\qquad c_1''=\operatorname{Im}c_1=-\operatorname{Im}c_2 \]

\[ c_1=c_1'+jc_1'' \quad \text{and} \quad c_2=c_1'-jc_1''. \]

If we make use of the trigonometric identities for the sine and cosine of the sum of two angles, we obtain the possibility of eliminating from relations (1-6) expressions of the form

\[ c'_1 \cos k\varphi_1 + c''_1 \sin k\varphi_1 \quad \text{and} \quad c''_1 \cos k\varphi_1 - c'_1 \sin k\varphi_1 . \]

As a result, after simplification, we shall have

\[ \frac{m_{k+1}}{k(k+1)} + \frac{m_k}{k}\frac{2|\alpha_1|}{|p_1|^2} + m_{k-1}\frac{1}{|p_1|^2} \approx 0, \]

and also

\[ \frac{m_{k+2}}{(k+1)(k+2)} + \frac{m_{k+1}}{k+1}\frac{2|\alpha_1|}{|p_1|^2} + m_k\frac{1}{|p_1|^2} \approx 0. \]

Solving jointly, we obtain

\[ \left. \begin{aligned} \alpha_1 &\approx \frac{k+1}{2} \frac{\left|k m_{k+2}m_{k-1}-(k+2)m_{k+1}m_k\right|} {\left|(k+2)m_{k+1}^2-(k+1)m_k m_{k+2}\right|},\\[6pt] |p_1|^2 &= \alpha_1^2+\omega_1^2 \approx \frac{(k+1)(k+2)\left|(k+1)m_k^2-km_{k+1}m_{k-1}\right|} {\left|(k+2)m_{k+1}^2-(k+1)m_k m_{k+2}\right|}. \end{aligned} \right\} \tag{1-7} \]

Thus, complex-conjugate roots with the smallest modulus are also expressed in terms of the moments of the weight function. Let us consider a procedure for determining the moments. Since the roots of equation (4) do not depend at all on the form of the polynomial \(H(p)\), we may set \(H(p)=1\) and consider a certain calculated transfer function

\[ W_1(p)=\frac{1}{D(p)}, \tag{1-8} \]

having the same poles as the transfer function (1). The transfer function (1-8) corresponds to the time original \(y(t)\), that is

\[ W_1(p)=\int_{0}^{\infty} y(t)e^{-pt}\,dt. \tag{1-9} \]

Differentiate (1-9) \(k\) times with respect to \(p\) and set \(p=0\); then, taking (1-3) into account, we may write

\[ \left. \frac{d^k W_1(p)}{dp^k} \right|_{p=0} = (-1)^k \int_{0}^{\infty} t^k y(t)\,dt = (-1)^k m_k, \]

or

\[ m_k=(-1)^k W_1^{(k)}(0). \tag{1-10} \]

Let us note that we have retained the previous notation for the moments of the function \(y(t)\). It is convenient for us to introduce a new characteristic \(M_k\) according to the formula

\[ W_1^{(k)}(0)=(-1)^k k! M_k, \tag{1-11} \]

then (1-10) can be written in the form

\[ m_k=k!M_k. \tag{1-12} \]

From (1-8) it follows that

\[ W_1(p)\cdot D(p)=1. \]

Applying Leibniz’s rule, we obtain

\[ \left. \frac{d^k}{dp^k}\left[W_1(p)\cdot D(p)\right] \right|_{p=0} = \sum_{\nu=0}^{k} \binom{k}{\nu} W_1^{(k-\nu)}(0)D^{(\nu)}(0) = \begin{cases} 1, & k=0,\\ 0, & k\ne 0, \end{cases} \tag{1-13} \]

If we take into account (1-11), as well as the equalities

\[ \binom{\nu}{k}=C_k^\nu=\frac{k!}{\nu!(k-\nu)!}, \]

\[ D^{(\nu)}(0)=\nu!\beta_{\nu+1}, \]

then we shall have the system

\[ \sum_{\nu=0}^{\kappa}(-1)^{\kappa-\nu}M_{\kappa-\nu}\beta_{\nu+1} = \begin{cases} 1 & \kappa=0,\\ 0 & \kappa=1,2,3,\ldots \end{cases} \tag{1-14} \]

or, in expanded form,

\[ \left. \begin{aligned} M_0\beta_1&=1\\ -M_1\beta_1+M_0\beta_2&=0\\ M_2\beta_1-M_1\beta_2+M_0\beta_3&=0\\ -M_3\beta_1+M_2\beta_2-M_1\beta_3+M_0\beta_4&=0\\ &\text{— — — — — — — — —} \end{aligned} \right\} \tag{1-15} \]

The system (1-14) is easily solved. Successively we can find

\[ \left. \begin{aligned} M_0&=\frac{1}{\beta_1},\\[4pt] M_1&=\frac{\beta_2}{\beta_1^2}=\frac{1}{\beta_1}M_0\beta_2,\\[4pt] M_2&=\frac{\beta_2^2}{\beta_1^3}-\frac{\beta_3}{\beta_1^2} =-\frac{1}{\beta_1}(M_0\beta_3-M_2\beta_2),\\[4pt] M_3&=\frac{1}{\beta_1}(M_0\beta_4-M_1\beta_3+M_0\beta_2),\\ &\text{— — — — — — — — — —} \end{aligned} \right\} \tag{1-16} \]

\[ M_\kappa=(-1)^{\kappa+1}\frac{1}{\beta_1} \sum_{\nu=0}^{\kappa-1}(-1)^\nu M_\nu\beta_{\kappa-\nu+1} \quad(\kappa=1,2,\ldots). \tag{1-17} \]

Let us note that for \(\kappa>n\) (\(n\) is the degree of the polynomial \(D(p)\)) the first terms in (1-17), for which \(\nu<\kappa-n\), vanish. In addition, it should be borne in mind that \(\beta_{n+1}=1\).

We have obtained a practically very convenient algorithm for the successive determination of the quantities \(M_k\), by which the moments \(m_k\) can easily be found. The latter, however, need not be computed, since the calculation formulas (1-4) and (1-7), as a result of substitution of (1-12), acquire the more convenient form:

\[ \left. \begin{aligned} |p_1|&\simeq(\kappa+1)^{-1}\frac{|m_\kappa|}{|m_{\kappa+1}|} =\frac{m_\kappa}{m_{\kappa+1}},\\[6pt] \alpha_1&\simeq \frac{1}{2}\, \frac{|M_{\kappa+2}\cdot M_{\kappa-1}-M_{\kappa+1}\cdot M_\kappa|} {|M_{\kappa+1}^2-M_\kappa\cdot M_{\kappa+2}|},\\[6pt] |p_1|^2&=\alpha_1^2+\omega_1^2\simeq \frac{|M_\kappa^2-M_{\kappa+1}\cdot M_{\kappa-1}|} {|M_{\kappa+1}^2-M_\kappa\cdot M_{\kappa+2}|}. \end{aligned} \right\} \tag{1-18} \]

Thus, the sequence of real quantities \(M_\kappa\) \((\kappa=0,1,2,\ldots)\), constructed by the recurrence formula (1-17), solves the problem of determining the simple root of smallest modulus as

real as well as complex. But before applying formula (1-18), it is necessary to know the type of the root. It turns out that the character of the sequence \(M_k\) \((k=0,1,2,\ldots)\) makes it possible to judge this. It can be rigorously shown that if the indicated sequence changes monotonically without a change of sign, then the root of smallest modulus is real. For a complex root the sequence has an oscillatory character, and sometimes the sign changes quite irregularly. One more question arises: what value of \(k\) should be taken in formulas (1-18)?

An approximate value of the root begins to stand out already at \(k=3\text{--}4\). As \(k\) increases, the accuracy grows if, of course, the subsequent terms are computed with high precision, for example on an arithmometer. In calculations on an ordinary slide rule it is inadvisable to take a value of \(k\) greater than \(9\text{--}10\), since the accumulated rounding errors do not make it possible to find a more accurate value of the root.

In practice one should proceed as follows. For three successive values of \(k\) (for example, for \(k=4,5,6\) or \(k=5,6,7\)) three values of the root are determined. If they differ little from one another and the difference between them (in the real and imaginary parts) changes sign, then one should take the value corresponding to the larger value of \(k\), or their arithmetic mean. If, however, the difference between them is large and of the same sign, then this is an indication of a multiple root. Its approximate value can be determined by further increasing \(k\).

In a completely analogous way one can find the root of equation (3) that is smallest in modulus, that is, the zero of the transfer function (1) that is smallest in modulus. To construct the sequence \(M_k\) \((k=1,2,\ldots)\) in this case, it is necessary in the recurrence formula (1-17) to replace the coefficients \(\beta\) by the corresponding coefficients \(\alpha\), and to take into account that \(\alpha_{m+1}\ne 1\).

§ 2. Determination of the root of largest modulus and of the other roots

In equation (4) we make the substitution

\[ p=\frac{1}{z} \tag{2-1} \]

and reduce it to the form

\[ D_1(z)=\beta_1 z^n+\beta_2 z^{n-1}+\ldots+\beta_{n-1}z^2+\beta_n z+1=0. \tag{2-2} \]

Let \(z_1\) be the root of smallest modulus of this equation; then the root \(p_n\) of largest modulus of the original equation (4) will obviously be equal to

\[ p_n=\frac{1}{z_1}. \tag{2-3} \]

Thus, the problem is reduced to determining the root \(z_1\) of smallest modulus of equation (2-2), which can be done by the method already considered.

Now we must keep in mind that the order in which the coefficients of equation (2-2) occur is reversed. In place of the coefficient \(\beta_m\) in equation (4) there now stands the coefficient \(\beta_{n-m+2}\); therefore in the recurrence formula (1-17) we must make the substitution \(\beta_{\kappa-\nu+1}\to\beta_{n-\kappa+\nu+1}\), as a result of which we shall have

\[ N_k=(-1)^{k+1}\sum_{\nu=0}^{k-1}(-1)^\nu N_\nu \beta_{n-k+\nu+1}\quad (k=1,2,\ldots). \tag{2-4} \]

To avoid confusion, instead of \(M_k\) we introduce the notation \(N_k\). Let us write (2-4) in expanded form:

\[ \begin{aligned} N_0&=1=\beta_{n+1}\\ N_1&=N_0\beta_n=\beta_n\\ N_2&=-N_0\beta_{n-1}+N_1\beta_n=\beta_n^2-\beta_{n-1}\\ N_3&=N_0\beta_{n-2}-N_1\beta_{n-1}+N_2\beta_n\\ &\quad -\ -\ -\ -\ -\ -\ -\ -\ - . \end{aligned} \tag{2-5} \]

If \(z_1\) is a real root, then

\[ |p_n|=\frac{1}{|z_1|}\simeq \frac{N_{k+1}}{N_k}. \tag{2-6} \]

If, however, \(z_1\) and \(z_2\) form a complex-conjugate pair \((z_1=-a_1+jb_1\) and \(z_2=-a_1-jb_1)\) with the smallest modulus, then the modulus of the root \(p_n\) and its real part are determined by the formulas:

\[ \begin{aligned} \sigma_n&=\frac{a_1}{|z_1|^2}\simeq \frac{1}{2}\, \frac{|N_{k+2}N_{k-1}-N_{k+1}N_k|} {|N_k^2-N_{k+1}N_{k-1}|},\\ |p_n|^2&=\sigma_n^2+\omega_n^2= \frac{1}{a_1^2+b_1^2}\simeq \frac{|N_{k+1}^2-N_kN_{k+2}|} {|N_k^2-N_{k+1}N_{k-1}|}. \end{aligned} \tag{2-7} \]

Thus, by constructing with the aid of the recurrence formula (2-4) the sequence of quantities \(N_k\) \((k=0,1,2,\ldots)\), we can find the root \(p_n\) of the original equation (4) that is largest in modulus; moreover, all the remarks made earlier remain fully valid.

If one of the roots of an algebraic equation of degree \(n\) has been determined, the problem reduces to solving an equation of degree \(n-1\). Therefore, from the standpoint of the labor involved in determining the remaining roots, the most favorable case will be the case in which the initially found root is complex. In this case we can lower the degree of the equation by two units, since for real coefficients of the equation complex roots always form conjugate pairs. If, in addition, the root largest in modulus also turns out to be complex, then the degree of the equation can be lowered by four units. The case of real roots is the least favorable.

Having constructed an equation of reduced degree, we can again apply our algorithm and find part of the roots; then, if necessary, construct an equation of still lower degree and again find roots by the method described, and so on. In this way we can determine all the roots. Naturally, the labor involved increases as the degree of the equation grows; moreover, the accuracy of determining the subsequent roots will depend on the accuracy of the roots found initially.

To implement the method described, it is expedient to use Newton’s formulas, which relate the coefficients of the equation to the power sums of the roots [1, 2]. Denote

\[ S_k=p_1^k+p_2^k+\ldots+p_n^k=\sum_{i=1}^{n}p_i^k \quad (k=1,2\ldots), \tag{2-8} \]

then

\[ S_k=-k\beta_{n-k+1}-\sum_{i=1}^{k-1}\beta_{n-i+1}S_{k-i} \quad (k=1,2\ldots). \tag{2-9} \]

The recurrence formulas (2-9) are also called Newton’s formulas. For the first several values of \(k\) we obtain

\[ \left. \begin{aligned} S_1 &= -\beta_n,\\ S_2 &= -2\beta_{n-1}-\beta_n S_1=-2\beta_{n-1}+\beta_n^2,\\ S_3 &= -3\beta_{n-2}-\beta_{n-1}S_1-\beta_n S_2 =3(\beta_n\beta_{n-1}-\beta_{n-2})-\beta_n^3,\\ S_4 &= -4\beta_{n-3}-\beta_{n-2}S_1-\beta_{n-1}S_2-\beta_n S_3,\\ &\quad -\ -\ -\ -\ -\ -\ -\ -\ -\ -\ -\ -\ -\ -\ -\ -\ - . \end{aligned} \right\} \tag{2-10} \]

In addition to (2-8), let us also note the formula relating the product of the roots and the constant term of equation (4) (Vieta’s theorem)

\[ p_1\cdot p_2\cdots p_n=\prod_{i=1}^{n}p_i=(-1)^n\beta_1 . \tag{2-11} \]

We shall now consider a procedure for determining the roots of equation (4) for \(n=3,4,5,6\).

1. Cubic equation

Two cases are possible:

a) The root of smallest modulus forms a complex-conjugate pair \((p_1=-\alpha_1+j\omega_1\) and \(p_2=-\alpha_1-j\omega_1)\).

From (2-8) it follows (\(k=1\)) that the only real root will be equal to

\[ p_3=s_1-(p_1+p_2)=-\beta_3+2\alpha_1 . \tag{2-12} \]

From (2-11) we can find

\[ p_3=-\frac{\beta_1}{\alpha_1^2+\omega_1^2}. \tag{2-13} \]

One of these formulas may serve to check the computations.

b) The root of smallest modulus has turned out to be real. From (2-8) we find

\[ p_2+p_3=s_1-p_1, \]

\[ p_2^2+p_3^2=s_2-p_1^2. \]

Squaring the first and subtracting the second from the result, we then obtain

\[ 2p_2p_3=(s_1-p_1)^2-(s_2-p_1^2). \]

Next, taking into account that

\[ (p_2-p_3)^2=p_2^2-2p_2p_3+p_3^2 =2(s_2-p_1^2)-(s_1-p_1)^2, \]

we shall have two linear equations

\[ p_2-p_3=\pm\sqrt{2(s_2-p_1^2)-(s_1-p_1)^2}, \]

\[ p_2+p_3=s_1-p_1. \]

Their simultaneous solution gives the desired formula

\[ p_{2,3}=\frac{1}{2}(s_1-p_1)\pm \frac{1}{2}\sqrt{2(s_2-p_1^2)-(s_1-p_1)^2}. \tag{2-14} \]

For checking we have the independent relation

\[ p_1\cdot p_2\cdot p_3=-\beta_1. \tag{2-15} \]

2. Fourth-degree equation

a) The root of smallest modulus forms a complex-conjugate pair
\((p_1=-\alpha_1+j\omega_1;\ p_2=-\alpha_1-j\omega_1)\).

If we carry out all the derivations, entirely analogous to those we performed in deriving formula (2-14), and take into account that

\[ p_1+p_2=-2\alpha_1 \quad \text{and} \quad p_1^2+p_2^2=2(\alpha_1^2-\omega_1^2), \]

then we obtain

\[ p_{3,4}=\frac{1}{2}(s_1+2\alpha_1)\pm \frac{1}{2}\sqrt{2[s_2-2(\alpha_1^2-\omega_1^2)]-(s_1+2\alpha_1)^2}. \tag{2-16} \]

For checking, we again use the formula

\[ p_3\cdot p_4=\frac{\beta_1}{\alpha_1^2+\omega_1^2}. \tag{2-17} \]

b) The root of smallest modulus \((p_1)\) is real. We construct the sequence \(N\) and find the root of largest modulus. Suppose that it forms a complex-conjugate pair
\((p_3=-\alpha_3+j\omega_3\) and \(p_4=-\alpha_3-j\omega_3)\).

The remaining root \(p_2\) can only be real, and its value will be equal to

\[ p_2=s_1-p_1+2\alpha_3=-\beta_4-p_1+2\alpha_3. \tag{2-18} \]

The check in this case can be carried out by the formula

\[ p_2=\frac{\beta_1}{p_1(\alpha_3^2+\omega_3^2)}. \tag{2-19} \]

c) This case is characterized by the fact that both roots \((p_1\) and \(p_4)\), the smallest and the largest in modulus, are real. The other two roots are found by a formula analogous to (2-14) and (2-16):

\[ p_{2,3}=\frac{1}{2}(s_1-p_1-p_4)\pm \frac{1}{2}\sqrt{2(s_2-p_1^2-p_4^2)-(s_1-p_1-p_4)^2}. \tag{2-20} \]

For checking we have

\[ p_2\cdot p_3=\frac{\beta_1}{p_1\cdot p_4}. \tag{2-21} \]

3. Fifth-degree equation

a) The root of smallest modulus forms a pair \(p_1=-\alpha_1+j\omega_1\) and \(p_2=-\alpha_1-j\omega_1\), while the root of largest modulus forms a pair \(p_4=-\alpha_4+j\omega_4\) and \(p_5=-\alpha_4-j\omega_4\). The root \(p_3\) will obviously be real and can be determined from the formula

\[ p_3=s_1+2\alpha_1+2\alpha_4=-\beta_5+2\alpha_1+2\alpha_4. \tag{2-22} \]

For checking we have

\[ p_3=-\frac{\beta_1}{(\alpha_1^2+\omega_1^2)\cdot(\alpha_4^2+\omega_4^2)}. \tag{2-23} \]

b) In this case we combine two variants:

1. The root of smallest modulus is real \((p_1)\), while the largest forms the pair
   \(p_4=-\alpha_4+j\omega_4;\ p_5=-\alpha_4-j\omega_4\).

2. Conversely, the root of smallest modulus forms the pair
   \(p_1=-\alpha_1+j\omega_1;\ p_2=-\alpha_1-j\omega_1\), while the largest is real.

In the first variant

\[ p_{2,3}=\frac{1}{2}(s_1-p_1+2\alpha_4)\pm \frac{1}{2}\sqrt{\,2\left[s_2-p_1^2-2(\alpha_4^2-\omega_4^2)\right]-(s_1-p_1+2\alpha_4)^2\,}. \tag{2-24} \]

For checking we shall use the formula

\[ p_2\cdot p_3=-\frac{\beta_1}{p_1(\alpha_4^2+\omega_4^2)}. \tag{2-25} \]

In the case of the second variant we obtain the formula

\[ p_{3,4}=\frac{1}{2}(s_1-p_5+2\alpha_1)\pm \frac{1}{2}\sqrt{\,2\left[s_2-p_5^2-2(\alpha_1^2-\omega_1^2)\right]-(s_1-p_5+2\alpha_1)^2\,}, \tag{2-26} \]

and for checking we have

\[ p_3\cdot p_4=-\frac{\beta_1}{p_5(\alpha_1^2+\omega_1^2)}. \tag{2-27} \]

c) Both roots turned out to be real. In this case, in order to determine the three other roots it is necessary to solve an auxiliary equation of the third degree. Suppose that it has the form

\[ p^3+\beta_3^1p^2+\beta_2^1p+\beta_1^1=0. \tag{2-28} \]

From (2-10) and (2-8) it follows that:

\[ \left. \begin{aligned} s_1^1&=-\beta_3^1=p_2+p_3+p_4,\\ s_2^1&=-2\beta_2^1+(\beta_3^1)^2=p_2^2+p_3^2+p_4^2,\\ s_3^1&=3(\beta_3^1\beta_2^1-\beta_1^1)-(\beta_3^1)^3=p_2^3+p_3^3+p_4^3. \end{aligned} \right\} \tag{2-29} \]

According to (2-8), for the original equation we have

\[ s_1=p_1+p_2+p_3+p_4+p_5 \]

\[ s_2=p_1^2+p_2^2+p_3^2+p_4^2+p_5^2 \]

\[ s_3=p_1^3+p_2^3+p_3^3+p_4^3+p_5^3. \]

Comparing, we arrive at the formulas:

\[ \left. \begin{aligned} s_1^1&=s_1-p_1-p_5,\\ s_2^1&=s_2-p_1^2-p_5^2,\\ s_3^1&=s_3-p_1^3-p_5^3. \end{aligned} \right\} \tag{2-30} \]

From (2-29) we easily find

\[ \beta_3^1=-s_1^1 \]

\[ \beta_2^1=\frac{1}{2}(s_1^1)^2-\frac{1}{2}s_2^1 \]

\[ \beta_1^1=\beta_3^1\beta_2^1-\frac{1}{3}\left[s_3^1+(\beta_3^1)^3\right]. \tag{2-31} \]

Solving equation (2-28), we find the missing roots \(p_2, p_3\), and \(p_4\). For checking the computations we have

\[ p_2\cdot p_3\cdot p_4=-\frac{\beta_1}{p_1\cdot p_5}. \tag{2-32} \]

4. Equation of the sixth degree

a) This case is characterized by the same conditions as case a) for an equation of the fifth degree, that is, we have two pairs of roots \(p_1=-\alpha_1+j\omega_1\) and \(p_2=-\alpha_1-j\omega_1\), as well as \(p_5=-\alpha_5+j\omega_5\) and \(p_6=-\alpha_5-j\omega_5\). The other two roots are found from the formula

\[ p_{3,4}=\frac{1}{2}(s_1+2\alpha_1+2\alpha_5)\pm \frac{1}{2}\sqrt{2\left[s_2-2(\alpha_1^2-\omega_1^2)-2(\alpha_5^2-\omega_5^2)\right]-(s_1+2\alpha_1+2\alpha_5)^2}. \tag{2-33} \]

For checking, we have the equality

\[ p_3\cdot p_4=\frac{\beta_1}{(\alpha_1^2+\omega_1^2)(\alpha_5^2+\omega_5^2)}. \tag{2-34} \]

b) See case b) for the equation of the fifth degree.

First variant: the real root \(p_1\) and the pair of complex-conjugate roots with the largest modulus \(p_5=-\alpha_5+j\omega_5\) and \(p_6=-\alpha_5-j\omega_5\) have been determined. To determine the other three roots, we construct an auxiliary equation of the third degree, having the form (2-28). The coefficients of this equation are determined by formulas (2-31); however, now

\[ \left. \begin{aligned} s_1^1&=s_1-p_1+2\alpha_5,\\ s_2^1&=s_2-p_1^2-2(\alpha_5^2-\omega_5^2),\\ s_3^1&=s_3-p_1^3+2\alpha_5(\alpha_5^2-3\omega_5^2). \end{aligned} \right\} \tag{2-35} \]

The check formula has the form

\[ p_2\cdot p_3\cdot p_4=\frac{\beta_1}{p_1(\alpha_5^2+\omega_5^2)}. \tag{2-36} \]

In the second variant, when the pair of roots \(p_1=-\alpha_1+j\omega_1\) and \(p_2=-\alpha_1-j\omega_1\) and the real root with the largest modulus \(p_6\) have been determined, instead of (2-33) we shall have

\[ \left. \begin{aligned} s_1^1&=s_1-p_6+2\alpha_1,\\ s_2^1&=s_2-p_6^2-2(\alpha_1^2-\omega_1^2),\\ s_3^1&=s_3-p_6^3+2\alpha_1(\alpha_1^2-3\omega_1^2). \end{aligned} \right\} \tag{2-37} \]

The coefficients of the auxiliary cubic equation (2-28) are still determined by formulas (2-31), and instead of (2-36) we may write

\[ p_3\cdot p_4\cdot p_5=\frac{\beta_1}{p_6(\alpha_1^2+\omega_1^2)}. \tag{2-38} \]

c) Both roots \(p_1\) and \(p_6\) are real.

We construct an auxiliary equation of the fourth degree

\[ p^4+\beta_4^1p^3+\beta_3^1p^2+\beta_2^1p+\beta_1^1=0. \tag{2-39} \]

According to the method considered earlier, we may write

\[ \left. \begin{aligned} s_1^1&=s_1-p_1-p_6\\ s_2^1&=s_2-p_1^2-p_6^2\\ s_3^1&=s_3-p_1^3-p_6^3\\ s_4^1&=s_4-p_1^4-p_6^4. \end{aligned} \right\} \tag{2-40} \]

From system (2-10) we find the coefficients of equation (2-39)

\[ \left. \begin{aligned} \beta_4^1&=-s_1^1,\\ \beta_3^1&=\frac{1}{2}(s_1^1)^2-\frac{1}{2}s_2^1,\\ \beta_2^1&=\beta_4^1\beta_3^1-\frac{1}{3}\left[s_3^1+(\beta_4^1)^3\right],\\ \beta_1^1&=-\frac{1}{4}\left[s_4^1+\beta_4^1s_3^1+\beta_3^1s_2^1+\beta_2^1s_1^1\right]. \end{aligned} \right\} \tag{2-41} \]

Thus, we have reduced the problem to solving an equation of the fourth degree. For verification we can again use formula (2-11). Equations of higher degrees can be solved in an entirely analogous way. The matter reduces to constructing auxiliary equations whose degrees will depend on the nature of the initially determined roots with the smallest and largest moduli. The method for determining the coefficients of the auxiliary equations is clear from the preceding. It is applicable to equations of any degree. We note that the formulas we have obtained for determining the other roots, as well as the formulas for determining the coefficients of the auxiliary equations, are exact. Naturally, they will be approximate if approximate values of the initially found roots with the smallest and largest moduli are substituted into them.

Let us consider several examples*, illustrating the computational technique by the proposed method. All computations will be carried out on an ordinary slide rule.

Example 1. The characteristic equation of a belt-drive mechanism has the form

\[ D(p)=p^3+8.2p^2+157p+231=0. \]

Thus, we have

\[ \beta_3=8.2;\quad \beta_2=157;\quad \beta_1=231\quad(\beta_4=1). \]

Let us find the roots of this equation. To determine the root smallest in modulus, we construct the sequence \(M\). Formulas (1-17) give

\[ M_0=\frac{1}{\beta_1}=43.3\cdot 10^{-4};\quad M_1=\frac{\beta_2}{\beta_1}M_0=29.2\cdot 10^{-4}; \]

\[ M_2=-\frac{1}{\beta_1}(M_0\beta_3-M_1\beta_2)=18.3\cdot 10^{-4}; \]

\[ M_3=\frac{1}{\beta_1}(M_0\beta_4-M_1\beta_3+M_0\beta_2)=11.4\cdot 10^{-4}; \]

\[ M_4=-\frac{1}{\beta_1}(-M_1\beta_4+M_2\beta_3-M_3\beta_2)=7.2\cdot 10^{-4}; \]

\[ M_5=\frac{1}{\beta_1}(M_2\beta_4-M_3\beta_3+M_4\beta_2)=4.53\cdot 10^{-4}; \]

\[ M_6=-\frac{1}{\beta_1}(-M_3\beta_4+M_4\beta_3-M_5\beta_2)=2.87\cdot 10^{-4}. \]

* All examples are borrowed from the book by V. V. Solodovnikov et al., The Frequency Method for Constructing Transient Processes, GITTL, Moscow, 1955. There, too, the root values are given to three significant digits. We shall call these values exact.

We shall stop at this. The terms of the sequence decrease monotonically; consequently, the root of smallest modulus is real.

Let us form three ratios:

\[ \frac{M_3}{M_4}=1,58;\qquad \frac{M_4}{M_5}=1,59;\qquad \frac{M_5}{M_6}=1,575; \]

their arithmetic mean, accurate to three significant figures, is equal to \(1,58\), which coincides with the exact value. Thus, the desired root is \(p_1=-1,58\).

Let us find the other two roots. By formulas (2-10) we find

\[ s_1=-\beta_3=-8,2;\qquad s_2=\beta_3^2-2\beta_2=-246,8. \]

Bearing in mind that

\[ s_1-p_1=-8,2+1,58=-6,62 \quad \text{and} \quad s_2-p_1^2=-249,3, \]

by formula (2-14) we obtain

\[ p_{2,3}=-3,31\pm j\,11,65. \]

The check formula (2-15) gives

\[ p_1=-\frac{\beta_1}{p_2\cdot p_3}=-1,575. \]

The exact value of the roots \(p_{2,3}\) is equal to \(-3,32\pm j\,11,6\).

Example 2. The characteristic equation of an electrohydraulic tracking system has the form

\[ D(p)=p^4+103p^3+3065p^2+149250p+1081500=0, \]

that is,

\[ \beta_5=1;\quad \beta_4=103;\quad \beta_3=3065;\quad \beta_2=149250;\quad \beta_1=1081500. \]

Let us find the roots of this equation. To determine the root of smallest modulus, we construct the sequence \(M\) by formulas (1-17)

\[ M_0=0,925\cdot 10^{-6};\quad M_1=12,7\cdot 10^{-6};\quad M_2=1,75\cdot 10^{-6};\quad M_3=0,205\cdot 10^{-6}; \]

\[ M_4=0,0245\cdot 10^{-6};\quad M_5=0,00295\cdot 10^{-6};\quad M_6=0,000355\cdot 10^{-6}. \]

The character of the sequence indicates that the desired root is real. Let us form the ratios:

\[ \frac{M_3}{M_4}=8,36;\qquad \frac{M_4}{M_5}=8,3;\qquad \frac{M_5}{M_6}=8,3. \]

We take the value \(p_1=-8,3\) as more reliable, since it corresponds to larger values of \(k\). This value of the root coincides with the exact one. We now determine the root with the greatest modulus; for this purpose we construct the sequence \(N\) by formulas (2-4):

\[ N_0=1;\qquad N_1=\beta_4=103;\qquad N_2=\beta_4^2-\beta_3=0,754\cdot 10^4; \]

\[ N_3=\beta_2-N_1\beta_3+N_2\beta_4=60,93\cdot 10^4; \]

\[ N_4=-(\beta_1-N_1\beta_2+N_2\beta_3-N_3\beta_4)=57,8\cdot 10^6; \]

\[ N_5=-N_1\beta_1+N_2\beta_2-N_3\beta_3+N_4\beta_4=51\cdot 10^8; \]

\[ N_6=-(N_2\beta_1-N_3\beta_2+N_4\beta_3-N_5\beta_4)=44,31\cdot 10^{10}; \]

\[ N_7=-N_3\beta_1+N_4\beta_2-N_5\beta_3+N_6\beta_4=38,84\cdot 10^{12}; \]

\[ N_8=-(N_4\beta_1-N_5\beta_2+N_6\beta_3-N_7\beta_4)=33,39\cdot 10^{14}. \]

The sequence changes monotonically without a change of sign; consequently, the sought root is real. We form the ratios according to (2-6):

\[ \frac{N_6}{N_5}=87;\qquad \frac{N_7}{N_6}=87.5;\qquad \frac{N_8}{N_7}=85.9. \]

We take \(p_4=-85.9\), which coincides with the exact value. Thus we have case c) for a fourth-degree equation. The two other roots are found from formula (2-20) \((S_1=-103;\ S_2=4480)\)

\[ p_{2,3}=-4.4\pm j38.7, \]

The exact value of the roots is \(-4.42\pm j38.7\).
The check formula (2-21) gives

\[ p_2\cdot p_3=\frac{\beta_1}{p_1\cdot p_4}=1480 \]

or

\[ \sqrt{p_2\cdot p_3}=|p_3|=38.5, \]

while the modulus of the root found is approximately \(38.7\).

Example 3. For the airplane–autopilot system we have a sixth-degree characteristic equation with coefficients:

\[ \beta_7=1;\quad \beta_6=16.4;\quad \beta_5=107.4;\quad \beta_4=364.2;\quad \beta_3=1146.5;\quad \beta_2=771.2;\quad \beta_1=292.1 \]

i.e.,

\[ D(p)=p^6+16.4p^5+107.4p^4+364.2p^3+1146.5p^2+771.2p+292.1=0. \]

Let us find the roots of this equation. We begin, as usual, with the root of smallest modulus. We construct the sequence \(M\)

\[ M_0=34.4\cdot10^{-4};\qquad M_1=91\cdot10^{-4};\qquad M_2=106\cdot10^{-4};\qquad M_3=-33\cdot10^{-4}; \]

\[ M_4=-403\cdot10^{-4};\qquad M_5=-492\cdot10^{-4};\qquad M_6=205\cdot10^{-4}; \]

\[ M_7=1985\cdot10^{-4};\qquad M_8=3960\cdot10^{-4};\qquad M_9=3050\cdot10^{-4}. \]

The nature of the sequence leaves no doubt that the sought root forms a complex-conjugate pair. By formula (1-18) we have:
for \(\kappa=5\)

\[ a_1\simeq \frac{1}{2}\frac{|M_7\cdot M_4-M_6\cdot M_5|}{|M_6^2-M_5\cdot M_7|} =\frac{1}{2}\cdot0.686=0.343; \]

for \(\kappa=6\)

\[ a_1\simeq \frac{1}{2}\frac{|M_8M_5-M_7\cdot M_6|}{|M_7^2-M_6\cdot M_8|} =\frac{1}{2}\cdot0.752=0.376; \]

for \(\kappa=7\)

\[ a_1\simeq \frac{1}{2}\frac{|M_9\cdot M_6-M_8\cdot M_7|}{|M_8^2-M_7\cdot M_9|} =\frac{1}{2}\cdot0.75=0.375. \]

We take the last value. We find the square of the modulus of the root \((\kappa=7)\)

\[ a_1^2+\omega_1^2\simeq \frac{|M_7^2-M_8\cdot M_6|}{|M_8^2-M_7\cdot M_9|} =0.324, \]

whence

\[ \omega_1=\sqrt{0.324-\alpha_1^2}=0.428. \]

Thus, \(p_{1,2}=-0.375\pm j0.428\). The exact value is \(-0.377\pm j0.428\).

Let us determine the root of greatest modulus. We construct the sequence \(N\)

\[ N_0 = 1; \qquad N_1 = 16,4; \qquad N_2 = 161,6; \qquad N_3 = 1,244 \cdot 10^4; \]

\[ N_4 = -7,804 \cdot 10^3; \qquad N_5 = -220,7 \cdot 10^3; \qquad N_6 = -2,617 \cdot 10^6; \]

\[ N_7 = -23,23 \cdot 10^6; \qquad N_8 = -169 \cdot 10^6; \qquad N_9 = -974 \cdot 10^6. \]

The sequence changes sign; consequently, the desired root is complex. By formula (2-7) we find, for \(k = 5, 6, 7\),

\[ \alpha_5 \approx \frac{1}{2}\, \frac{\left|N_7 \cdot N_4 - N_6 \cdot N_5\right|} {\left|N_5^2 - N_6 \cdot N_4\right|} = \frac{1}{2}\cdot 14 = 7 \quad (k=5); \]

\[ \alpha_5 \approx \frac{1}{2}\cdot 13,85 = 6,92 \quad (k=6); \]

\[ \alpha_5 \approx \frac{1}{2}\,14,2 = 7,1 \quad (k=7). \]

We take the value corresponding to \(k=7\); then

\[ \alpha_5^2 + \omega_5^2 \approx \frac{\left|N_8^2 - N_7 \cdot N_9\right|} {\left|N_7^2 - N_8 \cdot N_6\right|} = 61,7, \]

whence

\[ \omega_5 = \sqrt{61,7 - \alpha_5^2} = 3,36. \]

Thus, \(p_{5,6} = -7,1 \pm j\,3,36\). The exact value of the root is \(-7,15 \pm j\,3,35\). We have case a) for an equation of the sixth degree. The two remaining roots are found by formula (2-33):

\[ s_1 = -16,4; \qquad s_2 = 54,2; \]

\[ s_1 + 2\alpha_1 + 2\alpha_5 = -1,45; \qquad s_2 - 2(\alpha_1^2 - \omega_1^2) - 2(\alpha_5^2 - \omega_5^2) = -24. \]

Formula (2-33) gives

\[ p_{3,4} \approx -0,725 \pm j\,3,55. \]

The exact value of the root is \(-0,648 \pm j\,3,74\). The discrepancy is appreciable. As was already noted, it is associated with the error in the roots found initially.

Conclusion

The examples considered show that the proposed method is sufficiently simple and convenient for engineering calculations on an ordinary slide rule and arithmometer. It is also convenient for computations on a digital computer.

REFERENCES

1. V. A. Zaguskin, Handbook of Numerical Methods for Solving Algebraic and Transcendental Equations. Fizmatgiz, Moscow, 1960.
2. A. G. Kurosh, Course of Higher Algebra. Gostekhizdat, Moscow, 1952.