Thread: Find Area using trapozeoidal Rule

How many times does the loop in computeFa actually run?

The code doesn't appear to be calculating the area of each trapezoid. Using double deltax = 2.*a/num_steps, the area of the first trapezoid is delta_area = ((f(-a) + f(-a+deltax))/2.)*deltax;.

The total area using integration is 4.*pow(a,3)/e, where e = exp(1.) == 2.7182818284590452353602874713527 .

If a is 2, and numsteps is 100, how many times does your loop execute?

Originally Posted by Kamal Joub

that part was wrong i edited it out there shouldn't be /12.5 , the code above is wrong

I made an example program, that outputs the results instead of plotting them, using a conventional and non-optimized trapezoidal calculation. I'm curious what the differences are with your version.

Code:

// trp.c trapezoidal test program

#include <stdio.h>
#include <stdlib.h>
#include <math.h>

// Define symbolic constant
#define MAX_SIZE 500  // maximum size of arrays
// Structure definitions
typedef struct
{
    // INPUT from the user
    double aStart;  // Initial value of dimension a for plotting
    double aEnd;    // Final value of dimension a for plotting
    double inc;     // incrementation for value of a
    int num_steps;  //  number of steps to determine h (Trapezoidal Rule)
    // Calculated values
    int n;       // number of points to compute
    double a[MAX_SIZE];  // values  of dimension a
    double area[MAX_SIZE];  // area using analytical equation
    double areaTrap[MAX_SIZE]; // area using Trapezoidal Rule
}  SHAPE;

// function prototypes
void getUserInput(SHAPE *);
double getPositiveValue(char *);
void calculateAreaAnalytical(SHAPE *);
void calculateAreaTrapezoid(SHAPE *);
double computeFa(double, double );

/*--------------------------------------------
Function: main
Description:  Overall control of the program.
Gets the input from the user, calculates areas
using analytical solution and using Trapezoidal rule,
and plot the 2 curves.
----------------------------------------------*/
int main()
{
int i;
    SHAPE shape;  // Input and output data
    // Get the user input
    getUserInput(&shape);
    // Calculations
    calculateAreaAnalytical(&shape);
    calculateAreaTrapezoid(&shape);
    // output 
    for(i = 0; i < shape.n; i++)
        printf("%20.12lf %20.12lf\n", shape.area[i], shape.areaTrap[i]);
    return 0;
}

/*----------------------------------------------------------
Function: getUserInput
Parameters:
    pPtr - reference to SHAPE structure variable.  Members used
            num_steps - number of steps for determining h
            aStart - starting value for a
            aEnd - end value for a
            inc - incrementation value for dimension a
            n - number of elements used in time/area arrays
Description: Gets from the user values for range of values for
             a, an incremetation value for and the number of steps
             to be used with the Trapezoidal rune (determines h)
             and stores in appropriate variables.
             Ensures that aStart is less than aEnd. Computes n
             and ensures that it is less than MAX_SIZE.
-------------------------------------------------------------*/
void getUserInput(SHAPE *pPtr)
{
    // Get input
    pPtr->num_steps = (int)getPositiveValue("number of steps for defining h");
    do
    {
        pPtr->aStart = getPositiveValue("the start value for a");
        pPtr->aEnd = getPositiveValue("the end value for a");
        pPtr->inc = getPositiveValue("incrementation of a");
        pPtr->n = (int)(1+(pPtr->aEnd - pPtr->aStart)/pPtr->inc);
        if(pPtr->n > MAX_SIZE)
            printf("Incrementation too small for range of dimension a (%d)\n", pPtr->n);
        if(pPtr->n < 0)
            printf("The end value for a must be larger than its start value\n");
    }
    while(pPtr->n > MAX_SIZE || pPtr->n <= 0);
    printf("\n");
}

/*----------------------------------------------------------
Function: getPositiveValue
Parameters:
    prompt - reference to string to include in the user prompt
Returns
    value: positive value obtained from the user.
Description: Prompt the user for a value (using the prompt string)
    and checks that the value is positive.
-------------------------------------------------------------*/
double getPositiveValue(char *prompt)
{
    double value; // Value entered by the user.
    do
    {
        printf("Please enter a value for %s: ", prompt);
        scanf("%lf",&value);
        if(value <= 0.0)
            printf("The value must be greater than zero.\n");
    }
    while(value <= 0.0);
    return(value);
}

