2016-11-17 2 views
1

я получаю некоторые собственные значения матрицыSymPy: упростить большее выражение с биномиальной формулой и квадратичным дополнением

import sys 
import mpmath 
from sympy import * 

X,Y,Z = symbols("X,Y,Z") 
Rxy,Rxz, Ry,Ryx,Ryz, Rz,Rzy,Rzz = symbols("Rxy,Rxz, Ry,Ryx,Ryz, Rz,Rzy,Rzz") 


J = Matrix([ 
     [  -1,  0,  0], 
     [  0, -Ry*Y, Ry*Rzy*Y], 
     [Rxz*Rz*Z, Ryz*Rz*Z, -Rz*Z]]) 

, которые являются следующим:

{-Ry*Y/2 - Rz*Z/2 + sqrt(Ry**2*Y**2 + 4*Ry*Ryz*Rz*Rzy*Y*Z - 2*Ry*Rz*Y*Z + Rz**2*Z**2)/2: 1, 
-Ry*Y/2 - Rz*Z/2 - sqrt(Ry**2*Y**2 + 4*Ry*Ryz*Rz*Rzy*Y*Z - 2*Ry*Rz*Y*Z + Rz**2*Z**2)/2: 1, 
-1: 1} 

позволяет просто смотреть на собственном значении одного:

In [25]: J.eigenvals().keys()[0] 
Out[25]: -Ry*Y/2 - Rz*Z/2 + sqrt(Ry**2*Y**2 + 4*Ry*Ryz*Rz*Rzy*Y*Z - 2*Ry*Rz*Y*Z + Rz**2*Z**2)/2 

Я хочу упростить этот термин следующим образом: коэффициент из 1/2 и (это важно) radicant.

я могу преобразовать radicant следующего путем добавления квадратичного дополнения

Ry**2*Y**2 + 4*Ry*Ryz*Rz*Rzy*Y*Z - 2*Ry*Rz*Y*Z + Rz**2*Z**2  | + 4*Ry*Rz*Y*Z -4*Ry*Rz*Y*Z 

, что приводит к

Ry**2*Y**2 + Rz**2*Z**2 + 2*Ry*Rz*Y*Z - 4*Ry*Rz*Y*Z + 4*Ry*Ryz*Rz*Rzy*Y*Z 

, которые могут быть факторизуются в

(Ry*Y + Rz*Z)**2 - 4*Ry*Rz*Y*Z*(1 - Ryz*Rzy) 

с этими оценками полноцикловыми Собственное значение должно выглядеть так:

-1/2*(Ry*Y + Rz*Z - sqrt((Ry*Y + Rz*Z)**2 - 4*Ry*Rz*Y*Z*(1 - Ryz*Rzy))) 

Этот расчет очень важен для меня, потому что я должен оценить, является ли собственное значение < 0. Это намного проще в последней форме.

Позвольте мне показать вам, что я делал до сих пор.

In [24]: J.eigenvals().keys()[0] 
Out[24]: -Ry*Y/2 - Rz*Z/2 + sqrt(Ry**2*Y**2 + 4*Ry*Ryz*Rz*Rzy*Y*Z - 2*Ry*Rz*Y*Z + Rz**2*Z**2)/2 

In [25]: J.eigenvals().keys()[0].factor() 
Out[25]: -(Ry*Y + Rz*Z - sqrt(Ry**2*Y**2 + 4*Ry*Ryz*Rz*Rzy*Y*Z - 2*Ry*Rz*Y*Z + Rz**2*Z**2))/2 

In [26]: J.eigenvals().keys()[0].simplify() 
Out[26]: -Ry*Y/2 - Rz*Z/2 + sqrt(Ry**2*Y**2 + 4*Ry*Ryz*Rz*Rzy*Y*Z - 2*Ry*Rz*Y*Z + Rz**2*Z**2)/2 

So .simplify() не меняет результат вообще. .factor() просто исключает -1/2. Если я правильно помню, я могу передать аргумент в .factor() как Y или Z, какая переменная должна быть факторизована. Но я получаю много немного разных собственных значений в качестве вывода, и я не хочу указывать каждый аргумент фактора() вручную (если это решение даже работает).

Я также попытался вычислить собственные значения самостоятельно, вычислив определитель и решив определение == 0 ... Я также использовал метод детерминат.фактор() и решил его впоследствии, но лучший результат этого подхода был таким же, как J.eigenvals(). клавиши() [0] .factor().

У вас есть идеи, как решить эту проблему?

Спасибо заранее

Alex

+0

'factor' в SymPy означает« фактор в произведение многочленов ». Таким образом, это будет работать, только если вы дадите ему подмножество терминов, которые влияют на продукт. – asmeurer

ответ

1

Такого рода вещи просят много (см также, например, этот вопрос: Expression simplification in SymPy), но есть на самом деле не является хорошим способом в SymPy сделать это. Проблема в том, что такие «частичные» факторизации не уникальны (может быть несколько способов преобразования многочлена в сумму произведений).

Я открыл this issue об этом в журнале отслеживания проблем SymPy.Я показал там так, что вы можете получить близко (здесь a это термин, под квадратным корнем)

In [92]: collect(expand(a.subs(Ry*Y, x - Rz*Z)), x, func=factor).subs(x, Ry*Y + Rz*Z) 
Out[92]: 
     2 2                 2 
- 4⋅Rz ⋅Z ⋅(Ryz⋅Rzy - 1) + 4⋅Rz⋅Z⋅(Ry⋅Y + Rz⋅Z)⋅(Ryz⋅Rzy - 1) + (Ry⋅Y + Rz⋅Z) 

Здесь я временно заменяющий Ry*Y + Rz*Z с переменной x, так что я могу получить квадрат термин, который вы хотите.

Я не мог найти способ приблизиться к тому, что вы хотите (т. Е. Коэффициент Ryz*Rzy - 1 из оставшихся условий).

+0

Хорошо, спасибо за ваш ответ. К сожалению, я пишу сценарий, где результирующее выражение время от времени отличается. Ну, я думаю, мне нужно упростить свои выражения вручную. - Благодарю вас за помощь –

Смежные вопросы