Numerical Integration

// main.cpp
#include "integrate.h"
#include <iostream>
#include <cassert>
#include <iomanip>
#include <string>
#include <cmath>
using namespace std;

const double EPS = 1e-6;
const double R = 4.0;

const double a = -4.0, b = 4.0;

// Вычисление значения функций f(x) и ее первообразной F(x)

double F (double x) { return 0.5 * (R * R * asin(x / R) + x * sqrt(R * R - x * x)); }
double f (double x) { return sqrt(R * R - x*x); }

// Правило Рунге оценки погрешности

void Runge(int& n, double& I1, double& I2, Method method)
{
	if (fabs(I1 - I2) >= EPS)
	{
		I2 = I1;

		n *= 2;

		I1 = method(f, a, b, n);

		Runge(n, I1, I2, method);
	}
}

// вспомогательные функции для вывода таблицы на экран

Method choice(int i)
{
	switch (i)
	{
		case 1: return LeftRectangle;
		case 2: return RightRectangle;
		case 3: return Rectangle;
		case 4: return Trapezoid;
		case 5: return Simpson;
	}
}

string method(int i)
{
	switch (i)
	{
		case 1: return "LeftRectangle";
		case 2: return "RightRectangle";
		case 3: return "Rectangle";
		case 4: return "Trapezoid";
		case 5: return "Simpson";
	}
}

int main()
{
	cout << setprecision(12) << fixed;

	assert(-R <= a && b <= R);

	// шапка таблицы

	cout << left << setw(25) << "Method" << setw(12) << "n" << setw(20) << "I1" << setw(20) << "|I1 - RealAns|" << setw(15) << "|I1 - I2|" << "\n";

	// выбор методов вычисления интеграла:
	// 1 - метод левых прямоугольников
	// 2 - метод правых прямоугольников
	// 3 - метод средних прямоугольников
	// 4 - метод трапеций
	// 5 - метод Симпсона (парабол)

	int variants[] = { 2, 5 };
	int size = sizeof variants / sizeof(*variants);

	// реальный ответ

	double RealAns = F(b) - F(a);

	// вывод таблицы на экран

	for (int i = 0; i < size; ++i)
	{
		int n = 4, cur = variants[i];

		double I1 = choice(cur)(f, a, b, n);
		double I2 = choice(cur)(f, a, b, n / 2);

		Runge(n, I1, I2, choice(cur));

		cout << left << setw(25) << method(cur) << setw(12) << n << setw(20) << I1 << setw(20) << fabs(I1 - RealAns) << setw(15) << fabs(I1 - I2) << "\n";
	}

	return 0;
}

//integrate.cpp
#include <cassert>

typedef double (*Function)(double);

// Метод левых прямоугольников
double LeftRectangle(Function f, double a, double b, int n)
{
	assert(n >= 2);

	double res = 0.0;

	double h = (b - a) / (n - 1), x = a;

	for (int i = 1; i <= n - 1; ++i, x += h)
		res += f(x);

	return res * h;
}

// Метод правых прямоугольников
double RightRectangle(Function f, double a, double b, int n)
{
	assert(n >= 2);

	double res = 0.0;

	double h = (b - a) / (n - 1), x = a + h;

	for (int i = 2; i <= n; ++i, x += h)
		res += f(x);

	return res * h;
}

// Метод средних прямоугольников
double Rectangle(Function f, double a, double b, int n)
{
	assert(n >= 2);

	double res = 0.0;

	double h = (b - a) / (n - 1), x = a + h / 2;

	for (int i = 1; i <= n - 1; ++i, x += h)
		res += f(x);

	return res * h;
}

// Метод трапеций
double Trapezoid(Function f, double a, double b, int n)
{
	assert(n >= 2);

	double res = (f(a) + f (b)) / 2.0;

	double h = (b - a) / (n - 1), x = a + h;

	for (int i = 2; i <= n - 1; ++i, x += h)
		res += f(x);

	return res * h;
}

// Метод Симпсона (парабол)
double Simpson(Function f, double a, double b, int n)
{
	assert(n >= 2);

	double res1 = 0.0, res2 = 0.0;

	double h = (b - a) / 2.0 / (n - 1), x = a + h;

	for (int i = 2; i <= 2 * n - 2; ++i, x += h)
		if (i & 1) res1 += f(x);
		else res2 += f(x);

	double res = f(a) + f(b) + 2 * res1 + 4 * res2;

	return res * h / 3.0;
}

//integrate.h
#ifndef INTEGRATE
#define INTEGRATE

typedef double (*Function)(double);
typedef double (*Method)(Function, double, double, int);

double LeftRectangle(Function, double, double, int);
double RightRectangle(Function, double, double, int);
double Rectangle(Function, double, double, int);
double Trapezoid(Function, double, double, int);
double Simpson(Function, double, double, int);

#endif