1st-order linear ODE with tridiagonal matrix. Efficient solutions?

Question

I have a 1st-rder linear ODE system where the system is characterized by $A$. Given an initial state $x_0$, I want the state at some later time $t$, efficiently.

$A$ happens to be a symmetric tridiagonal matrix where the coefficients are the same in each diagonal line.

$$ A = \begin{pmatrix} a & b & 0 & 0 & ... \\ b & a & b & 0 & ... \\ 0 & b & a & b & ... \\ 0 & 0 & b & a & ... \\ \vdots &&\ddots&\ddots&\ddots\end{pmatrix} $$

I need to solve this for many different $A$'s and also many different $x_0$'s. The general solution to that is of course

$$x_t = e^{A t}x_0$$

One way to do this is to eigen-decompose $A$ into $A = F D F^\top$. Then the solution becomes

$$x_t = (F e^{D t} F^\top) x_0$$

Which is a little better because now the multiplication with $x_0$ is $O(n^2)$, but the eigen-decomposition is $O(n^3)$ still. But this doesn't take advantage of the constant values of $a$ and $b$ or the symmetric tridiagonal nature of $A$. So it seems like I should be able to do better.

Answer 1

I'm not aware of any special algorithms for computing matrix exponentials for matrices with tridiagonal structure, but you probably shouldn't be solving your ODE by directly computing the matrix exponential anyway. Computing matrix exponentials is a tricky business. There's a famous paper by Cleve Moler (inventor of MATLAB) and Charles van Loan that you should check out titled "Nineteen Dubious Ways to Compute the Matrix Exponential". It goes through a variety of schemes people have tried over the years and discusses the problems with each one. That's not to say it can't be done but that it's something that needs to be handled with care.

The good news is that there are plenty of ways to numerically solve your ODE without having to resort to the matrix exponential, and these are all almost certainly more efficient. As an example, consider Euler's method. You want the solution at time $t = t_f$ to $x'(t) = Ax(t)$ with initial condition $x(0) = x_0$. Pick $K + 1$ equally-spaced time-points $0 = t_0 < t_1 < \cdots < t_K = t_f$, and let $h = (t_f - t_0)/K$ be the spacing between them. Let $x_k$ denote the approximation to the solution at time $t_k$. Euler's method is just the iteration $x_{k + 1} = x_k + hAx_k$.

Each step requires one matrix-vector multiply and one addition of two vectors. Since your matrix is tridiagonal, the matrix-vector multiply will be $O(3n)$, and the vector addition will be $O(n)$. With $K + 1$ steps, you're looking at an overall computational complexity of $O(nK)$. As far as accuracy goes, Euler's method converges at a rate of $O(h) = O(K^{-1})$, so if you double the number of points you use, you cut the error in half.

If you want something with better convergence properties than Euler's method, there are many other possibilities. Some things to search for are "Runge-Kutta methods" and "one-step methods." More sophisticated methods require more function-evaluations (matrix-vector multiplies, in this case) at each time step but can attain better accuracy with fewer time steps.

At any rate, there's no need to do an $O(n^3)$ eigenvalue computation if you don't want to!