# OK to post homework # Aurora Hiveley, 2/6/25, Assignment 5 Help := proc(): print(`RandomState(n,K), GeneralPauli(n)`): end: with(linalg): with(combinat): ## random state vector # RandomState(n,K): inputs a pos. integer n and a pos. integer K # outputs a random state vector with n complex components whose real and imaginary parts # are from -K to K, which is then normalized RandomState:=proc(n,K) local ra,i,V,k: ra:=rand(-K..K): V := [seq(ra()+I*ra(), i=1..n)]: k:=sqrt(add(abs(V[i])^2,i=1..n)): [seq(V[i]/k,i=1..n)]: end: ## experimenting with eigenvalues & eigenvectors # Experiment with three different random Hermitian matrices using RHM(n,K) with n=3 and K=10 # check that the eigenvalues are real and that the eigenvectors (pure states) are all orthogonal to each other. # for i from 1 to 3 do # A := RHM(3,10): # eigen := EVC(A): # for v in eigen[1] do # v := simplify(v): # if v<>conjugate(v) then ## print(v); ## RETURN(FAIL): # if an eigenvalue isn't real, returns FAIL # fi: # od: # pairs := choose(eigen[2],2): # for p in pairs do # if IP(p[1],p[2])<>0 then # RETURN(FAIL): # if a pair of eigenvectors isn't orthogonal, returns FAIL # fi: # od: # od: # some eigenvalues appear to come out imaginary, still. i suspect this is a result of maple's difficulty with simplification? # however, orthogonality seems to behave as expected! ## experimenting with prob dist # Experiment with five different random Hermitian matrices using RHM(n,K) with n=3 and K=10 # and corresponding five random state vector and apply ProbDist(L,A) to them. # For each case check that the prob dist consists of non-negative numbers that add up to 1. # for i from 1 to 5 do # M := RHM(3,10): # V := RandomState(3,10): # P := ProbDist(M,V): # for p in P do # if not p>=0 then # RETURN(FAIL): # fi: # od: # if simplify(add(p,p in P))<>1 then # RETURN(FAIL): # fi: # od: # there's also some issues here with imaginary numbers, perhaps because of the simplify command? i'm not certain ## general pauli matrices # GeneralPauli(n): inputs a NORMALIZED vector n=[nx,ny,nz], then constructs the 2x2 matrix # at the bottom of p. 349 of Susskind-Frieman GeneralPauli := proc(n) local i : if nops(n)<>3 or add(abs(n[i])^2,i=1..3)<>1 then # verifies that input vector is normalized and of length 3 RETURN(FAIL): fi: [[n[3], n[1]-I*n[2]], [n[1]+I*n[2], -n[3]]]: end: ## Find its eigenvalues and eigenvectors. # example: # M := GeneralPauli(RandomState(3,10)): # EVC(N); ### copied from C5 #C5.txt, Feb. 6, 2025 Experimental Mathematics (Rutgers, Dr. Z.) Help5:=proc(): print(` PM(), RHM(n,K) , EV(L), EVC1(L,g), EVC(L), IP(A,B), ProbDist(L,A) `): end: with(linalg): #PM(): the three famous Pauli matrices sigma_x,sigma_y, sigma_z PM:=proc(): [ [[0,1],[1,0]], [[0,-I],[I,0]], [[1,0],[0,-1]] ]: end: #RHM(n,K): a random n by n Hermitian matrix with entriries from -K to K RHM:=proc(n,K) local ra,i,j,L: ra:=rand(-K..K): for i from 1 to n do for j from 1 to i-1 do L[i,j]:=ra()+I*ra(): L[j,i]:=conjugate(L[i,j]): od: L[i,i]:=ra(): od: [seq([seq(L[i,j],j=1..n)],i=1..n)]: end: #EV(L): inputs a Hermitian matrix (say n by n) and outputs the LIST of eigenvalues EV:=proc(L) local x,L1,i,n: n:=nops(L): L1:=[seq(L[i]-[0$(i-1),x,0$(n-i)],i=1..n)]: [solve(det(L1),x)]: end: #EVC1(L,g): inputs a matrix L and an eigenvalue g finds the NORMALIZED eigenvector EVC1:=proc(L,g) local v,i,j,var,eq,V,n,L1,V1,va,k: n:=nops(L): L1:=[seq(L[i]-[0$(i-1),g,0$(n-i)],i=1..n)]: if simplify(det(L1))<>0 then RETURN(FAIL): fi: V:=[seq(v[i],i=1..n)]: var:={seq(v[i],i=1..n)}: eq:={seq(add(L1[i,j]*V[j],j=1..n),i=1..n)}: var:=solve(eq,var): #V1 the set of unknowns that is arbitrary V1:={}: for va in var do if op(1,va)=op(2,va) then V1:=V1 union {op(1,va)}: fi: od: if nops(V1)<>1 then RETURN(FAIL): fi: V:=subs(var,V): V:=subs(V1[1]=1,V): k:=sqrt(add(abs(V[i])^2,i=1..n)): [seq(V[i]/k,i=1..n)]: end: #EVC(L): inputs a Hermitian matrix and outputs the pair ListOfEigenvalues, ListOfEigenvectors EVC:=proc(L) local n,ev,evc,i: n:=nops(L): ev:=EV(L): evc:=[seq(EVC1(L,ev[i]),i=1..n)]: [ev,evc]: end: #IP(A,B): the inner product of state-vectors A and B IP:=proc(A,B) local n,i: n:=nops(A): if nops(B)<>n then RETURN(FAIL): fi: simplify(add(A[i]*conjugate(B[i]),i=1..n)): end: #ProbDist(L,A): Given an obervalble represented by the Hermitian matrix L and a state vector A #outputs first all the eigenvectors (pure states) and the prob. of winding in the respective pure states #followed by the expected value of A (i.e. if you make the observation L and the state is A the average in the long run ProbDist:=proc(L,A) local evc,i,A1,k,PD,n: n:=nops(L): k:=sqrt(add(abs(A[i])^2,i=1..n)): A1:=[seq(A[i]/k,i=1..n)]: evc:=EVC(L): PD:=[seq(abs(IP(A1,evc[2][i]))^2,i=1..n)]: [PD,add(PD[i]*evc[1][i],i=1..n)]: end: