# Trouble with an integral

Hi,

Let f be the following function:

f(t, z0, z1) = (z0z1 - 1)(z0 - z1)/((tz0 - 1)(tz1 - 1)(t - z0)(t - z1)z0)

f is a nice holomorphic function and I would like to compute the contour integral

\int_{|z1|=1} \int_{|z0|=1} f(t, z0, z1) dz0/z0 dz1/z1

Assuming 0 < t < 1, the result should be -4 pi^2/(t^2-1).

In Sage:

var("z0, z1, theta0, theta1")

assume(0 < t < 1)

f = (z0z1 - 1)(z0 - z1)/((tz0 - 1)(tz1 - 1)(t - z0)(t - z1)z0)

g = f.subs({z0: exp(Itheta0), z1: exp(Itheta1)})

Now if I integrate with respect to theta0 first: g.integrate(theta0, 0, 2pi) Sage answers that the integral is zero. If I integrate with respect to theta1 first: factor(g.integrate(theta1, 0, 2pi).integrate(theta0, 0, 2pi)) Sage answers that the integral is -4pi^2/t^2 which is also clearly wrong...

What can I do to make Sage compute it right? (This is a test I need to compute much more complicated integrals after this so I have to make sure Sage gives me the right answer).

edit retag close merge delete

Sort by » oldest newest most voted

For the above rational function f and for the integration on curves that could without doubt detect if a pole is enclosed the following alternative solution worked for me:

var( 'z0,z1,t0,t1,t' );
assume( 0 < t < 1 );

E = (z0*z1 - 1) * (z0 - z1) / ( (t*z0 - 1)*(t*z1 - 1)*(t - z0)*(t - z1)*z0 )
assume( 0<t<1 )

DEN = E.denominator()
E0  = sum( [ E.residue( z0==pole ) for pole, mult in DEN.roots( z0 ) if abs(pole) < 1 and mult==1 ] )
E1  = sum( [ E.residue( z1==pole ) for pole, mult in DEN.roots( z1 ) if abs(pole) < 1 and mult==1 ] )

DEN0 = E0.denominator()
DEN1 = E1.denominator()

E01  = sum( [ E0.residue( z1==pole ) for pole, mult in DEN0.roots( z1 ) if abs(pole) < 1 ] )
E10  = sum( [ E1.residue( z0==pole ) for pole, mult in DEN1.roots( z0 ) if abs(pole) < 1 ] )

print "E   = %s\n" % E
print "E0  = %s"   % E0.factor()
print "E01 = %s\n" % E01
print "E1  = %s"   % E1.factor()
print "E10 = %s\n" % E10


(One really needs a machine provable True in abs(pole) < 1. This may fail in more complicated cases to be correctly recognized.) Results printed:

E   = (z0*z1 - 1)*(z0 - z1)/((t*z0 - 1)*(t*z1 - 1)*(t - z0)*(t - z1)*z0)

E0  = -(2*t*z1 - z1^2 - 1)/((t*z1 - 1)*(t - z1)*(t + 1)*(t - 1))
E01 = 1/(t^2 - 1)

E1  = 1/((t + 1)*(t - 1)*z0)
E10 = 1/(t^2 - 1)


The Residue Theorem is now only needed to conclude.

Unfortunately, i could not convince sage to correctly compute the integrals:

f = (z0*z1 - 1) * (z0 - z1) / ( (t*z0 - 1)*(t*z1 - 1)*(t - z0)*(t - z1)*z0 )
g = f.subs( { z0: exp( i*t0 ) ,
z1: exp( i*t1 ) } )

I10 = g . integrate( t1, 0, 2*pi ) . integrate( t0, 0, 2*pi ) . factor()
I01 = g . integrate( t0, 0, 2*pi ) . integrate( t1, 0, 2*pi ) . factor()

print "I10 is", I10
print "I01 is", I01


Here we get at least in both cases the same false result:

sage: version()
'SageMath version 7.6, Release Date: 2017-03-25'

I10 is -4*pi^2/t^2
I01 is -4*pi^2/t^2

more

After posting i realized that a similar post with similar notation (!) was done one hour before...

