This is an old version, view current version.

## 13.8 Solving a system of linear ODEs using a matrix exponential

Linear systems of ODEs can be solved using a matrix exponential. This can be considerably faster than using one of the ODE solvers.

The solution to $$\frac{d}{dt} y = ay$$ is $$y = y_0e^{at}$$, where the constant $$y_0$$ is determined by boundary conditions. We can extend this solution to the vector case: $\frac{d}{dt}y = A \, y$ where $$y$$ is now a vector of length $$n$$ and $$A$$ is an $$n$$ by $$n$$ matrix. The solution is then given by: $y = e^{tA} \, y_0$ where the matrix exponential is formally defined by the convergent power series: $e^{tA} = \sum_{n=0}^{\infty} \dfrac{tA^n}{n!} = I + tA + \frac{t^2A^2}{2!} + \dotsb$

We can apply this technique to the simple harmonic oscillator example, by setting $y = \begin{bmatrix} y_1 \\ y_2 \end{bmatrix} \qquad A = \begin{bmatrix} 0 & 1 \\ -1 & -\theta \end{bmatrix}$

The Stan model to simulate noisy observations using a matrix exponential function is given below.

In general, computing a matrix exponential will be more efficient than using a numerical solver. We can however only apply this technique to systems of linear ODEs.

data {
int<lower=1> T;
vector[2] y0;
array[T] real ts;
array[1] real theta;
}
model {
}
generated quantities {
array[T] vector[2] y_sim;
matrix[2, 2] A = [[ 0,  1],
[-1, -theta[1]]]
for (t in 1:T) {
y_sim[t] = matrix_exp((t - 1) * A) * y0;
}