Refactoring

CMakeLists.txt

@@ -97,4 +97,3 @@ add_mmps_plugin(
 # ПРИМЕР ЗАПУСКА ИЗ КАТАЛОГА build
 #   bin/Release/mmps_host.exe  (Visual Studio)
 #   bin/mmps_host              (Linux/Makefiles)
-# Каталог плагинов и значения k1, k2, h программа запросит через std::cin.

include/api.hpp

@@ -30,7 +30,7 @@ struct OutputParameters{
     size_t point_count;
 };
 
-inline constexpr std::uint32_t PLUGIN_API_VERSION = 2;
+inline constexpr std::uint32_t PLUGIN_API_VERSION = 3;
 
 #if defined(_WIN32) && defined(PLUGIN_BUILD)
 #define PLUGIN_EXPORT extern "C" __declspec(dllexport)
@@ -56,4 +56,3 @@ inline constexpr const char* SYMBOL_NAME = "plugin_name";
 inline constexpr const char* SYMBOL_DESCRIPTION = "plugin_description";
 inline constexpr const char* SYMBOL_SOLVER = "solver";
 inline constexpr const char* SYMBOL_WRITE_OUTPUT = "write_output";
-

include/state.hpp

@@ -28,8 +28,7 @@ inline State derivative(const State& state,
     };
 }
 
