Решение символьной матрицы Симпи для задачи на собственные значения возвращает пустой список - PullRequest
Новая
       27

Решение символьной матрицы Симпи для задачи на собственные значения возвращает пустой список

0 голосов
/ 23 ноября 2018

Я очень новый пользователь Python.Я буду признателен за любую помощь, которую вы можете предоставить.Я пытаюсь решить проблему собственных значений.У меня есть матрица, которую я называю «data3» с элементами в терминах неизвестной переменной «omega».Определитель матрицы даст мне характеристический многочлен в терминах «омега», из которого я могу определить нужные мне корни или собственные значения.

Ниже приведен мой код: 

import sympy as sp
import numpy as np
import scipy as sc

omega=sp.Symbol('omega');
data1=sp.Matrix(sc.genfromtxt('./Z_real.csv',dtype=float,delimiter=',')); 
data2=sp.Matrix(sc.genfromtxt('./Z_imag.csv',dtype=float,delimiter=','));
data2a=data2*1j*omega;
data3=sp.Matrix(data1+data2a);

eqn=sp.Eq(sp.det(data3),0);
print(sp.solve(eqn));

Как видно выше, я импортирую data1 и data2, которые затем объединяю, чтобы сформировать матрицу комплексных чисел (с «omega», который мне нужно решить).Например, моя матрица data3 выглядит следующим образом при выводе на консоль iPython: 

Matrix([
[             0.536,  -1.119e-8*I*omega, -1.3558e-8*I*omega, -3.8778e-9*I*omega],
[ -1.119e-8*I*omega,              0.536, -3.8778e-9*I*omega, -1.3558e-8*I*omega],
[-1.3558e-8*I*omega, -3.8778e-9*I*omega,              0.536,  -1.119e-8*I*omega],
[-3.8778e-9*I*omega, -1.3558e-8*I*omega,  -1.119e-8*I*omega,              0.536]])

Однако строка кода sp.solve (eqn) возвращает для меня пустой список.То, что я ожидаю увидеть, это корни матричного полинома (из детерминанта).Может кто-нибудь посоветовать, пожалуйста, что я делаю не так?Также, если кто-нибудь покажет мне альтернативные способы решения проблемы омега, это тоже будет здорово.Если вам нужна информация о данных1 и данных2, ниже приведены матрицы, которые я использую для своего теста:

data1 

Matrix([
[0.536,   0.0,   0.0,   0.0],
[  0.0, 0.536,   0.0,   0.0],
[  0.0,   0.0, 0.536,   0.0],
[  0.0,   0.0,   0.0, 0.536]])

и data2 

Matrix([
[       0.0,  -1.119e-8, -1.3558e-8, -3.8778e-9],
[ -1.119e-8,        0.0, -3.8778e-9, -1.3558e-8],
[-1.3558e-8, -3.8778e-9,        0.0,  -1.119e-8],
[-3.8778e-9, -1.3558e-8,  -1.119e-8,        0.0]])

Большое спасибоза ваше время.

1 Ответ

0 голосов
/ 04 декабря 2018

Спасибо всем, кто прочитал пост.Я сейчас исправил эту проблему.Вместо непосредственного решения характеристического полинома матрицы, я обнаружил, что намного проще получить определитель (теперь с помощью алгоритма Берковица - без особых причин) в форме уравнения полинома, извлечь коэффициенты полинома с помощью sympy Poly, затем подключитькорни внутри NumPy корни функционируют.Работает как шарм!Пожалуйста, найдите код ниже (данные, которые я импортирую, не имеют значения, метод работает для любой матрицы sympy): 

import numpy as np
import scipy as sc
import sympy as sp
import csv

omega=sp.Symbol('omega');
data_Za=sp.Matrix(sc.genfromtxt('./Z_matrix.csv',dtype=float,delimiter=','));
L_mata=sp.Matrix(sc.genfromtxt('./L_mat.csv',dtype=float,delimiter=','));
C_mata=sp.Matrix(sc.genfromtxt('./C_mat.csv',dtype=float,delimiter=','));
N_total=int(sc.genfromtxt('./N_total.csv',dtype=float,delimiter=','));

data_Z=data_Za*omega;
L_mat=L_mata*omega;
C_mat=C_mata*(1/omega);

data_matrix=data_Z+L_mat+C_mat;
det_eqn=data_matrix.berkowitz_det();

Вот так я импортирую матрицы, добавляю их, чтобы получить конечную матрицу, которую я хочу, и вычисляюопределитель - пока ничего нового. 

det_poly_coeff=sp.Poly(det_eqn*(omega**N_total),omega);
#print(det_poly_coeff.coeffs());

Матрица, которую я использую, немного сложна, потому что, как вы можете видеть, что я сделал с C_mat, я получу термины в моем детерминанте (характеристический полином) которые имеют переменную omega даже в знаменателе.Функция Поли Симпи не нравится.Поэтому я умножаю определитель на omega**N_total, где N_total - степень полиномиального уравнения (а также наивысшая степень омеги, которая встречается в знаменателе). 

root_find_sq=np.roots(det_poly_coeff.coeffs());
#print(root_find_sq);

Включение коэффициентов в корни-пустышкидает мне корни уравнения, которые представляют собой значения, которые я хотел решить. 

root_find=np.sqrt(abs(root_find_sq));

Наконец, квадратный корень, потому что корни, которые были решены, были для omega**2, а не для omega.

Надеюсь, кому-нибудь еще, у кого может быть эта проблема, будет полезен этот пост.Приветствия.:)

...