new_ode

#include <cmath>
#include <cstdio>
#include <iostream>
#include <vector>
#include <chrono>

class RungeKutta {
protected:
    double t;
    std::vector<double> X, XH, K1, K2, K3, K4, FX;

public:
    RungeKutta(uint32_t N) {
        Init(N);
    }

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

    void Init(uint32_t N) {
        X.resize(N);
        XH.resize(N);
        K1.resize(N);
        K2.resize(N);
        K3.resize(N);
        K4.resize(N);
        FX.resize(N);
    }

    void SetInit(double t0, std::vector<double> X0) {
        t = t0;
        if (X.empty()) Init(X0.size());
        for (size_t i = 0; i < X.size(); i++)
            X[i] = X0[i];
    }

    virtual std::vector<double> F(double t, std::vector<double> X) = 0;

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

        K1 = F(t, X);

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

        K2 = F(t + h / 2.0, XH);

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

        K3 = F(t + h / 2.0, XH);

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

        K4 = F(t + h, 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]);

        t = t + h;
    }
};

class DOPRI8 {
protected:
    double t;
    std::vector<double> X, XH, K1, K2, K3, K4, K5, K6,
        K7, K8, K9, K10, K11, K12, K13, FX;

public:
    DOPRI8(uint32_t N) {
        Init(N);
    }

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

    void Init(uint32_t N) {
        X.resize(N);
        XH.resize(N);
        K1.resize(N);
        K2.resize(N);
        K3.resize(N);
        K4.resize(N);
        K5.resize(N);
        K6.resize(N);
        K7.resize(N);
        K8.resize(N);
        K9.resize(N);
        K10.resize(N);
        K11.resize(N);
        K12.resize(N);
        K13.resize(N);
        FX.resize(N);
    }

    void SetInit(double t0, std::vector<double> X0) {
        t = t0;
        if (X.empty()) Init(X0.size());
        for (size_t i = 0; i < X.size(); i++)
            X[i] = X0[i];
    }

    virtual std::vector<double> F(double t, std::vector<double> X) = 0;

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

        K1 = F(t + h, X);

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

        K2 = F(t + h / 18.0, XH);

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

        K3 = F(t + h / 48.0 + h / 16.0, 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(t + h / 32.0 + 3 * h / 32.0, 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(t + 5 * h / 16.0 - 75.0 * h / 64.0 + 75.0 * h / 64.0, 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(t + 3 * h / 80.0 + 3.0 * h / 16.0 + 3.0 * h / 20.0, XH);

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

        K7 = F(t + 29443841.0 * h / 614563906.0 +
            77736538.0 * h / 77736538.0 - 28693883.0 * h / 1125000000.0 + 23124283.0 * h / 1800000000.0, XH);

        std::vector<double> 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] + K8[i] * t1;

        K8 = F(t + t1, XH);

        std::vector<double> 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] + K9[i] * t2;

        K9 = F(t + t2, XH);

        std::vector<double> 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] + K10[i] * t3;

        K10 = F(t + t3, XH);

        std::vector<double> 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] + K11[i] * t4;

        K11 = F(t + t4, XH);

        std::vector<double> 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] + K12[i] * t5;

        K12 = F(t + t5, XH);

        std::vector<double> 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] + K13[i] * t6;

        K13 = F(t + t6, XH);

        std::vector<double> 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 };

        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]);

        t = t + h;
    }
};

class TA : public DOPRI8 {
private:

    const double b = 0.15; // bifurcation parameter

public:
    TA(uint32_t N) : DOPRI8(N) {}

    std::vector<double> F(double h, std::vector<double> X) override {
        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;
    }

    static void Test() {
        double dt = 0.001;
        TA task(3);
        int n = 0;
        std::vector<double> X0 = { 0, 1, 0 };
        task.SetInit(0, X0);
        while (n < 20000) {
            n++;
            std::cout << n <<" Time = " << task.t << "; x = " << task.X[0] << "; y = " << task.X[1] << "; z = " << task.X[2] << '\n';
            task.NextStep(dt);
        }
    }
};



int main() {
    auto start = std::chrono::steady_clock::now();
    TA::Test();
    auto end = std::chrono::steady_clock::now();
    std::cout << "time " << std::chrono::duration_cast<std::chrono::seconds>(end - start).count() << " sec\n";
    return 0;
}