Thread: Exponentiating a Matrix

Code:

ii=i;

%Pauli ops
sx=[0,1;1,0];
sz=[0,0;0,1];

%Number of steps, define time slice
dt=T/N;

%Physical parameters
G=1/T1;

trajectory=[];

%Set up the part of the Hamiltonian that doesn't have detuning
drive=0.03;
H=2*pi*(drive*sx-ii*(G/2)*sz);
psi=[1;0];

for step=1:N
     Heff=H+noise(step)*sz;
     U=expm(-ii*dt*Heff);
     psi=U*psi;
     
     %Quantum Jump Step
     p=G*dt*abs(psi(2))^2
     jump=rand();
     if jump<p
          psi= [1,0];
     end
     %Record current state in the history
     trajectory(step,:)=psi;
end

The trajectory function is another small averaging function.

I can't speak for the internals of Matlab, but I would bet quite a bit that internally they put the matrix in Jordan form and use that for exponentiation. There are two cases for a two-by-two matrix (diagonalizable or non-diagonalizable) and each of them have a very straightforward exp formula.

Originally Posted by tabstop

I can't speak for the internals of Matlab, but I would bet quite a bit that internally they put the matrix in Jordan form and use that for exponentiation. There are two cases for a two-by-two matrix (diagonalizable or non-diagonalizable) and each of them have a very straightforward exp formula.

I'm sure there's a slick way to do it. I'm not there yet in my skills. Currently, I'm just working on developing a taylor series for the expm[A]. Once that works, I'll try to develop some structure for complex numbers...

I've ran into another wall. Ouch... I've created A^n for the taylor expansion, but the answer is in a func which is not in main. How do I take the result of A^n and map it to something usable in main?
The code below shows A^n, but the parts in red are what I'm trying to map to something in main. At this point, it just maps result in the func to mediate, which is a global variable. Anyways, it's printing a 2x2 zero's matrix instead of result.

Code:

#include <stdio.h>
#include <stdlib.h>
#include <math.h>

int mediate[2][2];

void multmtrx(int A[2][2], int *B, int *result)
{
	for (int i=0; i<2; i++){
		for (int j=0; j<2; j++)
		for (int k=0; k<2; k++)
			result[2*i+j]+=A[i][k]*B[2*k+j];}
}

void printmtrx(int *result)
{
	printf("\nThis is your new matrix.\n");
	for (int i=0; i<2; i++)
	{	for (int j=0; j<2; j++)
		printf(" %d", result[2*i+j]);
		printf("\n");
	}
}

void expm(int A[2][2], int n){
	int *temp;
	
	int *result;
	
	do {
		temp = (int*)malloc(sizeof(int) * 4);
		temp[0] = 0;
		temp[1] = 1;
		temp[2] = 1;
		temp[3] = 0; 
		
		for (int i = 1; i < n; i++) {
		result = (int*)malloc(sizeof(int) * 4);
		result[0] = 0;
		result[1] = 0;
		result[2] = 0;
		result[3] = 0;

		multmtrx(A, temp, result);
		int *tempPointer;
		tempPointer = result;
		free(temp);
		temp = tempPointer;
		}
		[COLOR="Red"]
                /* Here I am trying to save result from a function so I can use it in main for something               else.
		/* How do I map the result to something usable?

		for (int i=0; i<2; i++)
			for (int j=0; j<2; j++)
				mediate[i][j]=result[2*i+j];
		
		*/
		[/COLOR="Red"]
                printmtrx(temp);
			
	} 
	while (n==1);
	}

int main()
{
	int sigma_x[2][2]={{0,1},{1,0}};
	int n;	
	start:
	printf("\nHow many times do you want to multiply the matrix?\n");
		printf("This returns A^n\n");
		scanf("%d", &n);
		printf("\nThe matrix will be multiplied %d times.", n);
		expm(sigma_x, n);
		/*
		printmtrx(&mediate[2][2]);
		*/
		[/COLOR="Red"]
                printf("\nAgain?");
		scanf("%d", &n);
		if (n==1) {goto start;}	 
		else
	return 0;
}

Does anyone have some advice?

You can't return an array from a function. If you need that result, you must pass a matrix into your function to hold the result.

I am passing sigma_x[2][2] into expm()... i.e. expm(sigma_x).

then it's job after that is to compute sigma_x^n. then it prints out the result then returns back to main.....
so

Code:

for (int i = 1; i < n; i++) {
		result = (int*)malloc(sizeof(int) * 4);
		for (int i=0; i<4; i++)
		result[i] = 0;
		
		multmtrx(A, temp, result);
		int *tempPointer;
		tempPointer = result;
		 for (int i=0; i<2; i++) 
			for (int j=0; j<2; j++)
		free(temp);
		temp = tempPointer;
		}

do i need to map result back to the input argument A[2][2]?
Thanks for the response, it's just not clear on what to do.

You need to pass in the matrix to hold the result.

Code:

expm(matrix_that_goes_in, matrix_that_comes_back);

more specifically, you need to pass in a pointer to a matrix to hold the result, so that changes made to it will show up in main.

Originally Posted by KBriggs

more specifically, you need to pass in a pointer to a matrix to hold the result, so that changes made to it will show up in main.

Note to OP: don't do this.

I got it to work. Thanks. Right now, I'm having a problem divinding by a friggin' integer... the first section of code is my factorial func

Code:

int factorial(int k){
	if (k==0){return 1;}
	int j=1;
	for (int i=1; i<=k; i++)
		j=i*j;
	return j;
}

This is in the main()

Code:

x=factorial(i);
	for (int j=0; j<2; j++)
		for (int k=0; k<2; k++)
		term[j][k]=mediate[j][k]/x;

when i print out term, it is a zero matrix, but if I multiply by x it works...

if divide it's zero, if multiply it does the right multiplication..

Were you expecting to be able to hold 0.5 in an int variable? Why did you expect this?

...uh...oops.

Originally Posted by tabstop

Note to OP: don't do this.

Sorry, not sure what I was thinking when I posted that. Must not have read very carefully.

I got the code to work. Thanks everyone for your help.