Почему решение SymPy дает пустой набор для полинома высокого порядка? - PullRequest
Новая
       52

Почему решение SymPy дает пустой набор для полинома высокого порядка?

0 голосов
/ 24 апреля 2020

Я должен найти точки равновесия, в которых нульклины пересекаются . Мой код, как показано ниже. 

>>> from sympy import symbols, Eq, solve
>>> A,M = symbols('A M')
>>> dMdt = Eq(1.05 - (1/(1 + pow(A,5))) - M)
>>> dAdt = Eq(M*1 - 0.5*A - M*A/(2 + A))
>>> solve((dMdt,dAdt), (M,A))
[]

Почему он не дает решения?

1 Ответ

1 голос
/ 24 апреля 2020

Вы поймете, почему, когда я работаю, чтобы получить решение.

Я собираюсь написать уравнения как e1 и e2 - использование Eq без второго аргумента больше не работает (или делает это с предупреждение в последних версиях SymPy): 

>>> from sympy import solve, nsimplify, factor, real_roots
>>> from sympy.abc import A, M
>>> e1 = (1.05 - (1/(1 + pow(A,5))) - M)
>>> e2 = (M*1 - 0.5*A - M*A/(2 + A))

Решить для M, используя e1 

>>> eM = solve(e1, M)[0]

Заменить на e2 

>>> e22 = e2.subs(M, eM); e22
-0.5*A - 0.05*A*(21.0*A**5 + 1.0)/((A + 2)*(A**5 + 1.0)) + 0.05*(21.0*A**5 + 1.0)/(A**5 + 1.0)

Получить числитель и знаменатель 

>>> n,d=e22.as_numer_denom()

Найдите действительные корни для этого выражения (которое зависит только от A) 

>>> rA = real_roots(n)

Найдите соответствующие значения M, подставив каждое в eM: 

>>> [(a.n(2), eM.subs(A, a).n(2)) for a in rA]
[(-3.3, 1.1), (-1.0, zoo), (-0.74, -0.23), (0.095, 0.050)]

То, что root для A = -1 является ложным - если вы посмотрите на свой знаменатель e1, вы увидите, что такое значение вызывает деление на ноль. Так что root можно игнорировать. Остальные могут быть проверены графически .

Почему не решили дать решение? Он не мог дать решение для этого многочлена высокого порядка в замкнутой форме. Даже если вы вычислите числитель, описанный выше (и превратите числа в Rationals с nsimplify), у вас есть коэффициент степени 7: 

>>> factor(nsimplify(n))
-(A + 1)*(A**4 - A**3 + A**2 - A + 1)*(5*A**7 + 10*A**6 - 21*A**5 + 5*A**2 + 10*A - 1)/10
...