UPDATE: I've attempted to answer this question
I'm trying to invert a Laplace Transform using Cauchy's Residue Theorem. The transform is
$$\Theta(s)=\frac{F(s)}{G(s)}$$
where
$$ F(s)=[1+\gamma P(s)] [As^{-1}+\theta(0)] $$
and $$ G(s)=[1+\gamma P(s)] (s+A) + \gamma BP(s) $$
in which $\theta(0)$ is the initial condition of the inverse transform function in the time domain, $\gamma$, $A$, and $B$ are real constants, and $P(s)$ is given by the series
$$ P(s)=s\sum_{k=1}^{k=N} \frac{C_k}{s+\lambda_k^2} $$
where $N \ge 50$. The terms $C_k$ and $\lambda_k$ are the real Fourier-Bessel coefficients and real eigenvalues that are the outcome of another coupled problem.
Here is the approach that I used:
There is one pole at $s=0$ and the remaining $N$ poles $z_k, k=1\cdots N$ are the zeroes of $G(s)$, all of which lie on the negative Real Axis, and each of which I have calculated with residuals pretty close to machine roundoff. Following the Residue Theorem, the Inverse Laplace Transform can be written as
$$ \theta(t)=1+\sum_{k=1}^{k=N}e^{z_k t}F(z_k) \lim_{s\to z_k}\frac{s-z_k}{G(s)} $$
whereupon the indeterminate limit can be resolved with L'Hôpital's Rule yielding
$$ \theta(t)=1+\sum_{k=1}^{k=N} e^{z_k t}\frac{F(z_k)}{G'(s)} $$
where $G'(s)$ is the derivative w.r.t. $s$ of $G(s)$
$$ G'(s)=[1+\gamma P(s)] +\gamma (s+A+B)P'(s) $$
and $P'(s)$ is the derivative w.r.t. $s$ of $P(s)$
$$ P'(s)=\sum_{k=1}^{k=N} \frac{C_k}{s+\lambda_k^2}-s\sum_{k=1}^{k=N} \frac{C_k}{(s+\lambda_k^2)^2} $$
Issues:
I've observed that the evaluation of $\theta(t)$ doesn't seem quite right especially for small times and in particular for larger values of $\gamma$. Moreover, it demonstrates a sensitivity to the eigenvalues and Coefficients $C_k$ that is not consistent with expected observations. Rather than inundate the reader with more details at this stage, I'll be glad to update the question with more relevant information following any comments.
My specific questions:
Is the approach that I've used with the Residue Theorem fundamentally correct? Is there anything else that I should be aware of when solving problems of the type above?
Many thanks, as always.
Maybe incorporating a pre-processing to reduce the transfer function to its simplest expression, make the inversion process more accurate. For $N = 2$, the final transfer function is presented, adequate to apply the method
$$ \Theta(s) = \frac{A \left(C_1 \gamma s \left(\lambda _2^2+s\right)+\left(\lambda _1^2+s\right) \left(C_2 \gamma s+\lambda _2^2+s\right)\right)}{s \left(C_1 \gamma s \left(\lambda _2^2+s\right) (A+B+s)+\left(\lambda _1^2+s\right) \left(C_2 \gamma s (A+B+s)+(A+s) \left(\lambda _2^2+s\right)\right)\right)} $$
with denominator coeficients
$$ \left\{ \begin{array}{c} 0 \\ A \lambda _1^2 \lambda _2^2 \\ \left(\lambda _2^2+A+(A+B) \gamma C_2\right) \lambda _1^2+\left(A+(A+B) \gamma C_1\right) \lambda _2^2 \\ \lambda _1^2+\lambda _2^2+A+\gamma C_2 \left(\lambda _1^2+A+B\right)+\gamma C_1 \left(\lambda _2^2+A+B\right) \\ \gamma \left(C_1+C_2\right)+1 \\ \end{array} \right. $$
Follow a symbolic processing script in MATHEMATICA to implement this simplification
Attached a python script to obtain the inverse Laplace transform from a transfer function $G(s) = \frac{N(s)}{D(s)}$ with the $N(s) = a_n s^n + \cdots + a_0$ and $D(s) = b_m s^m+\cdots b_0$ The $N$ coefficients are given in num and the $D$ coefficients in den both in power decreasing style. $G(s)$ is assumed strictly proper. The script gives a sequence of response terms which should be added to compound the time inversion.
The inversion problem gives as result
That should be added to give the response to an impulse. In the figure are shown two curves. One the assumed exact and the other, the computed by the script.