使用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分解中获得Jacoian fsolve？

我能做的就像是

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

ind对角线条目在哪里r,但我怀疑这是最好的方式.

Answer 1

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))