使用fsolve检查微分方程的稳定性

发布时间：2020-12-16 22:28:19 所属栏目：Python 来源：网络整理

导读：我想找到微分方程的平衡点,并检查平衡点是否稳定. 这是一个最小的工作示例 import numpy as npfrom scipy.optimize import fsolvedim = 2A = np.random.uniform(size = (dim,dim))sol,infodict,ier,mesg = fsolve(lambda x: 1-np.dot(A,x),np.ones(dim),full

我想找到微分方程的平衡点,并检查平衡点是否稳定.

这是一个最小的工作示例

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是否稳定,雅可比行列的所有特征值必须具有负实部.然而,Jacobian没有保存在infodict中,而是QR分解被保存在infodict中.

如何从fsolve的QR分解中获得Jacoian？

我能做的就像是

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

ind是r的对角线条目,但我怀疑这是最好的方法.

最佳答案

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

（编辑：李大同）

【声明】本站内容均来自网络，其相关言论仅代表作者个人观点，不代表本站立场。若无意侵犯到您的权利，请及时与联系站长删除相关内容!