RK-4_ode

#include <cmath>
#include <iostream>
#include <fstream>
#include <array>

const int N = 3; // размерность системы
const int number_of_generated_points = 100000;
const int area_size = 10; // параметр, отвечающий за размер области генерации точек

class RungeKutta {
public:

    double b = 0.11;
    std::array<double, N> X, XH, K1, K2, K3, K4, FX;


    RungeKutta() {
        X.fill(0);
        XH.fill(0);
        K1.fill(0);
        K2.fill(0);
        K3.fill(0);
        K4.fill(0);
        FX.fill(0);
    };

    double getX(int index) {
        return X.at(index);
    }

    void SetInit(std::array<double, N> X0) {
        for (size_t i = 0; i < X.size(); i++)
            X[i] = X0[i];
    }

    std::array<double, N> F(std::array<double, N>& X) {
        FX[0] = sin(X[1]) - b * X[0];
        FX[1] = sin(X[2]) - b * X[1];
        FX[2] = sin(X[0]) - b * X[2];
        return FX;
    }

    void NextStep(double h) {
        if (h < 0) return;

        K1 = F(X);

        for (size_t i = 0; i < X.size(); i++)
            XH[i] = X[i] + K1[i] * (h / 2.0);

        K2 = F(XH);

        for (size_t i = 0; i < X.size(); i++)
            XH[i] = X[i] + K2[i] * (h / 2.0);

        K3 = F(XH);

        for (size_t i = 0; i < X.size(); i++)
            XH[i] = X[i] + K3[i] * h;

        K4 = F(XH);

        for (size_t i = 0; i < X.size(); i++)
            X[i] = X[i] + h / 6.0 * (K1[i] + 2.0 * K2[i] + 2.0 * K3[i] + K4[i]);
    }
};

class DOPRI8 {
public:
    double b = 0.11;
    std::array<double, N> X, XH, K1, K2, K3, K4, K5, K6,
        K7, K8, K9, K10, K11, K12, K13, FX;


    DOPRI8() {
        X.fill(0);
        XH.fill(0);
        K1.fill(0);
        K2.fill(0);
        K3.fill(0);
        K4.fill(0);
        K5.fill(0);
        K6.fill(0);
        K7.fill(0);
        K8.fill(0);
        K9.fill(0);
        K10.fill(0);
        K11.fill(0);
        K12.fill(0);
        K13.fill(0);
        FX.fill(0);
    }

    double getX(uint32_t index) {
        return X.at(index);
    }

    void SetInit(std::array<double, N> X0) {
        for (size_t i = 0; i < X.size(); i++)
            X[i] = X0[i];
    }

    std::array<double, N> F(std::array<double, N>& X) {
        FX[0] = sin(X[1]) - b * X[0];
        FX[1] = sin(X[2]) - b * X[1];
        FX[2] = sin(X[0]) - b * X[2];
        return FX;
    }

