Schrodinger's Equation differential

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:

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

• Sure, a scaled sine function (e.g. $5\sin(x)$) still qualifies as a "sinusoidal oscillation". Dec 20 '21 at 20:26
• @Tyberius thanks for that. I have added an additional part to the question, does that mean that the algorithm is wrong? Dec 20 '21 at 21:18
• 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. Dec 20 '21 at 21:29
• @Tyberius in the second plot, the only thing I changed was the step size from $dx = 0.01$ to $dx = 0.08$ Dec 20 '21 at 21:30
• 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). Dec 20 '21 at 22:02

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.

• 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 Dec 20 '21 at 23:08
• @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. Dec 21 '21 at 2:15
• 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. Dec 21 '21 at 2:32