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)