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

25

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

0sinc(xr)rE(r)dr

где , и . Здесь - модифицированная функция Бесселя второго рода. В моем конкретном случае у меня , и .xR+λ,κ,ν>0Kλ=0,00313κ=0,00825ν=0,33E(r)=r4(λκ2+r2)ν5/2Kν5/2(λκ2+r2)xR+λ,κ,ν>0Kλ=0.00313κ=0.00825ν=0.33

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

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

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

  • Можно ли использовать качестве элементов в описанном одномерном методе Левина?ψ kxkψ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

torbonde
источник
Что такое в вашем интеграле? x
Педро
Любое положительное, реальное число. Я только что обновил свой пост.
torbonde
Если бы вы могли показать некоторый код и ошибки, вероятно, не так уж и сложно решить большинство из них. Конечно, сначала попробуйте внимательно прочитать ошибку и посмотреть, сможете ли вы ее устранить самостоятельно.
Деннис Джаэруддин
Я сделаю комментарий позже сегодня с некоторым кодом и ошибками. Или завтра.
torbonde
Хорошо, я забыл. Но теперь я обновил свой пост примером (я разделил интеграл на две части, явно вычислив ). sinc
torbonde

Ответы:

12

Я написал свой собственный интегратор, 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.g

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

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

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

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

Pedro
источник
Большое спасибо. Я посмотрю на различные методы. Для моих целей ошибка не должна быть такой же маленькой, как стандартная в уравнении integral(я думаю, что 1e-10), но 1.7e + 07 по-прежнему действительно очень велика. Возможно, другое преобразование принесет пользу, как вы упоминаете.
torbonde
@ cimrg.joe: Обратите внимание, что оценка ошибки является оценкой абсолютной ошибки, основанной, среди прочего, на максимальных абсолютных значениях подынтегральной функции. В некоторых крайних случаях возвращаемое значение может быть вполне нормальным. Если вам нужна точность в десять цифр, тогда я настоятельно рекомендую использовать методы типа Левина, о которых я упоминал в конце своего поста.
Педро
Возможно, мне не нужны десять цифр точности, но я думаю, что мне нужно как минимум пять. Может ли ваш метод произвести это?
torbonde
Этот метод не может гарантировать такую ​​точность для вашего интеграла, поскольку значения в правом конце интервала на несколько порядков больше, чем сам интеграл.
Педро
11

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

У вас есть доступ к 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}]]

Участок от х = 0,5 до х = 10

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

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. Я просто обновлю свой вопрос несколькими добавленными вопросами, касающимися статьи @Pedro, на которую ссылаются.
torbonde
Хорошо, как вы говорите, вам придется обойтись с Matlab. Я добавлю еще один ответ об этом.
Андрей Мойлан
5

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

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

Андрей Мойлан
источник
Я думал о реализации метода Левина сам, но я не уверен, что готов принять вызов. Я думаю, что мне нужно немного лучше понять метод. Может быть, я мог бы поговорить с моим советником об этом. В любом случае, причина, по которой я спросил о методах Филона, заключается в том, что они кажутся более простыми в реализации. И так как мне не нужна чрезвычайно высокая точность, но это часть моей магистерской диссертации, трудности весят.
torbonde
Я взглянул на библиотеку chebfun (которая впечатляет) и пример интеграции Levin. Но я не могу заставить его работать. Я фактически отправил вопрос относительно этого здесь .
torbonde
0

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

Джоэл Врум
источник