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.

  • $\begingroup$ Yes, I have tried explicit schemes, but the time step requirements (stability conditions) are too restrictive. This implicit approach already performs better than an explicit scheme. $\endgroup$ – MPA Oct 5 '18 at 13:14
  • 1
    $\begingroup$ The general solution is $u(t) = A^{-1} f + \exp(-At) v$ where $v = u(0)-A^{-1} f$. You might check out Moler and van Loan, Nineteen Dubious Ways to Compute the Exponential of a Matrix, Twenty-Five Years Later. $\endgroup$ – Robert Israel Oct 5 '18 at 16:02

Your Answer

By clicking "Post Your Answer", you acknowledge that you have read our updated terms of service, privacy policy and cookie policy, and that your continued use of the website is subject to these policies.

Browse other questions tagged or ask your own question.