In [1]:
# ELECTRON IN CUBE - solution of Schrodinger equation (for single variable X)
# * EM lectures: derivation of Schr.equation + detailed comments
# * this script: demonstration how Schr.equation can be solved in Python/SymPy
#   1) we take Schr.equation for eln in cube (for variable X)
#      (based on our previous derivation that F(x,y,z)=F(x)*F(y)*F(z)
#   2) we solve the equation just for X to receive characteristic function F(x)
#      (supposing that the solutions for Y and Z are analogous: F(x) ~ F(y) ~ F(z)
#      (this is logical due to symmetry; full justification is given in EM lectures
#   3) solution below = characteristic function F(x), which describes state of the eln
#      reminder: |F(x)|**2 dx = probability that we find electron in location dx
#      analogously: |F(x,y,z)|**2 dxdydz = probability...electron in location dxdydz
In [2]:
# (0) Initialize sympy
In [3]:
import sympy as sp
In [4]:
sp.init_printing()
In [5]:
# (1) Define equation to solve
In [6]:
# (1a) Prepare symbols
# (additional information, such as [positive=True] is usualy very useful
# (it may simplify solution and results greatly; it is like in wxMaxima, where the program asks
(x) = sp.symbols('x')
(F) = sp.Function('F', real=True)
(E,m,h) = sp.symbols('E m h', real=True, positive=True)
In [7]:
# (1b) Define equation
# (recommendation for unevaluated derivatives: sp.Derivative
e1 = sp.Eq(-h**2/2*m * sp.Derivative(F(x),x,2), E * F(x))
In [8]:
# (1c) This is the equation that we are going to solve...
display(e1)
$\displaystyle - \frac{h^{2} m \frac{d^{2}}{d x^{2}} F{\left(x \right)}}{2} = E F{\left(x \right)}$
In [9]:
# (2) Manual modification/simplification
# (multiply both sides by -h/(m*h**2)
In [10]:
e2 = sp.Eq(e1.lhs*(-2/(m*h**2)), e1.rhs*((-2/(m*h**2))))
In [11]:
display(e2)
$\displaystyle \frac{d^{2}}{d x^{2}} F{\left(x \right)} = - \frac{2 E F{\left(x \right)}}{h^{2} m}$
In [12]:
# (3) Manual modification/simplification
# (in right-hand-side of the equation, include all constants into one
In [13]:
# (3a) Prepare the constant
# (not necessary, but useful: we add information that the constant is real and positive
# (unless we do this here, the final solution contains k in absolute values |k|, which is ugly
(k) = sp.symbols('k', real=True, positive=True)
In [14]:
# (3b) Make the substitution of rhs = right-hand-side
e3 = sp.Eq(e2.lhs, e2.rhs.subs(((2*E)/(m*h**2)),k**2))
In [15]:
 display(e3)
$\displaystyle \frac{d^{2}}{d x^{2}} F{\left(x \right)} = - k^{2} F{\left(x \right)}$
In [16]:
# (4) Solve the equation
# (here the final solution is simple
# (the equation was manualy prepared in convenient form
# (if the constants were not so well defined, the solution is MUCH MORE complex!
In [17]:
e4 = sp.dsolve(e3,F(x))
In [18]:
display(e4)
$\displaystyle F{\left(x \right)} = C_{1} \sin{\left(k x \right)} + C_{2} \cos{\left(k x \right)}$
In [19]:
# (5) Verification of all possible solutions
In [20]:
# (5a) Three solutions depending on zero/non-zero values of constants C1 and C2
(A,B) = sp.symbols('A B', real=True, positive=True)
f1 = A*sp.sin(k*x)
f2 = B*sp.cos(k*x)
f3 = f1+f2
# Important note:
# - from the point of view of mathematics:
#   solutions f1,f2 are just special cases of f3
# - from the point of view of physics behind the scenes:
#   solutions f1,f2 are basic solutions (which have their quantum numbers)
#   solution f3 is linear combination of solutions f1,f2 which is also a solution
#   analogy with orbitals in atom:
#   ...primary solution: complex atomic orbitals = complex AO
#   ...linear combination of complex AO => real AO (which are used for visualizations)
#   ...linear combination of real AO => hybrid/hybridized AO (which are used in bond theory)
In [21]:
# (5b) For completness, display the equation and the solutions
print('Final form of Schrodinger equation that we have solved:')
display(e3)
print('Solutions that we found and that we want to verify:')
display(sp.Eq(F(x),f1))
display(sp.Eq(F(x),f2))
display(sp.Eq(F(x),f3))

Final form of Schrodinger equation that we have solved:
$\displaystyle \frac{d^{2}}{d x^{2}} F{\left(x \right)} = - k^{2} F{\left(x \right)}$ 
Solutions that we found and that we want to verify:
$\displaystyle F{\left(x \right)} = A \sin{\left(k x \right)}$
$\displaystyle F{\left(x \right)} = B \cos{\left(k x \right)}$
$\displaystyle F{\left(x \right)} = A \sin{\left(k x \right)} + B \cos{\left(k x \right)}$
In [22]:
# (5c) Prepare function for testing equality
# (just for convenience, not to repeat the same code three times
# (this single-purpose function should insert [solution] into our final [equation]=e3
def verify_solution(equation,solution):
    # LHS: (1) take lhs (2) substitute F(x)->solution, (3) perform derivatives
    # RHS: (1) take rhs (2) substitute F(x)->solution, (no doit ~ no derivatives
    eq_lhs = equation.lhs.subs(F(x),solution).doit()
    eq_rhs = equation.rhs.subs(F(x),solution)
    # Print results in simple/brief form (print instead of display)
    print('Verification of solution:',solution)
    print('LHS:',eq_lhs)
    print('RHS:',eq_rhs)
    print('LHS = RHS:', eq_lhs - eq_rhs == 0)
In [23]:
# (5d) Verify the three solutions
In [24]:
verify_solution(equation=e3, solution=f1)

Verification of solution: A*sin(k*x)
LHS: -A*k**2*sin(k*x)
RHS: -A*k**2*sin(k*x)
LHS = RHS: True
In [25]:
verify_solution(equation=e3, solution=f2)

Verification of solution: B*cos(k*x)
LHS: -B*k**2*cos(k*x)
RHS: -B*k**2*cos(k*x)
LHS = RHS: True
In [26]:
verify_solution(equation=e3, solution=f3)

Verification of solution: A*sin(k*x) + B*cos(k*x)
LHS: -k**2*(A*sin(k*x) + B*cos(k*x))
RHS: -k**2*(A*sin(k*x) + B*cos(k*x))
LHS = RHS: True
In [27]:
# CONCLUSIONS
# * we solved special case of Schrodinger equation: Derivative(F(x),x,2) = -k^2*F(x)
#   (this Schr.equation describes states of electron in cube along x-axis
# * we found the general solution: F(x) = A*sin(kx) + B*cos(kx)
#   (and its two special cases: F(x) = A*sin(kx), F(x) = B*cos(kx)
#   (in fact the physically meaningful solutions are just the two special cases!
# * we verified that all solutions are valid
#   (see the last section of this script
#   (more details => EM lectures; this is just a verification of the solution in sympy
# * important note - what is the physical meaning of the solution F(x)?
#   F(x) = state of the electron
#   |F(x)|**2 = probability density function
#   |F(x)|**2 dx = probability that the eln is localized in space element dx
# * complete solution in this case: F(x,y,z) = F(x)*F(y)*F(z)
#   (derivation of this solution -> EM-lectures