Check the stability of differential equation using fsolve

user_7393694 0
Viewed 554

I want to find the equilibrium point of a differential equation and check whether the equilibrium point is stable.

Here's a minimal working example

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)

To find whether the solution sol is stable all the eigenvalues of the Jacobian must have negative real part. However the Jacobian is not saved in the infodict, rather the QR decomposition is saved in the infodict.

How do I get the Jacoian back from the QR decomposition of fsolve?

All I could do was something like

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

Where ind are the diagonal entries of r, however I doubt this is the best possible way.

scipy matrix python
2 Answers

QR decomposition is cheap: it takes a fixed number, around n**3, of operations, compared to finding eigenvalues which is an iterative process. Indeed, one of eigenvalue-finding algorithms involves the iteration of QR decomposition. So, knowing the QR factors does not really bring you much closer to having eigenvalues. And the cost of reconstructing the matrix by multiplication (also less than n**3 operations) is negligible compared to the cost of finding the eigenvalues.

The conclusion is that reconstructing the Jacobian by multiplication is the way to go here. What you are doing (finding the eigenvalues of Q factor alone?) is not correct. First, recover the R matrix from its given flat form using np.triu_indices; then multiply Q by R; then find the eigenvalues.

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

Need also to take transpose of infodict["fjac"] to get Q. Also I had a bug if I labeled the matrix as "r". I would also recommend directly checking the Jacobian So :

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