Definition of boundaries for the 1D quantum harmonic oscillator


202
views
0
5 months ago by
Hello,

I'm solving the 1D quantum harmonic oscillator equation \[ u'' + (x^2 - K) u = 0, \qquad u(-x_0) = u(x_0) = 0. \] When I define the boundaries of the problem using 'near(...)', I don't get the correct solutions to the equation. If instead I use the "x[0] - bound > tol ..." option as shown in the code, I get the correct answer. Can any one shine some light on how to use correctly the "near" method so I get the right solutions to the problem? The reason for my question is the following: when I'm solving the infinite square well problem \[ u'' - Ku = 0, \qquad u(0) = u(x_0) = 0, \] the boundary conditions using "near" are the ones that work this time, but not the others. I really appreciate any insight on what's wrong with the definition of boundary conditions in the programs.

This is the code for the harmonic oscillator:

'''1D Quantum Harmonic Oscillator'''

import matplotlib.pyplot as plt
from fenics import *

N = 500            # Number of finite elements
leftBound = -5.0   # Left bound
rightBound = 5.0   # Right bound
funcFamily = 'CG'  # Family of interpolation polynomial
funcDegree = 4     # Degree of interpolation polynomial
scale = 1E-7       # Scale factor
tol = DOLFIN_EPS   # Tolerance

K = 5.0  # Dimensionless energy (for n = 0, 1, 2, then K = 1, 3, 5, ...)

# Create a mesh and define function space
mesh = IntervalMesh(N, leftBound, rightBound)
V = FunctionSpace(mesh, funcFamily, funcDegree)

# Dirichlet conditions for left and right boundaries
u_D = Expression('x[0] - 0', degree=1)

class LeftBoundary(SubDomain):
    def inside(self, x, on_boundary):
        return on_boundary and (x[0] - leftBound) < tol   # This one works
        #return on_boundary and near(x[0], leftBound)     # This one doesn't

class RightBoundary(SubDomain):
    def inside(self, x, on_boundary):
        return on_boundary and (x[0] - rightBound) > tol   # This one works
        #return on_boundary and near(x[0], rightBound)     # This one doesn't

# Collecting boundary conditions
LB = DirichletBC(V, u_D, LeftBoundary())
RB = DirichletBC(V, u_D, RightBoundary())
bcs = [LB, RB]

# Define variational problem
u = TrialFunction(V)
v = TestFunction(V)
b = Expression('x[0]*x[0] - K', degree=2, K=K)
f = Constant(0.0)
a = dot(grad(u),grad(v))*dx + dot(b*u,v)*dx
L = f*v*dx

# Compute solution
u = Function(V)
solve(a == L, u, bcs)

# Plot solution and mesh
plt.subplot(211)
plot(u)

plt.subplot(212)
plot(dot(u,u))
​

This is the code for the infinite quantum well:

'''1D Infinite Square Well'''

import matplotlib.pyplot as plt
from fenics import *

N = 200            # Number of finite elements
leftBound = 0.0    # Left bound
rightBound = 1.0   # Right bound
funcFamily = 'CG'  # Family of interpolation polynomial
funcDegree = 5     # Degree of interpolation polynomial
tol = DOLFIN_EPS   # Tolerance

n = 3                                 # Quantum number
width = abs(rightBound - leftBound)   # Well's width
K = n*n*pi*pi/(width*width)           # Dimensionless energy

# Create a mesh and define function space
mesh = IntervalMesh(N, leftBound, rightBound)
V = FunctionSpace(mesh, funcFamily, funcDegree)

# Dirichlet conditions for left and right boundaries
u_D = Expression('x[0] - 0', degree=1)

class Left(SubDomain):
    def inside(self, x, on_boundary):
        #return on_boundary and (x[0] - leftBound) < tol   # not working
        return near(x[0], leftBound)                       # working

class Right(SubDomain):
    def inside(self, x, on_boundary):
        #return on_boundary and (x[0] - rightBound) > tol  # not working
        return near(x[0], rightBound)                      # working

left = Left()
right = Right()

boundaries = MeshFunction("size_t", mesh, mesh.topology().dim() - 1)
boundaries.set_all(0)
left.mark(boundaries, 1)
right.mark(boundaries, 2)

# Collecting boundary conditions
bcL = DirichletBC(V, u_D, boundaries, 1)
bcR = DirichletBC(V, u_D, boundaries, 2)
bcs = [bcL, bcR]

# Define variational problem
u = TrialFunction(V)
v = TestFunction(V)
b = Expression('K', degree=0, K=K)
f = Constant(0.0)
a = dot(grad(u),grad(v))*dx - dot(b*u,v)*dx
L = f*v*dx

# Compute solution
u = Function(V)
solve(a == L, u, bcs)

# Plot solution
plt.subplot(211)
plot(u)

plt.subplot(212)
plot(dot(u,u))
Community: FEniCS Project

2 Answers


0
5 months ago by
So for the harmonic oscillator the method using near works just fine for me but your other method doesn't.
Adding abs() obviates this problem for me:
class LeftBoundary(SubDomain):
    def inside(self, x, on_boundary):
        return on_boundary and abs(x[0] - leftBound) < tol

class RightBoundary(SubDomain):
    def inside(self, x, on_boundary):
        return on_boundary and abs(x[0] - rightBound) < tol​

Same goes for your second code. Just use the absolute value and require it to be less than your tolerance.

0
5 months ago by
Klunkean is right, in that the near() function is doing exactly what you tell it to, and applying the BCs, while the BC on one end is failing when you use inequalities without abs().  However, the inequality-based BCs without abs() only "work" by chance in some case. 

Based on hazy recollections of my undergraduate physics classes, I think the "real problem" here is that the solution you are looking for is the solution to the eigenproblem on the infinite domain.  Using the odd function \(x\) for a BC only reproduces the desired eigenmodes from the infinite-domain problem when the eigenmodes themselves are odd functions (i.e., \(K=3,7,\ldots\)).  Essentially, the magnitude of the BC selects the arbitrary scaling factor of the eigenfunction.  (Physically, this scaling factor is selected by wavefunction normalization.)  If you use a nonzero even function, e.g.,

u_D = Constant(1.0e-8)​

for \(K=1,5,\ldots\), you will see that it produces the corresponding eigenfunctions of the Hamiltonian (up to the arbitrary constant factor) when using near() to enforce both BCs.
Please login to add an answer/comment or follow this question.

Similar posts:
Search »