## Context

I'm trying to numerically solve the following differential equation: $\frac{\mathrm{d} u}{\mathrm{d} t} = -Au + f$, where $u$ and $f$ are vectors, and $A$ is an $N \times N$ matrix, with $N > 2^{10}$ typically. For a small step $t + \Delta t$, I have the following approximate solution:

$u(t + \Delta t) = \exp \left( - A \Delta t \right) \left[ A^{-1} \left(I - \exp \left( - A \Delta t \right) \right) f + u(t) \right]$

where $A^{-1}$ is the inverse of $A$ and $I$ is the identity matrix. This solution works well for me, except that computing $\exp \left(-A \Delta t\right)$ is incredibly expensive. I want to speed things up by exploiting the structure of $A$, which is constructed as:

$A = D_k F D_{\alpha} F^{-1} D_k$

with $D_k$ being a diagonal matrix with diagonal entries $(0, 1, \ldots , \frac{N}{2} - 1, -\frac{N}{2}, -\frac{N}{2}+1, \ldots -1)$, $D_{\alpha}$ being a diagonal matrix with no particular structure of its diagonal entries. $F$ represents the discrete Fourier transform matrix (i.e. $\mathrm{DFT}(x) = \tilde{x} = Fx$), and $F^{-1}$ represents the inverse discrete Fourier transform matrix (i.e. $x = F^{-1} \tilde{x}$).

## My question

All of the components of $A$ have nice properties, like having known eigenvalues/vectors, and being trivially invertible, etc. Can I use these sub-properties somehow to expedite the computation of $\exp \left(-A \Delta t\right)$ (or better yet, compute $\exp \left(-A \Delta t \right) b$, where $b$ is a vector)? Or are there other ways of computing the exponential of a matrix?

## Something to think about

In my search for a solution, I've come across this paper, which offers an iterative solution approach for $\exp \left(A \right) B$ as:

$B_{i+1} = r_m \left(s^{-1} A \right) B_i, \;\;\; i=0:s-1, \;\;\; B_0 = B$

with $r_m$ being the Taylor expansion of $\exp \left( s^{-1} A \right)$ given by:

$r_m = \sum_{j=0}^{m} \frac{\left(s^{-1} A \right)^j}{j!}$

If the eigenvectors/values of $A$ are known a-priori based on the eigenvalues of $D_k$, $D_{\alpha}$, and $F$, then the computation of $r_m$ is greatly simplified. Furthermore, $A$ exhibits a block structure, and is relatively sparse in practice. I'm willing to sacrifice accuracy in favour of performance by e.g. discarding near-zero values.