Численная оценка сильно колебательного интеграла

11

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

I(λ)=cos(λcosx)sinxxdx

должен быть аппроксимирован для больших значений с использованием метода седловой точки в комплексной плоскости.λ

Из-за своей высокой колебательной природы этот интеграл очень трудно оценить, используя большинство других методов. Это два фрагмента графика подынтегрального выражения для в разных масштабах:λ=10

cos (10 cos (x)) sin (x) / x

Асимптотическое приближение ведущего порядка

I1(λ)=cos(λ14π)2πλ

и дальнейшее (гораздо меньшее) уточнение добавляет термин

I2(λ)=18sin(λ14π)2πλ3

График аппроксимированных значений как функции выглядит следующим образом:λ

Я (лямбда) ок

Теперь возникает мой вопрос: чтобы наглядно увидеть, насколько хорошее приближение, я хотел бы сравнить его с «реальным значением» интеграла или, точнее, с хорошим приближением к тому же интегралу с использованием независимого алгоритма. Из-за малости коррекции субаренды, я ожидаю, что это будет очень близко.

Я попытался оценить интеграл для некоторого используя другие алгоритмы, но с очень небольшим успехом: Mathematica и Matlab, использующие числовой интегратор по умолчанию, не могут произвести значащее значение (и сообщают об этом явно), mpmath, используя оба в два раза экспоненциальный Подстановка и метод Гаусса-Лежандра дают очень шумные результаты, хотя он имеет небольшую тенденцию колебаться вокруг значений, которые дает метод седловой точки, как может показать этот график:λtanh(sinh)

приблизительно

Наконец, я попытал счастья с интегратором Монте-Карло, используя пример важности, который я реализовал, но мне также не удалось получить стабильных результатов.

Кто-нибудь имеет представление о том, как этот интеграл может быть независимо оценен для любого фиксированного значения или около того?λ>1

doetoe
источник
Функция четная?
Никогуаро
Да, это даже
Doetoe
Вы пытались превратить свой интеграл в ODE?
Никогуаро
1
Нет, дифференцирование по а затем численное решение дифференциального уравнения. x
Никогуаро
1
Ваш первый сюжет, кажется, показывает другую функцию, чем ваш интеграл. А именно, похоже, заменен на . Т.е. сюжет имеет функциюλ x x ( cos ( λ x cos x ) sinc x )λλxx(cos(λxcosx)sincx)
Руслан

Ответы:

12

Используйте теорему Планшереля для оценки этого интеграла.

Основная идея заключается в том, что для двух функций ,f,g

I=f(x)g(x)dx=F(k)G(k)dk

где - преобразования Фурье от . Ваши функции имеют относительно небольшую поддержку в спектральной области. Здесь и должны иметь аналитическое преобразование Фурье (или ряд), как разложение Якоби-Анжера . Вы можете обрезать бесконечные ряды примерно до членов из-за суперэкспоненциального убывания функции Бесселядля, Надеюсь это поможет.F,Gf,gsinx/xrect(k)cos(λcosx)λ | J н ( х ) | n > | х |λ|Jn(x)|n>|x|

Изменить : На самом деле, вы должны использовать здесь представления рядов Фурье вместо преобразований. Путь преобразования приводит к получению асимптотического представления, которое у вас уже есть (оказывается, это просто ). Приведенная выше теорема Планшереля также работает для рядов Фурье с областью интегрирования на последнем интеграле.πJ0(λ)[0,2π]

SMH
источник
Спасибо, это очень хорошая идея!
Doetoe
7

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

πN+π2

асимптотике

Нетрудно догадаться, что асимптотический ряд имеет вид Чтобы численно проверить, что достаточно построить разницу между интегральным и ведущим асимптотическим выражением.

I(λ)2πλ[cos(λπ4)+c1sin(λπ4)λ+c2cos(λπ4)λ2+c3sin(λπ4)λ3+]
c1=18

int := NIntegrate[Cos[l*Cos[x]]*Sinc[x], {x, 0, 20.5*Pi}]; 
Plot[{l*(Sqrt[2*l/Pi]*int - Cos[l-Pi/4]), Sin[l-Pi/4]/8}, {l, Pi/4, 20}]

В результате вы получите довольно хороший синус, который совпадает с тем, который вы получили выше.

18

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

J[l_?NumericQ] := Block[{n=500},
  f[k_] := NIntegrate[Cos[l*Cos[x]]*Sinc[x], {x,0,(n+k)*Pi+Pi/2},
  Method->{"DoubleExponential"}, AccuracyGoal->14, MaxRecursion->100];
  1/2*((f[0]+f[1])/2+(f[1]+f[2])/2)
]
t = Table[{l, l^2*(Sqrt[2*l/Pi]*J[l] - Cos[l-Pi/4] - 1/8*Sin[l-Pi/4]/l)}, 
    {l, 4*Pi+Pi/4, 12*Pi+Pi/4, Pi/36}];
Fit[t, Table[Cos[l-Pi/4+Pi/2*n]/l^n, {n, 0, 10}], l]

Это дает следующий ответ.

c2=9128,c3=751024,c4=367532768,

объяснение

Простой пример

Для иллюстрации я собираюсь использовать более простой пример синуса - интеграл Позвольте мне представить, что меня интересует значение , но я этого не знаю.

S(x)=0xsin(y)ydy.
S()=π2

синус

