Интеграция функции Goempertz в C++ glitch

Question

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

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

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

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

Для некурящего, продолжительность жизни может быть рассчитана с помощью: констант b = 0,0005, c = 1,07. И для курильщика продолжительность жизни может быть рассчитана с помощью констант b = 0,0010, c = 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 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)/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; 
} 

Answer 1

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

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

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

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

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

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

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

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

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

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

Answer 2

Как я понимаю, вы сделали ошибку с константами 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; 
}