8.4: Приклад - проблема частинок в коробці
- Page ID
- 79569
Для демонстрації використання розріджених матриць для розв'язання скінченно-різницевих рівнянь розглянуто одномерну задачу частинки в коробці. Це складається з 1D незалежного від часу хвильового рівняння Шредінгера,
\[-\frac{1}{2} \, \frac{d^2\psi}{dx^2} = E\psi(x), \quad 0 \le x \le L,\]
разом з граничними умовами Діріхле
\[\psi(0) = \psi(L) = 0.\]
Аналітичний розчин нам добре відомий; до нормалізаційного фактора власні стани і енергетичні власні значення
\[\psi_m(x) = \sin\Big(m\pi x / L\Big),\quad E_m = \frac{1}{2}\, \left(\frac{m\pi}{L}\right)^2, \quad m = 1, 2, 3, \dots\]
Ми напишемо програму, яка шукає числове рішення. Використовуючи триточкове правило для дискретизації другої похідної, скінченно-різницеві матричні рівняння стають
\[-\frac{1}{2h^2}\begin{bmatrix} -2 & 1 \\ 1 & -2 & \ddots \\ & \ddots & \ddots & 1 \\ & & 1 & -2\end{bmatrix} \begin{bmatrix}\psi_0 \\ \psi_1 \\ \vdots \\ \psi_{N-1}\end{bmatrix} = E \begin{bmatrix}\psi_0 \\ \psi_1 \\ \vdots \\ \psi_{N-1}\end{bmatrix},\]
де
\[h = \frac{L}{N + 1}, \quad \psi_n = \psi\Big(x = h(n+1)\Big)\]
Наступна програма будує матричне рівняння скінченної різниці та відображає перші три розв'язки:
from scipy import *
import scipy.sparse as sp
import scipy.sparse.linalg as spl
import matplotlib.pyplot as plt
## Solve the 1D particle-in-a-box problem for box length L,
## using N discretization points. The parameter nev is the
## number of eigenvalues/eigenvectors to find. Return three
## arrays E, psi, and x. E stores the energy eigenvalues;
## psi stores the (non-normalized) eigenstates; and x stores
## the discretization points.
def particle_in_a_box(L, N, nev=3):
dx = L/(N+1.0)
x = linspace(dx, L-dx, N)
I = ones(N)
## Set up the finite-difference matrix.
H = sp.dia_matrix(([I, -2*I, I], [-1,0,1]), shape=(N,N))
H *= -0.5/(dx*dx)
## Find the lowest eigenvalues and eigenvectors.
E, psi = spl.eigsh(H, k=nev, sigma=-1.0)
return E, psi, x
def particle_in_a_box_demo():
E, psi, x = particle_in_a_box(1.0, 1000)
## Print the energy eigenvalues.
print(E)
## Plot |psi(x)|^2 vs x for each solution found.
fig = plt.figure()
axs = plt.subplot(1,1,1)
for n in range(len(E)):
fig_label = "State #" + str(n)
plt.plot(x, abs(psi[:,n])**2, label=fig_label, linewidth=2)
plt.xlabel('x')
plt.ylabel('|psi(x)|^2')
## Shrink the axis by 20%, so the legend can fit.
box = axs.get_position()
axs.set_position([box.x0, box.y0, box.width * 0.8, box.height])
## Print the legend.
plt.legend(loc='center left', bbox_to_anchor=(1, 0.5))
plt.show()
particle_in_a_box_demo()
Гамільтонова матриця, яка є тридіагоналлю, побудована у форматі розрідженої матриці DIA. власні значення та власні вектори знайдені з eigsh, що може бути використано, оскільки гамільтонова матриця, як відомо, є Ермітієвою. Зверніть увагу, що ми викликаємо eigsh за допомогою параметра sigma, вказуючи власному розв'язнику знайти власні значення, найближчі за значенням до\(-1.0\):
E, psi = spl.eigsh(H, k=nev, sigma=-1.0)
Це дозволить знайти найнижчі власні значення енергії, оскільки в цьому випадку всі власні значення енергії є позитивними. (Ми використовуємо\(-1.0\) замість 0.0, тому що алгоритм не працює добре, коли сигма дорівнює рівному нулю.) Якщо присутній негативний потенціал, нам доведеться знайти іншу оцінку для нижньої межі власних значень енергії, і передати її сигмі.
Крім того, ми могли б назвати eigsh з входом whic='SA '. Це скаже власному вирішувачу знайти власне значення з найменшим значенням. Ми уникаємо цього, оскільки алгоритм ARPACK власний розв'язувач відносно неефективний при знаходженні невеликих власних значень, у якому режимі (і іноді він може навіть не сходитися, якщо k занадто малий). Як правило, якщо у вас є можливість вивести нижню межу для потрібних власних значень, то краще використовувати сигму.
Запуск програми виводить найменші власні значення енергії:
[ 4.93479815 19.73914399 44.41289171]
Це добре узгоджується з аналітичними результатами\(E_1 = \pi^2/2 = 4.934802\)\(E_2 = 2\pi^2 = 19.739208\), і\(E_3=9\pi^2/2 = 44.413219\). Він також виробляє сюжет, показаний нижче, що так само, як і очікувалося.
