Алгоритм диагонализации матрицы

[B3LYP]

Мне нужно научиться решать ядерное уравнение Шредингера в гармоническом приближении. Есть набор силовых постоянных - вторых производных энергии молекулы по декартовым координатам атомов (первые производные должны быть нулевыми), по этой информации можно рассчитать частоты колебаний. Мне сказали, что матрицу силовых постоянных нужно диагонализировать. Я пытался это сделать и столкнулся с проблемами:

topic158094.html

Мой алгоритм диагонализации выдал другие значения, чем алгоритм, который madschumacher-у предоставила программа NumPy.
Суть проблемы кажется я могу сформулировать. Я не понимаю, что такое "диагонализация матрицы", потому что не знаю, куда делся в моей задаче столбец справа от матрицы. Т.е. я со своим алгоритмом привык диагонализировать матрицу для решения системы линейных уравнений, у них есть этот самый столбец справа (подскажите как он называется по русски и по английски), который при диагонализации тоже меняется. А тут не то.
Мой алгоритм работает так. Сначала находим самый большой по модулю элемент в первом столбце, и меняем строку с этим элементом на первую строку. Далее прибавляем ко всем строкам, начиная со второй, первую, умноженную на коэффициент чтобы в первом столбце везде был ноль (кроме первой строки разумеется). Далее повторяем это же со втором столбцом: сначала находим самый большой элемент от второй строки до последней, меняем эту строку на вторую, и далее прибавляем вторую строку к первой и третьей, четвёртой,...последней, так чтобы везде во втором столбце кроме второй строки остались нули. И так далее.

[dgwuqtj]

Так это разные диагонализации. Есть приведение матрицы $A$ (не обязательно квадратной) к ступенчатой форме с помощью домножения с одной стороны (скажем, слева) на обратимую матрицу. Это, кстати, над произвольным полем. А есть приведение квадратной матрицы $A$ к жордановой нормальной форме преобразованиями вида $A \mapsto U A U^{-1}$ для обратимой $U$, уже только над алгебраически замкнутым полем. В общем случае обе эти стандартные формы не являются диагональными. Ступенчатая форма диагональна тогда и только тогда, когда $A$ обратима. А жорданова нормальная форма диагональна тогда и только тогда, когда $A$ полупростая (по определению полупростоты). Есть ещё другие "диагонализации", например, для симметричных (или эрмитовых) матриц...

[Missir]

B3LYP в сообщении #1647703 писал(а):

Я не понимаю, что такое "диагонализация матрицы"

и не нужно, просто решайте проблему собственных значений, правая часть для этого не нужна

[ozheredov]

B3LYP в сообщении #1647703 писал(а):

Мой алгоритм диагонализации

это здорово, но почему бы для начала не использовать нечто "из коробки" - типа lapack например (numpy видимо не подойдет, поскольку ТС пишет свою исполняемую прогу, а не скрипт на Питоне).

Посмотреть, как будет работать весь алгоритм в сборе (диагонализация это как я понимаю всего лишь маленький винтик у этого механизма), а потом перейти к написанию своего модуля диагонализации (если не нравится как он работает или просто для души)

[B3LYP]

dgwuqtj в сообщении #1647705 писал(а):

Так это разные диагонализации. Есть приведение матрицы $A$ (не обязательно квадратной) к ступенчатой форме с помощью домножения с одной стороны (скажем, слева) на обратимую матрицу. Это, кстати, над произвольным полем. А есть приведение квадратной матрицы $A$ к жордановой нормальной форме преобразованиями вида $A \mapsto U A U^{-1}$ для обратимой $U$, уже только над алгебраически замкнутым полем. В общем случае обе эти стандартные формы не являются диагональными. Ступенчатая форма диагональна тогда и только тогда, когда $A$ обратима. А жорданова нормальная форма диагональна тогда и только тогда, когда $A$ полупростая (по определению полупростоты). Есть ещё другие "диагонализации", например, для симметричных (или эрмитовых) матриц...

Извиняюсь, совсем не понимаю. Но я знаю что матрица силовых постоянных симметрична, т.е. элементы Uij и Uji равны. Какие бывают диагонализации для таких матриц?
И как называется диагонализация которую я описал?

