Python, solving Schrodinger equation using Runge-Kutta

Tags:
1. Mar 29, 2015

Milpool

1. The problem statement, all variables and given/known data
I'm currently working on a project in which I have to solve the energy eigenvalues of the Schrodinger equation to compute the mass of certain Mesons. We've been taught very little programming (so apologies that my understanding is very basic), and are therefore given any stock programs we will need, for example my Runge-Kutta. I need to define a function that guesses a value of energy, solves my equation and spits out the solution for r=0, so that I can minimise this to satisfy the boundary condition of the problem that u(r=0)=0. I cannot seem to work out how to return u(0) from my Runge-Kutta program. Below is the runkut program and my definition of my Schrodinger equation.

Code (Python):
def derivs(n, x, y):
"Defines Schrodinger equation"
dy=zeros(n+1,float)
if x==0: y[2]==0 #otherwise divide by zero error
else:
dy[1]=y[2]
dy[2]=(-E+x+((L)/x**2))*y[1]
return dy

def runkut(n, x, y, h):
"Advances the solution of diff eqn defined by derivs from x to x+h"
y0=y[:]
k1=derivs(n, x, y)
for i in range(1,n+1): y[i]=y0[i]+0.5*h*k1[i]
k2=derivs(n, x+0.5*h, y)
for i in range(1,n+1): y[i]=y0[i]+h*(0.2071067811*k1[i]+0.2928932188*k2[i])
k3=derivs(n, x+0.5*h, y)
for i in range(1,n+1): y[i]=y0[i]-h*(0.7071067811*k2[i]-1.7071067811*k3[i])
k4=derivs(n, x+h, y)
for i in range(1,n+1):
a=k1[i]+0.5857864376*k2[i]+3.4142135623*k3[i]+k4[i]
y[i]=y0[i]+0.16666666667*h*a

x+=h
return (x,y)

N=150
N0=34
L=2.0
E=1.0 #introduce loop
xmax=25.0
xmin=0.0
dx=(xmax-xmin)/N

2. Relevant equations
-

3. The attempt at a solution
I have attempted to print values of x and y[1] for a range of x as a test
Code (Python):

x=0.0; y=[0, 0.0, 0.0]

z=[0 for j in range(0,N)]
for j in range(0,N):
(x,y) = runkut(2,x,y,dx)
print x, y[1]
z[j]=y[1]

However all I am returned with is the value of x increasing with y[1]=0.0 each time. Also tried the method for x decreasing from x=10 down to zero, and was met with the same output. Also tried to print y, which just returned the same array defined above. I am somewhat confused as a library program which is available to us defines a different differential equation, uses the same runge kutta and with the above piece of code returns the numerical solutions of the equation in question.

As a note once I have a functioning code returning u(0), am then comfortable as to where to go from there, it is just wrapping my head around how to simply return u(0). Any help would be hugely appreciated, thanks.

Last edited: Mar 29, 2015
2. Mar 29, 2015

BvU

What is the equation you are trying to solve ?
My guess is that the library example is solving a rather different kind of equation

3. Mar 29, 2015

Milpool

My equation is

$$\frac{d^2 u}{d x^2} - \left[ \epsilon + \frac{l(l+1)}{x^2} - V(x) \right] u = 0 \ \ (V(x) = x)$$
Where I am trying to eventually recover values of epsilon. The library program has been edited to (hopefully) solve this in my code above.
The library program I used was solving
$$\frac{d^2 u}{d x^2} + \left[ \epsilon - V(x) \right] u = 0 \ \ \ \ (V(x) = e^{-2 \beta x} - e^{\beta x})$$
Where the original runkut is identical, and the derivs function was slightly different (as they are obviously different Schrodinger equations).
The library derivs is below
Code (Python):
def derivs(n, x, y):
"Defines Schrodinger equation"
dy=zeros(n+1,float)
dy[1]=y[2]
dy[2]=(-E+V(x))*y[1]
return dy
And V(x) is defined at a later point in the program. Thanks for the reply!

4. Mar 30, 2015

BvU

Yes. My issue is that your program is producing solutions of a differential equation that isn't yet defined (E is unknown), whilst what you really want is to find the discrete values of E for which reasonable solutions exist at all. So there must be some higher level loop around this thing you have now.

For the moment the thing doesn't get going, because you start off with y = 0 -- and it stays zero. So try nudging either y(1) or y(2) away from 0 at x = 0.

PS Since this is a radial problem: I miss a term ${\tfrac {2} {x} }{\partial u\over \partial x}$ when I compare with e.g. eqn (11) here

and: this thread anything for you ? It's not the same and it isn't the solution, but it might help thinking things over.

PPS here I found another goody:
and then writing that out:
Perhaps you can use some of that nifty approach too. And the pictures further down show the kind of results you hope to achieve too (with a different V(r) ).

Last edited: Mar 30, 2015
5. Mar 30, 2015

Milpool

Thanks, checking those examples out now. Yeah sorry I didn't make this clear in my explanation, for the time being I was simply trying to get something together that returns u(0) for a single guess value of E. After which I will be using the fmin function to minimise this so that u(0)->0, and then the value of E will be returned.

Last edited: Mar 30, 2015