Более надежный интегратор cpp

Я пытаюсь создать релятивистское распределение Войта, свертку релятивистского Breit- Распределение Вигнера и функция Гаусса.

MWE:

double relativisticBreitwigner_pdf(double energy, double width, double mass){
  double massSquare = pow(mass,2);
  double widthSquare = pow(width,2);
  double gamma = sqrt(massSquare*(massSquare+widthSquare));
  double k = (2*sqrt(2)*mass*width*gamma)/(M_PI*sqrt(massSquare + gamma) );
  return k/(pow((pow(energy,2)-massSquare),2) + massSquare*widthSquare);
}

double gaussian_pdf(double energy, double sigma, double mass){
  return (1.0/(sigma*sqrt(2*M_PI)))*exp(-(1.0/2.0)*pow((energy-mass)/sigma,2.0));
}

double relativisticVoigt_pdf(double energy, double width, double mass, double sigma, double range=100.0){

  auto f = [&](double dummy) { return ( relativisticBreitwigner_pdf(dummy+mass,width,mass)*gaussian_pdf(energy-dummy,sigma,mass) );};
  boost::math::quadrature::tanh_sinh<double> integrator;
  return integrator.integrate(f,-range,range);
  //  return boost::math::quadrature::trapezoidal(f,-range,range,sqrt(std::numeric_limits<double>::epsilon()),10000);    }

Эта функция relativisticVoigt_pdf(...) правильно создает релятивистское распределение Фойгта, однако в распределении есть много неправильных «провалов», которые связаны со значением, возвращаемым из integrator.integrate(f,-range,range); не быть правильным.

Если я уменьшу диапазон интегрирования, размер/количество этих провалов будет меньше, но тогда диапазон релятивистского распределения Фойгта будет обрезан.

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

Я предполагаю, что проблема связана с ошибками округления в методе интегрирования из-за большого количества задействованных небольших чисел, но это всего лишь предположение. Есть ли более мощный/надежный интегратор, который хорошо работает в этом режиме?

диапазон интеграции tanh_sinh = 100 tanh_sinh range=100

диапазон трапециевидной интеграции = 250 трапециевидный диапазон=250


c++ numerical-methods numerical-integration numerical-computing
person John Cunningham    schedule 11.02.2021    source источник
comment
1. int в двойное преобразование может быть проблемой в первой функции. 2. Есть и другие константы, кроме M_PI, которые вам, вероятно, следует использовать вместо того, чтобы вычислять их самостоятельно. 3. Пожалуйста, добавьте больше форматирования в свой код. Можно писать однострочники, но их сложно отлаживать. 4. Какие входные области падения? Вероятно, хорошей идеей будет вручную рассчитать/отладить, что там должно произойти.   -  person Sugar    schedule 11.02.2021
comment
Пробовали ли вы другие методы интеграции? Даже базовый? Я был бы удивлен, если бы эти провалы были вызваны ошибками округления.   -  person Damien    schedule 11.02.2021
comment
Также, где вы находите данные для розовой линии?   -  person Sugar    schedule 11.02.2021
comment
и, да, как сказал @Damien, вы пытались проверить накопление ошибок точности и значения переполнения / nan после определенных шагов? большие значения могут привести к большим ошибкам, которые вы видите на провалах   -  person Sugar    schedule 11.02.2021
comment
Привет, я не понимаю, о какой проблеме с преобразованием int в double вы говорите, или какие константы, кроме M_PI, есть? В любом случае, с 1-й или 2-й функцией проблем нет, они везде правильные.   -  person John Cunningham    schedule 11.02.2021
comment
Я пробовал оба метода интеграции, которые у меня есть в MWE (один закомментирован), и у них обоих есть похожие проблемы, я прикрепил изображение трапециевидного случая для диапазона = 250, но провалы находятся в разных местах в зависимости от выбранный интегратор и диапазон интегрирования Провал на уровне 544, например, с tanh_sinh и диапазоном = 100 является релятивистским Voigt_pdf(544,2,085,501,1) = 3,48874e-06   -  person John Cunningham    schedule 11.02.2021
comment
Розовая линия взята из TMath::Voigt((энергия-масса),сигма,ширина)   -  person John Cunningham    schedule 11.02.2021

Ответы (1)


arrow_upward
1
arrow_downward

Для тех, кто в будущем столкнется с подобной проблемой, убрав допуск boost::math::quadrature::trapezoidal(...) по умолчанию (sqrt(std::numeric_limits::epsilon()=1.48e- 8, который в моем коде выше я вставил явно, но это значение по умолчанию, если оно не вставлено) до большего значения 1e-6 мне удается заставить трапециевидный интегратор работать во всем диапазоне.

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

Обновлять:

Хотя описанное выше помогло, но не устранило проблему. Чтобы решить эту проблему, вы можете взять параметр L1, возвращаемый трапециевидной формой, как показано ниже:

double error; 
double L1; 
boost::math::quadrature::trapezoidal(f,-range,range,1e-6,10000,&error,&L1);

Затем проверьте, является ли L1 == 0 или очень маленьким, если это так, измените диапазон, пока это не так.

person John Cunningham    schedule 11.02.2021
comment
Создавая больше дистрибутивов, я вижу, что это решение просто делает проблему более редкой, она все еще возникает для определенных значений. Есть ли какое-либо решение, чем правильно решает это? - person John Cunningham; 14.02.2021