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?
0 个评论
采纳的回答
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'.
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.
更多回答(0 个)
另请参阅
类别
在 Help Center 和 File Exchange 中查找有关 Particle & Nuclear Physics 的更多信息
Community Treasure Hunt
Find the treasures in MATLAB Central and discover how the community can help you!
Start Hunting!