Summation Program

[vpascuzzi]

Hello:

I am trying to make a program that will calculate the integral of a function using the Trapezoidal Rule.

So far I have this:

/*
trap.cpp - Trapezoid rule.
    To evaluate a definite integral of a function from x = a to x = b using the Trapezoid rule.
    (ref. Washington, pages 749) Note: The function to be integrated is hard-coded into the program.
    User input: a, b, n.
    Output: Calculated definite integral (area).
*/
 
#include <iostream>
#include <math.h>
#include <iomanip>
using namespace std;
 
int main()
{
    double trapezoid( double a , double b , long n );    // prototype
 
    double a , b , area;
    long n;
 
    cout << "Lower bound, a> ";
    cin >> a;
    cout << "Upper bound, b> ";
    cin >> b;
 
    cout << "Intervals, n> ";
    cin >> n;
 
//    repeat with various values of n to see its effect.
 
    while ( n != 0 )
    {
        area = trapezoid( a , b , n );
        cout << fixed << setprecision( 12 );
        cout << "\nDefinite integral ( trapezoid ): " << area  << "\n\n";
 
        cout << "Intervals, n (0 to quit)> ";
        cin >> n;
    }
    return 0;
}
 
double myfunc( double y )
{
    return ( 1.0 / y );     // Example 1 on page 750
}
 
double trapezoid( double a , double b , long n )
{
    double myfunc( double y );        // prototype
 
    double h = ( ( b - a ) / n );
    double y = 0;
    double sum = 0;
 
    for ( int i = 0; i <= n; i ++ )
    {
        y += h;
        sum = sum + myfunc( y );
    }
    sum = sum * ( h / 2 );
 
    return sum;
}
 

This unfortunately does not give me the correct answer to the problem. The answer should be [ 67 / 60 ] or [ 1.1166666667 ].

Any help would be greatly appreciated. Thank you in advance.

[willakawill]

Hi. What result are you getting and what input are you using to get that result?

[Ganon11]

Is myfunc representing the original function or the integral function?

[DeMan]

1) what answer are you getting...could it be some cast that the compiler is doing without you knowing? (ie does the result look suspiciously close to an "int" say, when it should be a float.

2) I have never seen prototypes declared inside a method (this may not be incorrect, but I'm intrigued because it's something I've not seen - especially in the case where myfunc appears to be defined before it's prototype)

3) you are always adding the same value to y...I think (although integration was not my favourite math) that each value of h must be different.....that h should be the sum of areas along a given interval, not merely the entire inteval....

[Ganon11]

double trapezoid( double a , double b , long n )
{
    double myfunc( double y );        // prototype
 
    double h = ( ( b - a ) / n );
    double y = 0;
    double sum = 0;
 
    for ( int i = 0; i <= n; i ++ )
    {
        y += h;
        sum = sum + myfunc( y );
    }
    sum = sum * ( h / 2 );
 
    return sum;
}
 

This function appears to be computing the Trapezoidal Rule incorrectly. Specifically, I think your for... loop is incorrect. Currently, it increments y to the next position to be calculated, then adds the function's value at y to sum. But for the Trap. Rule, you need the first and second endpoint of the Trapezoid. You should be finding the area from y to y + h and adding this value to sum.

I'd suggest that, instead of looping i from 1 to n, use y as the index. You can use the header

for (y = a; y < b - h; y += h) {
   // Find the area
   // Add area to sum
}

[vpascuzzi]

The answer I get is 1.0416666666667.

I know there is something wrong in the trapezoid(double,double,long) function, where I am assigning y+=h.

I'm just not sure how to correct it. I'm missing something.

I believe the problem I am facing is that the first y (y sub-0) and the last y (y sub-n) are supposed to be added seperately, while the y values in between the first and last are to be multiplied by 2.

Equation:
Area of trapezoid = h / 2 [ ( y0 + yn ) + Σ yi ]

function to be integrated:
( 1 / x ) dx

