This is an old version, view current version.

13.3 Solving a System of Linear ODEs using a Matrix Exponential

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$$ id:ode.linODEs

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$$ id:ode.linOEs.sln

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!} + ...$$ id:ode.matrix_exp.def

We can apply this technique to the simple harmonic oscillator example, by setting

$$y = \left[\begin{array}{c} y_1 \\ y_2 \\ \end{array}\right] \ \ \ \ \ A = \left[\begin{array}{cc} 0 & 1 \\ -1 & -\theta \\ \end{array}\right]$$ id:ode.sho_matrix

The Stan model to simulate noisy observations using a matrix exponential function is given below. Because we are performing matrix operations, we declare y0 and y_hat as vectors, instead of using arrays, as in the previous example code.

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 ODEs.

data {
  int<lower=1> T;
  vector[2] y0;
  real ts[T];
  real theta[1];
}
model {
}
generated quantities {
  vector[2] y_hat[T];
  matrix[2, 2] A = [[ 0,  1],
                    [-1, -theta[1]]]
  for (t in 1:T)
    y_hat[t] = matrix_exp((t - 1) * A) * y0;
  // add measurement error
  for (t in 1:T) {
    y_hat[t, 1] += normal_rng(0, 0.1);
    y_hat[t, 2] += normal_rng(0, 0.1);
  }
}

id:sho-sim-me.figure

This Stan program simulates noisy measurements from a simple harmonic oscillator. The system of linear differential equations is coded as a matrix. The system parameters theta and initial state y0 are read in as data along observation times ts. The generated quantities block is used to solve the ODE for the specified times and then add random measurement error, producing observations y_hat. Because the ODEs are linear, we can use the matrix_exp function to solve the system.