Hello everyone, I have a question about numerical integration. The formula is shown below, where the integration path C_beta is a closed curve on the complex plane.

It’s important to note that the integration path C_\beta is not an analytical expression; it is a set of complex numbers obtained through a complicated code. I am familiar with simple MATLAB integration functions like integral, quadgk, and so on, but I’m unsure how to use these functions for my integration. I also strongly suspect that these functions may not work for my case, as the integration path C_\beta is not analytical. Therefore, I used a rather crude approach, expressing the integral as a series sum. This method works in many cases, but under certain parameters, it gives results that are unacceptable. Below is the code that produces the intolerable results for a specific set of parameter values(The data used is included in the attachment, the path C_\beta varies with different parameter sets. The C_\beta here only applies to the following set of parameter values.):

load list_beta.mat

t1=0;
t2=0.13;
t3=1/5;
gamma1=5/3;
gamma2=1/3;


R_jia=@(x)(t2-gamma2/2)*(x^-1)+(t1+gamma1/2)+t3*x;
R_jian=@(x)t3*(x^-1)+(t1-gamma1/2)+(t2+gamma2/2)*x;
q=@(x)R_jia(x)/sqrt(R_jia(x)*R_jian(x));
q_inv=@(x)R_jian(x)/sqrt(R_jia(x)*R_jian(x));
qSet=arrayfun(q,list_beta);
figure
scatter(real(list_beta),imag(list_beta),'o')

list_dq=qSet(2:end)- qSet(1:end-1);
list_q_inv=arrayfun(q_inv,list_beta);
list_q_inv_middle=0.5*(list_q_inv(2:end)+list_q_inv(1:end-1));

w=sum(list_q_inv_middle.*list_dq)*1j/(2*pi)

w = 0.1141 - 0.0012i

I have determined through a special method that the correct result of the integral for this set of parameters should be 0. However, when using the rough approach mentioned earlier, I obtained a value that deviates significantly from 0 for these parameters, which is unacceptable.

I’m a beginner in numerical integration. Could anyone suggest the appropriate method to correctly compute the integral for this set of parameters? (Note: My integral formula can indeed be simplified to obtain the correct result, which is 0, for this set of parameters. However, I’ve completely abandoned this simplification because it gives more incorrect results than the original formula for many other parameter values. Therefore, I prefer to use the original integral formula for the calculation.)

I would greatly appreciate any responses and thank you all in advance for your help.

[Torsten]

You say your integration path is a closed curve. If so, then close it (see above). This will definitely change the value of your integral.

[ma Jack]

Thank you very much for your comments. However, as mentioned in my question, the array list_beta forms a closed curve on the complex plane. I wrote the relevant plotting commands in my code, but it runs too slowly online, so I wasn't able to generate the plot. Perhaps you could try plotting it yourself.

Answer 1

[Torsten]

However, as mentioned in my question, the array list_beta forms a closed curve on the complex plane.

Isn't then dq = q' dbeta and you should integrate

q'(beta)*q^(-1)(beta) dbeta

over the list_beta values ?

I'm confused about your integral notation.

And are you sure that q^(-1) is just 1/q and not the inverse function of q with respect to beta ?

The different results in the code below might be related to the evaluation of the complex sqrt. If you replace q by the line that is commented out, the results (almost) agree.

load list_beta.mat
t1=0;
t2=0.13;
t3=1/5;
gamma1=5/3;
gamma2=1/3;
Rplus = @(x)(t2-gamma2/2)./x+(t1+gamma1/2)+t3*x;
Rminus = @(x)t3./x+(t1-gamma1/2)+(t2+gamma2/2)*x;
Rplusdot = @(x)-(t2-gamma2/2)./x.^2+t3;
Rminusdot = @(x)-t3./x.^2+(t2+gamma2/2);

q11 = @(x)Rplus(x)./sqrt(Rplus(x).*Rminus(x));
qinv11 = @(x) 1./q11(x);
qdot11 = @(x) 0.5*qinv11(x).*(Rplusdot(x).*Rminus(x)-Rplus(x).*Rminusdot(x))./Rminus(x).^2;
integrand11 = qinv11(list_beta).*qdot11(list_beta);
integral11 = 1i/(2*pi)*trapz(list_beta,integrand11)

integral11 = -2.6421e-05 + 6.6262e-17i

