I'm attempting to implement the trapezoidal method of calculating the integral of a curve. It's supposed to be the measure of water flow out of a tank's release valve.
The water flow is a function of time in seconds where Y is the rate of discharge and X is time.

Y = 200 * ( (2 / X + 1) - e(raised to the power of) -.693(X - 1)

The idea is to calculate total discharge and compare it to the tanks capacity. My code spits out an abominably huge number fr this, way over the purported 3107 gallon capacity.

My question is, have I made a beginner's mistake in my coding, or are my calculus skills really this rusty?

#include <cmath>
#include <iostream>
#include <iomanip>

using namespace std;

class Trap
{
 private:
 float sidea;//function of time
 float sideb;//function of time
 float width;//representation of time passed
 float area;//water discharged
 public:
 float Y(float);//calculates side lengths
 void Area(float);//durr
 void Width(float nwidth){width = nwidth;};//assings width for sep. runs
 float GetArea(){return area;};
};


int main()
{
 cout << "This program calculates the rate of flow from a 3107g tank\n";
 cout << "over 10000 seconds using the trapezoidal method.\n";
 cout << "The chart below displays the calculated total flow \n";
 cout << "based on 100, 1000, 10000, and 100000 trapezoids.\n\n";

 cout << " Trapezoids  Total Discharge(gallons) Tank Capacity(gallons)\n";

 Trap zoids[100000];
 float Traplimits[5] = {10, 100, 1000, 10000, 100000};

 int trapcount;
 float sum;

 //loops through the calculations with different numbers of trapezoids
 //being calculated for
 for(int i = 0; i < 5; i++)
   {
      trapcount = 0;
      sum = 0;

      //loops through array of trapezoids assigning new width
      //and then calculating area again
      for( int c = 0; c < Traplimits[i]; c++)
      {
           zoids[c].Width(10000/Traplimits[i]);
           zoids[c].Area(trapcount++);
           
      }

      //sums the areas of all calculated trapezoids in this run
      for( int p = 0; p < Traplimits[i]; p++)
           sum += zoids[p].GetArea();

      cout << fixed << showpoint << setprecision(2) << left;
      cout << setw(14) << Traplimits[i];
      cout << setw(27) << sum;
      cout << "3107\n";
      
   }

   cin.get();
   cin.get();
   return 0;

}

float Trap::Y(float time)
{
    return 200 * ( (2/(time + 1)) - exp(-0.693*(time + 1)) );
}

void Trap::Area(float tcount)
{
    // sidea  = Y(time passed until beginning of trapezoid
    sidea = Y(tcount++ * width);

    //sideb = Y(time passed up to end of trapezoid
    sideb = Y(tcount * width);
    
    //formula for area of a trapezoid
    area = 0.5 * (sidea + sideb) * width;
    cout << area << endl;
}

Debug it by throwing a bunch of output statements in the code to see where it deviates from what you expect or use the debugger that probably came with your compiler to monitor variable values as you step through the code.

Hi,

I am amazed at this code. The normal way to do a trapezium integral is simply to loop over the value like this

const int N(500);
const double startTime(0.0);
const double endTime(100.0);
double area(func(startTime)+func(endTime));
area/=2.0;
for(int i=1;i<N-1;i++) 
   {
      area+=func(i);
  }
area*=(endTime-startTime)/N;

What you are doing creating 100000 objects to do that is absolutely beyond me. However, I think you have the code mostly right ! So well done

BUT the problem is MATHEMATICALLY very ill-conditioned. I suspect that is why it was set. I guess that you will use Simpson next and then runge-kutta on the underlying differential equation. (I will comment just about ANY monotonically decreasing/increasing function is ill-conditioned under this integral method.

If you REALLY want to get the right answer with a trapiz. method then do the problem in log coordinates.

Edited 3 Years Ago by Reverend Jim: Fixed formatting

Thanks so much, and yes never done this kind of math before with c++ so i thought there might be some kind of misconception on my part. Thanks for the reply, i appreciate not telling me to use my debugger lol.

This question has already been answered. Start a new discussion instead.