#----------------------------------------# # Week 12: Die Funktion # # < a_ell^+ a_ell >(t) # #----------------------------------------# #--------------------# # Input Parameter: # #--------------------# L = 2 # Number of Lattice Sites N = 20 # Total Number of Particles #L = 5 #N = 10 eps = 1 # H = eps * H0 + u * Hint g = 0.5 u = g/N T = 500 # Time Horizon, t in [0,T] M = 20000 # M Zeitpunkte zum Plotten #psi0 = c(0,0,N,0,0) # Anfangszustand, darf einer der Basis-Vektoren sein, # also ein Vektor c(n1,n2,...,nL) mit Summe nj = N psi0 = c(N,0) # nj Teilchen am Gitterplatz j #-------------------------------# # set up many body basis # #-------------------------------# nB = choose(N+L-1,L-1) # = |B| nB # etwa 56 fuer L=4, N=5 B0 = t(combn(N+L-1,L-1)) # t = transponieren, dann nB Zeilen B = matrix(0,nrow=nB,ncol=L) # soll werden Basis B = {(n1,...,nL)| sum nj = N, nj>=0} for(j in 1:nB) { b = B0[j,] # das b sind jetzt 2 Trennstriche, die bei b[1] und b[2] bb = c( 0 , b , N+L ) # auf die Plus-Zeichen in (4) weiter oben gelegt werden m = diff( bb ) # das sind jetzt die mj aus (3) weiter oben n = m - rep(1,L) # und das dann also die nj aus (1) B[j,] = n } B # ist jetzt die gewuenschte Basis B mit dim(B) # B = {(n1,...,nL)| sum nj = N, nj>=0} #-------------------------------# # END set up many body basis # #-------------------------------# # L N nB nB^2 # 3 10 66 4'356 # 3 20 231 53'361 # 5 10 1001 1 Mio # 6 10 3003 9'018'009 # 6 12 6188 38'291'344 # 6 15 15504 240'374'016 # 6 20 53130 2'822'796'900 #-------------------------------# # set up Hamiltonian # #-------------------------------# # allocate some storage place for the result: nzmax = 10^7 # non zero matrix elements h0 = rep( 0 , nzmax ) # kinetic energy hint = rep( 0 , nzmax ) # diagonal part rows = rep( 0 , nzmax ) # location in full nB x nB sparse matrix cols = rep( 0 , nzmax ) # location in full nB x nB sparse matrix id = rep( 0 , nzmax ) # soll die Einheitsmatrix werden, # aber in dieser sparse-representation t0 = Sys.time() # monitor calc time nz = 1 # counter of non-zero matrix elements for(i in 1:nB) { m = B[i,] # Basisvektor m for(j in 1:nB) { n = B[j,] # Basisvektor n , wir berechnen # < m | H | n > # diagonal part: if(i==j) { res = 0 for(k in 1:L) { res = res + n[k]*(n[k]-1) # das ist (5) aus pdf-file } h0[nz] = 0 hint[nz] = res id[nz] = 1 rows[nz] = i cols[nz] = j nz = nz + 1 } # off-diagonal part: d = m - n k = which( d != 0 ) # nuetzlicher Befehl.. if( length(k)==2 ) { k1 = k[1] k2 = k[2] if( d[k1]*d[k2] == -1 & k2-k1==1 ) # we have a non zero contribution { if( m[k1] == n[k1]+1 ) { res = sqrt(n[k1]+1)*sqrt(n[k2]) # das ist die erste Zeile von (4).pdf h0[nz] = res hint[nz] = 0 id[nz] = 0 rows[nz] = i cols[nz] = j nz = nz + 1 } if( m[k1] == n[k1]-1 ) { res = sqrt(n[k1])*sqrt(n[k2]+1) # das ist die zweite Zeile von (4).pdf h0[nz] = res hint[nz] = 0 id[nz] = 0 rows[nz] = i cols[nz] = j nz = nz + 1 } } # END we have a non zero contribution } # END if( length(k)==2 ) } # NEXT j, basis element n = B[j,] # nur zur Info: if(i%%100==0) { print( c(i,nB) ) print( Sys.time() - t0 ) } } # NEXT i, basis element m = B[i,] # remove unused elements: if( nz < nzmax ) { h0 = h0[-(nz:nzmax)] hint = hint[-(nz:nzmax)] id = id[-(nz:nzmax)] rows = rows[-(nz:nzmax)] cols = cols[-(nz:nzmax)] } h = eps*h0 + u*hint nz # number of nonzero elements # etwa 23024 for L=6, N=10 Sys.time() - t0 # Finished setting up h #-------------------------------# # END set up Hamiltonian # #-------------------------------# #-----------------------------------# # store many-body hamiltonian # # as a sparse matrix # #-----------------------------------# # wenn man eine Matrix als duennbesetzte Matrix # speichern moechte, muss man das Matrix-package # laden: require(Matrix) sH = sparseMatrix(rows,cols,x=h,dims=c(nB,nB)) #---------------------------------# # Eigenwerte und Eigenvektoren: # #---------------------------------# calctime = Sys.time() res = eigen(sH, symmetric=TRUE) Sys.time() - calctime e = res$values V = res$vectors Vadj = t(V) # reelle Matrix test = Vadj %*% sH %*% V zapsmall(test) # ist Diagonalmatrix mit e auf Diagonale, tail(e) # also H V[,j] = e[j] V[,j] test2 = Vadj %*% V zapsmall(test2) # Einheitsmatrix #----------------------------------------# # # # Die Funktion < a_ell^+ a_ell >(t) # # # #----------------------------------------# t = seq(0,T,length=M) j0 = 0 for(j in 1:nB) { if( sum(B[j,]==psi0) == L ) { j0 = j } } if(j0==0) { stop("unzulaessiges psi0") } # U := = sum_j exp(-ite[j]) * * # = sum_j V[m,j] * exp(-ite[j]) * Vadj[j,j0] # =: sum_j V[m,j] * v0t[j,] # # U[m,] ist Vektor der Laenge M = Anzahl Zeitpunkte zum Plotten, # m indiziert Vielteilchen-Basisvektoren, ist kein Zeitindex U = matrix(0, nrow=nB, ncol=M ) Abs2U = matrix(0, nrow=nB, ncol=M ) # |U|^2 v0 = Vadj[,j0] v0t = matrix(0,nrow=nB,ncol=M) tcalc = Sys.time() for(j in 1:nB) { v0t[j,] = exp(-1i*t*e[j]) * v0[j] } U = V %*% v0t Abs2U = abs(U)^2 Sys.time() - tcalc En = matrix(0, nrow = L, ncol = M ) En[,1] = psi0 res = rep(0,M) for(ell in 1:L) { res = rep(0,M) for(j in 1:nB) { m = B[j,] res = res + m[ell] * Abs2U[j,] } En[ell,] = res } info = paste("g = ",g) plot(t,En[1,],ylim=c(0,N),main=info,ylab="En(t)") EN = En[1,] # EN = sum_ell < a_ell^+ a_ell >(t) for(j in 2:L) { points(t,En[j,],col=j) EN = EN + En[j,] readline("press enter to continue..") } points(t,EN,col="lightgrey") # check: muss gleich N sein