-inline State runge_kutta_step(const State& state,
-                              const SolverParameters& parameters){
+inline State runge_kutta_step(const State& state, const SolverParameters& parameters){
     double h = parameters.step;
 
     State k1 = derivative(state, parameters);

scripts/optimize.py

@@ -1,20 +1,32 @@
 import numpy as np
 from scipy.optimize import least_squares
 
+# Эталонные значения параметров
+K1_REFERENCE = 0.577
+K2_REFERENCE = 0.422
 
-REFERENCE_K1 = 0.577
-REFERENCE_K2 = 0.422
-TIMES = np.linspace(0.0, 10.0, 101)
-SUBSTEPS = 10
+K1_BOUNDS = (1e-5, 1e2)
+K2_BOUNDS = (1e-2, 1.0)
+
+T_BEGIN = 0.0
+T_END = 10.0
+H = 0.1
+
+
+def analytical_solution(t, k1, k2):
+    """
+    Аналитическое решение задачи Коши.
+    """
+    t = np.asarray(t, dtype=float)
 
-def analytical(t, k1, k2):
     x1 = np.exp(-k1 * t)
 
-    if k1 == k2:
+    if np.isclose(k1, k2):
         x2 = k1 * t * np.exp(-k1 * t)
     else:
-        x2 = k1 / (k2 - k1) * (
-            np.exp(-k1 * t) - np.exp(-k2 * t)
+        x2 = (
+            k1 / (k2 - k1)
+            * (np.exp(-k1 * t) - np.exp(-k2 * t))
         )
 
     x3 = 1.0 - x1 - x2
@@ -22,90 +34,97 @@ def analytical(t, k1, k2):
     return np.column_stack((x1, x2, x3))
 
 
-def solve_rk4(k1, k2):
+def derivative(state, k1, k2):
+    """
+    Правая часть системы дифференциальных уравнений.
+    """
+    x1, x2, x3 = state
 
-    def rhs(x):
-        x1, x2, x3 = x
+    return np.array([
+        -k1 * x1,
+        k1 * x1 - k2 * x2,
+        k2 * x2
+    ])
 
-        return np.array([
-            -k1 * x1,
-            k1 * x1 - k2 * x2,
-            k2 * x2,
-        ])
 
-    state = np.array([1.0, 0.0, 0.0])
+def rk4_step(state, h, k1, k2):
+    """
+    Один шаг метода Рунге--Кутты четвёртого порядка.
+    """
+    q1 = derivative(state, k1, k2)
+    q2 = derivative(state + h * q1 / 2.0, k1, k2)
+    q3 = derivative(state + h * q2 / 2.0, k1, k2)
+    q4 = derivative(state + h * q3, k1, k2)
 
-    solution = np.empty((len(TIMES), 3))
-    solution[0] = state
+    return state + h / 6.0 * (q1 + 2.0 * q2 + 2.0 * q3 + q4)
 
-    for j in range(1, len(TIMES)):
-        h = (TIMES[j] - TIMES[j - 1]) / SUBSTEPS
 
-        for _ in range(SUBSTEPS):
-            q1 = rhs(state)
-            q2 = rhs(state + h * q1 / 2.0)
-            q3 = rhs(state + h * q2 / 2.0)
-            q4 = rhs(state + h * q3)
+def rk4_solution(t, k1, k2):
+    """
+    Решение задачи Коши методом RK4.
+    """
+    solution = np.zeros((len(t), 3), dtype=float)
+    solution[0] = [1.0, 0.0, 0.0]
 
-            state = state + h / 6.0 * (
-                q1 + 2.0 * q2 + 2.0 * q3 + q4
-            )
+    for index in range(1, len(t)):
+        h = t[index] - t[index - 1]
 
-        solution[j] = state
+        solution[index] = rk4_step(
+            solution[index - 1],
+            h,
+            k1,
+            k2
+        )
 
     return solution
 
+time = np.linspace(T_BEGIN, T_END, 101)
 
-def residuals(parameters, reference):
+# Эталонные данные
+reference_solution = analytical_solution(time, K1_REFERENCE, K2_REFERENCE)
+
+def residuals(parameters):
+    """
+    Вектор невязок между моделью и эталонными данными.
+    """
     k1, k2 = parameters
 
-    calculated = solve_rk4(k1, k2)
+    numerical_solution = rk4_solution(time, k1, k2)
 
-    return (calculated - reference).ravel()
+    # Получается вектор из 101 * 3 = 303 невязок
+    return (numerical_solution - reference_solution).ravel()
 
 
-def main():
-    reference = analytical(TIMES, REFERENCE_K1, REFERENCE_K2)
+# Начальное приближение
+initial_parameters = np.array([0.3, 0.7])
 
-    result = least_squares(
-        fun=residuals,
-        x0=[0.3, 0.7],
-        args=(reference,),
-        bounds=(
-            [1e-5, 0.01],
-            [100.0, 1.0],
-        ),
-        method="trf",
-        max_nfev=200,
-        verbose=1,
+result = least_squares(
+    residuals,
+    initial_parameters,
+    bounds=(
+        [K1_BOUNDS[0], K2_BOUNDS[0]],
+        [K1_BOUNDS[1], K2_BOUNDS[1]]
     )
+)
 
-    k1_found, k2_found = result.x
+# result.x - solution
+found_k1, found_k2 = result.x
 
-    absolute_error_k1 = abs(k1_found - REFERENCE_K1)
-    absolute_error_k2 = abs(k2_found - REFERENCE_K2)
+absolute_error_k1 = abs(found_k1 - K1_REFERENCE)
+absolute_error_k2 = abs(found_k2 - K2_REFERENCE)
 
-    relative_error_k1 = absolute_error_k1 / REFERENCE_K1 * 100
-    relative_error_k2 = absolute_error_k2 / REFERENCE_K2 * 100
+relative_error_k1 = absolute_error_k1 / K1_REFERENCE * 100.0
+relative_error_k2 = absolute_error_k2 / K2_REFERENCE * 100.0
 
-    sum_squared_errors = np.sum(result.fun ** 2)
-    rmse = np.sqrt(np.mean(result.fun ** 2))
+print("Results")
+print(f"Found k1:{found_k1:.4f}")
+print(f"Theoretical k1:{K1_REFERENCE:.4f}")
+print(f"Absolute error of k1:{absolute_error_k1:.5e}")
+print(f"Relative error of k1:{relative_error_k1:.5e} %")
+print()
 
-    print(f"\nk1 = {k1_found:.4f}")
-    print(f"k2 = {k2_found:.4f}")
-
-    print(f"\nАбсолютное отклонение k1: {absolute_error_k1:.6e}")
-    print(f"Абсолютное отклонение k2: {absolute_error_k2:.6e}")
-
-    print(f"\nОтносительное отклонение k1: {relative_error_k1:.6e}%")
-    print(f"Относительное отклонение k2: {relative_error_k2:.6e}%")
-
-    print(f"\nСумма квадратов отклонений: {sum_squared_errors:.6e}")
-    print(f"RMSE: {rmse:.6e}")
-
-    if not result.success:
-        print("\nПоиск завершён без подтверждённой сходимости.")
-
-
-if __name__ == "__main__":
-    main()
+print(f"Found k2:{found_k2:.5f}")
+print(f"Theoretical k2:{K2_REFERENCE:.5f}")
+print(f"Absolute error of k2:{absolute_error_k2:.5e}")
+print(f"Relative error of k2:{relative_error_k2:.5e} %")
+print()

scripts/results.py

@@ -0,0 +1,50 @@
+import pandas as pd
+import matplotlib.pyplot as plt
+
+
+def export_component_comparison(
+    rk4_path: str,
+    ab2_path: str,
+    component: str,
+    ylabel: str,
+    out_name: str,
+):
+    """
+    Строит график зависимости одной компоненты x_i(t) от времени,
+    на котором результаты метода Рунге--Кутты и метода Адамса--Башфорта
+    представлены на общей координатной плоскости (п.4 базовой части задания).
+    """
+    df_rk4 = pd.read_csv(rk4_path)
+    df_ab2 = pd.read_csv(ab2_path)
+
+    t_rk4, y_rk4 = df_rk4["t"].values, df_rk4[component].values
+    t_ab2, y_ab2 = df_ab2["t"].values, df_ab2[component].values
+
+    plt.figure(figsize=(10, 6))
+
+    plt.plot(t_rk4, y_rk4, label="Рунге--Кутты 4-го порядка",
+              linewidth=2.2, color="#1f77b4")
+    plt.plot(t_ab2, y_ab2, label="Адамса--Башфорта 2-го порядка",
+              linewidth=1.8, color="#ff7f0e", linestyle="--")
+
+    plt.xlabel("t", fontsize=12)
+    plt.ylabel(ylabel, fontsize=12)
+    plt.grid(True, which="both", linestyle="--", alpha=0.6)
+    plt.legend(fontsize=11, loc="best", framealpha=0.9)
+
+    t_min = min(t_rk4.min(), t_ab2.min())
+    t_max = max(t_rk4.max(), t_ab2.max())
+    plt.xlim(t_min - 0.1, t_max + 0.1)
+    plt.ylim(-0.02, 1.05)
+
+    plt.tight_layout()
+    plt.savefig(f"figures/{out_name}.png", dpi=300)
+    plt.close()
+
+
+RK4_CSV = "../build/bin/results/runge_kutta_plugin.csv"
+AB2_CSV = "../build/bin/results/adams_bashforth_plugin.csv"
+
+export_component_comparison(RK4_CSV, AB2_CSV, "x1", r"$x_1(t)$", "x1_comparison")
+export_component_comparison(RK4_CSV, AB2_CSV, "x2", r"$x_2(t)$", "x2_comparison")
+export_component_comparison(RK4_CSV, AB2_CSV, "x3", r"$x_3(t)$", "x3_comparison")