Интеграция функции Гемперца в глюк C ++

Question

Интеграция функции Гемперца в глюк C ++

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

Функция Гемперца для человека в возрасте 50 лет может быть закодирована как:

exp((-b/log(c))*pow(c,50)*(pow(c,t)-1))

где b а также c являются константами, и нам нужно интегрировать его от 0 до бесконечности (очень большое число), чтобы получить ожидаемую продолжительность жизни.

Для некурящих, продолжительность жизни может быть рассчитана с помощью:
константы b = 0,0005, с = 1,07.
А для курильщика продолжительность жизни может быть рассчитана с
константы b = 0,0010, с = 1,07.

    const double A = 0; // lower limit of integration
const double B = 1000000000000; // Upper limit to represent infinity
const int N = 10000; //# number of steps of the approximationdouble g(double b, double c, double t)  //
{//b and c are constants, t is the variable of integration.
return exp((-b/log(c))*pow(c,50)*(pow(c,t)-1));
}

double trapezoidal(double Bconst, double Cconst)
{
double deltaX = (B-A)/N; //The "horizontal height" of each tiny trapezoid
double innerTrap = 0; //innerTrap is summation of terms inside Trapezoidal rule
for (int i = 0; i <= N; i++)
{
double xvalue;
if (i == 0) // at the beginning, evaluate function of innerTrap at x0=A
{
xvalue = A;
}
else if (i == N) //at the end, evaluate function at xN=B
{
xvalue = B;
}
else //in the middle terms, evaluate function at xi=x0+i(dX)
{
xvalue = A + i * deltaX;
}

if ((i == 0) || (i == N)) //coefficient is 1 at beginning and end
{
innerTrap = innerTrap + 1*g(Bconst, Cconst, xvalue);
}
else // for all other terms in the middle, has coefficient 2
{
innerTrap = innerTrap + 2*g(Bconst, Cconst, xvalue);
}
}
return (deltaX/2)*innerTrap;
}

int main()
{
cout << "years 50 year old nonsmoker lives: " << trapezoidal(0.0005,1.07) << endl;
cout << "years 50 year old smoker lives: " << trapezoidal(0.0010,1.07) << endl;
cout << "difference between life expectancies: " << trapezoidal(0.0005,1.07)-trapezoidal(0.0010,1.07) << endl;
return 0;
}

3

c++function loops numerical-integration

Решение

Другие решения

Как я понял, вы ошиблись с константами B а также N, B — сколько лет человек может прожить с определенной вероятностью и N это шаг интеграции. Следовательно B должно быть относительно маленьким (<100, потому что вероятность того, что человек проживет 50 + 100 лет или более, чрезвычайно мала) и N должно быть как можно больше. Вы можете использовать следующий код для решения вашей задачи

const double A = 0; // lower limit of integration
const double B = 100; // Upper limit to represent infinity
const int N = 1000000; //# number of steps of the approximation

double g(double b, double c, double t)  //
{//b and c are constants, t is the variable of integration.
return exp((-b/log(c))*pow(c,50)*(pow(c,t)-1));
}

double trapezoidal(double Bconst, double Cconst)
{
double deltaX = (B-A)/double(N); //The "horizontal height" of each tiny trapezoid
double innerTrap = 0; //innerTrap is summation of terms inside Trapezoidal rule
double xvalue = A + deltaX/2;
for (int i = 0; i < N; i++)
{
xvalue += deltaX;
innerTrap += g(Bconst, Cconst, xvalue);
}
return deltaX*innerTrap;
}

int main()
{
double smk = trapezoidal(0.0010,1.07);
double nonsmk = trapezoidal(0.0005,1.07);
cout << "years 50 year old nonsmoker lives: " << nonsmk  << endl;
cout << "years 50 year old smoker lives: " << smk << endl;
cout << "difference between life expectancies: " << nonsmk-smk << endl;
return 0;
}

1

Источник

Accepted Answer

Проблема заключается в выборе конечной x-координаты и количества срезов, по которым вы суммируете площадь:

const double A = 0;
const double B = 1000000000000;
const int N = 10000;

double deltaX = (B-A) / N;  //100 million!

Когда вы делаете дискретную интеграцию, как это, вы хотите, чтобы ваш deltaX быть маленьким по сравнению с тем, как меняется функция. Я предполагаю, что функция Гемпертца довольно сильно меняется от 0 до 100 миллионов.

Чтобы исправить это, просто внесите два изменения:

const double B = 100;
const int N = 10000000;

Это делает deltaX == 0.00001 и, кажется, дает хорошие результаты (21,2 и 14,8). Изготовление B больше не меняет окончательный ответ (если вообще), так как значение функции в этом диапазоне по существу 0.

Если вы хотите быть в состоянии выяснить, как выбрать хорошие значения B а также N процесс примерно такой:

За B найти значение x где результат функции достаточно мал (или изменение функции достаточно мало), чтобы игнорировать. Это может быть сложно для периодических или сложных функций.
Начните с малого N цените и вычисляйте ваш результат. Увеличение N в 2 раза (или около того), пока результат не достигнет желаемой точности.
Вы можете проверить, если ваш выбор B действителен, увеличив его и посмотрев, меньше ли изменение результата, чем вы хотите.

Например, мой выбор B а также N были очень консервативны. Они могут быть уменьшены до B = 50 а также N = 10 и до сих пор дают тот же результат 3 значимым цифрам.

1