spline_1

//19142, Аношин Сергей, Интерполирование сплайном
#include <iostream>
#include <fstream>
#include <cmath>
#include <iomanip>
#define A 0
#define B 1
#include "Прогонка.h"

using namespace std;

double func(double x) {
    return exp(x);
}

void bound_cond(char c, double* arr_spl, int N, double step) {
    if (c == 'a') {
        arr_spl[0] = func(A);
        arr_spl[N - 1] = func(B);
    }

    if (c == 'b') {
        arr_spl[0] = (1 - 2 * func(step) + func(2 * step)) / (pow(step, 2));
        arr_spl[N - 1] = (func(1) - 2 * func(1 - step) + func(1 - 2 * step)) / (pow(step, 2));
    }

    if (c == 'c') {
        arr_spl[0] = 0;
        arr_spl[N - 1] = 0;
    }
}

double* ind_h_spl(int N) {
    double step = (double(B - A) / (double(N) - 1));

    double A_sweep;
    double B_sweep;
    double C_sweep;

    A_sweep = B_sweep = (step / 6);
    C_sweep = double(2 * step) / 3;
    double* h_spline = new double[N];


    h_spline[1] = (func(A + 2 * step) - 2 * func(A + step) + func(A)) / (pow(step, 2));
    h_spline[N - 2] = (func(B) - 2 * func(B - step) + func(B - 2 * step)) / pow(step, 2);

    double* F_sweep_1 = new double[N - 4];

    F_sweep_1[0] = (func(A + 3 * step) - 2 * func(A + 2 * step) + func(A + step)) / step - ((step * h_spline[1]) / 6);
    F_sweep_1[N - 5] = (func(B - step) - 2 * func(B - 2 * step) + func(B - 3 * step)) / step - (step * h_spline[N - 2]) / 6;
    for (int i = 3; i < N - 3; i++)
        F_sweep_1[i - 2] = (func(A + (double(i) + 1) * step) - 2 * func(A + i * step) + func(A + (double(i) - 1) * step)) / step;

    double* ans = new double[N - 4];
    ans = sweep_method(N - 4, A_sweep, B_sweep, C_sweep, F_sweep_1);

    for (int i = 2; i < N - 2; i++)
        h_spline[i] = ans[i - 2];

    h_spline[0] = 2 * h_spline[1] - h_spline[2];
    h_spline[N - 1] = 2 * h_spline[N - 2] - h_spline[N - 3];

    delete[] ans;
    delete[] F_sweep_1;

    return h_spline;
}

double* help_spline(int N, char c) {
    double step = (double(B - A) / (double(N) - 1));

    double A_sweep;
    double B_sweep;
    double C_sweep;

    A_sweep = B_sweep = (step / 6);
    C_sweep = double(2 * step) / 3;
    double* h_spline = new double[N];

    double* F_sweep = new double[N - 2];
    bound_cond(c, h_spline, N, step);

    F_sweep[0] = (func(A + 2 * step) - 2 * func(A + step) + func(A)) / step - (step * h_spline[0]) / 6;
    F_sweep[N - 3] = (func(B) - 2 * func(B - step) + func(B - 2 * step)) / step - (step * h_spline[N - 1]) / 6;

    for (int i = 2; i < N - 2; i++)
        F_sweep[i - 1] = (func(A + (double(i) + 1) * step) - 2 * func(A + i * step) + func(A + (double(i) - 1) * step)) / step;

    double* ans = new double[N - 2];
    ans = sweep_method(N - 2, A_sweep, B_sweep, C_sweep, F_sweep);

    for (int i = 1; i < N - 1; i++)
        h_spline[i] = ans[i - 1];

    delete[] ans;
    delete[] F_sweep;


    return h_spline;
}

double spline(double* h_spl, int N, double x) {
    double step = (double(B - A) / double(double(N) - 1));
    int count;
    for (count = 1; count < N; count++)
        if (x < step * count) break;
    count = count - 1;
    return (((exp(count * step) * ((double(count) + 1) * step - x)) / step) + ((exp((double(count) + 1) * step) * (x - count * step)) / step) + h_spl[count] * (pow(((double(count) + 1) * step - x), 3) - pow(step, 2) * ((double(count) + 1) * step - x)) / (6 * step) + h_spl[count + 1] * (pow((x - count * step), 3) - pow(step, 2) * (x - count * step)) / (6 * step));
}

void error(char cas, ofstream& fin, double eps, int flag) {
    int count = 0;
    int N = 10;
    double step = (double(B - A) / ((double(N) - 1) * 10));
    double err = 1;
    double h_err;
    double p;


    while (err > eps) {
        err = 0;
        double* spl = new double[N];

        if (cas != 'i')
            spl = help_spline(N, cas);
        else {
            spl = ind_h_spl(N);
        }

        for (double i = A; i <= B; i = i + step) {
            if (abs(spline(spl, N, i) - func(i)) >= err)
                err = abs(spline(spl, N, i) - func(i));
        }

        if (count != 0) {
            p = log(h_err / err) / log(2);
        }

        h_err = err;

        if (flag == 1) {
            if (count != 0)
                fin << "number of nodes: " << N << "    " << "error: " << setprecision(30) << err << "  error's order: " << p << endl;
            else fin << "number of nodes: " << N << "    " << "error: " << setprecision(30) << err << " error's order: isn't defined " << endl;
            N = 2 * N;
        }

        if (flag == 0) {
            fin << N << "   " << setprecision(20) << log10(err) << endl;
            N = 1 + N;
        }

        step = (double(B - A) / ((double(N) - 1) * 10));

        count++;
        delete[] spl;
    }

}

int main() {
    double eps = 0.0000000001;

    ofstream error_1("C:\\Users\\Сергей\\Desktop\\C++ ВМЛА\\Сплайн\\error_spline_'a'.txt");
    error('a', error_1, eps, 1);

    ofstream error_2("C:\\Users\\Сергей\\Desktop\\C++ ВМЛА\\Сплайн\\error_spline_'b'.txt");
    error('b', error_2, eps, 1);

    ofstream error_3("C:\\Users\\Сергей\\Desktop\\C++ ВМЛА\\Сплайн\\error_spline_'c'.txt");
    error('c', error_3, eps, 1);

    ofstream error_4("C:\\Users\\Сергей\\Desktop\\C++ ВМЛА\\Сплайн\\error_spline_'i'.txt");
    error('i', error_4, eps, 1);

    ofstream error_5("C:\\Users\\Сергей\\Desktop\\C++ ВМЛА\\Сплайн\\diagram_error_spline_'i'.txt");
    error('i', error_5, eps, 0);

    ofstream error_6("C:\\Users\\Сергей\\Desktop\\C++ ВМЛА\\Сплайн\\diagram_error_spline_'a'.txt");
    error('a', error_6, eps, 0);

    return 0;
}