python_bio_sem_6/main.py

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)