Stability analysis of a time dependent Lyapunov equation

I want solution of coupled linear equations for which I am using time dependnet Lyapunov equation:

dP/dt=A*P + P*A'+ D, 

where A (A' is transpose) is a system matrix, P is covariance matrix to be found and D is diffusion matrix. I solve it staright-forward as coupled equations extracted from a matrix in terms of a col. vector but the solution blows up i.e. "instabilities", even though the real eignevalues of matrix A are negative. Whie seraching literatures for this, an options is given, using vectorization: A*P+P*A'=(A⊗I+I⊗A)*P,. using this, giving me an error of incompatibel sizes. I am not sure if this vectorization is enough to make sure the stability of the system but I dont know other options, please suggets if any. For that vectorization, I am using a function:

% showing just a stucture of the function file

function dpdt=(t,p0)

dpdt=zero(m,1)% m is number of rows

[m,n]=size(A);

P1=reshape(P0, m,1);

D1=reshape(D,m,1);% D is some matrix

b=eye(m);

a=real(eig(A))

if all(a<0)

dprdt=(tensorprod(A.',b)+tensorprod(b,A.')).*P1 + D1;% error appears here;

dpdt(1)=dprdt(1);

.

.

.

.

dpdt(m)=dprdt(m)

end

end

is the above right way of doing the vectorization? Any other suggestions for the stability of this time depedent Lyapunov equation, I shall be grateful. Thanks.

Answer 1

[Christine Tobler]

The tensorprod function doesn't do what you are looking for here, that would be kron (short for Kronecker product).

It also seems like P being a column vector of size mx1 would have the wrong size, are you looking for it to be a m x m matrix that you reshape to be m*m x 1?

In terms of performance, you should expect A*P + P*A' (and then possibly reshaping the result of this) to be much more efficient than forming the (A' kron I + I kron A') matrix explicitly.

[Anil]

Dear @Christine Tobler, Thanks for that. The vectorization transformation says that make that P into a vector and seems m*1. I am not sure about m*m*1. But using "kron" works for me. Thorugh the instability is still there.

As I mentioed, I solve the euqtaion dP/dt= A*P+P*A.' +D, and solution is not stable or blows up, even though the all eignevalues of matrix A are negative. It seems eigenvalue condition is not enough for the transient case as it Looks like off diagonal ellemnets are very sensitive for small perturbations.

[Christine Tobler]

The equation A*B + P*A.' implies that P is an m x m matrix. That means that if it is reshaped, it must still have m^2 elements, so mx1 would not be a possible result of that reshaping.

You might start out to check your code is doing the right thing by using diagonal matrices A - then eigenvalues being negative should be enough. For non-symmetric A, it's possible for the values to grow before starting to stabilize.

[Anil]

Dear @Christine Tobler, Yes, that's true, reshaped P is m*m x1. I missunderstood that. Thanks for the later part, I am checking it.