B

#### buzsh

##### Guest

`scipy.integrate.odeint()`

function to solve the Shroedinger equation and it works fine for the lower energy levels up to `n=6`

. I don't expect I'll have much need to go beyond that, but odeint returns `Excess work done on this call (perhaps wrong Dfun type).`

which only encourages me to extend what I can predict.The Shroedinger equation substitution I'm using is:

Code:

```
u'' - (l*(l+1)/r**2 - 2mu_e(E-V_emag(r))) * u = 0
=>
u' = v
v' = ((l*(l+1))/(r**2) - 2.0*mu_e*(E - V_emag(r)))*u
```

I'm then using

`scipy.integrate.odeint()`

on it and iterating through energies and using other functions I've defined that assess turning points and nodes in the result. The way I find the energy levels is finding the lowest possible E value where the number of turning points and nodes matches what it should; then incrementing `L`

by 1 and finding the new ground energy, e.g. if `L=0`

I'll find `n=1`

energy and if `L=3`

, I'll find the `n=2`

energy.Once the code increments to

`L=7`

it doesn't return anything useful. The range of `r`

has been extended but I've tried with keeping it the same to reduce the number of steps but to no avail. The code is self-taught so in my research I've read about Jacobians. I'm afraid I haven't worked out what they are or if I need one yet. Any ideas?
Code:

```
def v_emag(r):
v = -alpha/r
return v
def s_e(y,r,l,E): #Shroedinger equation for electromagntism
x = numpy.zeros_like(y)
x[0] = y[1]
x[1] = ((l*(l+1))/(r**2) - 2.0*mu_e*(E - V_emag(r)))*y[0]
return x
def i_s_e(l,E,start=0.001,stop=None,step=(0.005*bohr)):
if stop is None:
stop = ((l+1)*30-10)*bohr
r = numpy.arange(start,stop,step)
y = odeint(s_e,y0,r,args=(l,E))
return y
def inormalise_e(l,E,start=0.001,stop=None,step=(0.005*bohr)):
if stop is None:
stop = ((l+1)*30-10)*bohr
r = numpy.arange(start,stop,step)
f = i_s_e(l,E,start,stop,step)[:,0]
f2 = f**2
area = numpy.trapz(f2,x=r)
return f/(numpy.sqrt(area))
def inodes_e(l,E,start=0.001,stop=None,step=(0.005*bohr)):
if stop is None:
stop = ((l+1)*30-10)*bohr
x = i_s_e(l,E,start,stop,step)
r = numpy.arange(start,stop,step)
k=0
for i in range(len(r)-1):
if x[i,0]*x[i+1,0] < 0: #If u value times next u value <0,
k+=1 #crossing of u=0 has occured, therefore count node
return k
def iturns_e(l,E,start=0.001,stop=None,step=(0.005*bohr)):
if stop is None:
stop = ((l+1)*30-10)*bohr
x = i_s_e(l,E,start,stop,step)
r = numpy.arange(start,stop,step)
k = 0
for i in range(len(r)-1):
if x[i,1]*x[i+1,1] < 0: #If du/dr value times next du/dr value <0,
k=k+1 #crossing of du/dr=0, therefore a maximum/minimum
return k
l = 0
while l < 10: #The ground state for a system with a non-zero angular momentum will
E1 = -1.5e-08 #be the energy level of principle quantum number l-1, therefore
E3 = 0 #by changing l, we can find n by searching for the ground state
E2 = 0.5*(E1+E3)
i = 0
while i < 40:
N1 = inodes_e(l,E1)
N2 = inodes_e(l,E2)
N3 = inodes_e(l,E3)
T1 = iturns_e(l,E1)
T2 = iturns_e(l,E2)
T3 = iturns_e(l,E3)
if N1 != N2:# and T1 != T2: #Looks in lower half first, therefore will tend to ground state
E3 = E2
E2 = 0.5*(E1+E3)
elif N2 != N3:# and T2 != T3:
E1 = E2
E2 = 0.5*(E1+E3)
else:
print "Can't find satisfactory E in range"
break
i += 1
x = inormalise_e(l,E2)
if x[((l+1)**2)/0.005] > (x[2*((l+1)**2)/0.005]) and iturns_e(l,E2+1e-20)==1:
print 'Energy of state: n =',(l+1),'is: ',(E2*(10**9)),'eV'
l += 1
else:
E1 = E2+10e-20
```

