Filter löschen
Filter löschen

Problem with integration: there is a pole on a contour

15 Ansichten (letzte 30 Tage)
Shengfan Bi
Shengfan Bi am 14 Sep. 2024 um 3:20
Bearbeitet: Torsten am 14 Sep. 2024 um 13:49
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...

Akzeptierte Antwort

David Goodmanson
David Goodmanson am 14 Sep. 2024 um 7:29
Bearbeitet: David Goodmanson am 14 Sep. 2024 um 8:05
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.

Weitere Antworten (1)

Torsten
Torsten am 14 Sep. 2024 um 10:38
Bearbeitet: Torsten am 14 Sep. 2024 um 13:49
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)

Community Treasure Hunt

Find the treasures in MATLAB Central and discover how the community can help you!

Start Hunting!

Translated by