用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中。

如何从fsolve的QR分解中恢复雅各布

我所能做的就是

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

其中ind是r的对角线条目，但是我怀疑这是最好的方法。

Answer 1

QR分解的成本很低:与寻找特征值是一个迭代过程相比，它需要固定数量的操作，大约是n**3。实际上，特征值查找算法之一涉及QR分解的迭代。因此，知道QR因子并不能真正让你更接近特征值。与寻找特征值的成本相比，通过乘法(也小于n**3运算)重建矩阵的成本可以忽略不计。

结论是通过乘法重建雅可比矩阵是一种可行的方法。你在做什么(单独寻找Q因子的特征值？)是不正确的。首先，使用np.triu_indices从给定的平面形式恢复R矩阵；然后将Q乘以R；然后找到特征值。

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

Answer 2

我还需要转置infodict["fjac"]来获得Q。另外，如果我将矩阵标记为"r“，我也会有一个bug。我还建议直接检查雅可比矩阵，这样：

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)