How to decouple a first order ODE system by eliminating coupled functions

Question

Let's see an example of a first order ODE system : $$\begin{align*} y_1'&=a_{11}\cdot y_1+a_{12}\cdot y_2 \qquad(1)\\ y_2'&=a_{21}\cdot y_1+a_{22}\cdot y_2 \qquad(2) \end{align*}$$ where $y_1$ and $y_2$ are functions we are to solve, and $a_{11},\; a_{12},\; a_{21},\; a_{22}$ are constants. I know 2 methods to solve this, one is to eliminate $y_1$ or $y_2$ by replacing, and then we will get 2 decoupled second order ODEs

$$ y_1''-(a_{11}+a_{22})y_1'+(a_{11}a_{22}-a_{12}a_{21})y_1=0 $$ the other ODE is similar.

The other method: Try to calculate the eigenvalues and eigenvectors of the coefficient matrix $\{\{a_{11},a_{12}\},\{a_{21},a_{22}\}\} $, and then get the general solutions for system of $y_1$ and $y_2$ simultaneously. Though the second method is more popular in the textbooks, I like the first one more.

So my question is: Which function in Mathematica can decouple the ODE system by the first method automatically?

Answer 1

I believe going to the eigen-space precisely decouples the equations, so the two methods are strictly equivalent.

Lets define your matrix

 mat = {{a11,a12},{a21,a22}};

Now find the eigenvector matrix P (the so called transformation matrix)

 {\[Lambda]s, P} = Eigensystem[mat]; P = Transpose[P];

Let us check that going to the eigen-space makes the Matrix diagonal

 pmat= Inverse[P].mat.P // FullSimplify

 (* 
 1/2 (a11+a22-Sqrt[(a11-a22)^2+4 a12 a21])  0
 0  1/2 (a11+a22+Sqrt[(a11-a22)^2+4 a12 a21])
*)

Check that the diagonal elements are the eigenvalues

pmat[[1, 1]]/\[Lambda]s[[1]] // Simplify
pmat[[2, 2]]/\[Lambda]s[[2]] // Simplify

So if we define as new variables {z1[t],z2[t]} so that

Thread[{z1[t],z2[t]}=pmat.{y1[t],y2[t]} ]

The system in z is simply

z1'[t] ==  \[Lambda]s[[1]] z1[t]
z2'[t] ==  \[Lambda]s[[1]] z2[t]

which is decoupled.

Answer 2

Not really sure why you would want to prescribe the method Mathematica uses to solve the coupled ODEs. If you don't Mathematica solves them just fine:

eq1 = y1'[t] == a11*y1[t] + a12*y2[t];
eq2 = y2'[t] == a21*y1[t] + a22*y2[t];

FullSimplify[DSolve[{eq1, eq2}, {y1[t], y2[t]}, t]]

But if you want to go the eliminate & replace route you could do the following:

eq1D = Assuming[{a11 != 0, a12 != 0, a21 != 0, a22 != 0}, 
         FullSimplify[eq2 /. Solve[eq1, y2[t]] /. D[Solve[eq1, y2[t]], t]]][[1, 1]]

eq2D = Assuming[{a11 != 0, a12 != 0, a21 != 0, a22 != 0}, 
         FullSimplify[eq1 /. Solve[eq2, y1[t]] /. D[Solve[eq2, y1[t]], t]]][[1, 1]]

DSolve[eq1D, y1[t], t]

DSolve[eq2D, y2[t], t]