From 15e25474caae8ad3909fceaff2208213d7249ac3 Mon Sep 17 00:00:00 2001 From: Alexey-Sagaydak Date: Wed, 5 Feb 2025 22:39:37 +0700 Subject: [PATCH] =?UTF-8?q?=D0=BD=D0=B0=D1=87=D0=B0=D1=82=D0=B0=20=D1=80?= =?UTF-8?q?=D0=B0=D0=B1=D0=BE=D1=82=D0=B0=20=D0=BD=D0=B0=D0=B4=20DISPF(IJK?= =?UTF-8?q?),=20=D0=BF=D0=BE=D0=BA=D0=B0=20=D1=87=D1=82=D0=BE=20=D0=BF?= =?UTF-8?q?=D1=80=D0=B8=20=D0=B5=D0=B3=D0=BE=20=D0=B2=D1=8B=D0=B7=D0=BE?= =?UTF-8?q?=D0=B2=D0=B5=20=D0=BE=D1=82=D1=80=D0=B0=D0=B1=D0=B0=D1=82=D1=8B?= =?UTF-8?q?=D0=B2=D0=B0=D0=B5=D1=82=20DISPD?= MIME-Version: 1.0 Content-Type: text/plain; charset=UTF-8 Content-Transfer-Encoding: 8bit --- Client/Views/AboutWindow.xaml | 2 +- Server/src/Routers.cpp | 5 +- Server/src/Routers.hpp | 1 + .../odesolvers-lib/include/DISPFSolver.hpp | 78 ++++++++ Server/src/odesolvers-lib/src/DISPFSolver.cpp | 189 ++++++++++++++++++ 5 files changed, 273 insertions(+), 2 deletions(-) create mode 100644 Server/src/odesolvers-lib/include/DISPFSolver.hpp create mode 100644 Server/src/odesolvers-lib/src/DISPFSolver.cpp diff --git a/Client/Views/AboutWindow.xaml b/Client/Views/AboutWindow.xaml index 25ac094..fccb40e 100644 --- a/Client/Views/AboutWindow.xaml +++ b/Client/Views/AboutWindow.xaml @@ -48,7 +48,7 @@ - (), + parameters["J"].get(), parameters["K"].get()); + solver.Solve(t0, y0, tEnd, storage, tolerance); } else { diff --git a/Server/src/Routers.hpp b/Server/src/Routers.hpp index 05cf2c3..e3ff54e 100644 --- a/Server/src/Routers.hpp +++ b/Server/src/Routers.hpp @@ -9,6 +9,7 @@ #include "odesolvers-lib/include/RK23SSolver.hpp" #include "odesolvers-lib/include/STEKSSolver.hpp" #include "odesolvers-lib/include/DISPDSolver.hpp" +#include "odesolvers-lib/include/DISPFSolver.hpp" #include #include diff --git a/Server/src/odesolvers-lib/include/DISPFSolver.hpp b/Server/src/odesolvers-lib/include/DISPFSolver.hpp new file mode 100644 index 0000000..0bfca2c --- /dev/null +++ b/Server/src/odesolvers-lib/include/DISPFSolver.hpp @@ -0,0 +1,78 @@ +#pragma once +#include "Solver.hpp" + +class DISPFSolver : public Solver +{ +public: + DISPFSolver( + std::function( + double, + std::vector const &)> func, + double initialStep, + double I, + double J, + double K, + double q = 1.5 // Параметр адаптации шага (q > 1) + ) : Solver(func, initialStep), I(I), J(J), K(K), q(q) {} + + void Step( + double t, + std::vector &y, + double &h, + double tolerance) override; + + void Solve( + double t0, + std::vector const &y0, + double tEnd, + Storage &storage, + double tolerance) override; + + +private: + const int I, J, K; + const double q; // Параметр для формул (5.17-5.23) + const double SAFETY = 0.8; + const double MIN_SCALE = 0.2; + const double MAX_SCALE = 5.0; + + // Параметры для схем (5.2) и (5.12) + struct SchemeParams { + double p1, p2, p3; + double g; // Параметр g (зависит от схемы) + }; + + struct SchemeParams currentScheme; + + // Вычисление k1, k2, k3 для текущей схемы + void computeKTerms( + double t, + std::vector const &y, + double h, + std::vector &k1, + std::vector &k2, + std::vector &k3 + ); + + // Оценки ошибок и параметров адаптации + double computeAprime( + std::vector const &k1, + std::vector const &k2); + + double computeAddoublePrime( + double h, + std::vector const &k1, + std::vector const &yNext); + + double computeV( + std::vector const &k1, + std::vector const &k2, + std::vector const &k3); + + // Переключение между схемами + bool shouldSwitchToSchemeB( + double h, + double tolerance); + + void switchScheme(bool toSchemeB); +}; diff --git a/Server/src/odesolvers-lib/src/DISPFSolver.cpp b/Server/src/odesolvers-lib/src/DISPFSolver.cpp new file mode 100644 index 0000000..9f18f82 --- /dev/null +++ b/Server/src/odesolvers-lib/src/DISPFSolver.cpp @@ -0,0 +1,189 @@ +#include "../include/DISPFSolver.hpp" + +void DISPFSolver::Step( + double t, + std::vector &y, + double &h, + double tolerance) +{ + std::vector yTemp = y; + std::vector k1, k2, k3; + bool stepAccepted = false; + int attempts = 0; + + while (!stepAccepted && attempts < 10) + { + computeKTerms(t, yTemp, h, k1, k2, k3); + double Aprime = computeAprime(k1, k2); + double sn = log(tolerance / (pow(q, 2*attempts) * Aprime)) / (2 * log(q)); + + if (sn < 0) + { + h *= pow(q, sn); + attempts++; + continue; + } + + // Вычисление y_{n+1} и A'' + std::vector yNext = AddVectors( + yTemp, + AddVectors( + MultiplyScalarByVector(currentScheme.p1, k1), + AddVectors( + MultiplyScalarByVector(currentScheme.p2, k2), + MultiplyScalarByVector(currentScheme.p3, k3) + ) + ) + ); + std::vector kNext = MultiplyScalarByVector(h, f(t + h, yNext)); + double AddPrime = computeAddoublePrime(h, k1, kNext); + double vn = log(tolerance / (pow(q, 2*attempts) * AddPrime)) / (2 * log(q)); + + if (vn < 0) + { + h *= pow(q, vn); + attempts++; + continue; + } + + // Проверка устойчивости через V_n + double V = computeV(k1, k2, k3); + double rn = log(18.0 / (pow(q, attempts) * V)) / log(q); + + // Адаптация шага + double scale = SAFETY * std::min( + std::min(pow(q, sn), pow(q, vn)), + std::min(pow(q, rn), MAX_SCALE) + ); + scale = std::clamp(scale, MIN_SCALE, MAX_SCALE); + + if (V > 18.0 / pow(q, attempts) || Aprime > tolerance / pow(q, 2*attempts)) + { + h *= scale; + attempts++; + } + else + { + y = yNext; + h *= scale; + stepAccepted = true; + } + } + + if (!stepAccepted) + { + throw std::runtime_error("Step adaptation failed"); + } + + // Проверка переключения на схему Б + if (shouldSwitchToSchemeB(h, tolerance)) + { + switchScheme(true); + } +} + +void DISPFSolver::computeKTerms( + double t, + std::vector const &y, + double h, + std::vector &k1, + std::vector &k2, + std::vector &k3) +{ + k1 = MultiplyScalarByVector(h, f(t, y)); + + std::vector y2 = AddVectors(y, MultiplyScalarByVector(2.0/3.0, k1)); + k2 = MultiplyScalarByVector(h, f(t, y2)); + + std::vector y3 = AddVectors( + AddVectors(y, MultiplyScalarByVector(1.0/3.0, k1)), + MultiplyScalarByVector(1.0/3.0, k2) + ); + k3 = MultiplyScalarByVector(h, f(t, y3)); +} + +double DISPFSolver::computeAprime( + std::vector const &k1, + std::vector const &k2) +{ + std::vector diff = SubtractVectors(k2, k1); + return (std::abs(1 - 6 * currentScheme.g) / 4.0) * Norm(diff); +} + +double DISPFSolver::computeAddoublePrime( + double h, + std::vector const &k1, + std::vector const &yNext) +{ + std::vector diff = SubtractVectors(MultiplyScalarByVector(h, yNext), k1); + return (std::abs(1 - 6 * currentScheme.g) / 6.0) * Norm(diff); +} + +double DISPFSolver::computeV( + std::vector const &k1, + std::vector const &k2, + std::vector const &k3) +{ + double maxRatio = 0.0; + for (size_t i = 0; i < k1.size(); ++i) + { + double denominator = k2[i] - k1[i]; + denominator = (std::abs(denominator) < 1e-12) ? 1e-12 : denominator; + double ratio = std::abs((k3[i] - k2[i]) / denominator); + maxRatio = std::max(maxRatio, ratio); + } + + return 3.0 * maxRatio; +} + +bool DISPFSolver::shouldSwitchToSchemeB( + double h, + double tolerance) +{ + // Условие переключения: прогнозируемый шаг для схемы Б больше + return (h * MAX_SCALE > currentScheme.g * tolerance); +} + +void DISPFSolver::switchScheme(bool toSchemeB) +{ + if (toSchemeB) + { + currentScheme = {7.0/9.0, 16.0/81.0, 2.0/81.0, 0.0}; // Параметры схемы Б + } + else + { + currentScheme = {1.0/4.0, (3 - 18*currentScheme.g)/4, 9.0/2.0 * currentScheme.g, 0.15}; + } +} + +void DISPFSolver::Solve( + double t0, + std::vector const &y0, + double tEnd, + Storage &storage, + double tolerance) +{ + double t = t0; + std::vector y = y0; + double h = stepSize; + storage.Add(t, y); + switchScheme(false); // Начинаем с алгоритма А + + while (t < tEnd) + { + double hAttempt = std::min(h, tEnd - t); + if (hAttempt < 1e-12) break; + + try + { + Step(t, y, hAttempt, tolerance); + t += hAttempt; + storage.Add(t, y); + h = hAttempt; + } + catch (const std::runtime_error&) + { + h *= 0.5; + } + } +}