Вычисление интеграла функции на заданном отрезке

Ниже приведены методы вычисления интеграла заданной функции на заданном отрезке.
 Стоит отметить, что в "чистом" виде эти методы дают достаточно большие погрешности,
 либо выполняются слишком долго, поэтому в реальных задачах все методы используются
 вместе с Правилом Рунге, идея которого состоит в том, что сначала интеграл
 вычисляется для шага h, а затем для шага h/2. Есть разность между полученными
 результатами лежит за рамках допустимой погрешности, шаг вычислений уменьшается
 в 2 раза и расчеты повторяются.

/*///////////////////////////////////////////////////////////////
 
 f(x) = 1/4*x^2+2, x=[2,10], GNU C++, 2011-01-21
 
 Artix, master@7masterov.ru, icq:53666599, skype:artixmaster
 
 * Error in code? Nothing is perfect!
 * Free source for free Linux, use it for free!
 * Please, do not remove this comment!
 
///////////////////////////////////////////////////////////////*/
 
#include <stdio.h>
#include <math.h>
#include <stdlib.h>
#include <time.h>
 
#define STEP 500
 
typedef long double ldouble;
 
ldouble f(ldouble x)
{
    return 1./4.*x*x+2.;
}
 
// Метод Прямоугольников
ldouble integralRectangle(ldouble a=2, ldouble b=10, long steps=STEP)
{
    ldouble h = b-a;
    ldouble delta = h/steps;
    ldouble s = 0;
    long count = 0;
    ldouble x=a;
    for(long i=0;i<steps;i++)
    {
        s+=f(x+delta/2.);
        x+=delta;
        count+=1;
    }
    printf("integralRectangle() num iterations: %ld\n",count);
    return delta * s;
}
 
// Метод Трапеций
ldouble integralTrapeze(ldouble a=2, ldouble b=10, long steps=STEP)
{
    ldouble h = b-a;
    ldouble delta = h/steps;
    ldouble s = 0;
    long count=0;
    ldouble x=a;
    for(long i=0;i<steps;i++)
    {
        s+=f(x)+f(x+delta);
        x+=delta;
        count+=1;
    }
    printf("integralTrapeze() num iterations: %ld\n",count);
    return delta * s / 2.;
}
 
// Метод Парабол (метод Симпсона)
ldouble integralSimpson(ldouble a=2, ldouble b=10, long steps=STEP)
{
    ldouble h = b-a;
    ldouble delta = h/steps;
    ldouble s = 0;
    long count = 0;
    ldouble x=a;
    for(long i=0;i<steps;i++)
    {
        s+=delta/6 * ( f(x) + 4*f((x+x+delta)/2) + f(x+delta) );
        x+=delta;
        count+=1;
    }
    printf("integralSimpson() num iterations: %ld\n",count);
    return s;
}
 
// Метод Симпсона с оценкой ошибки (Правило Рунге)
ldouble integralSimpsonRunge(ldouble a=2, ldouble b=10, ldouble e=0.02)
{
 
    ldouble h = b - a;
    ldouble s = f(a) + f(b);
    ldouble x, r_2n, s_n;
    ldouble s_2n = s;
    long count = 0;
    do {
        s_n = s_2n;
        h /= 2;
        s_2n = 0;
        x = a + h;
 
        do {
            s_2n += f(x)*2;
            x += h*2;
            count+=1;
        }
        while (x<b);
 
        s += s_2n;
        s_2n = h/3 * (s + s_2n);
        r_2n = fabs(s_n - s_2n) / 15;
    }
    while (r_2n>e);
    printf("integralSimpsonRunge() num iterations: %ld\n",count);
    return s_2n;
}

// Метод Прямоугольников с оценкой ошибки (Правило Рунге)
ldouble integralRectangleRunge(ldouble a=2, ldouble b=10, ldouble e=0.02)
{
    ldouble h = b-a;
    ldouble s1 = h*f((a+b)/2.);
    long n = 1;
    long count = 0;
    ldouble s2;
    ldouble fe;
    do {
        n*=2; h/=2.;
        s2=s1; s1=0;
        for(long i=0;i<n;i++) {
            s1+=f(a+h*i+h/2);
            count+=1;
        }
        s1*=h;
        fe = fabs(s2-s1)*10.;
    } while (fe>=e);
    printf("integralRectangleRunge() num iterations: %ld\n",count);
    return s1;
}

// Метод Трапеций с оценкой ошибки (Правило Рунге)
ldouble integralTrapezeRunge(ldouble a=2, ldouble b=10, ldouble e=0.02)
{
    ldouble h = b-a;
    ldouble s1 = h*f((a+b)/2);
    long n = 1;
    long count = 0;
    ldouble s2;
    ldouble fe;
    do {
        n*=2; h/=2;
        s2=s1; s1=0;
        for(long i=0;i<n;i++) {
            s1+=f(a+h*i)+f(a+h*i+h);
            count += 1;
        }
        s1*=h/2;
        fe = fabs(s2-s1)*10;
    } while (fe>=e);
    printf("integralTrapezeRunge() num iterations: %ld\n",count);
    return s1;
}

// Метод Монте-Карло (стохастический)
ldouble integralMonteKarlo(ldouble a=2, ldouble b=10, long steps=STEP)
{
    // Detect min/max of f() to find rectangle square
    ldouble dx = b-a;
    ldouble delta = dx/steps;
    ldouble min=0,max=f(a+dx/2);
    for(ldouble x=a;x<=b;x+=delta)
    {
        if (f(x)>max) max=f(x);
    }
    ldouble dy = max-min;
    ldouble s = dy*dx;
    ldouble sdelta = dy/steps;
 
    // Random experiment
    srand(time(NULL));
    long success = 0;
    for(long i=0;i<steps*steps;i++)
    {
        ldouble x = a+delta*(rand() % steps) + 1./(rand()%steps);
        ldouble y = min+sdelta*(rand() % steps) + 1./(rand()%steps);
        if (y<=f(x)) success+=1;
    }
    return s*success*1./steps/steps;
}
 
int main()
{
    ldouble f;
    ldouble realResult = f = 248.0/3.0+16;
    ldouble f1 = integralRectangle();
    ldouble f2 = integralTrapeze();
    ldouble f3 = integralSimpson();
    ldouble f4 = integralSimpsonRunge();
    ldouble f5 = integralMonteKarlo();
    ldouble f6 = integralRectangleRunge();
    ldouble f7 = integralTrapezeRunge();
    printf("Real result, range x=[2,10], f=%Lf\n", realResult);
    printf("Integral rectangle result: %.8Lf (E=%.8Lf)\n", f1,f-f1);
    printf("Integral trapeze result: %.8Lf (E=%.8Lf)\n", f2,f-f2);
    printf("Integral simpson result: %.8Lf (E=%.8Lf)\n", f3,f-f3);
    printf("Integral simpson-runge result: %.8Lf (E=%.8Lf)\n", f4,f-f4);
    printf("Integral montekarlo result: %.8Lf (E=%.8Lf)\n", f5,f-f5);
    printf("Integral rectangle-runge result: %.8Lf (E=%.8Lf)\n", f6,f-f6);
    printf("Integral trapeze-runge result: %.8Lf (E=%.8Lf)\n", f7,f-f7);
    printf("Cycle step per unit: %d\n",STEP);
    return 1;
}