100 lines
2.7 KiB
Python
100 lines
2.7 KiB
Python
import numpy as np
|
|
from scipy.optimize import least_squares
|
|
|
|
|
|
# В эти моменты сравниваем численное решение с эталоном.
|
|
TIMES = np.linspace(0.0, 10.0, 101)
|
|
|
|
|
|
# 1. Создаём эталон ОДИН РАЗ.
|
|
def analytical(t, k1, k2):
|
|
x1 = np.exp(-k1 * t)
|
|
|
|
if k1 == k2:
|
|
x2 = k1 * t * np.exp(-k1 * t)
|
|
else:
|
|
x2 = k1 / (k2 - k1) * (
|
|
np.exp(-k1 * t) - np.exp(-k2 * t)
|
|
)
|
|
|
|
x3 = 1.0 - x1 - x2
|
|
|
|
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
|
|
|
|
return np.array([
|
|
-k1 * x1,
|
|
k1 * x1 - k2 * x2,
|
|
k2 * x2,
|
|
])
|
|
|
|
x = np.array([1.0, 0.0, 0.0])
|
|
solution = [x.copy()]
|
|
|
|
# Между соседними эталонными точками делаем 10 шагов.
|
|
# Внутренний шаг получается около 0.01.
|
|
for left, right in zip(TIMES[:-1], TIMES[1:]):
|
|
h = (right - left) / 10
|
|
|
|
for _ in range(10):
|
|
q1 = rhs(x)
|
|
q2 = rhs(x + h * q1 / 2)
|
|
q3 = rhs(x + h * q2 / 2)
|
|
q4 = rhs(x + h * q3)
|
|
|
|
x = x + h / 6 * (q1 + 2*q2 + 2*q3 + q4)
|
|
|
|
solution.append(x.copy())
|
|
|
|
return np.array(solution)
|
|
|
|
|
|
# 3. Сравниваем пробное решение с неизменным эталоном.
|
|
def residuals(parameters):
|
|
k1, k2 = parameters
|
|
|
|
calculated = solve_rk4(k1, k2)
|
|
|
|
# Обе таблицы имеют размер 101 × 3.
|
|
differences = calculated - reference
|
|
|
|
# Оптимизатор ожидает одномерный массив отклонений.
|
|
return differences.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,
|
|
)
|
|
|
|
|
|
# 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}")
|