/*----------------------------------------------------------
Function: calculateAreaAnalytical
Parameters:
    sPtr - reference to SHAPE structure variable.  Members used
            aStart - starting value for a
            inc - incrementation value for dimension a
            n - number of elements used in time/area arrays
            a - array for saving values of a
            area - array for saving area values
Description: Fills in the arrays with n points of
    a/fA(a) values using the analytical solution for the
    area fA(a).
-------------------------------------------------------------*/
void calculateAreaAnalytical(SHAPE *sPtr)
{
    double aValue;
    int ix;
    for(ix = 0; ix < sPtr->n; ix = ix + 1)
    {
        aValue = sPtr->aStart + ix * sPtr->inc;
        sPtr->a[ix] = aValue;
        sPtr->area[ix] = 4.0*pow(aValue,3)/exp(1);
        aValue = aValue + sPtr->inc;
    }
}

/*----------------------------------------------------------
Function: calculateAreaTrapezoid
Parameters:
    sPtr - reference to SHAPE structure variable.  Members used
            num_steps - number of steps for determining h
            aStart - starting value for a
            aEnd 0 ending value of a
            inc - incrementation value for dimension a
            n - number of elements used in time/area arrays
            a - array for saving values of a
            areaTrap - array for saving area values
Description: Fills in the arrays with n points of
    a/area values using the Trapezoidal rule for the
    distance.
-------------------------------------------------------------*/
void calculateAreaTrapezoid(SHAPE *sPtr)
{
int i, j;
double a;
double x0, x1;   // start, end x values for each trapezoid
    for(j = 0; j < sPtr->n; j++){
        a = sPtr->aStart + j * sPtr->inc;
        sPtr->areaTrap[j] = 0.;
        // non-optimized trapezoidal calculation
        for(i = 1; i < sPtr->num_steps; i++){
            x0 = -a + 2. * a *(i-1) / sPtr->num_steps;
            x1 = -a + 2. * a *(i  ) / sPtr->num_steps;
            sPtr->areaTrap[j] += (computeFa(a,x0) + computeFa(a,x1)) / 2. * (x1-x0);
        }
    }
}

/*----------------------------------------------------------
Function: computeFa
-------------------------------------------------------------*/
double computeFa(double a, double x)
{
    return((a*a - x*x)*exp(-x/a));
}

I needed a bit more clean up, and it's too late to edit prior post. So here's another example. I added Simpson rule to the code for comparison.

Code:

// cpsmptrp.c - simpson / trapezoid test program

#include <stdio.h>
#include <stdlib.h>
#include <math.h>

// Define symbolic constant
#define MAX_SIZE 500  // maximum size of arrays
// Structure definitions
typedef struct
{
    // INPUT from the user
    double aStart;  // Initial value of dimension a for plotting
    double aEnd;    // Final value of dimension a for plotting
    double aInc;    // incrementation for value of a
    int n;          // number of points to compute
    int num_steps;  // number of steps to integrate
    // Calculated values
    double a[MAX_SIZE];        // values of a
    double area[MAX_SIZE];     // area using analytical equation
    double areaSimp[MAX_SIZE]; // area using Simpson Rule
    double areaTrap[MAX_SIZE]; // area using Trapezoidal Rule
}SHAPE;

// function prototypes
void getUserInput(SHAPE *);
double getPositiveValue(char *);
void calculateAreaAnalytical(SHAPE *);
void calculateAreaSimpson(SHAPE *);
void calculateAreaTrapezoid(SHAPE *);
double computeF(double, double );

/*--------------------------------------------
Function: main
Description:  Overall control of the program.
Gets the input from the user, calculates areas
----------------------------------------------*/
int main()
{
int i;
    SHAPE shape;  // Input and output data
    // Get the user input
    getUserInput(&shape);
    // Calculations
    calculateAreaAnalytical(&shape);
    calculateAreaSimpson(&shape);
    calculateAreaTrapezoid(&shape);
    // output 
    for(i = 0; i < shape.n; i++)
        printf("%12.4lf %20.12lf %20.12lf %20.12lf\n",
               shape.a[i], shape.area[i], shape.areaSimp[i], shape.areaTrap[i]);
    return 0;
}

/*----------------------------------------------------------
Function: getUserInput
Parameters:
    pPtr - reference to SHAPE structure variable.  Members used
            num_steps - number of steps for integration
            aStart - starting value for a
            aEnd - end value for a
            aInc - increment for a
            n - number of elements of a[] used 
Description: Gets from the user values for range of values for
             a, an incremetation value for and the number of steps
             to be used with the Trapezoidal rune (determines h)
             and stores in appropriate variables.
             Ensures that aStart is less than aEnd. Computes n
             and ensures that it is less than MAX_SIZE.
-------------------------------------------------------------*/
void getUserInput(SHAPE *pPtr)
{
    // Get input
    pPtr->num_steps = (int)getPositiveValue("number of steps for defining h");
    do
    {
        pPtr->aStart = getPositiveValue("the start value for a");
        pPtr->aEnd = getPositiveValue("the end value for a");
        pPtr->aInc = getPositiveValue("incrementation of a");
        pPtr->n = (int)(1+(pPtr->aEnd - pPtr->aStart)/pPtr->aInc);
        if(pPtr->n > MAX_SIZE)
            printf("Incrementation too small for range of dimension a (%d)\n", pPtr->n);
        if(pPtr->n < 0)
            printf("The end value for a must be larger than its start value\n");
    }
    while(pPtr->n > MAX_SIZE || pPtr->n <= 0);
    printf("\n");
}

