Solution of Diffusion equation with exponential integrator

Question

It is well known that stiff differential equations can be well tackled with the exponential integrator. I was trying to solve a diffusion equation with an exponential integrator and got stuck at matrix exponentiation.

The equation that I am trying to solve is the following, $$\frac{\partial f}{\partial t}=\frac{\partial}{\partial x}\left(D(x)\frac{\partial f}{\partial x}\right)-\frac{f}{E(x)}$$. Using central differencing I got,

$$\frac{dF}{dt}=\mathcal{A}F$$ where $F$ is a vector $(f_{0},f_{1},...,f_{n})$ consisting the value of $f$ at $n+1$ discrete points. $\mathcal{A}$ is an usual tridagnal matrix arising from central differencing of the diffusion term. With exponential integrator, the evolution of $F$ is given by, $$F^{n+1}=e^{\mathcal{A}dt}F^{n}$$ where $dt$ is the timestepping.

I am stuck while exponentiation the tridiagonal matrix $\mathcal{A}$. How to exponentiate the matrix in C ? Is there any algorithm which can do it or do I need to approximate the exponential in the following way, $$e^{\mathcal{A}dt}\approx\sum_{i=0}^{k}\frac{(\mathcal{A}dt)^{i}}{i!}$$ approximating up to $k$.

Answer 1

There are several options:

1. Find a Similarity transform $S$ that diagonalizes $A$, and then $e^{A dt}=S^{-1} e^{Ddt}S$ with $D$ diagonal.
2. Realize that in any case you are dealing with an discrete time approximation ($A$ changes with time) and you can use a rational (Cayley transform) approximation $e^{dtA}\sim {(I+dt A)\over (I-dt A)}$ which is higher order, and implicitly solvable.