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