/*----------------------------------------------------------
Function: getPositiveValue
Parameters:
    prompt - reference to string to include in the user prompt
Returns
    value: positive value obtained from the user.
Description: Prompt the user for a value (using the prompt string)
    and checks that the value is positive.
-------------------------------------------------------------*/
double getPositiveValue(char *prompt)
{
    double value; // Value entered by the user.
    do
    {
        printf("Please enter a value for %s: ", prompt);
        scanf("%lf",&value);
        if(value <= 0.0)
            printf("The value must be greater than zero.\n");
    }
    while(value <= 0.0);
    return(value);
}

/*----------------------------------------------------------
Function: calculateAreaAnalytical
Parameters:
    sPtr - reference to SHAPE structure variable.  Members used
            aStart - starting value for a
            aInc - incrementation value for dimension a
            n - number of elements used in time/area arrays
            a - array for saving values of a
            area - array for saving area values
Description: Fills in area[] using analytical solution
-------------------------------------------------------------*/
void calculateAreaAnalytical(SHAPE *sPtr)
{
    double a;
    int i;
    for(i = 0; i < sPtr->n; i = i + 1)
    {
        a = sPtr->aStart + i * sPtr->aInc;
        sPtr->a[i] = a;
        sPtr->area[i] = 4.0 * pow(a,3) / exp(1);
    }
}

/*----------------------------------------------------------
Function: calculateAreaSimpson
Parameters:
    sPtr - reference to SHAPE structure variable.  Members used
            num_steps - number of steps for determining h
            aStart - starting value for a
            aEnd - ending value of a
            aInc - incrementation value for dimension a
            n - number of elements used in time/area arrays
            a - array for saving values of a
            areaSimp - array for saving area values
Description: Fills in areaSimp[] using Simpson rule
-------------------------------------------------------------*/
void calculateAreaSimpson(SHAPE *sPtr)
{
int i, j;
double a;
double x0, x1, xm;  // start, end, mid x values for each trapezoid
double f0, f1, fm;
    for(j = 0; j < sPtr->n; j++){
        a = sPtr->aStart + j * sPtr->aInc;
        sPtr->areaSimp[j] = 0.;
        x1 = -a;
        f1 = computeF(a,x1);
        for(i = 1; i <= sPtr->num_steps; i++){
            x0 = x1;
            f0 = f1;
            x1 = -a + 2. * a * i / sPtr->num_steps;
            xm = (x0 + x1) / 2.;
            f1 = computeF(a,x1);
            fm = computeF(a,xm);
            sPtr->areaSimp[j] += (f0 + 4. * fm + f1) / 6. * (x1-x0);
        }
    }
}

/*----------------------------------------------------------
Function: calculateAreaTrapezoid
Parameters:
    sPtr - reference to SHAPE structure variable.  Members used
            num_steps - number of steps for determining h
            aStart - starting value for a
            aEnd - ending value of a
            aInc - incrementation value for dimension a
            n - number of elements used in time/area arrays
            a - array for saving values of a
            areaTrap - array for saving area values
Description: Fills in the areaTrap[] using Trapezoidal rule
-------------------------------------------------------------*/
void calculateAreaTrapezoid(SHAPE *sPtr)
{
int i, j;
double a;
double x0, x1;   // start, end x values for each trapezoid
double f0, f1;
    for(j = 0; j < sPtr->n; j++){
        a = sPtr->aStart + j * sPtr->aInc;
        sPtr->areaTrap[j] = 0.;
        x1 = -a;
        f1 = computeF(a,x1);
        for(i = 1; i <= sPtr->num_steps; i++){
            x0 = x1;
            f0 = f1;
            x1 = -a + 2. * a * i / sPtr->num_steps;
            f1 = computeF(a,x1);
            sPtr->areaTrap[j] += (f0 + f1) / 2. * (x1-x0);
        }
    }
}

/*----------------------------------------------------------
Function: computeF
-------------------------------------------------------------*/
double computeF(double a, double x)
{
    return((a*a - x*x)*exp(-x/a));
}