How can I numerically solve the following integration in matlab?

Question

Luqman Saleem 2020-6-28

0
链接

此问题的直接链接

https://ww2.mathworks.cn/matlabcentral/answers/556030-how-can-i-numerically-solve-the-following-integration-in-matlab

评论： David Goodmanson 2020-7-2

在 MATLAB Online 中打开

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];

14 个评论
显示 12更早的评论隐藏 12更早的评论

Luqman Saleem 2020-6-28

编辑：Luqman Saleem 2020-6-28

在 MATLAB Online 中打开

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.

David Goodmanson 2020-6-29

在 MATLAB Online 中打开

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 2020-7-1

编辑：Luqman Saleem 2020-7-1

在 MATLAB Online 中打开

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

"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?
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 2020-7-1

在 MATLAB Online 中打开

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 2020-7-2

在 MATLAB Online 中打开

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.]

请先登录，再进行评论。

请先登录，再回答此问题。

How can I numerically solve the following integration in matlab?

14 个评论
显示 12更早的评论隐藏 12更早的评论

回答（0 个）

另请参阅

类别

标签

产品

版本

Community Treasure Hunt

How can I numerically solve the following integration in matlab?

14 个评论 显示 12更早的评论隐藏 12更早的评论

回答（0 个）

另请参阅

类别

标签

产品

版本

Community Treasure Hunt

14 个评论
显示 12更早的评论隐藏 12更早的评论