How can I numerically solve the following integration in matlab?

2 views (last 30 days)
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);
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];
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];
David Goodmanson
David Goodmanson on 2 Jul 2020
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;
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.]

Sign in to comment.

Answers (0)

Community Treasure Hunt

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

Start Hunting!

Translated by