Goempertz function integration in C ++ fails
I'm trying to find the Trapezoidal Rule score of the Goempertz function and use it to measure the difference between life expectancy for a 50-year-old smoker and a 50-year-old non-smoker, but my code gives me shit answers.
The Goempertz function for a person aged 50 can be coded as:
exp((-b/log(c))*pow(c,50)*(pow(c,t)-1))
where b
and c
are constants and we need to integrate it from 0 to infinity (a very large number) to get the life expectancy.
For a non-smoker, life expectancy can be calculated from: constants b = 0.0005, c = 1.07. And for a smoker, life expectancy can be calculated using the constant 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;
}
source to share
The problem lies in the choice of the end x-coordinate and the number of slices that you sum over the area:
const double A = 0;
const double B = 1000000000000;
const int N = 10000;
double deltaX = (B-A) / N; //100 million!
When you do discrete integration, you want yours to deltaX
be small compared to how the function changes. I would guess that the Goempertz function varies quite a bit from 0 to 100 million.
To fix this, just make two changes:
const double B = 100;
const int N = 10000000;
This does deltaX == 0.00001
and appears to give good results (21.2 and 14.8). Creation B
no longer changes the final answer (if at all), since the function value in that range is essentially 0.
If you want to know how to choose good values B
and N
, the process is something like this:
- For,
B
find a valuex
where the function result is small enough (or the function change is small enough) to be ignored. This can be tricky for periodic or complex functions. - Start with a small value
N
and calculate the result. ZoomN
in 2x (or something) until the result converges to the desired precision. - You can check if your selection is valid by
B
increasing it and see if the change in the result is less than your desired precision.
For example, my choice B
and N
was very conservative. They can be reduced to a level B = 50
and N = 10
still give the same result to 3 significant figures.
source to share
As I understand it, you made a mistake with constants B
and N
. B
- the number of years during which a person can live with a certain probability, and N
- the step of integration. Therefore, it B
should be relatively small (<100, since the probability that a person will live 50 + 100 years or more is extremely small), but N
should be as large as possible. You can use the following code to solve your problem.
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;
}
source to share