Update optimize.py

scripts/optimize.py

@@ -2,11 +2,11 @@ import numpy as np
 from scipy.optimize import least_squares
 
 
-# В эти моменты сравниваем численное решение с эталоном.
+REFERENCE_K1 = 0.577
+REFERENCE_K2 = 0.422
 TIMES = np.linspace(0.0, 10.0, 101)
+SUBSTEPS = 10
 
-
-# 1. Создаём эталон ОДИН РАЗ.
 def analytical(t, k1, k2):
     x1 = np.exp(-k1 * t)
 
@@ -22,11 +22,8 @@ def analytical(t, k1, k2):
     return np.column_stack((x1, x2, x3))
 
 
-reference = analytical(TIMES, 0.577, 0.422)
-
-
-# 2. Численный расчёт для ЛЮБОЙ пробной пары k1, k2.
 def solve_rk4(k1, k2):
+
     def rhs(x):
         x1, x2, x3 = x
 
@@ -36,64 +33,79 @@ def solve_rk4(k1, k2):
             k2 * x2,
         ])
 
-    x = np.array([1.0, 0.0, 0.0])
-    solution = [x.copy()]
+    state = np.array([1.0, 0.0, 0.0])
 
-    # Между соседними эталонными точками делаем 10 шагов.
-    # Внутренний шаг получается около 0.01.
-    for left, right in zip(TIMES[:-1], TIMES[1:]):
-        h = (right - left) / 10
+    solution = np.empty((len(TIMES), 3))
+    solution[0] = state
 
-        for _ in range(10):
-            q1 = rhs(x)
-            q2 = rhs(x + h * q1 / 2)
-            q3 = rhs(x + h * q2 / 2)
-            q4 = rhs(x + h * q3)
+    for j in range(1, len(TIMES)):
+        h = (TIMES[j] - TIMES[j - 1]) / SUBSTEPS
 
-            x = x + h / 6 * (q1 + 2*q2 + 2*q3 + q4)
+        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)
 
-        solution.append(x.copy())
+            state = state + h / 6.0 * (
+                q1 + 2.0 * q2 + 2.0 * q3 + q4
+            )
 
-    return np.array(solution)
+        solution[j] = state
+
+    return solution
 
 
-# 3. Сравниваем пробное решение с неизменным эталоном.
-def residuals(parameters):
+def residuals(parameters, reference):
     k1, k2 = parameters
 
     calculated = solve_rk4(k1, k2)
 
-    # Обе таблицы имеют размер 101 × 3.
-    differences = calculated - reference
-
-    # Оптимизатор ожидает одномерный массив отклонений.
-    return differences.ravel()
+    return (calculated - reference).ravel()
 
 
-# 4. Подбираем параметры.
-result = least_squares(
-    fun=residuals,
-    x0=[0.3, 0.7],           # Начальная пробная пара.
-    bounds=(
-        [1e-5, 0.01],       # Нижние границы.
-        [100.0, 1.0],       # Верхние границы.
-    ),
-    method="trf",
-    max_nfev=200,
-    verbose=1,
-)
+def main():
+    reference = analytical(TIMES, REFERENCE_K1, REFERENCE_K2)
+
+    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,
+    )
+
+    k1_found, k2_found = result.x
+
+    absolute_error_k1 = abs(k1_found - REFERENCE_K1)
+    absolute_error_k2 = abs(k2_found - REFERENCE_K2)
+
+    relative_error_k1 = absolute_error_k1 / REFERENCE_K1 * 100
+    relative_error_k2 = absolute_error_k2 / REFERENCE_K2 * 100
+
+    sum_squared_errors = np.sum(result.fun ** 2)
+    rmse = np.sqrt(np.mean(result.fun ** 2))
+
+    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Поиск завершён без подтверждённой сходимости.")
 
 
-# 5. Проверяем найденные параметры.
-k1_found, k2_found = result.x
-
-print(f"\nУспешное завершение: {result.success}")
-print(f"Причина остановки: {result.message}")
-
-print(f"k1 = {k1_found:.10f}")
-print(f"k2 = {k2_found:.10f}")
-
-print(f"Отклонение k1: {abs(k1_found - 0.577):.6e}")
-print(f"Отклонение k2: {abs(k2_found - 0.422):.6e}")
-
-print(f"Сумма квадратов отклонений: {np.sum(result.fun**2):.6e}")
+if __name__ == "__main__":
+    main()