Fredholm's determinant
1 Ansicht (letzte 30 Tage)
Ältere Kommentare anzeigen
Bogdan MP
am 19 Sep. 2021
Kommentiert: Bogdan MP
am 24 Sep. 2021
I have to calculate some Fredholm's determinant (FD) for that I found one work (arXiv:0904.1581), where the algorithm for calculating the FD was described. I was trying to repeat the following code from the mentioned work (see page 39 from arXiv:0904.1581):
m = 64;
[w, x] = ClenshawCurtis(0, inf, m);
w2 = sqrt(w);
[xi, xj] = ndgrid(x, x);
K1 = @(x,y) airy((x + y) / 2) / 2;
F10 = det(eye(m) - (w2' * w2) .* K1(xi, xj))
>> F10 = 0.831908066202953
and
KAi = @AiryKernel;
F20 = det(eye(m) - (w2' * w2) .* KAi(x, x))
>> F20 = 0.969372828355262
The function ClenshawCurtis() has the following code:
function [w,c] = ClenshawCurtis(a, b, m)
m = m - 1;
c = cos((0 : m) * pi / m);
M = [1 : 2 : m-1]'; l = length(M); n = m - l;
v0 = [2 ./ M ./ (M-2); 1 / M(end); zeros(n, 1)];
v2 = -v0(1 : end - 1) - v0(end : -1 : 2);
g0 = -ones(m, 1); g0(1 + l) = g0(1 + l) + m; g0(1 + n) = g0(1 + n) + m;
g = g0 / (m ^ 2 + mod(m, 2));
w = ifft(v2 + g); w(m + 1) = w(1);
c = ((1 - c) / 2 * a + (1 + c) / 2 * b)';
w = ((b - a) * w / 2)';
end
And the expression for the @AiryKernel is given in arXiv:0904.1581 on page 2, formula (1.6). In could be defined as
function fun = test(x, y)
eps = 0.00001;
if(x ~= y)
fun = ( airy(x) .* airy(1, y) - airy(y) .* airy(1, x) ) ./ (x-y);
else
x = y + eps;
fun = ( airy(x) .* airy(1, y) - airy(y) .* airy(1, x) ) ./ (x - y);
end
end
When I try to run the same code, I get F10 = Nan and F20 = Nan.
What could be the problem?
0 Kommentare
Akzeptierte Antwort
Ashutosh Singh Baghel
am 24 Sep. 2021
Hi Bogdan,
I understand that the problem is related to the inputs you provided to the function 'clenshawCurtis' as 'a' = 0 and 'b' = 'inf'. Please try to provide a scalar number to 'b'.
Please refer to the file on the MATLAB Central File Exchange which performs adaptive Clenshaw-Curtis quadrature. Download the file located at the following URL:
Note that MathWorks does not guarantee or warrant the use or content of these submissions. Any questions, issues, or complaints should be directed to the contributing author.
Weitere Antworten (0)
Siehe auch
Kategorien
Mehr zu Particle & Nuclear Physics finden Sie in Help Center und File Exchange
Community Treasure Hunt
Find the treasures in MATLAB Central and discover how the community can help you!
Start Hunting!