CODEWITHSUNDEEP
pythonscilab
user169808
user169808

Reputation: 523

using solve_bvp to solve Schrödinger equation

Hi I'm trying to sove the Schrödinger equation for a infinite square well with a barrior using solve_bvp but I don't understand how to do two things with solve_bvp. one is setting the range, the x axis I am trying to give it goes from -b to b but it is solving for 0 to b, I had previously tried solving with solve_ivp but I was running into the same problem. the other problem I have is how to set the boundary condition, I need to take into account the first and second derivative, but I'm a bit lost because I cannot understand the documentation very well. I was following an example found here, but it wasn't giving me the correct results so I'm trying to rewrite it. I'm also discussing the problem on another site but that one is a bit more general than what I'm trying to solve here.

this is the code that I'm trying to get to work

from scipy.integrate import solve_bvp
from pylab import *
import numpy as np

## dimentions         _
## |       :       |   | Vmax
## |    ___:___    |   |        _
## |___|       |___|  _| Vmin   _|- Vb
##-B  -A   0   A   B*Alpha
## \______/|\______/
##     Po       Pf

Vmax = 55
Vmin = 0
Vb = 50
A = 1
B = 4
alpha = 1
Po = -A-B
Pf = A+B*alpha
psi_b = array([0,1])   # Wave function boundry [first dir, secound dir]

N = 1000 #steps
psi = np.zeros([N,2])
x = linspace(Po, Pf, N)    # x-axis


def V(z):
    '''
    #Potential function in the finite square well.
    '''
    val=[]
    for i in z:
        if -A <=i <= A:
            val.append( Vb)
        elif i<=Po:
            val.append( Vmax)
        elif i>=Pf:
            val.append(Vmax)
        else:
            val.append(Vmin)
    return val

def boundary_conditions(psi_Po, psi_pf):
    return (psi_b)

def SE(z, p, E):
    state0 = p[1]
    state1 = 1.0*(V(z) - E)*p[0]
    return np.array([state0, state1])

def Wave_function(energy,psi):
    y_ = np.zeros((2, x.size))
    psi = solve_bvp(SE(x, psi, y_), boundary_conditions, x, y_)
    print(len(psi.y))
    return psi.y

def pltPotentail():
    plot(x,V(x))    

def main():
    en = linspace(0, Vb, 10)   # vector of energies where we look for the stable states
    figure(1)
    pltPotentail()
    Wave_function(0,psi)
    show()
if __name__ == "__main__":
    main()

Upvotes: 0

Views: 1045

Answers (1)

Russell Anderson
Russell Anderson

Reputation: 21

I've solved the finite square well potential using scipy.integrate.solve_bvp here.

Upvotes: 2

Related Questions