( 2017-04-24 23:17:34 +0200 )edit

In the other post, there is a further division by $z_0 z_1$. After reruning the code with the new E...

E = (z0*z1 - 1) * (z0 - z1) / ( (t*z0 - 1)*(t*z1 - 1)*(t - z0)*(t - z1)*z0 ) / (z0*z1)

E0  = sum( [ E.residue( z0==pole ) for pole, mult in DEN.roots( z0 ) if abs(pole) < 1 and mult==1 ] )
E1  = sum( [ E.residue( z1==pole ) for pole, mult in DEN.roots( z1 ) if abs(pole) < 1 and mult==1 ] )

E01  = sum( [ E0.residue( z1==pole ) for pole, mult in E0.denominator().roots( z1 ) if abs(pole) < 1 ] )
E10  = sum( [ E1.residue( z0==pole ) for pole, mult in E1.denominator().roots( z0 ) if abs(pole) < 1 ] )


we "also" get...

sage: print "E01 = %s\n" % E01.factor()
E01 = -1/((t + 1)*(t - 1))

sage: print "E10 = %s\n" % E10 # .factor()
E10 = 0

( 2017-04-24 23:27:30 +0200 )edit

The above was recomputed in a new sage session. The results agree:

sage: E10.factor()
-1/((t + 1)*(t - 1))
sage: E01.factor()
-1/((t + 1)*(t - 1))

( 2017-04-24 23:47:12 +0200 )edit

Thank you very much for your answer and also for providing me something almost automated.

You need to remove the mult == 1 in the first computation to get the correct result.

The initial problem is to compute the Hilbert series of the ring of invariants of a certain representation of the orthogonal group. Here the test is for the standard representation of SO(4).

( 2017-04-25 01:12:57 +0200 )edit

Applying Cauchy's residue theorem gives the result you expected (modulo a sign):

sage: t = var('t', domain='real'); assume(0<t); assume(t<1)
sage: (z0, z1) = var('z0, z1', domain='complex')
sage: f(t, z0, z1) = (z0*z1 - 1)*(z0 - z1)/((t*z0 - 1)*(t*z1 - 1)*(t - z0)*(t - z1)*z0)
sage: fi = f/(z0*z1)
sage: g0 = fi.residue(z0==t) + fi.residue(z0==0)
sage: g10 = g0.residue(z1==t) + g0.residue(z1==0)
sage: g10 = (g10 * (2*pi*I)^2).simplify_full()
sage: g10
4*pi^2/(t^2 - 1)


i think the sign is the good one, because the change of variables is rather $\frac{dz_j}{i z_j}$ for $j=0, 1$, isn't it?

moreover, for the integrals in $\theta$ i'm getting $-4\pi^2/t^2$, irrespective of the order of integration. but this is clearly contradictory with the result above and also with the symbolic integrals pre-evaluated at $t=0$:

sage: (theta0, theta1) = var('theta0, theta1', domain='real')
sage: g = -f.subs({t:0, z0: exp(I*theta0), z1: exp(I*theta1)})
sage: g.integrate(theta0, 0, 2*pi).integrate(theta1, 0, 2*pi).simplify_full()
-4*pi^2


for this, i would suspect some bug is going on with the trigonometric integrals, although i cannot confirm the issue.. maybe try with other interfaces? (algorithm keyword).

more

for the integrals in $\theta$, i'm getting $-4\pi^2/t^2$ irrespective of the order of integration.

( 2017-04-24 21:32:45 +0200 )edit

The Cauchy residue theorem was actually my first guess but since I will work with more complicated things that might take me out of the workd of rational functions, I prefered to leave Sage choose the best option.

There is indeed something fishy with the integrate function. You can actually set t equal to any value in the range (0, 1) before evaluation the integral and you'll get the corresponding value $-4\pi^2/t^2$...

Once again thank you very much!

( 2017-04-25 00:46:03 +0200 )edit
1

@mforets - Can you report that bit on Trac? Usually there is something interesting going on with branch choices when this kind of thing happens.

( 2017-04-25 15:27:14 +0200 )edit