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.