Fredholm's determinant

3 次查看(过去 30 天)
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?

采纳的回答

Ashutosh Singh Baghel
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.
  1 个评论
Bogdan MP
Bogdan MP 2021-9-24
Yes, the problem was with the limits.
Thank you for the answer!

请先登录,再进行评论。

更多回答(0 个)

类别

Help CenterFile Exchange 中查找有关 Particle & Nuclear Physics 的更多信息

产品


版本

R2018a

Community Treasure Hunt

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

Start Hunting!

Translated by