Numerical computation of Laurent coefficients for rational functions

Question

In terms of a project I am trying to implement the Carathéodory-Fejér algorithm, if interested see [1] or [2]. However, in one of the steps I need to compute the Laurent coefficients of a rational function $f$, which does only has single poly, on the unit circle. The usual formula $$c_k=\frac {1}{2\pi i} \int_\Gamma \frac{f(z)}{z^{k+1}} dz=\frac {1}{2\pi} \int_0^{2\pi} f(e^{i \varphi})e^{-ik\varphi}d\varphi,$$ is due to oscillation not well-suited for numerical integration. Does any one have some ideas, what I could try in order to determine those coefficients for an arbitrary rational function, with just single pols (no pols on the unit circle itself)?

Answer 1

In case someone is interest, this is what I was coming around. Basically it is the same way you would do by hand, but it is well suited for numerical calculation.

Since $f$ is a rational function we can apply the polynomial to obtain the representation $$f(z)=w(z)+\frac{p(z)}{q(z)},$$ where $m=\mathrm{deg}(p)<\mathrm{deg}(q)=n$. Further, we compute the partial fraction decomposition to obtain the representation $$\frac{p(z)}{q(z)}= \sum_{k=1}^m\frac{b_k}{z-z_k},$$ while $z_k$ denote the roots of the denominator $q$. The coefficients $b_k$ can be efficiently determined by the Barycentric representation of $p$, which states $$p(z)= \left(\prod_{k=1}^m (z-z_k)^{-1} \right)\sum_{k=1}^m\frac{\omega_k}{z-z_k}p(z_k)\qquad \text{for} \quad \omega_k =\prod_{0=j\neq k}^m(z_k-z_j)^{-1}.$$ Thus, we have $b_k= \omega_k p(z_k)$. After appropriate ordering we choose $0\leq \nu \leq m$ such that, $|z_k|<1$ for $k\leq \nu$ and $|z_k|>$ for $k>\nu$. We put this new representation of $f$ into the Laurent coefficient formula and obtain \begin{align*} c_l&=\frac{1}{2 \pi i} \int_\Gamma \frac{w(z)}{z^{l+1}}dz+\frac{1}{2 \pi i}\sum_{k=1}^\nu \int_\Gamma\frac{b_k}{(z-z_k)z^{l+1}}dz+ \frac{1}{2\pi i} \sum_{k=\nu+1}^m\int_\Gamma\frac{b_k}{(z-z_k)z^{l+1}}dz\newline &=:I_1+I_2+I_3. \end{align*} Then, the first integral is given by the generalized Chauchy integral theorem $$I_1 = \frac{w^{(l)}(0)}{l!}.$$ In a similar way we proceed for $I_3$ and get $$I_3 = \sum_{k=\nu+1}^m \frac{h_k^{(l)}(0)}{l!}$$ for $h_k=\frac{b_k}{z-z_k}$. Finally for $I_2$ we apply the residue theorem and obtain \begin{align*}I_2&= \sum_{k=1}^\nu \mathrm{Res}\left( \frac{b_k}{(z-z_k)z^{l+1}},0 \right)+\mathrm{Res}\left( \frac{b_k}{(z-z_k)z^{l+1}},z_k \right)\newline &=\sum_{k=1}^\nu\frac{h_k^{(l)}(0)}{l!}-(l+1)\frac{b_k}{z_k^{l+2}}, \end{align*} again for $h_k=\frac{b_k}{z-z_k}$.