q12 = @(x)sqrt(Rplus(x)./Rminus(x));
qinv12 = @(x) 1./q12(x);
qdot12 = @(x) 0.5*qinv12(x).*(Rplusdot(x).*Rminus(x)-Rplus(x).*Rminusdot(x))./Rminus(x).^2;
integrand12 = qinv12(list_beta).*qdot12(list_beta);
integral12 = 1i/(2*pi)*trapz(list_beta,integrand12)

integral12 = -2.6421e-05 + 8.3931e-17i

integrand21 = qinv11(list_beta);
integral21 = 1i/(2*pi)*trapz(q11(list_beta),integrand21)

integral21 = 0.1142 - 0.0012i

integrand22 = qinv12(list_beta);
integral22 = 1i/(2*pi)*trapz(q12(list_beta),integrand22)

integral22 = 4.7550e-06 + 3.5339e-17i

curvelength = cumsum(sqrt( (real(list_beta(2:end))-real(list_beta(1:end-1))).^2 + ...
                           (imag(list_beta(2:end))-imag(list_beta(1:end-1))).^2 ));
curvelength = [0;curvelength];
plot(curvelength,[real(integrand11),imag(integrand11)])
hold on
plot(curvelength,[real(integrand12),imag(integrand12)])
grid on

●

[ma Jack]

Sir, thank you very much for your insightful response. The first doubt in your response was correct. Indeed, what I meant to express is q′(β)⋅q−1(β) dβ.

However, there are still a few doubts I have not resolved. First, why does the complex square root have such a significant impact on the results, as you mentioned at the beginning?

Secondly, aren’t q11 and q12 mathematically equivalent? The latter is simply a simplification of the former. But why does this difference in expression cause such a significant difference in the integral results (i.e., integral21 and integral22) in MATLAB or computational environments? Apologies, I’m a novice in numerical computation, but I’m very interested in understanding the fundamental reason behind this noticeable discrepancy. Perhaps this also ties back to the first doubt.

Finally, what is the purpose of including the following code in the response?

curvelength = cumsum(sqrt( (real(list_beta(2:end))-real(list_beta(1:end-1))).^2 + ...
                           (imag(list_beta(2:end))-imag(list_beta(1:end-1))).^2 ));
curvelength = [0;curvelength];
plot(curvelength,[real(integrand11),imag(integrand11)])
hold on
plot(curvelength,[real(integrand12),imag(integrand12)])
grid on

Since integral11 and integral12 are already close to the correct result, why analyze their computation process?

In any case, I sincerely appreciate your thoughtful responses. They have given me many valuable insights to consider.

[Torsten]

Secondly, aren’t q11 and q12 mathematically equivalent? The latter is simply a simplification of the former.

No.

load list_beta.mat
t1=0;
t2=0.13;
t3=1/5;
gamma1=5/3;
gamma2=1/3;
Rplus = @(x)(t2-gamma2/2)./x+(t1+gamma1/2)+t3*x;
Rminus = @(x)t3./x+(t1-gamma1/2)+(t2+gamma2/2)*x;
Rplusdot = @(x)-(t2-gamma2/2)./x.^2+t3;
Rminusdot = @(x)-t3./x.^2+(t2+gamma2/2);
q11 = @(x)Rplus(x)./sqrt(Rplus(x).*Rminus(x));
q12 = @(x)sqrt(Rplus(x)./Rminus(x));
q11(list_beta(1))

ans = -0.0001 - 0.6280i

q12(list_beta(1))

ans = 0.0001 + 0.6280i

You can best see the problem with the complex sqrt with this example. Although z1 is almost equal to z2, their respective sqrt is taken from two different branches:

z1 = exp(-pi*1i)

z1 = -1.0000 - 0.0000i

z2 = exp(pi*1i)

z2 = -1.0000 + 0.0000i

sqrt(z1)

ans = 0.0000 - 1.0000i

sqrt(z2)

ans = 0.0000 + 1.0000i

Finally, what is the purpose of including the following code in the response?

curvelength = cumsum(sqrt( (real(list_beta(2:end))-real(list_beta(1:end-1))).^2 + ...
                           (imag(list_beta(2:end))-imag(list_beta(1:end-1))).^2 ));
curvelength = [0;curvelength];
plot(curvelength,[real(integrand11),imag(integrand11)])
hold on
plot(curvelength,[real(integrand12),imag(integrand12)])
grid on

I was just interested to see the function you try to integrate over C_beta and whether the curves agree despite the difference after evaluating q as q11 and q12, respectively (see above).

