# Exponentiate the master equation matrix

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?

## 1 Answer

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.

• Please can you have a look at the edit I have made to the question. Thanks Nov 6, 2016 at 17:23
• To understand your code read what I wrote about calculating the exponential of a diagonalisable matrix. I don't know what GammaTot is supposed to be, but from the way the code is written it looks like they are a calculating $\exp(D)$ and then multiplying by $B^{-1}$ and $B$ to recover $\exp(A)$. Nov 6, 2016 at 17:29
• Although they are setting the diagonal entries to something that definitely does not evaluate to $\exp(D)$. Nov 6, 2016 at 17:32
• I understood that. What is don't understand is when are they calculating $exp(D)$ in the code? It is almost clear that if they calculated $exp(D)$ and then pre and post multiply by eigenvector matrix we can recover $exp(A)$, but where are they doing so? Nov 6, 2016 at 17:34
• Exactly, my point. What are they setting it to? Nov 6, 2016 at 17:35