Как сделать символьные вычислениями в питоне для поиска стационарных точек системы диффуров?

Question

[chips23]

Мне надо сделать задание: напишите код на Python с использованием символьных вычислений библиотеки SymPy для того, чтобы автоматизировать поиск стационарных точек системы дифференциальных уравнений. На вход подается система из двух (любых) дифференциальных уравнений, одно из которых содержит вспомогательную функцию U(x1, x2), которая, в свою очередь, содержит вспомогательные функции psi, psi1, phi, phi' (производная от phi). Все вспомогательные функции задаются пользователем, т.е. вводятся с клавиатуры. Ваша задача - найти стационарные точки системы с использованием символьных вычислений.
То есть нужен код который будет работать с диффурами без управления и с управление (метод АКАР)
Ниже я привел код который может находить стационарные точки для диффуров, но с примером из статьи он не справляется(https://scm.etu.ru/assets/files/2025/sbornik/113-1...)
Пример такой: есть система dx1/dt = a1*x1 - b1*x1*x2
dx2/dt = -a2*x2 + b2*x1*x2 + U
где U = phi1 - a1*x1 + b1*x1*x2 - (x2 - phi)/T1 , производная от phi обозначается как phi1(она вычисляется автоматически в коде по phi)
phi(...) = a1/b1 + 1/(T2*b1) - c/(T2*b1*x1)
в коде вместо x*(целевая константа) ввожу c

import sympy as sp
import sys
# Символьные переменные
x1, x2 = sp.symbols('x1 x2')
U = None

# --- Этап 1: Ввод уравнений ---
while True:
    print("Введите правые части системы дифференциальных уравнений (возможно с управлением U):")
    f1_str = input("dx1/dt = ").strip()
    f2_str = input("dx2/dt = ").strip()

    if not f1_str or not f2_str:
        print("Ошибка: уравнения не должны быть пустыми.")
        continue

    break

# --- Этап 2: Проверка на наличие управления ---
u_in_f1 = 'U' in f1_str
u_in_f2 = 'U' in f2_str

# --- Этап 3: Ввод U и вспомогательных функций ---
if u_in_f1 or u_in_f2:
    while True:
        try:
            print("\nОбнаружено управление U. Введите выражение для U(x1, x2):")
            u_str = input("U = ")

            local_dict = {'x1': x1, 'x2': x2}
            required_helpers = []

            for helper in ['psi', 'psi1', 'phi']:
                if helper in u_str:
                    required_helpers.append(helper)

            # Ввод вспомогательных функций с проверками
            for helper in required_helpers:
                while True:
                    try:
                        print(f"\nВведите выражение для {helper}(x):")
                        helper_expr = input(f"{helper}(...) = ").strip()
                        expr = sp.sympify(helper_expr, locals=local_dict)

                        # Проверка переменных
                        if not (x1 in expr.free_symbols or x2 in expr.free_symbols):
                          raise ValueError(f"{helper} должна зависеть хотя бы от x1 или x2.")

                        local_dict[helper] = expr
                        break
                    except Exception as e:
                        print(f"Ошибка в {helper}: {e}. Повторите ввод.")

            # Если используется phi1 — вычисляем автоматически
            if 'phi1' in u_str:
                if 'phi' in local_dict:
                    phi_expr = local_dict['phi']
                    vars_in_phi = phi_expr.free_symbols & {x1, x2}
                    local_dict['phi1'] = sum(sp.diff(phi_expr, v) for v in vars_in_phi)
                    print(f"\nАвтоматически вычислена phi1 = {local_dict['phi1']}")
                else:
                    print("Ошибка: В U используется phi1, но phi не задана.")
                    continue  # вернуть пользователя назад

            # Пробуем построить U
            U = sp.sympify(u_str, locals=local_dict)
            break

        except Exception as e:
            print(f"Ошибка в выражении U: {e}. Повторите ввод.")
else:
    print("\nУправление U не используется.")

# --- Этап 4: Преобразование уравнений ---
try:
    f1 = sp.sympify(f1_str, locals={'U': U} if U else {})
    f2 = sp.sympify(f2_str, locals={'U': U} if U else {})
except Exception as e:
    print(f"Ошибка при обработке уравнений: {e}")
    sys.exit(1)

# --- Этап 5: Составление системы ---
eq1 = sp.Eq(f1, 0)
eq2 = sp.Eq(f2, 0)

# --- Этап 6: Решение ---
print("\nУравнения стационарности:")
sp.pprint(eq1)
sp.pprint(eq2)

stationary_points = sp.solve((eq1, eq2), (x1, x2), dict=True)

# --- Этап 7: Вывод ---
if stationary_points:
    print("\nСтационарные точки:")
    for i, sol in enumerate(stationary_points, 1):
        print(f"{i}) ")
        sp.pprint(sol)
else:
    print("\nСтационарные точки не найдены (возможно, нет аналитического решения).")