Интеграция функции Бесселя 0-го порядка с использованием MATLAB

Question

У меня есть простой вопрос. Я пытаюсь оценить несобственный интеграл от 0-го порядка функции Бесселя с использованием Matlab R2012a:

v = integral(@(x)(besselj(0, x), 0, Inf) 

, который дает мне v = 3.7573e + 09. Однако теоретически это должно быть v = 1. Когда я пытаюсь сделать

v = integral(@(l)besselj(0,l), 0, 1000) 

результат: v = 1.0047. Не могли бы вы кратко объяснить мне, что происходит с интеграцией? И как правильно интегрировать функции типа Бесселя?

Answer 1

Сначала я скептически считал, что интеграция по функции Бесселя приведет к конечным результатам. Mathematica/Wofram Alpha однако показал, что результат конечен, но это not for the faint of heart.

Однако, тогда я указал на this site где объясняется, как сделать это правильно, и что значение интеграла должно быть 1.

Я экспериментировал немного, чтобы проверить правильность своих утверждений:

F = @(z) arrayfun(@(y) quadgk(@(x)besselj(0,x), 0, y), z); 

    z = 10:100:1e4; 

    plot(z, F(z)) 

который дал:

так ясно, интеграл действительно кажется сходится к 1. Позор на Альфах Вольфрама!

(Обратите внимание, что это своего рода вводящий в заблуждение сюжет, попробуйте сделать это с помощью z = 10:1e4;, и вы поймете, почему. Но, ну, в любом случае, принцип все равно).

Эта цифра также показывает, какова проблема, с которой вы сталкиваетесь в Matlab; значение интеграла как затухающее колебание около 1 для увеличения x. Проблема в том, что увлажнение очень слабое - как вы можете видеть, мне нужно было пройти весь путь до 10,000, чтобы произвести этот график, тогда как амплитуда колебаний уменьшилась только на 0,5.

Когда вы пытаетесь сделать несобственный интеграл от возни с установкой «MaxIntervalCount», вы получите это:

>> quadgk(@(x)besselj(0,x), 0, inf, 'maxintervalcount', 1e4) 

    Warning: Reached the limit on the maximum number of intervals in use. 
    Approximate bound on error is 1.2e+009. The integral may not exist, or 
    it may be difficult to approximate numerically. 
    Increase MaxIntervalCount to 10396 to enable QUADGK to continue for 
    another iteration. 
    > In quadgk>vadapt at 317 
     In quadgk at 216 

Это не имеет значения, как высоко вы установите MaxIntervalCount; вы будете сталкиваться с этой ошибкой. Подобные вещи также случаются при использовании quad, quadl или аналогичных (эти функции R2012 integral).

Как показано в этом предупреждении и графике, интеграл просто не подходит для точного приближения любым квадратурным методом, реализованным в стандартном MATLAB (по крайней мере, я знаю).

Я считаю, что правильный аналитический вывод, как это делается на физическом форуме, действительно является единственным способом получить результат без необходимости прибегать к специализированным квадратурным методам.

Answer 2

From the docs сделать несобственный интеграл колебательных функций:

q = integral(fun,0,Inf,'RelTol',1e-8,'AbsTol',1e-13) 

в документации пример является

fun = @(x)x.^5.*exp(-x).*sin(x); 

, но я думаю, в вашем случае попробуйте:

q = integral(@(x)(besselj(0, x),0,Inf,'RelTol',1e-8,'AbsTol',1e-13)