How can I numerically solve the following integration in matlab?

I have an equation of the following form:

where A,B,C, and q are 3-by-3 matrix and Tr[...] represent trace. And

here b=1e3. The explicit form of A,B,C, q matrix is:

A = [1 0 0; 0 0 0; 0 0 0];
B = 1/(E-h+1i*1e-3);
C = conj(B);

%where
h = [0          -D*f(-x)       -conj(D)*f(-y)
  -conj(D)*f(x)   0           -D*f(x-y)
  -D*f(y)     -conj(D)*f(y-x)  0];
%and
D = 1 - 1i*2/sqrt(3);
f = @(k) 1 + exp(1i*k);

q = [0              D*1i*exp(-1i*x)         0;
  -conj(D)*1i*exp(1i*x) 0           -D*1i*exp(1i*(x-y));
     0            1i*conj(D)*exp(1i*(y-x))  0];

[Walter Roberson]

If B is a 3 x 3 matrix, then is it correct that in (E-h+1i*1e-3)^-1 that the ^-1 indicates inverse?

What is DD? You mention it twice but do not define it.

[Walter Roberson]

If I assume that DD is a scalar, then the expression to be integrated appears to itself have a closed form, but it is so long that I would expect that integration would be very expensive, and the integral seems unlikely to have a closed form.

[Luqman Saleem]

DD is actually D. I forget to change it. And ^-1 indicates inverse of matrix.

[Luqman Saleem]

@Walter Roberson, so it does not have any solution?

[Walter Roberson]

Still calculating.

I found reasonable ways to make the expression simpler. Unfortunately one of the main ways to speed up the calculation for numeric integration encounters a bug in R2020a that would lead to wrong answers.

I have my computer running a level of symbolic integration. I am not expecting that it will be able to handle it; I am expecting that numeric integration will be needed.

[Luqman Saleem]

function out = fun(x,y,E)

b = 0.1;

D = 1 - 1i*2/sqrt(3);
f = @(k) 1 + exp(1i*k);


h = [0          -D*f(-x)       -conj(D)*f(-y);
    -conj(D)*f(x)   0           -D*f(x-y);
    -D*f(y)     -conj(D)*f(y-x)  0];

q = [0              D*1i*exp(-1i*x)         0;
    -conj(D)*1i*exp(1i*x) 0           -D*1i*exp(1i*(x-y));
    0            1i*conj(D)*exp(1i*(y-x))  0];

A = zeros(3,3); A(1,1)=1;

N = -b*exp(b*E)/(exp(b*E)-1)^2;

B = ((E+1i*1e-3)*eye(3)-h)^(-1);
C = conj(B);
adf
term1 = C-B;
tr = trace(A*term1*q*term1);
out = tr*E*N;
end

then I run the following script:

integral3(@fun2,-pi,pi,-pi,pi,-inf,inf)

it gives me the following error:

Error using vertcat
Dimensions of arrays being concatenated are not consistent.