    void NextStep(double h) {
        if (h < 0) return;

        K1 = F(X);

        for (size_t i = 0; i < X.size(); i++)
            XH[i] = X[i] + K1[i] * (h / 18.0);

        K2 = F(XH);

        for (size_t i = 0; i < X.size(); i++)
            XH[i] = X[i] + K2[i] * (h / 48.0 + h / 16.0);

        K3 = F(XH);

        for (size_t i = 0; i < X.size(); i++)
            XH[i] = X[i] + K3[i] * (h / 32.0 + 3 * h / 32.0);

        K4 = F(XH);

        for (size_t i = 0; i < X.size(); i++)
            XH[i] = X[i] + K4[i] * (5 * h / 16.0 - 75.0 * h / 64.0 + 75.0 * h / 64.0);

        K5 = F(XH);

        for (size_t i = 0; i < X.size(); i++)
            XH[i] = X[i] + K5[i] * (3 * h / 80.0 + 3.0 * h / 16.0 + 3.0 * h / 20.0);

        K6 = F(XH);

        for (size_t i = 0; i < X.size(); i++)
            XH[i] = X[i] + K6[i] * (29443841.0 * h / 614563906.0 +
                77736538.0 * h / 77736538.0 - 28693883.0 * h / 1125000000.0 + 23124283.0 * h / 1800000000.0);

        K7 = F(XH);

        std::array<double, 7> T1 = { 16016141.0 / 946692911.0, 0.0, 0.0, 61564180.0 / 158732637.0, 22789713.0 / 633445777.0, 545815736.0 / 2771057229.0, -180193667.0 / 1043307555.0 };

        double t1 = 0.0;
        for (const auto& num : T1)
            t1 += h * num;

        for (size_t i = 0; i < X.size(); i++)
            XH[i] = X[i] + K7[i] * t1;

        K8 = F(XH);

        std::array<double, 8> T2 = { 39632708.0 / 573591083.0, 0.0, 0.0, -433636366.0 / 683701615.0, -421739975.0 / 2616292301.0, 100302831.0 / 723423059.0, 790204164.0 / 839813087.0, 800635310.0 / 3783071287.0 };

        double t2 = 0.0;
        for (const auto& num : T2)
            t2 += h * num;

        for (size_t i = 0; i < X.size(); i++)
            XH[i] = X[i] + K8[i] * t2;

        K9 = F(XH);

        std::array<double, 9> T3 = { 246121993.0 / 1340847787.0, 0.0, 0.0, -37695042795.0 / 15268766246.0, -309121744.0 / 1061227803.0, -12992083.0 / 490766935.0, 6005943493.0 / 2108947869.0, 393006217.0 / 1396673457.0, 123872331.0 / 1001029789.0 };

        double t3 = 0.0;
        for (const auto& num : T3)
            t3 += h * num;

        for (size_t i = 0; i < X.size(); i++)
            XH[i] = X[i] + K9[i] * t3;

        K10 = F(XH);

        std::array<double, 10> T4 = { -1028468189.0 / 846180014.0, 0.0, 0.0, 8478235783.0 / 508512852.0, 1311729495.0 / 1432422823.0, -10304129995.0 / 1701304382.0, -48777925059.0 / 3047939560.0, 15336726248.0 / 1032824649.0, -45442868181.0 / 3398467696.0, 3065993473.0 / 597172653.0 };

        double t4 = 0.0;
        for (const auto& num : T4)
            t4 += h * num;

        for (size_t i = 0; i < X.size(); i++)
            XH[i] = X[i] + K10[i] * t4;

        K11 = F(XH);

        std::array<double, 11> T5 = { 185892177.0 / 718116043.0, 0.0, 0.0, -3185094517.0 / 667107341.0, -477755414.0 / 1098053517.0, -703635378.0 / 230739211.0, 5731566787.0 / 1027545527.0, 5232866602.0 / 850066563.0, -4093664535.0 / 808688257.0, 3962137247.0 / 1805957418.0, 65686358.0 / 487910083.0 };

        double t5 = 0.0;
        for (const auto& num : T5)
            t5 += h * num;

        for (size_t i = 0; i < X.size(); i++)
            XH[i] = X[i] + K11[i] * t5;

        K12 = F(XH);

        std::array<double, 12> T6 = { 403863854.0 / 491063109.0, 0.0, 0.0, -5068492393.0 / 434740067.0, -411421997.0 / 543043805.0, 652783627.0 / 914296604.0, 11173962825.0 / 925320556.0, -13158990841.0 / 6184727034.0, 3936647629.0 / 1978049680.0, -160528059.0 / 685178525.0, 248638103.0 / 1413531060.0, 0.0 };

        double t6 = 0.0;
        for (const auto& num : T6)
            t6 += h * num;

        for (size_t i = 0; i < X.size(); i++)
            XH[i] = X[i] + K12[i] * t6;

        K13 = F(XH);

        std::array<double, 13> Tb = { 14005451.0 / 335480064.0, 0.0, 0.0, 0.0, 0.0, -59238493.0 / 1068277825.0, 181606767.0 / 758867731.0, 561292985.0 / 797845732.0, -1041891430.0 / 1371343529.0, 760417239.0 / 1151165299.0, 118820643.0 / 751138087.0, -528747749.0 / 2220607170.0, 1.0 / 4.0 };

        double sum = 0.0;
        for (const auto& b : Tb) {
            sum += b;
        }

        for (size_t i = 0; i < X.size(); i++)
            X[i] = X[i] + h * (K1[i] * Tb[0] + K2[i] * Tb[1] + K3[i] * Tb[2] + K4[i] * Tb[3] + K5[i] * Tb[4] + K6[i] * Tb[5] + K7[i] * Tb[6] + K8[i] * Tb[7] + K9[i] * Tb[8] + K10[i] * Tb[9] + K11[i] * Tb[10] + K12[i] * Tb[11] + K13[i] * Tb[12]);

    }

};

// ф-ия main в данном файле используется для генерации значений точек

//int main() {
//    std::ofstream out;
//    out.open("prep.txt");
//    if (out.is_open()) {
//        for (int i = 0; i < number_of_generated_points; ++i) {
//            double x = (rand() % (int)(2 * area_size)) - area_size;
//            double y = (rand() % (int)(2 * area_size)) - area_size;
//            double z = (rand() % (int)(2 * area_size)) - area_size;
//            out << x << " " << y << " " << z << '\n';
//        }
//    }
//    out.close();
//    return 0;
//}