Вы видите, что колеблется вокруг своего предельного значения аналогично тому, как частичные суммы чередующихся в знаковых рядах колеблются с верхним отсечением. Если вы хотите оценить такую ​​сумму, в соответствии с методом ускорения ряда Эйлера вы должны взять Или с точки зрения синусоидальной функции, вы должны интегрировать ее до точки между максимумом и минимумом колебаний. Как ясно видно из графика, такая точка приблизительно определяется как для больших значений аргумента. В более общем смысле такой точкой является точка, гдепроисходит.S(x)

SN=n=1N(1)nn.
SSN+12(1)N+1N+1.
S(x)0πN+π2sinxxdx
max|S(x)|

Твоя проблема

Возвращаясь к интегралу из курса Константина и Ярослава, вы можете видеть, что он ведет себя точно так же, как синус - интеграл как функция верхнего предела. Это означает, что вам нужно только вычислить значения с . Ниже приведен график нескольких таких значений с .

Ix0(λ)=20x0cos(λcos(x))sinc(x)dx
x0=πN+π2λ=12π

tab = Table[{x0, 2*NIntegrate[Cos[12*Pi*Cos[x]]*Sinc[x], {x, 0, x0}, 
     Method->{"DoubleExponential"}, AccuracyGoal->12, MaxRecursion->100]},
    {x0, 10*Pi+Pi/2, 30*Pi+Pi/2, Pi}];
tab1 = Table[(tab[[i]] + tab[[i+1]])/2, {i,1,Length[tab]-1}];
ListPlot[{tab, tab1}]

акк

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

SN=12(SN+SN+1)
SN

Давид Сайкин
источник
Ницца! Являются ли преподаватели курса вашими настоящими профессорами? Их курс является фантастическим, хотя и очень жестким и быстрым темпом
doetoe
@ doetoe да, я ученик Константина. Он поделился со мной ссылкой на ваш вопрос.
Давид Сайкин
6

Метод Ура для синусоидальных интегралов Фурье работает здесь, см .:

Ура, Такуя и Масатаке Мори. Надежная двойная экспоненциальная формула для интегралов типа Фурье. Журнал вычислительной и прикладной математики 112.1-2 (1999): 229-241.

Я написал реализацию этого алгоритма, но никогда не включал его в работу, чтобы добиться его быстрого (скажем, кэширования узлов / весов), но, тем не менее, я получаю последовательные результаты во всем, кроме точности с плавающей точкой:

float     = 0.0154244
double    = 0.943661538060268
long dbl  = 0.943661538066058702
quad      = 0.943661538066060288748574485677942
oct       = 0.943661538066060288748574485680878906503533004997613278231689169604876
asymptotic= 0.944029734

Вот код:

#include <iostream>
#include <boost/math/quadrature/ooura_fourier_integrals.hpp>
#include <boost/math/constants/constants.hpp>
#include <boost/multiprecision/float128.hpp>
#include <boost/multiprecision/cpp_bin_float.hpp>

template<class Real>
Real asymptotic(Real lambda) {
    using std::sin;
    using std::cos;
    using boost::math::constants::pi;
    Real I1 = cos(lambda - pi<Real>()/4)*sqrt(2*pi<Real>()/lambda);
    Real I2 = sin(lambda - pi<Real>()/4)*sqrt(2*pi<Real>()/(lambda*lambda*lambda))/8;
    return I1 + I2;
}

template<class Real>
Real osc(Real lambda) {
    using std::cos;
    using boost::math::quadrature::ooura_fourier_sin;
    auto f = [&lambda](Real x)->Real { return cos(lambda*cos(x))/x; };
    Real omega = 1;
    Real Is = ooura_fourier_sin<decltype(f), Real>(f, omega);
    return 2*Is;
}

template<class Real>
void print(Real val) {
   std::cout << std::defaultfloat;
   std::cout << std::setprecision(std::numeric_limits<Real>::digits10);
   std::cout <<  val <<  " = " << std::hexfloat << val;
}

int main() {
    using boost::multiprecision::float128;
    float128  s = 7;
    std::cout << "Asymptotic = " << std::setprecision(std::numeric_limits<float128>::digits10) << asymptotic(s) << "\n";
    std::cout << "float precision = ";
    print(osc(7.0f));
    std::cout << "\n";
    std::cout << "double precision= ";
    print(osc(7.0));
    std::cout << "\n";
    std::cout << "long double     = ";
    print(osc(7.0L));
    std::cout << "\n";
    print(osc(s));

    print(osc(boost::multiprecision::cpp_bin_float_oct(7)));
    std::cout << "\n";
}

Вы не можете видеть разницу между квадратурой и асимптотикой, потому что они лежат друг на друге, кроме как : λ0введите описание изображения здесь

user14717
источник
Спасибо, это действительно приятно! У меня еще не получилось, моя буст-установка не совместима, но сейчас я загружаю последнюю версию.
Doetoe
Просто чтобы быть уверенным: в 23 у вас есть cos (лямбда * cos (x)) / x без фактора sin (x) от подынтегрального выражения. Является ли ooura_fourier_sin этим фактором sin (x) для умножения подынтегрального выражения, переданного ему?
Doetoe
Я получил это работает. Он и его зависимости кажутся только заголовочными, поэтому мне даже не пришлось устанавливать или компилировать (кроме исполняемого файла). Я надеюсь, что это будет включено в повышение!
Doetoe
@doetoe: Да, фактор неявен и включен в квадратурные веса. Я рассматриваю возможность ускорить его, сделав рефакторинг для кэширования узлов и весов; посмотрим! sin(x)
user14717
@doetoe: это объединено в Boost 1.71. API немного отличается от этого ответа.
user14717