Error in fun2 (line 15)
q = [0              D*1i*exp(-1i*x)         0; (...much more stuff)

Then instead of writing matrices. I symbolically evaluated the trace term tr. And then re-run the integeral3(...) script. This time I get the following error:

Error using  ^ 
Incorrect dimensions for raising a matrix to a power. Check that the matrix is square and the power is a scalar. To perform elementwise matrix powers, use '.^'.

Error in fun2 (line 10)
N = -b*exp(b*E)/(exp(b*E)-1)^2; (... much more)

I tried to put "." in front of ^. But then it gave error again in "tr" term. Can you please help me with it.

[Walter Roberson]

The integral*() series evaluate at multiple locations simultaneously, expecting a result that is the same size. Your code expects to evaluate only at a scalar location.

[Walter Roberson]

Look at the definition of N, with its denominator of (exp(b*E)-1)^2 . That would be 0 if exp(b*E) == 1. log of both sides, and we get b*E == 0 . b is a known non-zero constant, so the solution is E = 0. Which is in range of -inf to +inf .

Therefore unless the situation can be removed by careful consideration of limits, you have a singularity in the integration.

[Walter Roberson]

Your b is 1000. b*E (in N) will be above 710 when E increases above 0.71 . exp(b*E) would then be exp(710) or more. But that overflows to infinity, and that leads to NaN.

Therefore for practical reasons, you cannot calculate the integral when E > 0.709782712893384 unless you use symbolic methods -- and even if you use symbolic methods, there are limits to the symbolic floating point that are probably going to give you problems.

[Walter Roberson]

The discontinuity at E = 0 is not removable: the function goes to +infinity from the right and -infinity from the left.

[David Goodmanson]

Hello Walter/Luqman,

Although there are practical issues because the trace function might vary quickly with E, still I believe that the integral exists. For the function N, the numerator gets large as E -> inf, but the denominator gets large faster. The quantity EN can be written in a better way as

EN = -(bE/4) / sinh(bE/2)^2

and EN --> 0 as E --> +-inf. Eventually sinh will overflow to inifinty and EN will underflow to 0, but by then the function is so small that the undeflow does not meaningfully affect the integration.

The other question is what happens at E= 0. Since h is hermitian, it is not hard to show that (E - h +i*1e-3) is nonsingular and B = (E - h +i*1e-3)^(-1) exists for all E. Same for C. That means that

tr(E) = trace(A*(C-B)*q*(C-B)) is well behaved at E = 0 and has a finite value tr(0) there.

As E--> 0, EN --> -1/(bE). So the total function tr(E) EN goes to --> -tr(0)/(bE). Add and subtract this function, so

tr(E) EN = [ tr(E)EN  +tr(0)/(bE) ]  -tr(0)/(bE).

The function in square brackets now has no singularity at the origin and can be integrated as usual. The second function contains the singularity, but it is an odd function of E, so Integral{-inf,inf} -tr(0)/(bE) = 0. For the E dependence, integrating the function in square brackets gives the answer.

The integration of -tr(0)/(bE) to give zero sounds like a trick, but If you want to do this more formally, ignoring the constants for the moment,

Integral{-X,-eps} 1/E dE + Integral{eps,X} 1/E dE  = -log(X/eps) + log(X/eps) = 0  

and the limit as X--> inf and eps--> 0 is still 0.

[Luqman Saleem]

@David Goodmanson/ Walter Robertson. Sir, Thank you very much for such detailed answers. I have a few more questions, please help.

1. "tr(E) = trace(A*(C-B)*q*(C-B)) is well behaved at E = 0 and has a finite value tr(0) there." In this statement, are you defining a new variable E = A*(C-B)*q*(C-B)? What does tr(0) mean here?
2. If tr(0) = trace (0) then is not it equal to zero? And hence, tr(0)/(bE) won't make any difference on the overall function.

Apart from the above two queries from Sir David, I have a few question about the things that I have been tried in last 2 days.

Instead of integrating over dE from -inf to +inf, I tried to integrate over only E=1 to E=2, just to check if it work or not. (There are no sigularities now). I defined a function "fun1(x,y,E)" (The explicit form of fun1(x,y,E) is given at the very end of this comment). And run the following script:

integral3(@(x,y,E)arrayfun(@fun1,x,y,E),-pi,pi,-pi,pi,1,2);

I get the following error:

Warning: Reached the maximum number of function
evaluations (10000). The result fails the global
error test. 
Warning: The integration was unsuccessful. 

Then the next that I tried was to take positive constant E. And do integration over x and y only using integral2() function. I found when E is in region [4 5.9], I get the same error given above. When E is not in range of [4 5.9], the integration is successful.

To increase Max Fun Evals, the next thing that I tried was to use quad2d() with 'MaxFunEvals' = 1e6. Unfortunetly, I am still getting the same error, Warning: The integration was unsuccessful. It brings me to my 3rd and last question:

3. Is there anything else that I can try?

Thank you very much for your time.

fun1(x,y,E):

function out = fun1(x,y,E)

%=============just some constants==========

DbyJ = 2/sqrt(3);
T = 1e-2;
eta = 1e-3;
b = 1/T;
D = 1+1i*DbyJ;


fk1 = 1+exp(1i*x);
fk2 = 1+exp(1i*y);
fk1k2 = 1+exp(1i*(x-y));

%=============Matrices==========
A = eye(3); A(1,1) = 1;

q = [0, 1i*D*exp(-1i*x), 0 ;
    -1i*conj(D)*exp(1i*x), 0,-1i*D*exp(1i*(x-y));
    0,1i*conj(D)*exp(1i*(y-x)),0];

h = [0 -D*conj(fk1) -conj(D)*conj(fk2);
    -conj(D)*fk1 0 -D*fk1k2;
    -D*fk2 -conj(D)*conj(fk1k2) 0];

B = inv((E-1i*eta)*eye(3) - h);

C = conj(B);

Term1 = A*(B-C)*q*(B-C);

trc = trace(Term1);

EN = -(b*(E)/4)/(sinh(b*(E)/2)^2);

out = trc*EN + trace(0)/(b*E);

[Walter Roberson]

As you might have noticed, integral3 cannot have the maximum function calls adjusted, but quad2d can. I ended up using

rrr = integral(@(E) quad2d(@(x,y) II(E,x,y), -pi, pi, -pi, pi,'failureplot',true,'maxfunevals', 20000),eps,0.1, 'arrayvalued',true)

This was pretty slow, and the integration still failed. The failure plots show problems mostly at the edges, at +/- pi for x and y. I examined the values numerically, and found the numeric problems I posted about earlier.

That is, regardless of what the theory says about the original formula, the practice was that it generates unavoidable singularities.

I did not try David's formulas.

The function in square brackets now has no singularity at the origin and can be integrated as usual.

The function in square brackets still contains something / bE and as E goes to 0, that goes to infinity, leading to a singularity in practice. Perhaps the left term happens to contain an offsettting subtraction of the same expression as E goes to 0, but that would just mean that you are subtracting two infinities, and even if you could show that the hypothetical terms canceled for any finite value, they would not cancel for the infinite value implied by E = 0.

I am not convinced that the approach can be used in any numeric sense, except perhaps by special-casing 0 and emitting a special theoretical value just for that case.

[David Goodmanson]

Hello Luqman and Walter,

[1] I should have used different notation here. Trace(A*(C-B)*q*(C-B)) is a function of E. So if we define

g(E) = trace(A*(C-B)*q*(C-B))

then g(E) is a well-behaved function of E.

[2] g(E) has a finite value g(0) when evaluated at E = 0. g(0) need not equal 0 and is probably not zero.

[3] As to Walter's point, it's true that a straight numerical integration of [ g(E)EN +g(0)/(bE) ], which has two canceling terms proportional to 1/E, is not going to work. You do need an analytic expression near the origin. This idea is not just theory, it definitely works in practice. Here is a simple example with the exponential integral. We have a smooth function f(x) (which is e^x) and wish to integrate f(x) / x through a singularity at x = 0. Very similar to the case at hand.

The idea is to pick a value x0 and define a segment -x0 to x0 that includes the origin. Do numerical integration of f(x) / x on either side of the segment. Within the segment, expand the smooth function in a taylor series,

f(x) = (c0 + c1*x + c2*x^2 + c3*x^3 ...)
% then
f(x)/x = (c0/x + c1 + c2*x + c3*x^2 ...)

[the c coefficents here are in reverse order from those provided by polyval.]

Upon integrating from -x0 to x0 all the odd functions, including c0/x disappear, leaving as the integrand

c1 + c3*x^2 ... 

which is integrated numerically (or could be done analytically in this case).

In practice the idea is to pick x0 large enough so that integrating the 1/x behavior outside the interval works, but small enough so that not too many terms are required in the taylor expansion.

% calculate Ei(5) = integral{-inf,5) exp(x)/x dx
% this integrates through the singularity at x = 0

xup = 5;
x0 = .01;         % +-x0 defines the small interval about x = 0  
x = linspace(-x0,x0,1000);
c = polyfit(x,exp(x),3)           % exp(x) is the smooth function
% eliminate const and x^2 terms and slide the coefficients
%    one place to the right to account for division by x
c(4) = 0; c(2) = 0;  c(end) = [];    
I2 = integral(@(x) polyval(c,x),-x0,x0);
I1 = integral(@f,-inf,-x0);      % lhs
I3 = integral(@f,x0,xup);        % rhs
I = I1+I2+I3
Icompare = expint(-xup)
function y = f(x)
y = exp(x)./x;
end
I =         40.185275355802887
Icompare = -40.185275355803164 - 3.141592653589793i

Matlab defines the exponential integral slightly differently but there is agreement to 11 decimal places.

Now for the case in question, I will take a slightly different approach from before. The function to integrate is

trace(A*(C-B)*q*(C-B))*EN  =  g(E)*(-bE/4)/sinh(bE/2)^2

To simpify matters, let -b/4 = a so we want to integrate

g(E) * (aE)/sinh(2aE)^2

As stated before, this falls off quickly as E --> +- inf. So one can numerically integrate everything but a small segment from -Es to Es that contains the singularity at the origin. Rewrite the integrand as

(1/a) * (g(E)/E) * [(aE)^2 / sinh(2aE)^2]

The function in brackets is well behaved at the origin where its value is 1/4, and can be calculated for all E by trapping for E=0 and setting the value to 1/4. It is also even in E, so it can be carried along and does not affect the basic approach from before. After using polyfit on g(E), the result is

(1/a) * integral {-Es,Es} (c1 + c3*E^2 ...) * [(aE)^2 / sinh(aE)^2]  dE

[again the c coefficents are in reverse order form those provided by polyval.]