ПВ численное интегрирование

#include <iostream>
#include <cmath>


/// Значение функции в точке x
double function(double x) {
    return x*x*x*x + x*x + 2.0;
}

/// Значение второй производной в точке x
double seconddif(double x, double h) {
    //u''(x) = 1 / hi ((ui+1-ui) / hi + 1 - (ui-ui-1) / hi)
    double f1 = function(x);
    double f2 = function(x + h);
    double f3 = function(x - h);
    return (f2 - 2.0 * f1 + f3) / (h * h);
}

double CalcRectangle(double a, double b, int n) {
    double result = 0;
    double h = (b - a) / n; // Шаг сетки

    for(int i = 0; i < n; i++) {
        // a + h * (i + 0.5) - X координата середины узла
        result += function(a + h * (i + 0.5));
        // std::cout << result << std::endl;
    }

    result *= h;

    result += seconddif((a + b) / 2.0, h) * h * h * (b - a) / 24.0;

    return result;
}


double CalcTrapeze(double a, double b, int n) {
    double result = 0;
    double h = (b - a) / n; // Шаг сетки


    for (int i = 1; i < n; i++) {
        // a + h * i - X координата узла
        result += function(a + h * i);
    }
    result += (function(a) + function(b)) / 2.0;
    result *= h;
    result -= h * h * (b - a) * seconddif((a + b) / 2.0, h) / 12.0 ;

    return result;
}


int main() {
    double integral = CalcRectangle(0, 5, 1);
    std::cout << integral << std::endl;

    integral = CalcRectangle(0, 5, 2);
    std::cout << integral << std::endl;

    integral = CalcRectangle(0, 5, 4);
    std::cout << integral << std::endl;

    std::cout << std::endl;

    integral = CalcTrapeze(0, 5, 1);
    std::cout << integral << std::endl;

    integral = CalcTrapeze(0, 5, 2);
    std::cout << integral << std::endl;

    integral = CalcTrapeze(0, 5, 4);
    std::cout << integral << std::endl;

    return 0;
}