I'm trying to learn how to exponentiate a matrix in C. In quantum mechanics, the time evolution is given by exp[-iHt/hbar] where H is, say in my case, a 2x2 matrix. I know that mathematically you taylor expand the exp, but this is where I'm running into problems. Any ideas?
Raising matrices to powers generally involves diagonalizing the matrix first, so you would want to start there. (EDIT: Or of course, you do the normal-people thing and use Matlab/Octave instead.)
Last edited by tabstop; 05-12-2010 at 10:22 AM.
It's not raising a matrix to a power -- the matrix is in the exponent.Originally Posted by tabstop
The OP is correct that it involves an expansion. What is giving you difficulty?
Code://try //{ if (a) do { f( b); } while(1); else do { f(!b); } while(1); //}
This is probably the silliest thing you've posted in a long time. What do you think "Taylor expansion" means?
Eh?
What we're talking about is something like exp(H) where H is a matrix. The OP said he understands it involves a Taylor expansion. I'm trying to figure out what part of it he's having trouble with.
EDIT: You said "Raising matrices to powers". That means H^x, not x^H.
Last edited by brewbuck; 05-12-2010 at 12:45 PM.
Code://try //{ if (a) do { f( b); } while(1); else do { f(!b); } while(1); //}
It's either the loop to infinity or the raising of a matrix to different powers. Both problems are solved by diagonalizing the matrix (because (b) that gives you a formula for raising a matrix to a power which leads to (a) being able to sum those powers in the usual Taylor series for e^x in a single variable to avoid having to actually compute an infinite number of terms).
While the matrix may be in the exponent, in practice the calculation is actually done by summing powers of the matrix:
exp(H)=I+H+(1/2!)H^2+(1/3!)H^3+...
So diagonalizing the matrix is a very good idea, because if H=P(^-1)DP, then exp(H)=P(^-1)exp(D)P, and exp(D) is very simple to calculate without having an infinite loop ^_^
Just to clear that argument up...
Last edited by KBriggs; 05-12-2010 at 02:42 PM.
Wait, why would I need to diagonalize my matrices? I wrote out a taylor expansion of e^(x) (this code werks well) but it's e^(H) I want. So I started to write code for my taylor expansion, but I need A^n, where n is an int and A is a matrix. I can't get it to work.
Code:#include <stdio.h> #include <math.h> int main(){ int A[2][2]={{0,1},{1,0}}; int M[2][2]={{0,1},{1,0}}; int I[2][2]={{1,0},{0,1}}; int n; printf("A is the sigma_x matrix.\n"); getchar(); printf("How many times do you want to multiply the matrix?\n"); scanf("%d", &n); for (int l=0; l<n; l++){ for (int i=0; i<2; i++) for (int j=0; j<2; j++){ for (int k=0; k<2; k++) M[i][j]=A[i][k]*A[k][j];}} printf("\nThis is your new matrix.\n"); for ( int i=0; i<2; i++) { for ( int j=0; j<2; j++) printf(" %d", M[i][j]); printf("\n");} getchar(); scanf("%d", &n); if (n==1) {main();} else return 0; }
My real program will be in taylor expanding pauli spin matrices.... so I know what my answer is... i just need to code it.
By diagonalizing the matrix you can get an analytic expression for the exponential of the matrix without ever having to multiply it n times. The exponential of a diagonal matrix is found simply by exponentiating all of the numbers on the diagonal.
You don't have to, it'll just make your program run orders of magnitude faster for large matrices. it will also give you an exact solution (within floating point error, I suppose) rather than the approximation that you would get by taking the first n terms in the taylor expansion.
Last edited by KBriggs; 05-19-2010 at 09:19 AM.
Ok. Sure, a diagonalized matrix will compute faster. What about is the matrix in the exponent is just randomized?
There are algorithms out there to diagonalize a general matrix, designed for this very reason.
The vast majority of random matrices can be diagonalized, the chances of randomly generating a singular matrix are pretty small. However, in the rare case that it happens, the Jordan Form of the matrix also is very easy to exponentiate, and Jordan Form exists always (diagonal form is just a special case of jordan form). It is also studied very extensively for computing applications, you should be able to find an algorithm that does it quite easily.
Last edited by KBriggs; 05-19-2010 at 01:05 PM.
What is the point of a random matrix in quantum mechanics?
Code://try //{ if (a) do { f( b); } while(1); else do { f(!b); } while(1); //}
I'm trying to convert some matlab code into C for a qubit simulation. So, I'm trying to create the taylor series for the time evolution operator. How could I get a function like:
to keep multiplying?Code:printf("How many times do you want to multiply the matrix?\n"); scanf("%d", &n); for (int l=0; l<n; l++){ for (int i=0; i<2; i++) for (int j=0; j<2; j++){ for (int k=0; k<2; k++) M[i][j]=A[i][k]*A[k][j];}}
If i can compute A^1,A^2, shouldn't I be able to get A^3, A^4, A^5,...,A^n?