[ma Jack]

Once again, thank you, sir, for your responses. They have been very enlightening and helpful.

I believe I’ve found the true cause of the error in the integration: the complex square root. This leads to two values with opposite signs, and we need to choose the correct root to ensure that the q-curve is continuous on the complex plane. Ensuring the continuity of the q-curve is key to making the integral correct.

Answer 2

[David Goodmanson]

Hi MJ,

The result depends on what goes on in the complex plane, so as an intro here are some plots of the angle of functions f(z) in the complex plane, where -pi < angle <= pi, as shown in the colorbar. The colors are cyclic, so going around a zero or a pole gives a plot with no discontinuities (code for this kind of plot is below). Note that for z, with a zero at the origin, the angles go from green to blue etc. as you go counterclockwise. For 1/z, with a pole at the origin, the angles go in the other direction, from green to yellow etc. as you go counterclockwise.

The lower plots of +sqrt(z) and -sqrt(z) each have a branch point at 0, and a branch cut along the negative x axis. Half the angles of z are contained in the sqrt plot and the other half in the -sqrt plot. These cuts are actually not a concern when q'/q is the integrand, since in the expression (2) below, the integrand is just a sum of contributions from poles.

For the problem at hand, writing z in place of beta, you have

(1/(2*pi*i)) Int dq/q = (1/(2*pi*i)) Int (dq/dz)/q dz

for a closed prescribed path in the z plane. (Here the i is in the denominator rather than the numerator as in your expression, to agree with convention for the argument theorem. Later the result can be corrrected by a factor of -1 to make up for this change). Due to the form of the integrand you can make use of a variant of the argument theorem to show that for this square root situation, the value of the integral runs from -1 to 1 in steps of 1/2, depending on the path. Let

tm = (t2-gamma2/2),      tp = (t1+gamma1/2)
[t3 tp tm]    coefficients of (R+)       (reversed order follows Matlab convention)
[tp tm t3]    coefficients of (R-)       (same)

and

q = sqrt(R+/R-) = sqrt ( (t3*z + tp + tm/z)/(tp*z + tm + t3/z) ) 
                = sqrt ( (t3*z^2 + tp*z + tm)/(tp*z^2 + tm*z + t3) )

Using

R+:  rp1,rp2 = roots( [t3 tp tm] )   R-:  rm1,rm2 = roots( [tp tm t3] )                 (1)
             = -4.2102, 0.0435                    = 0.0220 + 0.4894i,  0.0220 - 0.4894i
then
q = sqrt(t3/tp) (z-rp1)^n (z-rp2)^n (z-rm1)^-n (z-rm2)^-n             n =1/2

Lots of sqrt-type factors, but when q’/q is calculated, the constant in front goes away and the result is the sum of four terms, the first of which is

( d/dz (z-rp1)^n ) / (z-rp1)^n =  n / (z-rp1)

and the entire result is the sum of four poles in the complex plane,

q’/q =  n/(z-rp1) +n/(z-rp2) -n/(z-rm1)  -n/(z-rm2)                 n = 1/2          (2)

Each enclosed pole in the integration path gives a contribution (2*pi*i)/(2*pi*i)(+- n) = (+-n), so enclosing each R+ pole gives a contribution of 1/2 and enclosing each R- pole gives a contribution of -1/2 (pending the sign correction)

With the values given, the next two plots are angle and log magnitude of q'/q as in (2). Actually I did not use (2) but to ensure that everything works I used the full expression (denoted as method 1 in the code below) and compared the results with (2) to make sure they agreed, which they did. You can see the poles at the locations of the roots, and also a couple of zeros. Your path of integration is overlaid. It encloses three poles and is not close to enclosing the one on the left, where n = 1/2. The integral will be 1/2-1/2-1/2 = -1/2. (or 1/2 after the sign correction is made). The integral would be zero if all four poles were enclosed but they are not.

angle plot:

log magnitude plot:

Here is a closeup of the angle plot. The colors go green to yellow around the poles and green to blue around the two zeros at approximately (0,-0.4) and (0,0.6)

Explicit numerical integration of q'/q around your path gives a value of (before the sign correction)

I = -0.499999997124216 - 0.000000000000000i'

Since the function is analytic except for poles, any continuous path such as a circle can be used as well and will give the same integral as long as it encloses just the three poles on the right.

