# compute the limit of a rational function

Suppose I have a rational function defined by ($s$ complex) $$f(s)=w^T s(sI-Q)^{-1} v$$ for nonzero column vectors $w,v$ and a (large) square matrix $Q$. Further assume that $Q$ is singular and that it is known that $f(0)=\lim_{s\rightarrow 0} f(s)$ exists finite. Clearly, for any $s$ that is not an eigenvalue of $Q$, one can compute $f(s)$ by $$f(s)=s \cdot w^T y \text{ where y is the unique solution of } (sI-Q)y=v$$ which is quite efficient and numerically stable.

Question: is there a similar way of computing $f(0)$, that is (ideally) by solving a single linear system?

[Context. In my problem, $f(s)$ arises as $sF(s)$, where $F(s)$ is the Laplace transform of a certain (exponential) function $g(t)$. The result is needed as a method to compute $f(0)=\lim_{t\rightarrow \infty} g(t)$, which is the final value theorem for Laplace transform.]

• I'm not sure if it helps, but if we consider diagonal $Q$ (e.g. if we diagonalize), then $\lim_{s \to 0} s(sI-Q)^{-1}$ sends the zero diagonal entries to $1$ and the nonzero diagonal entries to $0$.
– user13113
Mar 11, 2015 at 23:28
• It may be simpler to write $sI = (sI-Q) + Q$ and use the alternative form $$f(s) = w^T v + (w^T Q) (sI-Q)^{-1} v$$
– user13113
Mar 14, 2015 at 19:46
• @Hurkyl. Hello, could you explain why this form is simpler? Mar 15, 2015 at 22:38

Let $\Theta$ be the direct sum of two copies of $$\begin{bmatrix}0&b+1\\b&c\end{bmatrix},$$ and let $\kappa:=[2b+1,c]^t\oplus[0,1]^t$ (for some integers $b\ne0,-1$ and $c\ne0$). Consider $H(s):=(1-s)(\Theta-s^{-1}\Theta^*)$. Then the function $f(s):=\langle H(s)^{-1}\kappa,\kappa\rangle$ is identically $0$ whenever defined (and, hence, all limits are $0$), whereas $\kappa$ is not in the image of $H(s)$ for $s=(1+1/b)^{\pm1}$, so these values (that should also be $0$) cannot be "computed directly".
Consider the Jordan canonical form $Q = U J U^{-1}$ for $Q$. Thus $w^T Q w = \sum_B w_B^T B v_B$ for Jordan blocks $B$, where $w_B$ and $v_B$ are the entries of $U^T w$ and $U^{-1} v$ corresponding to the block. For a Jordan block $B$ of size $m$ corresponding to the eigenvalue $0$, $s (sI - B)^{-1}$ is upper triangular with $1$ on the main diagonal and $1/s^j$ on the $j$'th diagonal above that. If it is known (for particular $w$ and $v$) that a finite limit exists, all the terms in $1/s^j$ in $w_B^T s (sI-B)^{-1} v_B$ must cancel, leaving the same result that you would have for $B = 0$. On the other hand, for a Jordan block corresponding to a nonzero eigenvalue, the limit is $0$. Thus if $P$ is the projection on the generalized eigenspace for eigenvalue $0$, the limit, if it exists, is $w^T P v$.
$P$ can be obtained as a polynomial $g(A)$, where if the largest Jordan block for eigenvalue $0$ has size $m_0$, $g(0) = 1$ and $g^{(j)}(0) = 0$ for $j < m_0$, while for every other eigenvalue $\lambda$ with largest Jordan block of size $m_\lambda$, $g(\lambda) = \ldots = g^{(m_\lambda-1)}(\lambda) = 0$.
Of course, this can't be numerically stable: eigenvalues that are exactly $0$ and those that are close to $0$ will produce very different results.