Blog Entry © Sunday, October 4, 2026, by James Pate Williams, Jr. Riemannian Integration Tests

// Integration1dTests.cpp
// Copyright Saturday, October 3, 2026
// James Pate Williams, Jr. BA, BS, MSwE, PhD
// One-dimensional integration tests
// Our canonical integrands are x^3 exp(x) for
// x = 0 to x = 1 and x^3 sin(x) for x = 0 to x = pi 

#include <iomanip>
#include <iostream>
using namespace std;

const int steps = 1024;

static double f1(double x)
{
	return x * x * x * exp(x);
}

static double f2(double x)
{
	return x * x * x * sin(x);
}

static double g1()
{
	return -2.0 * exp(1.0) + 6;
}

static double g2()
{
	double pi = 4.0 * atan(1.0);

	return -pi * pi * pi * cos(pi) + 6.0 * pi * cos(pi);
}

double TrapezoidalRule(
	int n, double a, double b, double(*fx)(double))
{
	double h = (b - a) / n;
	double s = 0.5 * (fx(a) + fx(b));
	double x = a + h;

	for (int i = 1; i < n; i++)
	{
		s += fx(x);
		x += h;
	}

	return h * s;
}

double SimpsonsRule(
	int n, double a, double b, double(*fx)(double))
{
	double h = (b - a) / n;
	double h2 = 2.0 * h;
	double s = 0.0;
	double t = 0.0;
	double x = a + h;

	for (int i = 1; i < n; i += 2)
	{
		s += fx(x);
		x += h2;
	}

	x = a + h2;

	for (int i = 2; i < n; i += 2)
	{
		t += fx(x);
		x += h2;
	}

	return h * (fx(a) + 4 * s + 2 * t + fx(b)) / 3.0;
}

static void PrintResult(double x)
{
	cout << setw(15) << fixed;
	cout << right << setprecision(12);
	cout << setfill(' ');
	cout << x;
}

static void PrintError(double x, double exact)
{
	cout << '\t';
	cout << setw(15) << fixed;
	cout << right << setprecision(12);
	cout << setfill(' ');
	cout << 100.0 * fabs(x - exact) /
		fabs(exact) << endl;
}

int main()
{
	double gAct1 = g1();
	double gAct2 = g2();
	double pi = 4.0 * atan(1.0);
	double trap1 = TrapezoidalRule(steps, 0.0, 1.0, f1);
	double simp1 = SimpsonsRule(steps, 0.0, 1.0, f1);
	double trap2 = TrapezoidalRule(steps, 0.0,pi, f2);
	double simp2 = SimpsonsRule(steps, 0.0, pi, f2);

	cout << "Functions are: " << endl;
	cout << "f1(x) = x ^ 3 * exp(x)" << endl;
	cout << "f2(x) = x ^ 3 * sin(x)" << endl;
	cout << "Trapezoidal and Simpson Use 1024 Steps" << endl;
	cout << "         Value           % Error";
	cout << endl;
	cout << "trap1 = ";
	PrintResult(trap1);
	PrintError(trap1, gAct1);
	cout << "simp1 = ";
	PrintResult(simp1);
	PrintError(simp1, gAct1);
	cout << "true1 = ";
	PrintResult(gAct1);
	cout << endl;
	cout << "trap2 = ";
	PrintResult(trap2);
	PrintError(trap2, gAct2);
	cout << "simp2 = ";
	PrintResult(simp2);
	PrintError(simp2, gAct2);
	cout << "true2 = ";
	PrintResult(gAct2);
	cout << endl;
	return 0;
}