начата работа над DISPF(IJK), пока что при его вызове отрабатывает DISPD

Alexey-Sagaydak · Feb 5, 2025 · 15e2547 · 15e2547
1 parent c9eed79
commit 15e2547
Show file tree

Hide file tree

Showing 5 changed files with 273 additions and 2 deletions.
diff --git a/Client/Views/AboutWindow.xaml b/Client/Views/AboutWindow.xaml
@@ -48,7 +48,7 @@
                 </StackPanel>
             </Grid>
 
-            <TextBlock Text="v1.0.2&#10;2025" 
+            <TextBlock Text="v1.1.0&#10;2025" 
                        FontSize="16" 
                        FontWeight="Bold"
                        TextAlignment="Center"

diff --git a/Server/src/Routers.cpp b/Server/src/Routers.cpp
@@ -61,7 +61,10 @@ void route::RegisterResources(hv::HttpService& router)
                     }
                     else if (method == "DISPF")
                     {
-                        throw std::runtime_error("Unsupported yet: " + method);
+                        DISPFSolver solver(
+                            odeFunction, initialStep, parameters["I"].get<int>(),
+                            parameters["J"].get<int>(), parameters["K"].get<int>());
+                        solver.Solve(t0, y0, tEnd, storage, tolerance);
                     }
                     else
                     {

diff --git 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 <sstream>
 #include <unordered_map>

diff --git a/Server/src/odesolvers-lib/include/DISPFSolver.hpp 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<std::vector<double>(
+            double,
+            std::vector<double> 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<double> &y,
+        double              &h,
+        double               tolerance) override;
+
+    void Solve(
+        double                     t0,
+        std::vector<double> 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<double> const &y,
+        double                     h,
+        std::vector<double>       &k1,
+        std::vector<double>       &k2,
+        std::vector<double>       &k3
+    );
+
+    // Оценки ошибок и параметров адаптации
+    double computeAprime(
+        std::vector<double> const &k1,
+        std::vector<double> const &k2);
+
+    double computeAddoublePrime(
+        double                     h,
+        std::vector<double> const &k1,
+        std::vector<double> const &yNext);
+
+    double computeV(
+        std::vector<double> const &k1,
+        std::vector<double> const &k2,
+        std::vector<double> 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
@@ -0,0 +1,189 @@
+#include "../include/DISPFSolver.hpp"
+
+void DISPFSolver::Step(
+    double               t,
+    std::vector<double> &y,
+    double              &h,
+    double               tolerance) 
+{
+    std::vector<double> yTemp = y;
+    std::vector<double> 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<double> yNext = AddVectors(
+            yTemp, 
+            AddVectors(
+                MultiplyScalarByVector(currentScheme.p1, k1),
+                AddVectors(
+                    MultiplyScalarByVector(currentScheme.p2, k2),
+                    MultiplyScalarByVector(currentScheme.p3, k3)
+                )
+            )
+        );
+        std::vector<double> 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<double> const &y,
+    double                     h,
+    std::vector<double>       &k1,
+    std::vector<double>       &k2,
+    std::vector<double>       &k3) 
+{
+    k1 = MultiplyScalarByVector(h, f(t, y));
+
+    std::vector<double> y2 = AddVectors(y, MultiplyScalarByVector(2.0/3.0, k1));
+    k2 = MultiplyScalarByVector(h, f(t, y2));
+
+    std::vector<double> 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<double> const &k1,
+    std::vector<double> const &k2) 
+{
+    std::vector<double> diff = SubtractVectors(k2, k1);
+    return (std::abs(1 - 6 * currentScheme.g) / 4.0) * Norm(diff);
+}
+
+double DISPFSolver::computeAddoublePrime(
+    double                     h,
+    std::vector<double> const &k1,
+    std::vector<double> const &yNext) 
+{
+    std::vector<double> diff = SubtractVectors(MultiplyScalarByVector(h, yNext), k1);
+    return (std::abs(1 - 6 * currentScheme.g) / 6.0) * Norm(diff);
+}
+
+double DISPFSolver::computeV(
+    std::vector<double> const &k1,
+    std::vector<double> const &k2,
+    std::vector<double> 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<double> const &y0,
+    double                     tEnd,
+    Storage                   &storage,
+    double                     tolerance) 
+{
+    double t = t0;
+    std::vector<double> 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;
+        }
+    }
+}