lab6.cpp

using namespace std;

#include <cmath>
#include <vector>
#include <cstdio>

double f(double x) {
    return 1 / x - 0.1 * x * x * sin(2 * x);
}

double *coef_a(double a, double b, int N) {
    double *A = new double[N];
    double h = (b - a) / N;
    double x = a;
    for (int i = 0; i < N; i++) {
        A[i] = f(x);
        x += h;
    }
    return A;
}

double *coef_d(double a, double b, double *C, int N) {
    double *D = new double[N];
    double h = (b - a) / N;

    D[0] = C[0];
    for (int i = 1; i < N; i++) {
        D[i] = (C[i] - C[i - 1]) / h;
    }
    return D;
}

double *coef_b(double a, double b, double *C, double *D, int N) {
    double *B = new double[N];
    double h = (b - a) / N;
    B[0] = 0;

    for (int i = 1; i < N; i++) {
        B[i] = (h * C[i] / 2) - (h * h * D[i] / 2) + (f(a + h * i) - f(a + h * (i - 1))) / h;
    }
    return B;
}

double *slar(double a, double b, int N) {
    double *matrix = new double[N * (N + 1)];

    double h = (b - a) / N;
    for (int i = 0; i < N; i++) {
        for (int j = 0; j < N + 1; j++) {
            matrix[i * (N + 1) + j] = 0;
        }
    }

    for (int i = 0; i < N; i++) {
        matrix[i * (N + 1) + i] = 4 * h;
        if (i > 0) {
            matrix[i * (N + 1) + i - 1] = h;
        }
        matrix[i * (N + 1) + i + 1] = h;
        double x_1 = a + h * (i - 1);
        double x_i = a + h * i;
        double x_n = a + h * (i + 1);
        matrix[i * (N + 1) + N] = 6 * (f(x_n) - 2 * f(x_i) + f(x_1)) / h;
    }
    return matrix;
}


double *single_division(double *A, int N) {
    double *x = new double[N];

    for (int j = 0; j < N; j++) {
        double m = 1.0 / A[j * (N + 1) + j];
        for (int i = j; i < N + 1; i++) {
            A[j * (N + 1) + i] *= m;
        }
        for (int k = j + 1; k < N; k++) {
            double a = A[k * (N + 1) + j];
            for (int l = j; l < N + 1; l++) {
                A[k * (N + 1) + l] -= A[j * (N + 1) + l] * a;
            }
        }
    }
    for (int j = N - 1; j >= 0; j--) {
        double s = 0;
        for (int i = N - 1; i > j; i--) {
            s += A[j * (N + 1) + i] * x[i];
        }
        x[j] = A[j * (N + 1) + N] - s;
    }
    return x;
}

double splain(double x, double a, double b, double *A, double *B, double *C, double *D, int N) {
    int i = 0;
    double h = (b - a) / N, xi = a + i * h;

    while (x >= xi) {
        xi += h;
        i++;
    }
    i--;
    xi = a + i * h;

    return A[i] + B[i] * (x - xi) + C[i] * (x - xi) * (x - xi) / 2 + D[i] * (x - xi) * (x - xi) * (x - xi) / 6;
}


int main() {
    double a = 2;
    double b = 11;
    int N = 55;

    double *M = slar(a, b, N);
    double *C = single_division(M, N);
    double *A = coef_a(a, b, N);
    double *D = coef_d(a, b, C, N);
    double *B = coef_b(a, b, C, D, N);

    double x = a;
    int n = (int((b - a) * N) >> 1) << 1;
    double step = (b - a) / (double)n;
    printf("|      x      |       y       |\n");
    for (int i = 0; i <= n; i++, x += step) {
        printf("| %.9f | %13.9f |\n", x, splain(x, a, b, A, B, C, D, N));
    }

    delete[] M;
    delete[] A;
    delete[] B;
    delete[] C;
    delete[] D;
}