<p>The Shroedinger equation substitution I'm using is:</p>

<pre><code>u'' - (l*(l+1)/r**2 - 2mu_e(E-V_emag(r))) * u = 0

=>

u' = v

v' = ((l*(l+1))/(r**2) - 2.0*mu_e*(E - V_emag(r)))*u

</code></pre>

<p>I'm then using <code>scipy.integrate.odeint()</code> on it and iterating through energies and using other functions I've defined that assess turning points and nodes in the result. The way I find the energy levels is finding the lowest possible E value where the number of turning points and nodes matches what it should; then incrementing <code>L</code> by 1 and finding the new ground energy, e.g. if <code>L=0</code> I'll find <code>n=1</code> energy and if <code>L=3</code>, I'll find the <code>n=2</code> energy.</p>

<p>Once the code increments to <code>L=7</code> it doesn't return anything useful. The range of <code>r</code> has been extended but I've tried with keeping it the same to reduce the number of steps but to no avail. The code is self-taught so in my research I've read about Jacobians. I'm afraid I haven't worked out what they are or if I need one yet. Any ideas?</p>

<pre><code>def v_emag(r):

v = -alpha/r

return v

def s_e(y,r,l,E): #Shroedinger equation for electromagntism

x = numpy.zeros_like

x[0] = y[1]

x[1] = ((l*(l+1))/(r**2) - 2.0*mu_e*(E - V_emag(r)))*y[0]

return x

def i_s_e(l,E,start=0.001,stop=None,step=(0.005*bohr)):

if stop is None:

stop = ((l+1)*30-10)*bohr

r = numpy.arange(start,stop,step)

y = odeint(s_e,y0,r,args=(l,E))

return y

def inormalise_e(l,E,start=0.001,stop=None,step=(0.005*bohr)):

if stop is None:

stop = ((l+1)*30-10)*bohr

r = numpy.arange(start,stop,step)

f = i_s_e(l,E,start,stop,step)[:,0]

f2 = f**2

area = numpy.trapz(f2,x=r)

return f/(numpy.sqrt(area))

def inodes_e(l,E,start=0.001,stop=None,step=(0.005*bohr)):

if stop is None:

stop = ((l+1)*30-10)*bohr

x = i_s_e(l,E,start,stop,step)

r = numpy.arange(start,stop,step)

k=0

for i in range(len(r)-1):

if x[i,0]*x[i+1,0] < 0: #If u value times next u value <0,

k+=1 #crossing of u=0 has occured, therefore count node

return k

def iturns_e(l,E,start=0.001,stop=None,step=(0.005*bohr)):

if stop is None:

stop = ((l+1)*30-10)*bohr

x = i_s_e(l,E,start,stop,step)

r = numpy.arange(start,stop,step)

k = 0

for i in range(len(r)-1):

if x[i,1]*x[i+1,1] < 0: #If du/dr value times next du/dr value <0,

k=k+1 #crossing of du/dr=0, therefore a maximum/minimum

return k

l = 0

while l < 10: #The ground state for a system with a non-zero angular momentum will

E1 = -1.5e-08 #be the energy level of principle quantum number l-1, therefore

E3 = 0 #by changing l, we can find n by searching for the ground state

E2 = 0.5*(E1+E3)

i = 0

while i < 40:

N1 = inodes_e(l,E1)

N2 = inodes_e(l,E2)

N3 = inodes_e(l,E3)

T1 = iturns_e(l,E1)

T2 = iturns_e(l,E2)

T3 = iturns_e(l,E3)

if N1 != N2:# and T1 != T2: #Looks in lower half first, therefore will tend to ground state

E3 = E2

E2 = 0.5*(E1+E3)

elif N2 != N3:# and T2 != T3:

E1 = E2

E2 = 0.5*(E1+E3)

else:

print "Can't find satisfactory E in range"

break

i += 1

x = inormalise_e(l,E2)

if x[((l+1)**2)/0.005] > (x[2*((l+1)**2)/0.005]) and iturns_e(l,E2+1e-20)==1:

print 'Energy of state: n =',(l+1),'is: ',(E2*(10**9)),'eV'

l += 1

else:

E1 = E2+10e-20

</code></pre>