ВУЗ: Не указан
Категория: Не указан
Дисциплина: Не указана
Добавлен: 07.04.2021
Просмотров: 794
Скачиваний: 1
функции, полученные интегрированием «вперед» и «назад», также су-
щественно отличаются. Используя таким образом вычисленные волно-
вые функции, можно сформулировать критерий, позволяющий уточнить
собственное значение с наперед заданной точностью и, следовательно,
найти для него и соответствующую волновую функцию. Рассмотрим
алгоритм, реализующий эту идею.
Зададим на интервале
[
a, b
]
сетку из
N
узлов с постоянным шагом
h
= (
b
−
a
)
/
(
N
−
1)
:
x
n
=
a
+ (
n
−
1)
h,
n
= 1
,
2
,
3
, . . . , N.
(14)
Граничные условия (10) приобретают вид
ψ
1
=
ψ
N
= 0
.
(15)
где использовано обозначение
ψ
n
≡
ψ
(
x
n
)
. Задачи Коши для дифферен-
циального уравнения (11) будем решать методом Нумерова (см. При-
ложение 2). В рамках этого метода значение функции в узле сетки
находят, используя ее значения либо в двух предшествующих узлах
(интегрирование «вперед»):
ψ
n
+1
= [2(1
−
5
c q
n
)
ψ
n
−
(1 +
c q
n
−
1
)
ψ
n
−
1
] (1 +
c q
n
+1
)
−
1
,
(16)
либо в двух последующих узлах (интегрирование «назад»):
ψ
n
−
1
= [2(1
−
5
c q
n
)
ψ
n
−
(1 +
c q
n
+1
)
ψ
n
+1
] (1 +
c q
n
−
1
)
−
1
,
(17)
Здесь
c
=
h
2
/
12
,
q
n
=
q
(
E, x
n
)
. При использовании формулы (16) необ-
ходимо знать
ψ
1
и
ψ
2
, а формулы (17) -
ψ
N
−
1
и
ψ
N
. Значения
ψ
1
и
ψ
N
нам известны (15), а
ψ
2
и
ψ
N
−
1
- нет. Однако, если
N
достаточно ве-
лико, то для простоты можно считать, что
ψ
2
=
d
1
,
ψ
N
−
1
=
d
2
, где,
d
1
,
d
2
- малые числа (в силу непрерывности волновой функции). Как по-
казывает практика, величины
d
1
,
d
2
могут варьироваться в довольно
широких пределах. Заметим, что абсолютные значения этих величин
не имеют принципиального значения (в разумных пределах), так как
волновая функция в дальнейшем будет нормироваться.
Как уже упоминалось выше, для точного собственного значения опе-
ратора Гамильтона процедуры интегрирования «вперед» и «назад» при
11
подходящем выборе
d
1
,
d
2
должны приводить к одинаковым резуль-
татам (с точностью до постоянного множителя и ошибок округления).
Также очевидно, что при неравенстве
E
собственному значению, эти
процедуры дают различные функции, отношение которых в различных
узлах сетки не равно постоянной величине. Для оценки близости
E
к собственному значению будем вычислять разность производных вол-
новых функций, полученных интегрированием «вперед» и «назад», в
некотором внутреннем узле сетки
x
m
f
=
dψ
>
(
x
)
dx
x
m
−
dψ
<
(
x
)
dx
x
m
.
(18)
Здесь
ψ
>
,
ψ
<
– волновые функции, полученные интегрированием «впе-
ред» и «назад» соответственно;
x
m
– узел сшивки производных. Есте-
ственно, перед вычислением (18) необходимо масштабировать функции
ψ
>
и
ψ
<
так, чтобы
ψ
>
(
x
m
) =
ψ
<
(
x
m
)
. Такое масштабирование будем
называть
математической нормировкой
. Вычисление
f
позволяет ор-
ганизовать процедуру поиска собственного значения методом бисекции
(см. блок-схему на рис. (5)). Допустим, что необходимо найти энергию
и волновую функцию основного состояния. Выберем нулевое приближе-
ние к энергии основного состояния следующим образом:
E
(0)
=
U
min
+
δ
,
где
δ
малая величина (
δ >
0
), и вычислим соответствующее
f
(0)
по
формуле (18). Затем будем увеличивать энергию с шагом
∆
E
до тех
пор, пока величины
f
на двух соседних шагах
i
и
i
−
1
не будут иметь
различные знаки. Если шаг
∆
E
был меньше, чем разность энергий со-
седних уровней (собственных значений) в актуальном диапазоне изме-
нения
E
, то можно быть уверенным, что искомое собственное значение
∈
[
E
(
i
−
1)
, E
(
i
)
]
. Далее для уточнения собственного значения с наперед
заданной точностью
ǫ
используется метод бисекции.
Для расчета производных в формуле (18) можно воспользоваться
какой-либо подходящей численной процедурой; например, по формуле
dψ
dx
x
m
=
ψ
(
x
m
−
2
h
)
−
ψ
(
x
m
+ 2
h
) + 8[
ψ
(
x
m
+
h
)]
−
ψ
(
x
m
−
h
)]
12
h
.
(19)
12
да
нет
E
=
E
(0)
;
f
(0)
;
i
= 0
i
=
i
+ 1
E
(
i
)
=
E
(
i
−
1)
+ ∆
E
f
(
i
)
·
f
(
i
−
1)
<
0
МЕТОД
БИСЕКЦИИ
ВЫВОД
ВВОД
Рис. 5. Блок-схема процедуры поиска собственного значения.
1.1.3. Компьютерная программа
Ниже представлена программа численного решения одномерного
стационарного уравнения Шредингера для электрона в прямоугольной
потенциальной яме с шириной 4 а. е. длины и с бесконечными стен-
ками; дно ямы находится на оси ординат при -1 а. е. энергии (см. рис.
3) на языке СКМ Mathematica. Сравнение полученных по этой про-
грамме данных с результатами, найденными по аналитическим форму-
лам (8) и (9), позволяет оценить корректность алгоритма, описанного в
предыдущем параграфе. Использовались атомные единицы Хартри (см.
13
Приложение 1).
1
Clear["Global‘*"];
2
Off[General::spell];
3
Off[General::spell1];
4
L = 2.;
5
A = -L; B = +L;
6
n = 161;
7
h = N[(B - A)/(n - 1)];
8
c = h^2/12.;
9
Umin = -1.;
10
W = 3.;
11
U[x_Real] := If[Abs[x] < L, Umin, W];
12
q[e_Real, x_Real] := 2*(e - U[x]);
13
Psi = Fi = X = F = Psi2 = Table[0., {n}];
14
r = (n - 1)/2 + 15;
15
d1 = d2 = 10.^(-9);
16
Deriv[Y_List, h_Real, m_Integer] :=
17
(Y[[m - 2]] - Y[[m + 2]] +
18
8*(Y[[m + 1]] - Y[[m - 1]]))/(12.*h);
19
Num[e_Real, d1_Real, d2_Real, r_Integer, n_Integer] :=
20
Module[{i, k, big, coef, p1, p2, f1, f2, f, x},
21
If[n != Length[Psi] || n != Length[Fi] ||
22
n != Length[X], {Print["ERROR in Num: n=", n];
23
Return[{13., Psi, Fi, X}]}];
24
For[i = 1, i <= n, i = i + 1,
25
X[[i]] = x = A + h*(i - 1); F[[i]] = c*q[e, N[x]]];
26
Psi[[1]] = Fi[[n]] = 0.;
27
Psi[[2]] = d1; Fi[[n - 1]] = d2;
28
For[i = 2, i < n, i = i + 1,
29
p1 = 2*(1 - 5*F[[i]])*Psi[[i]];
30
p2 = (1 + F[[i - 1]])*Psi[[i - 1]];
31
Psi[[i + 1]] = (p1 - p2)/(1 + F[[i + 1]])];
32
For[i = n - 1, i > 1, i = i - 1,
33
f1 = 2*(1 - 5*F[[i]])*Fi[[i]];
34
f2 = (1 + F[[i + 1]])*Fi[[i + 1]];
35
Fi[[i - 1]] = (f1 - f2)/(1 + F[[i - 1]])];
36
p1 = Abs[Max[Psi]]; p2 = Abs[Min[Psi]];
37
big = If[p1 > p2, p1, p2];
38
For[i = 1, i <= n,
39
i = i + 1, Psi[[i]] = Psi[[i]]/big];
14
40
coef = Psi[[r]]/Fi[[r]];
41
For[i = 1, i <= n, i = i + 1,
42
Fi[[i]] = coef*Fi[[i]]];
43
f = Deriv[Psi, h, r] - Deriv[Fi, h, r];
44
Return[N[f]]];
45
Energy = Input["Energy=?"];
46
e = N[Energy];
47
f = Num[e, d1, d2, r, n];
48
AA = A - 1; BB = B + 1;
49
Gr0 = Plot[U[z], {z, A - 0.1, B + 0.1},
50
PlotStyle -> {AbsoluteThickness[2], RGBColor[0, 0, 1]},
51
PlotRange -> {-1.1, 2.9}, DisplayFunction -> Identity];
52
Point1[i_Integer] := {X[[i]], Psi[[i]]};
53
H1 = Array[Point1, n];
54
Style1 = {AbsoluteThickness[2], RGBColor[1, 0, 0]};
55
Gr1 = ListPlot[H1, {Frame -> True, PlotStyle -> Style1,
56
PlotJoined -> True}, DisplayFunction -> Identity];
57
Point2[i_Integer] := {X[[i]], Fi[[i]]};
58
H2 = Array[Point2, n];
59
Style2 = {AbsoluteThickness[2], RGBColor[0, 0.7, 1]};
60
Gr2 = ListPlot[H2, {Frame -> True, PlotStyle -> Style2,
61
PlotJoined -> True}, DisplayFunction -> Identity];
62
Gr3 = Graphics[{RGBColor[0, 1, 0], Disk[{X[[r]], Psi[[r]]},
63
0.07]}, DisplayFunction -> Identity];
64
Show[{Gr0, Gr1, Gr2, Gr3}, PlotRange -> {-1.1, 2.9},
65
AspectRatio -> Automatic, Frame -> True,
66
GridLines -> Automatic, FrameTicks -> Automatic,
67
PlotLabel -> "Wave functions Psi, Fi and potential",
68
FrameLabel -> {"X", "U(X),Psi(X),Fi(X)"},
69
DefaultFont -> {"Arial", 14}, Background ->
70
RGBColor[1, 1, 1], DisplayFunction -> $DisplayFunction];
71
t = "-----------------------------------------";
72
t1 = "
";
73
ue = PaddedForm[e, {10, 9}];
74
uf = PaddedForm[f, {10, 9}];
75
Print[t];
76
Print[t1, "e=", ue];
77
Print[t1, "f=", uf];
78
Print[t];
15