Solving First-order simultaneous ODEs

Solving a system of ordinary differential equations (ODEs) has always been a cakewalk. Given a pen and paper, it takes almost no time to solve a simple system such as this:

$\frac{dx}{dt} + ax + by = 0$

$\frac{dy}{dt} + cx + dy = 0$

One can introduce the differential operator $D$ and solve as if it were a set of algebraic equations. But when the system of equations gets bigger, things become messy. They can be solved only by an in-built solver routine in a high-level language like MATLAB. However, feeding such a huge system can still be troublesome. The purpose of this post is to expose a cute little trick in mathematics and exploit it to our advantage.

To explain the trick, we consider the above set of differential equations. We rewrite them in matrix form as:

$\begin{bmatrix} x'\\ y' \end{bmatrix} + \begin{bmatrix} a &b\\ c &d \end{bmatrix} \begin{bmatrix} x\\ y \end{bmatrix} = \begin{bmatrix} 0\\ 0 \end{bmatrix}$

or simply,

$X' + AX = 0$

where, $X$ is the unknown vector whose elements are $x$ and $y$.
Now, here is the trick:

Split the matrix $A$ into 3 matrices as $P^{-1} \Lambda P$, where $\Lambda$ is a diagonal matrix, with the eigenvalues of $A$ as the diagonal elements. $P$ is a matrix, which has the eigenvectors of $A$ as its columns.

This is called as Eigenvalue decomposition. This is actually an one-liner in MATLAB. Now, the matrix equation becomes:

$X' + P^{-1}\Lambda P X = 0$

Pre-multiply by $P$ to get

$PX' + \Lambda PX = 0$

Now, take $PX$ to be $U$, and owing to the fact that $\Lambda$ is diagonal, we get a system of the form:

$\frac{du_1}{dt} + \lambda_1 u_1 = 0$

$\frac{du_2}{dt} + \lambda_2 u_2 = 0$

Now that we have separated the dependent variables, solving this should be a no-brainer. Finally, apply the reverse transformation $P^{-1}U$ to get the solution.

This works even when the right-hand side is non-zero. We just have to pre-multiply it by $P^{-1}$. Also, remember to transform the boundary conditions as well.

A variation to this problem may be something like this:

$AX' + BX = c$

where $A$ and $B$ are square matrices. In such a case, first pre-multiply by $A^{-1}$ and convert it to the 'standard form' discussed above. The rest should be easy.

Always remember: The trick behind the trick is knowing when to use the trick. Else there is no point in knowing the trick.