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.