#include <stdio.h>
#include <math.h>
#include <iostream>
#include <string>
#include <time.h>
#define EPS 0.01
using std::cout;
using std::cin;
using std::endl;
class NumericalIntegration {
private:
double (*func)(double); //функция
double left, right; //границы
double real_value; //эталон
double min, max; //минимальное и максимальное
//значения функции
double (NumericalIntegration::*tek_method)(int);
double fRand(double fMin, double fMax);
//функции численного интегрирования
double MiddRect(int count);
double Trapeze(int count);
double Simpson(int count);
double MonteKarlo(int count);
public:
enum methods{MIDDRECT, TRAPEZE, SIMPSON, MONTEKARLO};
NumericalIntegration(double (*f)(double),
double l, double r, double val, double min, double max);
~NumericalIntegration() {};
//аксессоры
void SetFunction(double (*f)(double)) { func = f; };
void SetBoards(double l, double r) { left = l; right = r; };
void SetRealValue(double v) { real_value = v; };
void SetMinValue(double v) { min = v; }
void SetMaxValue(double v) { max = v; }
void Calc(methods method, double eps);
};
NumericalIntegration::NumericalIntegration(double (*f)(double),
double l, double r, double val, double min, double max)
{
func = f;
left = l;
right = r;
real_value = val;
this->min = min;
this->max = max;
}
double NumericalIntegration::fRand(double fMin, double fMax)
{
double f = (double)rand() / RAND_MAX;
return fMin + f * (fMax - fMin);
}
double NumericalIntegration::MiddRect(int count)
{
double step = (right - left) / count;
double res = 0;
for (int i = 0; i < count; i++) {
res += func(left + step * (2 * i + 1) / 2);
}
return res * step;
}
double NumericalIntegration::Trapeze(int count)
{
double step = (right - left) / count;
double res = 0;
for (int i = 1; i <= count; i++) {
res += func(left + step * (i - 1)) + func(left + step * i);
}
return res * step / 2;
}
double NumericalIntegration::Simpson(int count)
{
double step = (right - left) / count;
double res = 0;
for (int i = 1; i < count; i += 2) {
res += func(left + step * (i - 1)) + 4 * func(left + step * i) + func(left + step * (i + 1));
}
return res * step / 3;
}
double NumericalIntegration::MonteKarlo(int count)
{
int k = 0;
for (int i = 0; i < count; i++) {
double x = fRand(left, right);
double y = fRand(min, max);
if (func(x) > y) k++;
}
return (right - left) * (max - min) * k / count;
}
void NumericalIntegration::Calc(methods method, double eps)
{
int err, k, count = 1;
std::string method_name;
switch (method) {
case MIDDRECT:
tek_method = &NumericalIntegration::MiddRect;
err = 3;
k = 2;
method_name = "MIDDRECT";
break;
case TRAPEZE:
tek_method = &NumericalIntegration::Trapeze;
err = 3;
k = 2;
method_name = "TRAPEZE ";
break;
case SIMPSON:
tek_method = &NumericalIntegration::Simpson;
err = 15;
k = 4;
method_name = "SIMPSON ";
break;
case MONTEKARLO:
tek_method = &NumericalIntegration::MonteKarlo;
err = 1;
count = 16;
method_name = "MONTEKARLO";
break;
default:
cout << "Unknown method" << endl;
return;
}
double I1;
double I2 = (this->*tek_method)(count);
do {
count *= 2;
I1 = I2;
I2 = (this->*tek_method)(count);
} while (abs(I1 - I2) / err > eps);
double richardson = (method == MONTEKARLO) ? I2 :
((1 << k) * I2 - I1) / ((1 << k) - 1);
cout << method_name << '\t'
<< count << '\t';
if (method != MONTEKARLO) cout << I2 << '\t'
<< richardson << '\t';
cout << abs(real_value - richardson) << endl;
}
double Function(double x) { return exp(x); }
int main()
{
srand(time(NULL));
NumericalIntegration integrate(Function, 0, 1,
exp(1.0) - 1, 0, Function(1));
integrate.Calc(NumericalIntegration::MIDDRECT, EPS);
integrate.Calc(NumericalIntegration::TRAPEZE, EPS);
integrate.Calc(NumericalIntegration::SIMPSON, EPS);
integrate.Calc(NumericalIntegration::MONTEKARLO, EPS);
getchar();
return 0;
}