GSL Неверные значения с численным интегрированием

У меня есть подинтегральном, который читаетGSL Неверные значения с численным интегрированием

static double Integrand(double k , void * params) 
{ 
long double *p = (long double *)params; 
    long double MassSquared = p[0]; 
    long double Temp = p[1]; 
    long double sign = p[2]; 

    long double EXP = std::exp(-std::sqrt(k*k+MassSquared/(Temp*Temp))); 
    double res = 0; 
    if(sign == -1) 
    { 
    res = k*k*std::log(1-exp(-std::sqrt(k*k+MassSquared/(Temp*Temp)))); 
    } 
    else if(sign == +1) 
    { 
    res = k*k*std::log(1+exp(-std::sqrt(k*k+MassSquared/(Temp*Temp)))); 
    } 
    return res; 
};

Это один находится в классе под названием VEffClass.

Моя интеграция Рутинное бы тогда

long double VEffClass::NumIntInf(long double low, double sign, long double 
MassSquared, long double Temp) 
{ 
//long double res = 0; 
double up = 
std::sqrt(std::log(C_threshold)*std::log(C_threshold)- MassSquared/(Temp*Temp)); 
long double StepSize = (up-low)/C_IntSteps; 
long double par[3] = {MassSquared, Temp, sign}; 


if(Temp == 0) return 0; 

gsl_integration_workspace *w = gsl_integration_workspace_alloc(5e+5); 
double result, error; 
gsl_function F; 
F.function = &VEffClass::Integrand; 
F.params = &par; 

gsl_integration_qags(&F, C_threshold, up, 1e-07, 0, 5e+5, w, &result, 
&error); 
gsl_integration_workspace_free(w); 

return result; 

}

с C_threshold = 1e-3;

Если я запустил NumIntInf (10^-3, +1, 100, 500), результат равен 1,83038. Результат , заданный Maple, составляет ~ 6 * 10^9. Если я использую простой метод Симпсонов 1/3, то он отличается примерно на 1% от результата Maple.

Может ли кто-нибудь указать мне на ошибку?

Спасибо за любые предложения

источник

2015-07-23 MrBaer

Вам необходимо упростить свое объяснение. Запишите точно, какой интеграл вы хотите вычислить и каковы значения параметров, используемые в вашем примере. Кажется, вы интегрируете экспоненциальный распад (std :: log (1 - exp (- ..)), и если это так, вам нужно быть осторожным при использовании адаптивного интегратора. Если интегральный x-диапазон слишком широк по сравнению с диапазон, в котором подынтегральное выражение отличное от нуля, тогда интегратор будет считать, что подынтегральное выражение равно нулю! –

Здравствуйте, Интеграл задается интегралом по функции Integrand, теоретически от 0 до бесконечности, но я вычисляю значение вверх так, чтобы Integrand там становится 10^-3 и все еще падает оттуда, поэтому я могу в принципе поставить его в 0 с точки вверх, и я могу запустить интеграл от 0 (или здесь значение C_threshold = 10^-3) вверх. Значения параметров, используемые в моем примере, даются MassSquared = 100, Temp = 500 и Sign = +1. Есть ли информация, которая вам нужна? – MrBaer

Всякий раз, когда вы задаете вопрос, который включает в себя численные результаты (а также жуки), нужно дать краткий код, который мы можем анализировать и работать, чтобы воспроизвести исходную проблему. Так что да, информация отсутствует. Нам нужно

(1) краткий код (нам не нужен VEffClass!) С интеграцией GSL (вы частично предоставили это, но ваш код по-прежнему довольно сложный для переполнения стека).

(2) краткий код с вашей реализацией правила Симпсона.

(3) Снимок вашего расчета клена.

(2) и (3) являются весьма важными, потому что ваш вопрос возникает из-за того, что методы (2) и (3) не согласуются с (1), которые вы предположили, являются неправильным ответом. Почему вы предположили, что (1) неверна, а не наоборот? У вас есть аналитическое решение для проверки номеров?

Когда я проанализировал в Wolfram Mathematica подынтегральное выражение с указанными в комментарии значениями параметров, я получил ответ, который согласуется с вычислением GSL. Поэтому я не могу воспроизвести вашу проблему и без (2) + (3) Я могу только догадываться, что пошло не так ...

Я приложил снимок моего ноутбука Mathematica.

источник

2015-07-24 14:15:38

Благодарим вас за подробный ответ! Это была не ошибка Programm, а через ваш ответ и после пытаясь воспроизвести его, я обнаружил, что существовал предварительный фактор Temp^4/(2 * Pi^2) отсутствует между двумя подпрограммами. – MrBaer

GSL Неверные значения с численным интегрированием

ответ

Смежные вопросы