An ode to ODEs

Earlier today, Master-o-Forests asked me for help on solving a system of coupled ordinary differential equations (ODEs). I regurgitated the standard method that I'd learned in college, then, as I was driving home, it started to seem sort of unsatisfactory to me, for reasons I won't go into. Anyway, I thought for a bit, and, to keep my brain otherwise occupied while I wait for my interminable simulations to run, I scribbled down an explanation that (to me, at least) seems a little more motivated.

Here's the basic problem statement. Suppose you've got two ODEs which are coupled in the following way:
\[ \frac{dx_1}{dt} = m_1 x_1 + m_2 x_2 \]\[ \frac{dx_2}{dt} = m_3 x_1 + m_4 x_2 \]Where the $m$'s are specified (they're just numbers). How can you find an exact solution for $x_1$ and $x_2$ as a function of $t$?

First, let's identify what the complication is here. The trouble is that the time-derivatives (on the left-hand side) are expressions of both $x_1$ and $x_2$. Life would be a lot easier if this wasn't the case -- if $dx_1/dt$ just had $x_1$ on the right-hand side, then we could separate the variables and integrate, and we'd be done! Fortunately, there's a way to uncouple the system so that we do have two very simple equations. To see this, the first thing to do is gather the $m$'s into a matrix,
\[ \frac{d}{dt}\begin{bmatrix}x_1 \\ x_2 \end{bmatrix} = \begin{bmatrix} m_1 & m_2 \\ m_3 & m_4 \end{bmatrix} \begin{bmatrix} x_1 \\ x_2 \end{bmatrix} \]which is written in more compact matrix notation:
\[ \frac{d\mathbf{x}}{dt} = \mathbf{M} \mathbf{x} \]Notice that we'd have an uncoupled system if M was a diagonal matrix (that is, if $m_2$ and $m_3$ were both equal to 0). Well, M isn't diagonal, but it can be diagonalized by finding its eigenvalues and eigenvectors:
\[ \mathbf{M} = \mathbf{S D S}^{-1} \]where D is diagonal (with its elements being the eigenvalues of M, which we'll call $\lambda_1$ and $\lambda_2$):
\[ \mathbf{D} = \begin{bmatrix} \lambda_1 & 0 \\ 0 & \lambda_2 \end{bmatrix} \]and S has the eigenvectors of M (we'll denote these as a and b) as its columns,
\[ \mathbf{S} = \begin{bmatrix} a_1 & b_1 \\ a_2 & b_2 \end{bmatrix} \]Plugging the diagonalized M into our equation for $d\mathbf{x}/dt$, we find
\[ \frac{d\mathbf{x}}{dt} = \mathbf{SDS}^{-1} \mathbf{x} \]and we can multiply both sides on the left by $\mathbf{S}^{-1}$:
\[ \mathbf{S}^{-1} \frac{d\mathbf{x}}{dt} = \mathbf{DS}^{-1} \mathbf{x} \]To clean up the notation, define a new vector $\mathbf{u} \equiv \mathbf{S}^{-1} \mathbf{x}$:
\[ \frac{d\mathbf{u}}{dt} = \mathbf{D} \mathbf{u} \] Which looks pretty much like our initial system of equations, with one important exception: D is diagonal! Therefore, these equations are uncoupled:
\[ \frac{du_1}{dt} = \lambda_1 u_1 \]\[ \frac{du_2}{dt} = \lambda_2 u_2 \]These are separable, so it's easy to solve for $u_1$ and $u_2$:
\[ u_1 = C_1 \mathrm{e}^{\lambda_1 t} \]\[ u_2 = C_2 \mathrm{e}^{\lambda_2 t} \]where $C_1$ and $C_2$ are as-yet-unknown constants of integration (they'll be specified by the initial conditions). Now that we know u, it's simple to also calculate x, since, by definition:
\[ \mathbf{x} = \mathbf{S u} \]So, that's our answer. It looks a little different than the usual result, but, multiplying it out explicitly, we find the general solution is written
\[ \begin{bmatrix}x_1 \\ x_2 \end{bmatrix} = C_1 \begin{bmatrix} a_1 \\ a_2 \end{bmatrix} \mathrm{e}^{\lambda_1 t} + C_2 \begin{bmatrix}b_1 \\ b_2 \end{bmatrix} \mathrm{e}^{\lambda_2 t} \]which is the expected result.

Anyway, your mileage may vary, but I thought that this makes it clearer where the exponentials in the solution come from, and why we care about diagonalizing M in the first place.