Python, solving Schrodinger equation using Runge-Kutta

[Milpool]

Homework Statement

I'm currently working on a project in which I have to solve the energy eigenvalues of the Schrödinger 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 Schrödinger equation.

def derivs(n, x, y):
      "Defines Schrödinger 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

Homework Equations

-

The Attempt at a Solution

I have attempted to print values of x and y[1] for a range of x as a test

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.

 

[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

 

[Milpool]

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

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 Schrödinger equations).
The library derivs is below

def derivs(n, x, y):
      "Defines Schrödinger 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!

 

[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:

regular solutions at the origin are of the form $R(r) = r^l \rho(r)$

and then writing that out:

the first term cancels the centrifugal potential term, and we find $$\rho '' (r)+ {2(l+1)\over r} \rho '(r) + \left (\epsilon - V(r)\rho(r) \right ) = 0$$
If the potential is regular at the origin, the singular term ${2(l+1)\over r}$ needs to be canceled by $\rho'(0) = 0$ .

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) ).

 

[Milpool]

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) ).

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.