Let's consider vector V = {X;Y}, therefore, we can combine the two equations into a single one as follows: {V'} = [A]{V}, where matrix [A] is now 2n*2n and is defined as: [A] = [A1, A2; -A2^T, -A1^T].

The exact solution to the Initial value problem is given as: {V} = exp(At)*V_0, where V_0 is the initial condition, i.e. V_0 = V(t = 0).

The computation of exp(At) is the major issue to solve the problem, by definition, we have:

exp(At) = I + At + (At)^2/(2!) + (At)^3/(3!) + ...+ (At)^n/(n!) +...

The trace of [A] = trace of [A1] + trace of [-A1^T] = 0. Thus the eigenvalues of [A] cannot be all of the same sign! (this is in case they were all real).

Nevertheless, the eigenvalues of [A] need not to be symmetric about the origin as is the case for A1 and A2. You can easily find a counter-example using simply 2*2 matrices.

In general, there is no special property that the eigenvalues of [A] satisfy.

As such, If you want to solve the system exactly, you either need to solve for exp(At) using the Taylor series as given before or you have to perform the eigenvalue decomposition of [A] using numerical technique (e.g. Krylov subspace method). and then you can say

exp(At) = [L]*exp(St)*[L]^-1 , where S is the diagonal matrix whose entries are the eigenvalues of [A], and [L] is the matrix combining all the column eigenvectors.

If you are interested in numerical methods to solve your system, There are plenty of methods which compute X(t) and Y(t) numerically.

By the way what is the size of your system (n = ?)

Errata to my previous answer.

My apologies, the eigenvalues of [A] are going to be symmetric symmetric with respect to the origin as Peter Breuer suggested. but indeed they can be complex.

the characteristic equation of A1 and A2 contains only even powers and so will be the one of [A].

Dear Peter and Hady:

apparently you did not notice that in the first equation matrices X and Y are postmultiplied by the coefficient matrices while in the second one they are premultiplied. The respective matrices do not commute so I cannot write the DE in the form you suggested.

I used the Kronecker product to vectorize the equations. I'll modify the question to include the details.

It seems X and Y are not vectors, they are matrices !!!

I answered the questions thinking that X and Y are column vectors... Ooops

Please clarify that so we know how to proceed.

On top, can you please identify the physical problem that these two coupled ODEs simulate?

Regards

Hady, surely, all these are [nxn] matrices. That's why I wrote *matrix* DE. Sorry if it there was a misunderstanding.

The problem relates to stability of time-delay DEs. It's an auxiliary system which has lto be solved to get a Lyapunov-Krasovsky functional of a special form. See, for instance, http://link.springer.com/chapter/10.1007/978-0-8176-8367-2_2

Thanks Dmitry for your clarification.

Well in this case, You have p = 2*n^2 total unknowns, So the vectorized form of the DE is

{W}' = [A]{W} where A is p by p matrix.

There is nothing that stops you from solving this system. I am not sure if the eigenvalues of A have some special property! I really cannot see at a glance. Are you looking for some special property you want to investigate?

Is there a good estimate for n = ? If n is large, then p is even larger. So you are somehow obliged to adopt a numerical scheme to solve for X and Y. (Runge-Kutta) is an option.

It would become really expensive to solve for exp(At) to solve for X and Y.

Hope that helps

Peter, sorry, I didn't get your point.

Consider again

(*)

X'= X*A1 + Y*A2

Y'= -A2^T*X - A1^T*Y

Changing coordinates to (X,Y)=(-Y,X) we get

X'= - A1^T*X + A2^T*Y

Y'= - X*A2 + Y*A1

which is obviously not equal to (*).

The only transformation which transforms the system to itself is (X, Y)-> (-Y^T, -X^T), but it seems to be of little help. The only conclusion which can be derived from it is that the system's eigenvalues are symmetric w.r.t. 0.

Peter, it is remarkable that the result you described leads to the same conclusion about the eigenvalues of the system while being apparently not applicable to the matrix case. I mean that when describing the system in the vector form (with Kronecker product etc) I cannot get the same result.

Hmm... there must be some relation between these two cases, I guess...

Well, it turns out that if we introduce new variables

V=X^T-Y and W=X^T+Y,

the system can be rewritten in the vector form as

q'=Kq,

where q=(vec(V)^T, vec(W)^T)^T,

K=[0, kron(I,A1)+kron(I,A2)*M; kron(I,A1)-kron(I,A2)*M, 0], and M is a permutation matrix. 

This form seems to be simpler than the initial one. It is also interesting that K and L (as defined in the question) have the same eigenvalues. However it is not clear to me why it is so.