# Can sage find the poles of a complex function, and the residues at those poles?

Hi all,

I searched the SAGE reference and didn't find the word residue in realtion to functions of a complex variable. The function I'm working with is not the ratio of two polynomials in Z, which can be factored out to find the poles and zeros. It's actually like this:

f(z) = A* sinh(z)/ (B* z^2* sinh(z)+ C* z* cosh(z)-D* sin(z))

got this from laplace transforming a partial differential equation I want to solve. I'm following the technique used in a paper I'm trying to apply to a different geometry/coordinate system, and if I can get the residues at the poles, I can get a series solution to the PDE. However, I can't fill in the omitted steps in the paper due to my never having taken a course in complex variables. Hence my interest in finding the residues at the poles programmatically.

I would greatly appreciate a clue or two. Thank you.

edit retag close merge delete

Although not directly Sage-related, perhaps you could use the trick described here [http://en.wikipedia.org/wiki/Residue_%28complex_analysis%29#Series_methods]. I know absolutely nothing about complex analysis. :-) On another note, the Sage reference manual says that there is a way to compute Laurent series (which would give you your residue) [http://www.sagemath.org/doc/reference/sage/symbolic/expression.html#sage.symbolic.expression.Expression.taylor]. But when I fiddled with it, I could only get taylor series. You could also take a look at this [http://www.sagemath.org/doc/reference/sage/rings/laurent_series_ring_element.html], but I couldn't figure anything useful out of that, either.

( 2011-04-11 16:51:27 -0500 )edit

You can find the zeroes of the denominator: "solve(f.denominator() == 0)". I don't know if there is a built-in way to compute residues from there.

( 2011-04-13 09:12:56 -0500 )edit

Just as a data point, Mma apparently does this: http://reference.wolfram.com/mathematica/ref/Residue.html. Unsurprisingly, rjf has also written on this topic.

( 2011-04-16 17:25:26 -0500 )edit

@John Palmieri - unfortunately, solve(f.denominator() == 0,x) doesn't help in this case, and using to_poly_solve=True gives a weird indexing error on the symbolic expression that comes out of the original non-solution that it's too late in this time zone for me to track down tonight. But I couldn't get W|A to do it, either; this probably doesn't have a closed-form solution, n'est pas?

( 2011-04-16 18:00:34 -0500 )edit

Sort by » oldest newest most voted

This won't work with a complicated example such as yours, but perhaps it will give you some ideas:

sage: f = sinh(x)/(sin(x)-1/2)
sage: (1/f == 0).solve(x)
[x == 1/6*pi]
sage: pole = _[0].rhs()
sage: (f*(x-pole)).limit(x=pole)
1/3*(e^(1/3*pi) - 1)*sqrt(3)*e^(-1/6*pi)
sage: f.taylor(x,pole,-1)*(x-pole)
2/3*sqrt(3)*sinh(1/6*pi)


the second line solves 1/f(x)=0, which for meromorphic functions is like solving f(x)=infinity. We then store the pole found in "pole". The last two lines calculates the residue, assuming the pole is simple (the two answers are identical).

Here is an example with a non-simple pole:

sage: f = sinh(x)/(sin(x)-1/2)-1/(x-pi/6)^2
sage: (1/f == 0).solve(x)
[x == 1/6*pi]
sage: pole = _[0].rhs()
sage: f(x=x+pole).taylor(x,0,-1).coeff(x,-1)
2/3*sqrt(3)*sinh(1/6*pi)


The substitution "x=x+pole" is needed - f.taylor(x,pole,-1).coeff(x-pole,-1), which should have given the same answer returns zero, since some automatic simplification takes place and messes things up:

sage: 1/(x-pole)
-6/(pi - 6*x)


hope this helps!

more

Apparently Maxima includes this, and we didn't even know it! See for instance here and the official docs.

(%i1) residue (s/(s**2+a**2), s, a*%i);
(%o1)                           1/2


In Sage, that would be doable.

sage: f(x)=1/x
sage: f.maxima_methods().residue(x,0)
1
f(x)=1/(x^2+1)
f.maxima_methods().residue(x,i)
-1/2*I


Well, we should wrap this!

This does not solve the problem of finding poles. Maple apparently has this, though perhaps not for such an example as the original one.

more

( 2011-04-16 17:58:52 -0500 )edit