-- 29.07.2024, 16:31 --

Missir в сообщении #1647706 писал(а):

и не нужно, просто решайте проблему собственных значений, правая часть для этого не нужна

Поясните? Собственные значения это элементы на диагонали, которые будут после диагонализации?

-- 29.07.2024, 16:32 --

ozheredov в сообщении #1647712 писал(а):

это здорово, но почему бы для начала не использовать нечто "из коробки" - типа lapack например (numpy видимо не подойдет, поскольку ТС пишет свою исполняемую прогу, а не скрипт на Питоне).

Посмотреть, как будет работать весь алгоритм в сборе (диагонализация это как я понимаю всего лишь маленький винтик у этого механизма), а потом перейти к написанию своего модуля диагонализации (если не нравится как он работает или просто для души)

Я "велосипедист", мне проще и интереснее самому закодить алгоритм, чем искать готовый.

[Missir]

B3LYP в сообщении #1647742 писал(а):

Поясните? Собственные значения это элементы на диагонали, которые будут после диагонализации?

в посте madschumacher подробно описано, что Вам нудно делать. Он там пишет "диагонализовать матрицу", но приведённая ниже формула обозначает проблему собственных значений. Вам нужно просто решить эту проблему и всё. Действительно, в случае симметричной положительно определённой матрицы, её можно привести к диагональному виду умножив слева и справа на матрицы собственных векторов. Но Вам то это всё зачем?, зачем усложнять задачу? Вам же нужно просто найти собственные числа матрицы. Имейте в виду, не любая диагонализация даст Вам на диагонали собственные числа. Никакое шаманство тут не пройдёт, нужно использовать корректные методы решения проблемы собственных значений.

-- 29.07.2024, 16:48 --

B3LYP, если хотите сами запрограммировать, то смотрите в сторону QR алгоритма

[Евгений Машеров]

Диагонализация не сводится к оставлению у матрицы только диагональных элементов, это было бы слишком просто и слишком бесполезно. Диагонализирующие операторы должны сохранять какое-то важное для определённой задачи свойство матрицы.
Для системы линейных уравнений $Ax=b$ действует оператор, одинаково прилагаемый к A и к b, так что равенство сохраняется.
В Вашей задаче отыскиваются экстремумы некоторой квадратичной формы, и допустимые операторы те, которые сохраняют эти экстремумы. Например, домножение справа и слева на ортогональную матрицу. На этом основан один из методов такого преобразования - Якоби. На каждом шаге домножают на матрицу вращения так, чтобы определённый элемент занулялся. После некоторого числа шагов диагональные элементы становятся малы, и матрица "диагонализирована". Если Вы желаете сами такую процедуру запрограммировать, рекомендовал бы начать с него, он достаточно прост. Хотя QR-метод быстрее, особенно если нам нужны только собственные значения, а не вектора.

[Alex Krylov]

Ваша "матрица силовых постоянных" - симметричная. Для нее по методике madschumacher нужно решить задачу на собственные значения, т.е. диагонализовать. В данном случае это синонимы.

Специально только что глянул... У NumPy имеются готовые функции для решения задачи на собственные значения как для матриц общего вида, так и для вещественных симметричных и комплексных эрмитовых матриц:

Цитата:

eigh
eigenvalues and eigenvectors of a real symmetric or complex Hermitian (conjugate symmetric) array.

eigvalsh
eigenvalues of a real symmetric or complex Hermitian (conjugate symmetric) array.

Вот именно эти функции и стоит использовать!

[B3LYP]

Alex Krylov в сообщении #1647862 писал(а):

Ваша "матрица силовых постоянных" - симметричная. Для нее по методике madschumacher нужно решить задачу на собственные значения, т.е. диагонализовать. В данном случае это синонимы.

https://en.wikipedia.org/wiki/Gaussian_elimination

Это оно?
Тут есть анимированная картинка, мне трудно её распознавать, могу сделать скриншоты и показать тут - так наверно будет понятнее. Если это нужно конечно.

