From Education

Tom's numerical integration of x^2 in C++

// complile with a command like "g++ integrate.cpp"
// run with the command "./a.out 1000 0 1"
//     to integrate y=x^2 using 1000 steps from 0 to 1
//     the theoretical answer is 1/3
	   
#include <iostream>
#include <iomanip>
#include <sstream>
using namespace std;

// This is the function "f" being integrated
double f(double x) {return x*x;}

// This provides a C++ method comparable to the intrinsic C routine atoi
int atoi(char *c) {
	stringstream param;
	int i;
	param << c;
	param >> i;
	return i;
}

// This provides a C++ method comparable to the intrinsic C routine atof
float atof(char *c) {
	stringstream param;
	float f;
	param << c;
	param >> f;
	return f;
}

// Function "init" calculates deltaX and initializes yVals
void init(int numSeg, double startX, double endX, double &deltaX, double yVals[]) {
	int i;
	double x, rangeX = endX - startX;
	deltaX = rangeX/numSeg;
	for(i=0; i <= numSeg; ++i) {
		x = startX + i*rangeX/numSeg;
		yVals[i] = f(x);
	}
}

// Add the area of all the trapazoids together, to estimate the integral
double integrate(double yVals[], double numSeg, double deltaX) {
	int i;
	double area=0;
	for (i=0; i<numSeg; ++i)
		area += (yVals[i] + yVals[i+1])*.5*deltaX;
	return area;
}

int main(int argc, char **argv) {
	int numSeg;
	double startX, endX, deltaX, *yVals;

// Get number of segments, starting and ending values for integral
	if (argc < 4) exit(0);
	numSeg =	atoi(argv[1]);
	startX =	atof(argv[2]);
	endX =		atof(argv[3]);
	if ((numSeg<1) || (startX >= endX)) exit(0);

// Allocate space for function values and calculate them
	yVals = new double[numSeg+1];
	init(numSeg, startX, endX, deltaX, yVals);

// Display result
	cout << "Area under curve y=x^2 from " << startX << " to " << 
		endX << " is " << setprecision (18) << 
		integrate(yVals, numSeg, deltaX) << endl;
	return 0;
}