|
@@ -0,0 +1,101 @@
|
|
|
+import numpy as np
|
|
|
+import matplotlib.pyplot as plt
|
|
|
+
|
|
|
+
|
|
|
+# Решатель ОДУ методом Хёйна
|
|
|
+def solve_ode(start, end, step, max_calls, tolerance, fs, initial_conditions):
|
|
|
+ """
|
|
|
+ Решить систему ОДУ методом Хёйна.
|
|
|
+
|
|
|
+ Параметры:
|
|
|
+ start: начальное время
|
|
|
+ end: конечное время
|
|
|
+ step: начальный шаг интегрирования
|
|
|
+ max_calls: максимальное количество вызовов функции fs
|
|
|
+ tolerance: точность решения
|
|
|
+ fs: функция-правая часть системы ОДУ
|
|
|
+ initial_conditions: начальные условия
|
|
|
+ """
|
|
|
+
|
|
|
+ # Переменные для хранения времени и вектора решения
|
|
|
+ t = start
|
|
|
+ v = np.array(initial_conditions)
|
|
|
+
|
|
|
+ # Счетчик вызовов функции fs
|
|
|
+ call_counter = [0]
|
|
|
+
|
|
|
+ # Вывести начальные значения
|
|
|
+ print(f"{t_0:13.6f}{step:13.6f}{0:13d}{0:13d}", *[f"{x:12.6f}" for x in v])
|
|
|
+
|
|
|
+ # Определение шага Хёйна
|
|
|
+ def heun_step(t, v, h):
|
|
|
+ """
|
|
|
+ Вычислить шаг Хёйна.
|
|
|
+
|
|
|
+ Параметры:
|
|
|
+ t: текущее время
|
|
|
+ v: текущее значение вектора решения
|
|
|
+ h: шаг интегрирования
|
|
|
+ """
|
|
|
+
|
|
|
+ k1 = fs(t, v, call_counter)
|
|
|
+ k2 = fs(t + h, v + h * k1, call_counter)
|
|
|
+ return v + (h / 2) * (k1 + k2)
|
|
|
+
|
|
|
+ # Цикл интегрирования
|
|
|
+ while t < end and call_counter[0] < max_calls:
|
|
|
+ # Вычислить коэффициенты Хёйна
|
|
|
+ k1 = fs(t, v, call_counter)
|
|
|
+ k2 = fs(t + step, v + step * k1, call_counter)
|
|
|
+ v1 = v + (step / 2) * (k1 + k2)
|
|
|
+
|
|
|
+ k2 = fs(t + step / 2, v + step / 2 * k1, call_counter)
|
|
|
+ v2 = v + (step / 4) * (k1 + k2)
|
|
|
+
|
|
|
+ # Вычислить следующий шаг Хёйна
|
|
|
+ v2 = heun_step(t + step / 2, v2, step / 2)
|
|
|
+
|
|
|
+ # Оценка ошибки
|
|
|
+ error = np.linalg.norm(v2 - v1) / 3
|
|
|
+ # Адаптация шага
|
|
|
+ if error > tolerance:
|
|
|
+ step /= 2
|
|
|
+
|
|
|
+ elif error < tolerance / 64:
|
|
|
+ step *= 2
|
|
|
+
|
|
|
+ # Если ошибка приемлема, то сделать шаг и вывести решение
|
|
|
+ if error < tolerance:
|
|
|
+ if t + step > end:
|
|
|
+ step = end - t
|
|
|
+
|
|
|
+ t += step
|
|
|
+ v = v1
|
|
|
+ print(f"{t:13.6f} {step:13.6f} {error:13.5e} {call_counter[0]:13d}", *[f"{x:12.6f}" for x in v])
|
|
|
+
|
|
|
+
|
|
|
+if __name__ == "__main__":
|
|
|
+ # Считать входные данные
|
|
|
+ t_0 = float(input())
|
|
|
+ T = float(input())
|
|
|
+ h_0 = float(input())
|
|
|
+ N_x = int(input())
|
|
|
+ eps = float(input())
|
|
|
+ n = int(input())
|
|
|
+
|
|
|
+ # Считать код функции-правой части системы ОДУ
|
|
|
+ function_code = []
|
|
|
+ for _ in range(n + 3):
|
|
|
+ line = input()
|
|
|
+ function_code.append(line)
|
|
|
+
|
|
|
+ # Определить функцию-правую часть системы ОДУ
|
|
|
+ function_definition = "\n".join(function_code)
|
|
|
+ exec(function_definition)
|
|
|
+
|
|
|
+ # Считать начальные условия
|
|
|
+ initial_conditions_str = input()
|
|
|
+ initial_conditions = [float(x) for x in initial_conditions_str.split()]
|
|
|
+
|
|
|
+ # Решить систему ОДУ
|
|
|
+ solve_ode(t_0, T, h_0, N_x, eps, fs, initial_conditions)
|