[мат-ламер]

B3LYP в сообщении #1648087 писал(а):

Это оно?

Вы всё же различайте решение системы линейных уравнений и нахождение собственных значений матрицы. По ссылке - решение системы. Если вам нужны частоты колебаний - то вам нужны собственные значения.

[Alex Krylov]

Нет!!!! Вот это: https://ru.wikipedia.org/wiki/%D0%A1%D0%BF%D0%B5%D0%BA%D1%82%D1%80%D0%B0%D0%BB%D1%8C%D0%BD%D0%BE%D0%B5_%D1%80%D0%B0%D0%B7%D0%BB%D0%BE%D0%B6%D0%B5%D0%BD%D0%B8%D0%B5_%D0%BC%D0%B0%D1%82%D1%80%D0%B8%D1%86%D1%8B

И там найдите раздел "4.2 Вещественные симметричные матрицы", относящийся к тому типу матриц, который имеется у вас.

Да... И не надо писать никаких своих программ решения задачи на собственные значения. Всё это давно написано и имеется хоть в NumPy, хоть в Matlab.

[B3LYP]

Я всё-таки ещё раз хотел бы спросить, эта статья - то что мне надо?

https://en.wikipedia.org/wiki/Gaussian_elimination

Тут на первой странице анимированная картинка, скопировал первые 5 кадров:

(Оффтоп)

Перед четвёртым и пятым кадром был ещё обмен третьей строки на четвёртую.
Я пишу на Delphi и мне там трудно находить библиотеки, привык всё делать сам. Я "велосипедист со стажем".

[Евгений Машеров]

Нет.
Точно не то, что Вам надо.
Это решение системы линейных уравнений.
Повторяю - диагонализировать матрицу можно разными способами, но чтобы от диагонализации была польза, диагонализирующие преобразования должны сохранять некоторое свойство. И в разных задачах это свойство различно. При решении линейных уравнений требуется, чтобы преобразованная к диагональному виду матрица, вместе с преобразованным вектором свободных членов, имела бы то же решение, что и исходная задача. В Вашей проблеме требуется, чтобы диагонализированная матрица имела бы те же собственные значения, что и исходная. Соответственно, преобразования иные.
Поскольку матрица у Вас симметричная, задача сильно упрощается (матрица всегда диагонализируется).
Алгоритмы расчёта можно узнать, скажем в Парлетт, "Симметрическая проблема собственных значений", в монографии Уилкинсона (но он рассматривает и несимметричный случай), в справочнике Уилкинсона и Райнша (если Алгол-60 не пугает). Но писать свои имеет смысл либо в порядке упражнения, либо, глубоко изучив существующие, совершенствовать их.

[B3LYP]

Значит задачу можно сформулировать как нахождение собственных значений симметричной матрицы? Собственные значения по английски называются eigenvalues?
Как я уже написал, я "профессиональный велосипедист", и мне кроме прочего было бы интересно самому разобраться и написать наконец правильный алгоритм. Раньше я не сталкивался с такими проблемами, потому что мне требовалось только решать системы линейных уравнений, и для этой задачи я написал свой код с алгоритмом изложенным в оп.
Если есть система линейных уравнений, мы можем прибавлять разные строки друг к другу, так чтобы в итоге остались числа только на диагонали; при этом меняются и числа в столбце справа от матрицы. Если умножить отдельную строку на число, ничего страшного не случится поскольку поменяется и числа в этой строке в матрице, и число справа. С "просто матрицей" же (не очень понимаю что это значит) просто так умножать строку на что-то нельзя. Какими свойствами должен обладать столбец справа, если на диагонали после диагонализации остаются только собственные значения? И что такое вообще собственные значения матрицы?

[dgwuqtj]

Если уж велосипедить профессионально, то надо прочитать теорию хотя бы на уровне университетского курса алгебры, где всё это определяется. Ну и если вам нужны не эффективные алгоритмы, а просто что-то работающее, то можно всё считать как учат студентов: посчитать характеристический многочлен матрицы, найти его корни (а это численное решение полиномиальных уравнений!), потом собственно решение систем линейных уравнений.