Advertisement
Not a member of Pastebin yet?
Sign Up,
it unlocks many cool features!
- #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;
- }
Advertisement
Add Comment
Please, Sign In to add comment
Advertisement