# OK to post homework # Aurora Hiveley, 2/11/25, Assignment 6 Help := proc(): print(`AllHM2(K), RealSV2(K), UncertaintyContest(K)`): end: with(combinat): with(linalg): # AllHM2(K): inputs a positive integer K and outputs ALL 2 by 2 matrices of the form [[a,b+I*c],[b-I*c,d]] # where a,b,c,d, are integers such that 1 ≤ a,b,c,d ≤ K (there are K^4 of them) AllHM2 := proc(K) local a,b,c,d,M: M := {}: for a from 1 to K do for b from 1 to K do for c from 1 to K do for d from 1 to K do M := M union {[[a,b+I*c],[b-I*c,d]]}: od: od: od: od: M: end: # RealSV2(K): inputs a positive integer K and outputs all normalizations of the vectors (i.e. make them unit vectors) # [a,b] where a and b are integers 1 ≤ a,b ≤ K (there are K2 of them RealSV2 := proc(K) local a,b,k,U: U := {}: for a from 1 to K do for b from 1 to K do k := sqrt(a^2 + b^2): # normalization factor U := U union {[a/k,b/k]}: od: od: U: end: # UncertaintyContest(K): inputs a positive integer K and outpts the triplet [A,B,V] (whose entries are in 1..K, V normalized) # that make CheckUP(A,B,V) as small as possible and as big as possible (if A and B commute, do not count them). UncertaintyContest := proc(K) local P,U,p,u,m,M,l: P := choose(AllHM2(K),2): U := RealSV2(K): for p in P do if Com(p[1],p[2])=[[0,0],[0,0]] then P := P minus {p}: # get rid of pairs which commute fi: od: m := [CheckUP(P[1][1], P[1][2],U[1]),[P[1][1],P[1][2],U[1]]]: # initialize max and min trackers M := [CheckUP(P[1][1], P[1][2],U[1]),[P[1][1],P[1][2],U[1]]]: for p in P do for u in U do l := CheckUP(p[1],p[2],u): if l < m[1] then m := [l, [p[1],p[2],u]]: fi: if l > M[1] then M := [l, [p[1],p[2],u]]: fi: od: od: [m[2],M[2]]: end: # What are UncertaintyContest(2)?, UncertaintyContest(3)? # UncertaintyContest(2) outputs: # smallest: [[[1, 2 + I], [2 - I, 2]], [[2, 2 + I], [2 - I, 1]], [sqrt(2)/2, sqrt(2)/2]] # biggest: [[[1, 1 + I], [1 - I, 1]], [[1, 1 + 2*I], [1 - 2*I, 2]], [(2*sqrt(5))/5, sqrt(5)/5]] # UncertaintyContest(3) outputs: # smallest: [[[1, 2 + I], [2 - I, 3]], [[3, 3 + I], [3 - I, 1]], [sqrt(2)/2, sqrt(2)/2]] # biggest: [[[1, 2 + 3*I], [2 - 3*I, 3]], [[2, 3 + 2*I], [3 - 2*I, 2]], [(3*sqrt(13))/13, (2*sqrt(13))/13]] ### copied from C6.txt # C6.txt, 2/10/25 #Feb. 10, 2025 C6.txt Help6:=proc(): print(`Mul(A,B), Com(A,B), RSV(n,K), Del(L,A), CheckUP(A,B,V) `): end: #Mul(A,B): the product of matrix A and B (assuming that it exists) Mul:=proc(A,B) local i,j,k,n: n:=nops(A): [seq([seq( add(A[i][k]*B[k][j] , k=1..n) , j=1..n)],i=1..n)]: end: #Com(A,B: The commutator of A and B Com:=proc(A,B): Mul(A,B)-Mul(B,A):end: #RSV(n,K): RSV:=proc(n,K) local ra,i,v,a: ra:=rand(-K..K): v:=[seq(ra()+I*ra(),i=1..n)]: a:=sqrt(add(v[i]*conjugate(v[i]),i=1..n)): [seq(v[i]/a,i=1..n)]: end: #Del(L,A): the "avergae error", i.e. standard deviation pf the prob. distribution of #measuring the observable L when the state is A Del:=proc(L,A) local n,pd,ev,i,av: n:=nops(L): pd:=ProbDist(L,A): ev:=pd[1]: av:=pd[3]: pd:=pd[2]: sqrt(add(ev[i]^2*pd[i],i=1..n)-av^2): end: #CheckUP(A,B,V): the ratio of Del(A,V)*Del(B,V)/Expectation(Com(A,B),V) (it should be >=1/2) CheckUP:=proc(A,B,V) local C: C:=Com(A,B): if abs(ProbDist(C,V)[3])=0 then RETURN(FAIL): fi: evalf(Del(A,V)*Del(B,V)/abs(ProbDist(C,V)[3])): end: #old stuff #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)]: [evc[1],PD,add(PD[i]*evc[1][i],i=1..n)]: end: