Метод численного интегрирования сложного колебательного интеграла

Мне нужно численно оценить интеграл ниже:

где , и . Здесь - модифицированная функция Бесселя второго рода. В моем конкретном случае у меня , и .x∈R+λ,κ,ν>0Kλ=0,00313κ=0,00825ν=$E(r) = r^4 (\lambda\sqrt{\kappa^2+r^2})^{-\nu-5/2} K_{-\nu-5/2}(\lambda\sqrt{\kappa^2+r^2})$$x \in \mathbb{R}_+$$\lambda, \kappa, \nu >0$$K$$\lambda = 0.00313$$\kappa = 0.00825$$\nu = 0.33$

Я использую MATLAB, и я попробовал встроенные функции integralи quadgk, что дает мне много ошибок (см. Ниже). Естественно, я пробовал и множество других вещей, таких как интегрирование по частям и суммирование интегралов от до .( k + 1 ) x $kx\pi$$(k+1)x\pi$

Итак, есть ли у вас какие-либо предложения относительно следующего метода?

ОБНОВЛЕНИЕ (добавленные вопросы)
Я прочитал статью @Pedro, на которую ссылается, и я не думаю, что это было слишком сложно понять. Однако у меня есть несколько вопросов:

• Можно ли использовать качестве элементов в описанном одномерном методе Левина?ψ $x^k$$\psi_k$
• Могу ли я вместо этого просто использовать метод Филона, поскольку частота колебаний фиксирована?

Пример кода
>> integral(@(r) sin(x*r).*sqrt(E(r)),0,Inf)
Warning: Reached the limit on the maximum number of intervals in use. Approximate
bound on error is 1.6e+07. The integral may not exist, or it may be difficult to
approximate numerically to the requested accuracy.
> In funfun\private\integralCalc>iterateScalarValued at 372
In funfun\private\integralCalc>vadapt at 133
In funfun\private\integralCalc at 84
In integral at 89

ans =

3.3197e+06

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

Чтобы использовать это для вашей проблемы, я сделал следующее:

>> lambda = 0.00313; kappa = 0.00825; nu = 0.33;
>> x = 10;
>> E = @(r) r.^4.*(lambda*sqrt(kappa^2 + r.^2)).^(-nu-5/2) .* besselk(-nu-5/2,lambda*sqrt(kappa^2 + r.^2));
>> sincp = @(x) cos(x)./x - sin(x)./x.^2;
>> f = @(r) sincp(x*r) .* r .* sqrt( E(r) );

Функция fтеперь ваша интеграция. Обратите внимание, что я только что присвоил любое старое значение x.

Чтобы интегрировать в бесконечную область, я применяю подстановку переменных:

>> g = @(x) f ( tan ( pi / 2 * x ) ) .* ( 1 + tan ( pi * x / 2 ).^2 ) * pi / 2;

т.е. интегрирование gот 0 до 1 должно быть таким же, как интегрирование fот 0 до . Разные преобразования могут давать разные результаты качества: математически все преобразования должны давать один и тот же результат, но разные преобразования могут давать более гладкие или более легко интегрируемые s.$\infty$g

Затем я вызываю свой собственный интегратор, quadccкоторый может работать с NaNs на обоих концах:

>> [ int , err , npoints ] = quadcc( g , 0 , 1 , 1e-6 )
int =
  -1.9552e+06
err =
   1.6933e+07
npoints =
       20761

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

Вы также можете посмотреть на более конкретные методы, такие как этот . Это немного сложнее, но определенно правильный метод для этого типа проблемы.

Как указывает Педро, методы типа Левина являются наилучшими методами для решения подобных проблем.

У вас есть доступ к Mathematica? Для этой проблемы Mathematica обнаружит и использует их по умолчанию:

In[1]:= e[r_] := 
 r^4 (l Sqrt[k^2 + r^2])^(-v - 5/2) BesselK[-v - 5/2, l Sqrt[k^2 + r^2]]

In[2]:= {l, k, v} = {0.00313, 0.00825, 0.33};

In[3]:= Block[{x = 10}, 
 NIntegrate[Sinc'[x r] r Sqrt[e[r]], {r, 0, \[Infinity]}, 
  PrecisionGoal -> 3]]

Out[3]= -112494.

Вот график в диапазоне значений х:

In[4]:= ListLinePlot[
 Table[NIntegrate[Sinc'[x r] r Sqrt[e[r]], {r, 0, \[Infinity]}, 
   PrecisionGoal -> 3], {x, .5, 10, 0.1}]]

Вы также можете вручную указать конкретный метод типа Левина, который будет применяться, что в этом случае может привести к небольшому улучшению производительности:

In[5]:= method = {"LevinRule", "Kernel" -> {Cos[r x], Sin[r x]}, 
   "DifferentialMatrix" -> {{0, -x}, {x, 0}}, 
   "Amplitude" -> {(
     3497.878840962873` Sqrt[(
      r^4 BesselK[-2.17`, 
        0.00313` Sqrt[
         0.00006806250000000001` + r^2]])/(0.00006806250000000001` + 
        r^2)^1.415`])/
     x, -((3497.878840962873` Sqrt[(
       r^4 BesselK[-2.17`, 
         0.00313` Sqrt[
          0.00006806250000000001` + r^2]])/(0.00006806250000000001` + 
         r^2)^1.415`])/(r x^2))}, "AdditiveTerm" -> 0};

In[6]:= Block[{x = 10}, 
 NIntegrate[Sinc'[x r] r Sqrt[e[r]], {r, 0, \[Infinity]}, 
  PrecisionGoal -> 3, Method -> method]]

Out[6]= -112495.

Смотрите документацию для подробностей о методах типа Левина в Mathematica .

Если у вас нет доступа к Mathematica, вы можете написать в Matlab метод Левина (или другой специализированный колебательный), как предлагает Педро.

Вы используете библиотеку chebfun для Matlab? Я только что узнал , что содержит реализацию основного метода Levin типа здесь . Реализация написана Олвером (одним из экспертов в области колебательных квадратур). Это не касается особенностей, адаптивного подразделения и т. Д., Но может быть именно тем, что вам нужно для начала.

Рекомендованное Педро преобразование - отличная идея. Вы пытались поиграться с параметрами в функции quadgk в Matlab? Например, используя преобразование Педро, вы можете сделать следующее:
quadgk(f, 0.0+eps, 1.0-eps, 'AbsTol', eps, 'MaxIntervalCount', 100000)
Использование этого дает мне решение:
-2184689.50220729
и занимает всего 0,8 секунды (с использованием значений, упомянутых выше: x = 10). У
Уолтера Гандера и Уолтера Гаутски есть статья по адаптивной квадратуре с Matlab. код, который вы также можете использовать (ссылка здесь )