# Fredholm's determinant

1 ビュー (過去 30 日間)
Bogdan MP 2021 年 9 月 19 日
コメント済み: Bogdan MP 2021 年 9 月 24 日
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 2021 年 9 月 24 日
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'.
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 2021 年 9 月 24 日
Yes, the problem was with the limits.

サインインしてコメントする。

R2018a

### Community Treasure Hunt

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

Start Hunting!

Translated by