OiO.lk Community platform!

Oio.lk is an excellent forum for developers, providing a wide range of resources, discussions, and support for those in the developer community. Join oio.lk today to connect with like-minded professionals, share insights, and stay updated on the latest trends and technologies in the development field.
  You need to log in or register to access the solved answers to this problem.
  • You have reached the maximum number of guest views allowed
  • Please register below to remove this limitation

Excess work done on call (perhaps wrong Dfun type) with scipy.integrate.odeint

  • Thread starter Thread starter buzsh
  • Start date Start date
B

buzsh

Guest
I'm writing a code in Python that predicts the energy levels of Hydrogen which I will use as a template for research into quarkonium energy levels. I'm using the 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>I'm writing a code in Python that predicts the energy levels of Hydrogen which I will use as a template for research into quarkonium energy levels. I'm using the <code>scipy.integrate.odeint()</code> function to solve the Shroedinger equation and it works fine for the lower energy levels up to <code>n=6</code>. I don't expect I'll have much need to go beyond that, but odeint returns <code>Excess work done on this call (perhaps wrong Dfun type).</code> which only encourages me to extend what I can predict.</p>

<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(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
</code></pre>
 

Latest posts

G
Replies
0
Views
1
Gamal Othman
G
Top