Проверьте стабильность дифференциального уравнения с помощью fsolve

Question

Я хочу найти точку равновесия дифференциального уравнения и проверить, устойчива ли точка равновесия.

Вот минимальный рабочий пример

import numpy as np
from scipy.optimize import fsolve

dim = 2
A = np.random.uniform(size = (dim,dim))
sol, infodict, ier, mesg = fsolve(lambda x: 1-np.dot(A,x),
        np.ones(dim), full_output = True)

Чтобы определить, устойчиво ли решение sol, все собственные значения якобиана должны иметь отрицательную вещественную часть. Однако якобиан не сохраняется в infodict, а QR-разложение сохраняется в infodict.

Как мне вернуть якояна из разложения QR fsolve?

Все, что я мог сделать, было что-то вроде

eigenvalues = np.linalg.eigvals(infodict["fjac"])*infodict["r"][ind]

Где ind - это диагональные элементы r, однако я сомневаюсь, что это наилучший из возможных способов.

Answer 1

Нужно также взять транспонирование infodict["fjac"], чтобы получить Q. Также у меня была ошибка, если я пометил матрицу как "r".Я бы также рекомендовал непосредственно проверить Якобиана. Так:

import numpy as np

matrixr=np.zeros((dim,dim))                    # dim= number of linear equations
matrixr[np.triu_indices(dim)]=infodict["r"]    # to unpack "r" into matrix
Jacobian=(infodict["fjac"].T).dot(matrixr)     # needs .T to get Q
eigenvalues=np.linalg.eigvals(Jacobian)

Answer 2

QR-разложение дешево: требуется фиксированное число, около n**3 операций, по сравнению с поиском собственных значений, что является итеративным процессом. Действительно, один из алгоритмов нахождения собственных значений включает в себя итерацию QR-декомпозиции. Таким образом, знание QR-факторов не слишком приближает вас к собственным значениям. И стоимость восстановления матрицы умножением (также менее чем n**3 операций) незначительна по сравнению со стоимостью нахождения собственных значений.

Вывод состоит в том, что реконструкция якобиана путем умножения - это путь, по которому можно пойти сюда. То, что вы делаете (находя только собственные значения Q-фактора?), Неверно. Сначала восстановите матрицу R из заданной плоской формы, используя np.triu_indices; затем умножьте Q на R; затем найдите собственные значения.

r = np.zeros((dim, dim))
r[np.triu_indices(dim)] = infodict["r"]
eigenvalues = np.linalg.eigvals(infodict["fjac"].dot(r))