Problem with integration: there is a pole on a contour

I have a problem with calculation of the complex integral

The integrand has a pole in z=1 and z=-1 lying on the unit circle.

My code are as following:

t = @(theta) exp(1i*theta);
f = @(t) (t^2+3*t+2*1i)/(t+4)/(t-1)/(t+1);
integrand = @(theta) f(t(theta))*1i*exp(1i*theta);
integral_result = integral(integrand, 0, 2*pi, 'ArrayValued', true, 'AbsTol', 1e-12, 'RelTol', 1e-12);

Warning: Minimum step size reached near x = 6.28319. There may be a singularity, or the tolerances may be too tight for this problem.

result = (1/(2*pi*1i)) * integral_result;
I_num=result;
disp(['numerical solution=',num2str(I_num)])

numerical solution=0.81432-0.96198i

How do i fix it...

Answer 1

[David Goodmanson]

Hi SB

You can use the principal value for the integration,which is defined as

(1/2) * [ (integral with path on one side of the pole) + (integral with path on other side of the pole)]

You have

pole 3     pole 2     pole 1
  -4         -1         +1

In this case there are two poles that get crossed, so there are 4 paths to consider

                  pole 2   pole 1    case              for example:
path passes to       L        L        a       path of radius 1 shifted to left by 1/2
                     L        R        b       path of radius 2      
                     R        L        c       path of radius .8
                     R        R        d       path of radius 1 shifted to right by 1/2

and the answer is 1/2 x 1/2 = 1/4 of the sum of the integrals along these paths. By explicit integration on those paths:

f = @(z) (z.^2+3*z+2i)./((z+4).*(z-1).*(z+1));
% case a
Ra = 1; shift = -1/2;
z = @(theta) Ra*exp(i*theta) + shift;
integranda = @(theta) f(z(theta)).*(Ra*exp(i*theta))*i;
Ia = integral(integranda,0,2*pi);
% case b
Rb = 2;
z = @(theta) Rb*exp(i*theta);
integrandb = @(theta) f(z(theta)).*(Rb*exp(i*theta))*i;
Ib = integral(integrandb,0,2*pi);
% case c
Rc = .8;                     % this integral will be zero since no poles are enclosed
z = @(theta) Rc*exp(i*theta);
integrandc = @(theta) f(z(theta)).*(Rc*exp(i*theta))*i;
Ic = integral(integrandc,0,2*pi);
% case d
Rd = 1;  shift = 1/2
z = @(theta) Rd*exp(i*theta) + shift;
integrandd = @(theta) f(z(theta)).*(Rd*exp(i*theta))*i;
Id = integral(integrandd,0,2*pi);
I = (Ia+Ib+Ic+Id)/4

    I =  0.4189 + 2.3038i

which agrees with the result from residue theory, no explicit integrartion.

A simpler way, still with numerical integration, is to draw two closely spaced vertical lines up the y axis, which splits the unit circle into a D and a backwards D. The ccw path around the unit circle is the same as the sum of ccw paths around each of the two D's since the integrals along the straight sides of the D's are in opposite directions and cancel each other. For the right hand D, the path going to just the left of the pole results in 0 since the new D encloses no poles. The path going just to the right of the pole does enclose the pole. You can alter the path to be a small circle enclosing the pole and keep the factor of (1/2) from the principal value definition. The backwards D works similarly which results in

    % case a
    Ra = 1/3;     % for example
    z = @(theta) Ra*exp(i*theta) + 1;
    integranda = @(theta) f(z(theta)).*(Ra*exp(i*theta))*i;
    Ia = (1/2)*integral(integranda,0,2*pi);
    % case b
    Rb = 1/3;     % for example
    z = @(theta) Rb*exp(i*theta) - 1;
    integrandb = @(theta) f(z(theta)).*(Rb*exp(i*theta))*i;
    Ib = (1/2)*integral(integrandb,0,2*pi);
    I = Ia+Ib

and gives the same result.

[Shengfan Bi]

Many thanks for your detailed response, David.

Answer 2

[Torsten]

In the usual sense, your integral does not exist. But you can compute its principal value:

syms z 
syms t real
f = (z^2+3*z+2*1i)/((z+4)*(z-1)*(z+1));
f = partfrac(f)
f = subs(f,z,exp(1i*t))*1i*exp(1i*t)
1/(2*sym(pi)*1i)*int(f,0,2*sym(pi),PrincipalValue=true)