Calculate zeros of rational function directly from residues of partial fraction expansion without calculating polynomial coefficients

Question

Imagine the following rational function in the Laplace-domain with $s = \mathrm{j}\omega$

$$G(s) = \sum_{i=1}^{m} \dfrac{c_i}{s-p_i}= h\,\dfrac{\prod\limits_{i=1}^{n}(s-z_i)}{\prod\limits_{i=1}^{m}(s-p_i)} = \dfrac{A_0 + A_1 s + \ldots A_n s^n}{B_0 + B_1 s + \ldots B_m s^m}$$

for which I know the poles $p_i$ and the residues $c_i$.

I would like to calculate the zeros $z_i$. This is actually trivial: one can calculate the coefficients $A_1 \ldots A_n$ of the numerator polynomial in various ways (most simple would be the Matlab function residue) and calculate the roots of this polynomial. However, in case floating-point-precision is an issue, the polynomial coefficients are the bottleneck and it would be advantageous to avoid that intermediate step.

Question

Do you know any approach, algortihm or source, which explains a possibility to directly calculate the zeros $z_i$ from the residues $c_i$ and poles $p_i$ without calculating the numerator polynomial as intermediate step?

---

At the moment I use Matlab's vpa to increase the precision, but it appears to be a waste of computational effort, if I am not interested in the polynomial but just the zeros and poles. So any hint would be appreciated.

Thank you very much!

Answer 1

$\color{brown}{\textbf{Alternative approach.}}$

From the given equation $$G(s)=\sum\limits_{i=0}^m\dfrac{c_i}{s-p_i} = h\;\dfrac{\prod\limits_{j=1}^n(s-z_j)} {\prod\limits_{i=1}^m(s-p_i)},\quad n\le m-1,\tag1$$

immediately should

• $$\text{If}\quad h_1= \lim\limits_{\large s\to \infty}sG(s) = \sum\limits_{i=1}^mc_i\not=0,\quad\text{then}\quad \dbinom nh = \dbinom {m-1}{h_1},\tag{2.1}$$
• $$\text{Else if}\quad h_2= \lim\limits_{\large s\to \infty}s^2G(s) \not=0,\quad\text{then}\quad \dbinom nh = \dbinom {m-2}{h_2},\dots\tag{2.2}$$ Also, easily to check the expression for the residue $\;c_k\;$ of $k-$th pole $p_k$ in the form of $$c_k =\lim\limits_{\large s\to p_k}(s-p_k)G(s) = h\;\dfrac{\prod\limits_{j=1}^n(p_k-z_j)} {\prod\limits_{i=1}^m(p_k-p_i)^{\large1-\delta_{\large ki}}},\quad (k=1\dots m),\tag3$$ or $$\prod\limits_{j=1}^n(p_k-z_j) = b_k,\quad (k=1\dots m),\tag4$$ where $$b_k = \dfrac{c_k}{h}\prod\limits_{i=1}^m(p_k-p_i)^{\large1-\delta_{\large ki}}\tag5$$ are the constants, and expression with the Kronekker symbol $\delta$ in $(3),(5)$ means the production of non-zero differences.

Solving of the system $(4)$ does not contain calculating the numerator polynomial and looks as the required alternative approach.

If some of the roots are known, then the model order decreases.

$\color{brown}{\textbf{Simple example.}}$

Let us consider the case $\quad\mathbf{ m=3,\; p_1=-1,\; p_2=0,\; p_3=1,\; c_i =1.}\quad$

$$G(s)=\dfrac1{s+1} + \dfrac1s+\dfrac1{s-1}\quad \left(=\dfrac{3s^2-1}{s^3-s}\right),$$ Then $\;h_1=\lim\limits_{s\to\infty} sG(s) =3\;\Rightarrow\; \mathbf{h=3,\;n=2},$ and from $(1)$ should $$G(s) = 3\dfrac{(s-z_1)(s-z_2)}{(s+1)s(s-1)},$$ $$\lim\limits_{s \to-1} (s+1)\left(\dfrac{c_1}{s+1} + \dfrac{c_2}s+\dfrac{c_3}{s-1}\right) = \lim\limits_{s \to-1} 3\;\dfrac{(s-z_1)(s-z_2)}{s(s-1)},$$ $$(-1-z_1)(-1-z_2) = b_1,$$ $$\mathbf{b_1 = \dfrac{(-1)(-1-1)}3 = \dfrac23},$$ $$b_2 = \dfrac{(0+1)(0-1)}3 = -\dfrac13,\quad b_3 = \dfrac{(1+1)(1+0)}3 = \dfrac23,$$ $$\begin{cases} (-1-z_1)(-1-z_2) = \dfrac23\\ (-z_1)(-z_2) = -\dfrac13\\ (1-z_1)(1-z_2) = \dfrac23, \end{cases}\;\Rightarrow\; \begin{cases} z_1z_2+(z_1+z_2) + 1 = \dfrac23\\ z_1z_2 = -\dfrac13\\ z_1z_2 -(z_1+z_2)+1 = \dfrac23, \end{cases}\tag6 $$ $$z_1=\dfrac1{\sqrt3},\quad z_2 = -\dfrac1{\sqrt3}.$$

Therefore, proposed approach works.

$\color{brown}{\textbf{Solving.}}$

Solution of the system $(6)$ can be obtained analytically.

Also, this system can be solved as the optimization task in the form of $$\vec z = \text{ argmin } \{((-1-z_1)(-1-z_2) - \,^2/_3)^2 + ((-z_1)(-z_2) + \,^1\!/_3)^2 +((1-z_1)(1-z_2) - \,^2/_3)^2 \}.$$

If $n\ge3,$ then attempts of the analytical solving of $(4)$ possibly lead to the same numerator polynomial. So the numeric solutions (such as the gradient descend method) can be recommended. Multivariable model allows to avoid the roots mixing via their resorting after each iteration step.

The alternative approach is the searching of the roots in the intervals between the poles via the dihotomy algorithm (without derivatives using). In this way, algorithm should use sufficient set of the initial points to obtain exactly $n$ roots.

At last, can be used combined algorithm, which uses

• simple dihotomy to obtain some real roots between the poles,
• the simplified model $(4)$ for the rest of roots in the complex presentation.

Good luck!

Answer 2

Usually when one is concerned with the numerical evaluation of roots of a given function, manipulating said function can cause many issues. Wilkinson's polynomials, which are of the form $\prod_i(s-z_i)+\epsilon s^k$ for very small $\epsilon$, are a case similar to yours. Indeed it turns out that trying to write your numerator in standard form may result in coefficients which are only marginally off, but the resulting roots are then changed drastically.

Finding the roots of rational functions directly, however, may be done by standard root-finding techniques. This doesn't require you to rewrite the given function in any way (though some methods, such as the Newton-Raphson method, requires derivatives).

If you are searching for real roots only and your function turns out to be $\mathbb R\to\mathbb R$, then there are bracketing methods which guarantee convergence if you can provide bounds where $G(s_p)>0$ and $G(s_m)<0$. Bisection, in particular, will not care even if you have $|G(s)|=\infty$, and is hence easy to implement if you wish to do so yourself, using the poles themselves as initial bounds. If you wish for faster convergence using bracketing methods, you may wish to look at this, which has things built-in to handle such extreme cases as well.

If the above is not what you want, you can of course just try the Newton-Raphson method or similar. I recall seeing somewhere that perhaps the roots of such rational functions can be estimated better somehow, but they are not my particular specialty.