# Simulating 1D time-independent Bose Einstein Condensation

Tags:
1. Jun 26, 2017

### svletana

Hello! I'm trying to simulate a one dimensional time independent BEC, I hope this is the right place to ask for help.

First of all, here's my code in Python.

Code (Python):
import sys
import numpy as np
import matplotlib.pyplot as plt

if len(sys.argv) == 1:
niter = 100
elif len(sys.argv) == 2:
niter = int(sys.argv[1]) # number of iterations
else:
print('Please insert only the iteration number')
exit(1)

# useful functions

# Laplacian operator for the second derivative
def Laplacian(x):
h = x[1]-x[0] # assume uniformly spaced points
n = len(x)
M = -2 * np.diag([1]*n)
for i in range(1,n):
M[i,i-1] = 1
M[i-1,i] = 1
return M/h ** 2

# for normalizing eigenstates
def Normalizate(U, x):
h = x[1] - x[0] # assume uniformly spaced points
n = len(x)
for j in range(0, n):
suma = 0.0
for i in range(1, n):
suma = suma + U[i, j]**2

suma = suma * h
rnorm = 1 / np.sqrt(suma)
# print j,’ integral (sin normalizar) =’,rnorm
# Normalization
rsign = 1
if U[1, j] < 0:
rsign = -1

rnorm = rnorm * rsign
for i in range(0, n):
U[i, j] = U[i, j] * rnorm

# parameters for Rubidium 87 (a.u.)
m = 87 # mass
a = 110 # scattering length
N = 500 # particle number
g = 4 * np.pi * a * N / m # coefficient for gross pitaevskii equation
w = 1 # frequancy of trapping potential
L = 1.5
dx = 600

# initialize wavefunction
x = np.linspace(-L, L, dx)
phi = [20 * np.exp(- (i ** 2)) for i in x]

W = 0.5 * m * (w ** 2) * (x ** 2)

# main loop
for i in range(niter):

# calculate density
density = [abs(p * p) for p in phi]

# schrodinger's eq in matrix form
T = -0.5 / m * Laplacian(x) # kinetic energy
P = g * np.diag(density)  #GPE non linear term

H = T + np.diag(W) + P #hamiltonian
E, U = np.linalg.eigh(H) #eigenvalues and eigenstates

Normalizate(U, x)

phi = U[:, 0]

print('chemical potential: ' + str(E[0]))

3. Jul 3, 2017

### svletana

Thanks for your help! :)
Instead of setting the real g value my teacher told me to try any value, so I did that. There was a problem with normalization as well, that was messing up the whole calculations. Thanks again!