1
$\begingroup$

I am working on a modified version of Schrodinger's equation (time-independent) where $\frac{d^2ψ}{dx^2}=-2(E-V)ψ$, where I have to consider $V = 0$ at all times.

I have been asked to use Python in order to solve this equation. I have used the following code:

e = 200

dx = 0.01   
xlim = 100

p = np.zeros(xlim+1)
x = np.zeros(xlim+1)

x[0] = 0
p[0] = 1 

for i in range(1, xlim):
   p[i+1] = 2 * p[i] - p[i-1] + (dx * dx) * (-2 * e * p[i])
   x[i+1] = x[i] + dx


plt.plot(x[1:], p[1:])
plt.show()

and have gotten the following output:

Wavefunction

When changing dx (the step size) to a slightly larger number, the oscillation doesn't look constant. Is this normal?:

enter image description here

$\endgroup$
8
  • $\begingroup$ Sure, a scaled sine function (e.g. $5\sin(x)$) still qualifies as a "sinusoidal oscillation". $\endgroup$
    – Tyberius
    Dec 20 '21 at 20:26
  • $\begingroup$ @Tyberius thanks for that. I have added an additional part to the question, does that mean that the algorithm is wrong? $\endgroup$
    – mkdirBed
    Dec 20 '21 at 21:18
  • 2
    $\begingroup$ How much larger is "slightly larger"? Too large and you might not be sampling finely enough, leading to sharper looking peaks and possibly skipping past the real peak. Also in your second plot, it looks like the amplitude are much smaller than in the first one. $\endgroup$
    – Tyberius
    Dec 20 '21 at 21:29
  • $\begingroup$ @Tyberius in the second plot, the only thing I changed was the step size from $dx = 0.01$ to $dx = 0.08$ $\endgroup$
    – mkdirBed
    Dec 20 '21 at 21:30
  • $\begingroup$ That's quite a bit larger. Looking at your first plot as an example, if you use $dx=0.01$ you get roughly 20 points to model each peak. Switching to $dx=0.08$, you will only get 2 points per peak. This could drastically affect the accuracy of the result (as seems to be the case from your second plot). $\endgroup$
    – Tyberius
    Dec 20 '21 at 22:02
2
$\begingroup$

Your code has a small issue that is distorting the amplitude. You are not initializing p[1], which throws off the Verlet steps right from the beginning. If you have this set, the amplitude will stay roughly between 1 and -1.

import matplotlib.pyplot as plt    
import numpy as np    
     
e = 200    
                
dx = 0.08       
xlim = 100    
     
p = np.zeros(xlim+1)    
x = np.zeros(xlim+1)    
     
x[0] = 0     
p[0] = 1     
p[1]=p[0]+ (dx * dx) * (-e * p[0])    
     
     
for i in range(1, xlim):    
   p[i+1] = 2 * p[i] - p[i-1] + (dx * dx) * (-e * p[i])    
   x[i+1] = x[i] + dx    
     
     
plt.plot(x[1:], p[1:])    
plt.show()    

You still can't use a step size as large as $dx=0.08$ as this will lose too much of the detail of the curve, but you should be able to change the step size slightly now without fundamentally changing the result.

$\endgroup$
3
  • $\begingroup$ While this fixes the amplitude issue, it seems to be having an issue with the wave type. As it should be sinusoidal, y should be 0 when x = 0. In this case y seems to be approaching 1 $\endgroup$
    – mkdirBed
    Dec 20 '21 at 23:08
  • 1
    $\begingroup$ @mkdirBed You set the initial value p[0]=1. You could choose any initial condition including p[0]=0. cosine is just a shifted sine function and could also be considered sinusoidal. One thing I would recommend before solving an equation numerical is try to get a rough sense of what the solution should look like. In your case, the differential equation $y''=Ay$ is very well known and has an analytical solution. You should try to see if you can solve this to determine what the exact solution should be. $\endgroup$
    – Tyberius
    Dec 21 '21 at 2:15
  • 2
    $\begingroup$ I should clarify that you would also need to set a non-zero initial condition for the velocity if you were going to set p[0]=0. Your code implicitly has the initial velocity set to zero. $\endgroup$
    – Tyberius
    Dec 21 '21 at 2:32

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service, privacy policy and cookie policy

Not the answer you're looking for? Browse other questions tagged or ask your own question.