I am using the Python library Sympy to calculate the partial fraction decomposition of the following rational function:
Z = (6.43369157032015e-9*s^3 + 1.35203404799555e-5*s^2 + 0.00357538393743079*s + 0.085)/(4.74334912634438e-11*s^4 + 4.09576274286244e-6*s^3 + 0.00334241812250921*s^2 + 0.15406018058983*s + 1.0)
I am calling the method apart() like this:
Z_f = apart(Z, full=True).doit()
And this is the result:
-1.30434068814287e+23/(s + 85524.0054884464) + 19145168.0/(s + 774.88576677949) + 91.9375/(s + 40.7977016133126) + 2.59375/(s + 7.79746609204661)
As you can see from the result, these are the residuals:
Note: I identify them like this:
Z_f_terms = Z_f.as_ordered_terms()
and then extract them in a for cycle.
The problem is that by using other tools I get different residuals.
By using GNU Octave's residue() I get these residuals:
Via Wolphram Alpha I get the same residuals as Octave see here the results (please wait until you see the partial fraction expansion) .
Why isn't the Sympy library not working as expected?
The poles are always correct. Only the residuals are not as I expect.
Thanks for any suggestion.
The problem is to do with floats. Presumably something in the algorithm used by apart(full=True) or RootSum.doit() is numerically unstable for approximate coefficients.
If you convert the coefficients to rational numbers before computing the RootSum then evalf gives the expected result:
In [27]: apart(nsimplify(Z), full=True).evalf()
Out[27]:
133.599202650992 1.07757928431867 0.395006955518971 0.564264854137341
──────────────────── + ─────────────────── + ──────────────────── + ────────────────────
s + 85524.0054884464 s + 774.88576677949 s + 40.7977016133126 s + 7.79746609204661
In [29]: apart(Z, domain="QQ", full=True).evalf() # rational field
Out[29]:
133.599202650994 1.07757928431867 0.395006955518971 0.564264854137341
──────────────────── + ─────────────────── + ──────────────────── + ────────────────────
s + 85524.0054884464 s + 774.88576677949 s + 40.7977016133126 s + 7.79746609204661
Note that the general algorithm and approach used for full=True is really intended for finding an exact representation of the decomposition by avoiding radicals. It is not intended for floats. Probably though normal apart without full=True should compute this decomposition.