123456789101112131415161718192021222324252627282930313233343536373839404142434445464748495051525354555657585960616263646566676869707172737475767778798081828384858687888990919293949596979899100101 |
- 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)
|