Advertisement
Not a member of Pastebin yet?
Sign Up,
it unlocks many cool features!
- #ifndef ODE
- #define ODE
- #include <cmath>
- #include <iostream>
- #include <fstream>
- #include <array>
- const int N = 3; // размерность системы
- const int number_of_generated_points = 1000000;
- 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;
- //}
- #endif // !ODE
Advertisement
Add Comment
Please, Sign In to add comment
Advertisement