lab6

#include <iostream>
#include <stdio.h>
#include <math.h>

#define N 2000

double a = 1;
double b = 6.5;

double F(double x) {
    return exp(pow(x + 5.0, 0.25)) * sin(3.0 * x) * cos(x);
}

int main() {
    double A[N], B[N], C[N], D[N];

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

    double *matrix[N];
    for (int i = 0; i < N; i++) {
        matrix[i] = new double[N + 1]{ 0 };
    }

    matrix[0][0] = 1;
    matrix[N - 1][N - 1] = 1;

    for (int i = 1; i < N - 1; i++) {
        for (int j = 0; j <= N; j++) {
            matrix[i][j] = 0;
            if (j == (i - 1)) matrix[i][j] = h;
            if (j == i) matrix[i][j] = 4.0 * h;
            if (j == (i + 1)) matrix[i][j] = h;
            if (j == N) matrix[i][j] = 6.0 * (F(a + h * (i - 1)) - 2 * F(a + h * i) + F(a + h * (i + 1))) / h;
        }
    }

    for (int k = 0; k < N - 1; k++) {
        double m = 1.0 / matrix[k][k];
        for (int j = 0; j <= N; j++) {
            matrix[k][j] = matrix[k][j] * m;
        }

        for (int i = k + 1; i < N; i++) {
            double m = matrix[i][k];
            for (int j = 0; j <= N; j++) {
                matrix[i][j] = matrix[i][j] - matrix[k][j] * m;
            }
        }
    }

    for (int i = 0; i < N; i++) {
        C[N - i - 1] = matrix[N - i - 1][N];
        for (int j = 0; j < i; j++) {
            C[N - i - 1] -= matrix[N - i - 1][N - 1 - j] * C[N - i];
        }
    }

    for (int i = 0; i < N; i++) {
        delete[] matrix[i];
    }

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

    B[0] = 0;

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

    int k = (b - a) * 10;
    double step = (b - a) / (double)k;
    double x = a;
    for (int i = 0; i <= k; i++) {
        int j = 0;

        while (x >= (a + j * h)) j++;
        --j;

        double xi = a + j * h;
        double y = A[j] + B[j] * (x - xi) + C[j] * (x - xi) * (x - xi) / 2.0 + D[j] * (x - xi) * (x - xi) * (x - xi) / 6.0;

        printf("%4.2f %7f\n", x, y);

        x += step;
    }
    return 0;
}