CODEWITHSUNDEEP
c++functionloopsnumerical-integration
shoestringfries
shoestringfries

Reputation: 279

Integrating Goempertz Function in C++ glitch

I'm trying to find the Trapezoidal Rule estimation of the Goempertz function and use it to measure the difference between the life expectancy for a 50 year old smoker and a 50 year old non smoker but my code has been giving me crap answers.

The Goempertz function for a person at age 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, the life expectancy can be calculated with: constants b = 0.0005, c = 1.07. And for a smoker, the life expectancy can be calculated with constants 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;
}

Upvotes: 3

Views: 126

Answers (2)

kvorobiev
kvorobiev

Reputation: 5070

As I understand you made a mistake with the constants B and N. B - the number of years that a person can live with a certain probability and N is integration step. Therefore B should be relatively small (<100, because the probability that a person will live 50+100 years or more is extremely small) and N should be as large as possible. You can use following code to solve your task

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;
}

Upvotes: 1

uesp
uesp

Reputation: 6204

The problem lies in your choice of end x-coordinate and the number of slices you sum the area over:

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

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

When you do a discrete integration like this you want your deltaX to be small compared to how the function changes. I would guess the Goempertz function changes quite a lot between 0 and 100 million.

To fix it simply make the two changes:

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

This makes deltaX == 0.00001 and seems to give good results (21.2 and 14.8). Making B larger doesn't change the final answer much (if at all) as the function value at this range is essentially 0.

If you'd like to be able to figure out how to choose good values of B and N the process is roughly like this:

  1. For B find the value of x where the function result is small enough (or the function change in is small enough) to ignore. This can be tricky for periodic or complex functions.
  2. Start with a small N value and compute your result. Increase N by a factor of 2 (or something) until the result converges to the desired accuracy.
  3. You can check if your choice of B is valid by increasing it and see if the change in result is less than your desired accuracy.

For example, my choices of B and N were very conservative. These could be reduced as low as B = 50 and N = 10 and still give the same result to 3 significant figures.

Upvotes: 1

Related Questions