oksana lab4

using namespace std;
#include <vector>
#include <cmath>
#include <cstdio>

const double a = 1;
const double b = 27;

double integral_function(double x) {
    return pow(sin(x), 2);
}

double atiderivative(double x) {
    return -sin(2 * x) / 4 + x / 2;
}

double second_derivative(double x) {
    return 2 * (sin(2 * x) + x * cos(2 * x));
}

double get_step(double eps) {
    double max = 0;
    for (double i = a; i <= b; i += 0.01) {
        double cur = abs(second_derivative(i));
        if (cur > max) max = cur;
    }

    return sqrt((12 * eps) / ((b - a) * max));
}

double calculate_integral(double a, double b) {
    return atiderivative(b) - atiderivative(a);
}

double trapezoid(double a, double b, double h) {
    double psum = integral_function(a) + integral_function(b);
    int n = (b - a) / h;
    for (int i = 1; i < n; ++i) {
        psum = psum + 2.0 * integral_function(a + i * h);
    }
    return (h / 2) * psum;
}

void runge(double a, double b, double error) {
    double n = 2 * ceil(1 / sqrt(error) * 0.5);
    double val_2n = trapezoid(a, b, (b - a) / n), R;
    do {
        double val_n = val_2n;
        val_2n = trapezoid(a, b, (b - a) / (n *= 2));
        R = fabs(val_n - val_2n) / 3;
    } while (R > error);
    printf("| %.7e | %.5e | %.8e |\n", error, (b - a) / n, abs(val_2n - calculate_integral(a, b)));
}

int main() {
    vector<double> errors;

    printf("Trapezoid method\n|  Eps  |     Step    | integral value | Prodused error |\n");
    for (double eps = 0.1; eps > 1E-10; eps *= 0.01) {
        int n = 2 * ceil((b - a) / get_step(eps) * 0.5);
        double h = (b - a) / n;

        double value = trapezoid(a, b, h);
        double error = fabs(value - calculate_integral(a, b));
        printf("| %.0e | %.5e | %.11f | %.8e |\n", eps, h, value, error);
        errors.push_back(error);
    }

    printf("\nRunge principle\n|  Given error  |     Step    | Prodused error |\n");
    for (int i = 0; i < errors.size(); ++i) {
        runge(a, b, errors[i]);
    }
}