load('list_beta')
lisplot = list_beta(1:100:end);
% q'/q
t1 = 0;
t2 = 0.13;
t3 = 1/5;
gamma1 = 5/3;
gamma2 = 1/3;
tm = (t2-gamma2/2);
tp = (t1+gamma1/2);
polyan = [tp tm t3];
rooan = roots(polyan)
polyap  = [t3 tp tm];
rooap = roots(polyap)

xx = linspace(-6,3,501);  % to make the expanded version change the
yy = linspace(-3,3,501);  % limits to -1,1 here 
[x y] = meshgrid(xx,yy);
z = x+i*y;
% method 1   full expression
rap =  t3*z + tp + tm./z;
ran =  tp*z + tm + t3./z;
q = rap.^(1/2).*ran.^(-1/2);
dqdz = ( 1/2)*rap.^(-1/2).*ran.^(-1/2).*(t3 -tm./z.^2) ...
    +  (-1/2)*rap.^( 1/2).*ran.^(-3/2).*(tp -t3./z.^2);
f1 = dqdz./q;
% method 2
f2 =   (1/2)*(1./(z-rooap(1)) + 1./(z-rooap(2))) ...
    + (-1/2)*(1./(z-rooan(1)) + 1./(z-rooan(2)));

figure(2)
zaplot(z,f2)
hold on
plot3(real(lisplot),imag(lisplot),5*ones(size(real(lisplot))))
hold off
fig(3)
zmplot(z,f2)

% contour integration
z = list_beta;
z = [z;z(1)];
rap =  t3*z + tp + tm./z;
ran =  tp*z + tm + t3./z;
q = rap.^(1/2).*ran.^(-1/2);
dqdz = ( 1/2)*rap.^(-1/2).*ran.^(-1/2).*(t3 -tm./z.^2) ...
    +  (-1/2)*rap.^( 1/2).*ran.^(-3/2).*(tp -t3./z.^2);
f1 = dqdz./q;
f2 =   (1/2)*(1./(z-rooap(1)) + 1./(z-rooap(2))) ...
    + (-1/2)*(1./(z-rooan(1)) + 1./(z-rooan(2)));
format long
I1 = (1/(2*pi*i))*contourint(z,f1)
I2 = (1/(2*pi*i))*contourint(z,f2)  % same answer
format

%-----------------------------
function zaplot(z,f)
% plot the angle of a complex function f(z) in the z plane
% z = x + i*y is a complex matrix, usually formed by meshgrid, e.g.
% xx = -2:.01:2;
% yy = -2:.01:2;
% [x y] = meshgrid(xx,yy);
% z = x+i*y;

x = real(z);
y = imag(z);
surf(x,y,angle(f),'edgecolor','none')
view(2)
colormap(wheelmap1)
% keep color from autoscaling, so that branch cuts show up
caxis([-pi,pi])
colorbar
end
% -------------------------------------
function zmplot(z,f)
% plot log magnitude of a complex function f(z) in the z plane
% z = x + i*y is a complex matrix, usually formed by meshgrid, e.g.
% xx = -2:.01:2;
% yy = -2:.01:2;
% [x y] = meshgrid(xx,yy);
% z = x+i*y;

x = real(z);
y = imag(z);
surf(x,y,log(abs(f)),'edgecolor','none')
view(2)
colormap(parula)
colorbar
end
% ------------------------------------------
function map = wheelmap1
% simplified version of wheelmap.
% 96x3 colormap in color wheel order.  Uses hsv2rgb.
% bottom to top of colorbar (map row 1 to map row 96) goes
% magenta -> red -> yellow -> green -> cyan -> blue -> magenta
% When plotting f(z) = z in the complex plane, the 2pi discontinuity does
% not appear and the colors go around z=0 in counterclockwise order. 
% When plotting the angle of a complex function, use
% caxis([-pi pi])
% to bypass autoscaling, so that branch cuts show up.

sat = 6/10;
h = linspace(0,1,97)';
h(end) = [];
s = sat*ones(size(h));
v = ones(size(h));
map = hsv2rgb([h s v]);
map = circshift(map,16,1);
end
% ---------------------------------
function I = contourint(z,f)
% complex contour integral Int f(z) dz by simple trapezoidal method.
% path does not have to be closed.
% z is a vector of points for the path in the complex plane
% f(z) is the integrand
%
% I = contourint(z,f)

I = (1/2)*sum((f(1:end-1)+f(2:end)).*diff(z));
end