Exponentiate the master equation matrix

Question

I am working on a problem where I need to solve the master equation for a system and I need to find the correlation function of the form $\langle x(0)x(t)\rangle$ for any quantity (e.g. position etc.) $x$ that takes on discrete values $x = x_1, x_2,\cdots, x_N$ when the system is, respectively, in states $i = 1, 2, \cdots, N$. A detailed explanation and notation can be found in this pdf -http://makarov.cm.utexas.edu/resources/Lecture-notes,-tutorials-etc./master_equations.pdf

They derive the correlator function to be $x^Te^{-Kt}\bar{x}$ where K is the probability transition (master) matrix.How can I solve this numerically? Can anyone provide me with a matlab code for the same?

Edit1-From a source I got a matlab code for a system which does the same. Assume Wmat is the master equation matrix.

sizeWmat=size(Wmat,1);
[Umat,Dmat]=eig(Wmat,'nobalance');
for m=1:sizeWmat
currLambda=Dmat(m,m);
if abs(currLambda)>1e-9;
Dmat(m,m)=1/(currLambda*GammaTot);
else
Dmat(m,m)=0;
end
end
invUmat=inv(Umat);

And then they used Umat*Dmat*invUmat in place of $exp(Wmat*t)$ Can somebody explain this?

Answer 1

For an $n\times n$ matrix $A$ the matrix exponential is defined as $$\exp(A)=\sum_{n=0}^\infty \frac{A^n}{n!}$$ If you want to calculate this numerically the obvious thing to do is do the sum and break off after some finite amount of steps. Bounding the error works in the same way as for the regular exponential map.

This is problematic because each time you calculate a new term in the sum you need to multiply a $n\times n$ matrix. If $n$ is big this is bad numerically.

What you could also do, if $A$ is diagonalisable, is to diagonalise $A$ so that $A=B^{-1}DB$ with $D$ a diagonal matrix: $$D=\begin{pmatrix}d_1& 0&\cdots&0\\ 0&d_2& \cdots &0\\ \vdots&\vdots&\ddots&\vdots\\ 0&0&\dots&d_n\end{pmatrix}$$ Diagonalising is also bad numerically, but you only need to do it once as opposed to $N$ times from above. So $A^n=(B^{-1}DB)(B^{-1}DB)\cdot ...\cdot (B^{-1}DB)=B^{-1}D^n B$. And $$\exp(A)=\sum_{n=0}^\infty \frac{(B^{-1}AB)^n}{n!}=B^{-1}\sum_{n=0}^\infty \frac{D^n}{n!}B=B^{-1}\exp(D)B$$ but $\exp(D)$ is easy to calculate: $$\exp(D)=\begin{pmatrix}e^{d_1}& 0&\cdots&0\\ 0&e^{d_2}& \cdots &0\\ \vdots&\vdots&\ddots&\vdots\\ 0&0&\dots&e^{d_n}\end{pmatrix}$$ It is not always possible to diagonalise a matrix, the best you can do is to find invertible $B$ so that $A=B^{-1}(D+N)B$ with $DN=ND$ and where $D$ is diagonal and $N$ is nilpotent, meaning there exists a $k$ with $N^k=0$. From $DN=ND$ you get: $$\exp(D+N)=\exp(D)\exp(N)$$ And from $N^k=0$ you know that to find $\exp(N)$ you only need to calculate a finite number of terms.

One example of such a procedure is decomposition into the Jordan normal form of a matrix.