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

partial fractionspolynomialsrational-functionsroots

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!

Best Answer

$\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 typical plot

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!

Related Question