where:
limit a = 1;
limit b = 3;
n (# intervals) = 4;
h = ( b - a ) / n

[Ganon11]

From what I remember (learned the trap. rule all the way back in like November XD),

The area under a function is the sum of n trapezoids, conputed with x values of x(sub)0, x(sub)1, x(sub)2...x(sub)n. The area of a trapezoid was:

area = (1 / 2) * h * (b(sub)1 - b(sub)2)

Now, h will be the change in any x(sub)i to x(sub)i+1, or delta x. The b(sub)1 and b(sub)2 values are the values of the function at x(sub)i and x(sub)i+1.

Now, to do this in programming, there are a few things that must be realized. Obviously, you are using a loop of some sort to add the areas together. How many times should the loop execute? If n = 4, then the area between a and b is cut into 4 trapezoids (for your example, the trapezoids are from x = 1 to x = 1.5, then 1.5 to 2, then 2 to 2.5, and 2.5 to 3). The loop is executing the area of 1 trapezoid, so it must execute 4 times - thus, from 1 to <= n will work (however, I think the for... loop I suggested will make things easier for you in the long run).

See if this makes any sense, and come back with any further questions.

[Ganon11]

By the way, I think the answer you are getting is a result of immediately executing "y+=h;" as the first line in the loop - I think this is skipping the trapezoid from 1 to 1.5, and instead calculating the trapezoid from 3 to 3.5, which is a smaller area.

[vpascuzzi]

double trapezoid( double a , double b , long n )
{
    double myfunc( double y );        // prototype
 
    double h = ( ( b - a ) / n );
    double y = h;
    double sum = 0;
 
    for ( int i = 1; i <= n; i ++ )
    {        
        if ( i == a || i == n)
        {
            y += h;
            sum += myfunc( y );
        }
        else
        {
            y = ( y + h );
            sum += myfunc( y )*2;
        }
 
        cout<<sum<<endl;
    }
 
    return sum;
}

This is what I am coming up with. I am trying to distinguish the first and last y-values (y(sub)0 and y(sub)n) from the other y-values because in order to find the summation, I have to add each y-value between y(sub)0 and y(sub)n twice.

So it would be something like:

A = h/2 [ ( y(sub)0 + y(sub)n ) + 2( y(sub)1 + y(sub)2 + y(sub)n-1 ) ]

I am still having problems figuring this out.

[Ganon11]

double trapezoid( double a , double b , long n )
{
    double myfunc( double y );        // prototype
 
    double h = ( ( b - a ) / n );
    double y = h;
    double sum = 0;
 
    for ( int i = 1; i <= n; i ++ )
    {        
        if ( i == a || i == n)
        {
            y += h;
            sum += myfunc( y );
        }
        else
        {
            y = ( y + h );
            sum += myfunc( y )*2;
        }
 
        cout<<sum<<endl;
    }
 
    return sum;
}

OK, I think I see what you're trying to do. Two immediate problems arise:

1) y+= h; and y = (y + h); do the exact same thing, just with different syntax. If you were trying to do something different with the second option, you might try something else.

2) You are checking if i == a in the if...condition, but it should be if i == 1. Remember, i is independent of a and b - it only works with n. You could choose 13 for a and 15 for b with n = 4, and i would never equal a.

Your for... loop should do the following:

1) Find the length of each base of the trapezoid. These will be your myfunc() values. Remember, you will have to explicitly call myfunc() twice - once with y, and again with y + some value.
2) Find the area of this single trapezoid using the formula I provided using the change in y and the two lengths found in (1).
3) Add this area to the total sum.
4) Increment y by some value.

[vpascuzzi]

OK, I think I see what you're trying to do. Two immediate problems arise:

1) y+= h; and y = (y + h); do the exact same thing, just with different syntax. If you were trying to do something different with the second option, you might try something else.

2) You are checking if i == a in the if...condition, but it should be if i == 1. Remember, i is independent of a and b - it only works with n. You could choose 13 for a and 15 for b with n = 4, and i would never equal a.

Your for... loop should do the following:

1) Find the length of each base of the trapezoid. These will be your myfunc() values. Remember, you will have to explicitly call myfunc() twice - once with y, and again with y + some value.
2) Find the area of this single trapezoid using the formula I provided using the change in y and the two lengths found in (1).
3) Add this area to the total sum.
4) Increment y by some value.

Problems:
1) Got it, I knew they did the same thing, I was just tweaking here and there and just wrote it out differently.

2) By setting i == 1 in the for..loop will that work then if the 'a' value is '3'? I was thinking that this would be the filter to determine y(sub)0 and y(sub)n from the others. Maybe I should rethink this?

for..loop:
1) I'm not sure what you mean by calling myfunc() twice EXPLICITY, once with 'y' and once with 'y + <some value>'.

Other than that, I should be able to get it. Just those parts need a little more clarification.


Thanks!

[Ganon11]

Problems:
2) By setting i == 1 in the for..loop will that work then if the 'a' value is '3'? I was thinking that this would be the filter to determine y(sub)0 and y(sub)n from the others. Maybe I should rethink this?

As it stands right now, your i vqariable is indicating which trapezoid you are calculating. When you start, i is 1, and you are calculating the 1st trapezoid. After the first execution, i is 2, and you are calculating the 2nd trapezoid. This continues until i == n, when you are calculating the nth trapezoid. To set apart the y(sub)0 and y(sub)n values, you need to set apart the first and last trapezoids, with happen when i == 1 and i == n. Remember, i is independent of a and b.

I'll reply to your other question when I get the chance!

[Ganon11]

for..loop:
1) I'm not sure what you mean by calling myfunc() twice EXPLICITY, once with 'y' and once with 'y + <some value>'.

I just meant you would have two calls, like:

double part1 = myfunc(y);
double part2 = myfunc(y + h);

inside the for loop.

[vpascuzzi]

I finally got it!

double trapezoid( double a , double b , long n )
{
                double myfunc( double x );                         // prototype
 
                double h , area = 0;
                h = ( b - a ) / n;
 
                for ( int i = a ; i < n; i ++ )
                   area += myfunc( a + i*h );
 
                area = ( h / 2.0 )*( myfunc( a ) + 2.0*area + myfunc( b ) );
 
                return area;
}

For the records, the only part I actually had to program is the italics - our instructor gave us the template. I figured out the best way to calculate the problem was to first summate all of the values between (sub)0 and (sub)n, then double that answer.

Next, I took the individual (sub)0 and (sub)1 values and added them in independantly. Works for me!

Thanks all for your help, and especially Ganon11.

[Ganon11]

Glad that it's giving you the correct results. I wonder how it will work with different functions?