By using this site, you agree to our updated Privacy Policy and our Terms of Use. Manage your Cookies Settings.
446,194 Members | 817 Online
Bytes IT Community
+ Ask a Question
Need help? Post your question and get tips & solutions from a community of 446,194 IT Pros & Developers. It's quick & easy.

Summation Program

P: 23
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:

Expand|Select|Wrap|Line Numbers
  1. /*
  2. trap.cpp - Trapezoid rule.
  3.     To evaluate a definite integral of a function from x = a to x = b using the Trapezoid rule.
  4.     (ref. Washington, pages 749) Note: The function to be integrated is hard-coded into the program.
  5.     User input: a, b, n.
  6.     Output: Calculated definite integral (area).
  7. */
  8.  
  9. #include <iostream>
  10. #include <math.h>
  11. #include <iomanip>
  12. using namespace std;
  13.  
  14. int main()
  15. {
  16.     double trapezoid( double a , double b , long n );    // prototype
  17.  
  18.     double a , b , area;
  19.     long n;
  20.  
  21.     cout << "Lower bound, a> ";
  22.     cin >> a;
  23.     cout << "Upper bound, b> ";
  24.     cin >> b;
  25.  
  26.     cout << "Intervals, n> ";
  27.     cin >> n;
  28.  
  29. //    repeat with various values of n to see its effect.
  30.  
  31.     while ( n != 0 )
  32.     {
  33.         area = trapezoid( a , b , n );
  34.         cout << fixed << setprecision( 12 );
  35.         cout << "\nDefinite integral ( trapezoid ): " << area  << "\n\n";
  36.  
  37.         cout << "Intervals, n (0 to quit)> ";
  38.         cin >> n;
  39.     }
  40.     return 0;
  41. }
  42.  
  43. double myfunc( double y )
  44. {
  45.     return ( 1.0 / y );     // Example 1 on page 750
  46. }
  47.  
  48. double trapezoid( double a , double b , long n )
  49. {
  50.     double myfunc( double y );        // prototype
  51.  
  52.     double h = ( ( b - a ) / n );
  53.     double y = 0;
  54.     double sum = 0;
  55.  
  56.     for ( int i = 0; i <= n; i ++ )
  57.     {
  58.         y += h;
  59.         sum = sum + myfunc( y );
  60.     }
  61.     sum = sum * ( h / 2 );
  62.  
  63.     return sum;
  64. }
  65.  
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.
Mar 14 '07 #1
Share this Question
Share on Google+
14 Replies


100+
P: 1,646
Hi. What result are you getting and what input are you using to get that result?
Mar 14 '07 #2

Ganon11
Expert 2.5K+
P: 3,652
Is myfunc representing the original function or the integral function?
Mar 14 '07 #3

DeMan
100+
P: 1,806
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....
Mar 14 '07 #4

Ganon11
Expert 2.5K+
P: 3,652
Expand|Select|Wrap|Line Numbers
  1. double trapezoid( double a , double b , long n )
  2. {
  3.     double myfunc( double y );        // prototype
  4.  
  5.     double h = ( ( b - a ) / n );
  6.     double y = 0;
  7.     double sum = 0;
  8.  
  9.     for ( int i = 0; i <= n; i ++ )
  10.     {
  11.         y += h;
  12.         sum = sum + myfunc( y );
  13.     }
  14.     sum = sum * ( h / 2 );
  15.  
  16.     return sum;
  17. }
  18.  
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

Expand|Select|Wrap|Line Numbers
  1. for (y = a; y < b - h; y += h) {
  2.    // Find the area
  3.    // Add area to sum
  4. }
Mar 14 '07 #5

P: 23
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
Mar 14 '07 #6

Ganon11
Expert 2.5K+
P: 3,652
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.
Mar 14 '07 #7

Ganon11
Expert 2.5K+
P: 3,652
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.
Mar 14 '07 #8

P: 23
Expand|Select|Wrap|Line Numbers
  1. double trapezoid( double a , double b , long n )
  2. {
  3.     double myfunc( double y );        // prototype
  4.  
  5.     double h = ( ( b - a ) / n );
  6.     double y = h;
  7.     double sum = 0;
  8.  
  9.     for ( int i = 1; i <= n; i ++ )
  10.     {        
  11.         if ( i == a || i == n)
  12.         {
  13.             y += h;
  14.             sum += myfunc( y );
  15.         }
  16.         else
  17.         {
  18.             y = ( y + h );
  19.             sum += myfunc( y )*2;
  20.         }
  21.  
  22.         cout<<sum<<endl;
  23.     }
  24.  
  25.     return sum;
  26. }
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.
Mar 14 '07 #9

Ganon11
Expert 2.5K+
P: 3,652
Expand|Select|Wrap|Line Numbers
  1. double trapezoid( double a , double b , long n )
  2. {
  3.     double myfunc( double y );        // prototype
  4.  
  5.     double h = ( ( b - a ) / n );
  6.     double y = h;
  7.     double sum = 0;
  8.  
  9.     for ( int i = 1; i <= n; i ++ )
  10.     {        
  11.         if ( i == a || i == n)
  12.         {
  13.             y += h;
  14.             sum += myfunc( y );
  15.         }
  16.         else
  17.         {
  18.             y = ( y + h );
  19.             sum += myfunc( y )*2;
  20.         }
  21.  
  22.         cout<<sum<<endl;
  23.     }
  24.  
  25.     return sum;
  26. }
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.
Mar 15 '07 #10

P: 23
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!
Mar 15 '07 #11

Ganon11
Expert 2.5K+
P: 3,652
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!
Mar 15 '07 #12

Ganon11
Expert 2.5K+
P: 3,652
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:

Expand|Select|Wrap|Line Numbers
  1. double part1 = myfunc(y);
  2. double part2 = myfunc(y + h);
inside the for loop.
Mar 15 '07 #13

P: 23
I finally got it!

Expand|Select|Wrap|Line Numbers
  1. double trapezoid( double a , double b , long n )
  2. {
  3.                 double myfunc( double x );                         // prototype
  4.  
  5.                 double h , area = 0;
  6.                 h = ( b - a ) / n;
  7.  
  8.                 for ( int i = a ; i < n; i ++ )
  9.                    area += myfunc( a + i*h );
  10.  
  11.                 area = ( h / 2.0 )*( myfunc( a ) + 2.0*area + myfunc( b ) );
  12.  
  13.                 return area;
  14. }
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.
Mar 15 '07 #14

Ganon11
Expert 2.5K+
P: 3,652
Glad that it's giving you the correct results. I wonder how it will work with different functions?
Mar 15 '07 #15

Post your reply

Sign in to post your reply or Sign up for a free account.