Лабораторная работа № 3-4 / Lab3_4new.doc

Лабораторная работа № 3-4

«Пакет MATHEMATICA»

Работа с таблицами, 2D -графикой, 3D -графикой, со специальными функциями. Изучение проходит на примере задачи рассеяния электрона на рассеивателе в виде сферического потенциального барьера.

Вдоль оси z слева направо на сферический рассеиватель падает поток электронов. Волновая функция электрона имеет вид плоской волны

(1)

где k -волновое число. Потенциальная энергия взаимодействия электрона с рассеивателем имеет вид

(2)

где a -радиус потенциального барьера (потенциального шара).

x ψ_s = f(θ)/r e^{i k r}

r

θ

ψ₀ = e^{i k z} z

V(r)

y

Уравнение Шредингера для стационарных состояний имеет вид

После преобразований это уравнение принимает вид

(3)

Решение уравнения (3) ищем с асимптотикой Зомерфельда

(4)

Дифференциальное уравнение (3) с асимптотикой (4) можно заменить эквивалентным интегральным уравнением

(5)

Здесь G(R) -функция Грина учитывает выбранную асимптотику.

Учитывая цилиндрическую симметрию задачи, решение интегрального уравнения (5) можно искать с помощью разложения по парциальным волнам.

(6)

где P_l(z) -полиномы Лежандра. В пакете Mathematica эта функция обозначается как LegendreP[l,z].

Неизвестная функция χ_l находится как решение следующего интегрального уравнения.

(7)

Здесь j_l(z) -сферическая функция Бесселя

(8)

где J_l(z) -функция Бесселя. В пакете Mathematica эта функция обозначается как BesselJ[l,z]. Интегральное ядро Γ(r, r′) в уравнении (7) определяется соотношениями

(9)

Здесь h_l(z) -сферическая функция Ханкеля

(8)

где H_l(z) -функция Ханкеля.

Лабораторная работа № 3

Разложение плоской волны по сферическим гармоникам.

Если в (7) положить потенциальную энергию равной нулю, т.е. убрать рассеиватель, то решением уравнения (7) будет

(9)

С другой стороны, в отсутствии рассеивателя, решением уравнения (3) является плоская волна (1). Поэтому, подставив (9) в разложение (6), получаем разложение плоской волны в ряд

(10)

Здесь θ - угол в сферической системе координат, поэтому . Сделаем замену

(11)

В формуле (11) есть две особенности, на которые стоит обратить внимание. Во-первых, ряд бесконечный, нам же для расчетов надо выбирать конечное число членов ряда L. Во-вторых, ряд (11) содержит произвольный параметр v, от значения которого величина экспоненты не должна зависеть. Такая независимость связана с тем, что ряд бесконечный, для конечного ряда это уже неверно. Так если взять конечное число членов ряда, то из (11) получим приближенное представление экспоненты.

(11a)

Левая часть соотношения (11a) зависит от одного параметра u , а правая от трех параметров u, v, L.

Заметим следующее, в задаче рассеяния бесконечный ряд (6) нам придется заменить конечной суммой аналогичной (11a). Конечная сумма, возникающая из ряда (6) будет также зависеть от этих трех параметров u, v, L. Возникает проблема - при заданных переменных u и v, сколько членов L суммы необходимо взять, чтобы бесконечный ряд вычислялся с заданной точностью. Если относительную ошибку вычисления ряда положить равной, например, δ = 0.01, то математически задача сведется к нахождению функции

(11b)

Если теперь в конечной сумме взять число членов ряда равным L ≥ Lr , то относительная ошибка вычисления бесконечного ряда будет не больше δ.

Так как в общем случае, эта задача требует много усилий и времени счета, сделаем разумное предположение, что если мы получим зависимость (11b) для плоской волны (11a), то этот результат можно перенести и на задачу рассеяния.

Задание 1.

Вычисление плоской волны с заданной точностью.

Создадим функцию

(12)

Эта функция является приближением экспоненты (11), и в пределе должна совпадать с ней

(13)

Напишем эту функцию в пакете Mathematica

Необходимо, для заданных значений u, v, |u| ≤ v, подобрать такое число членов ряда L, чтобы значение функции Fun[u, v, L] отличалось от точного значения не более чем на δ = 0.01. Заметим, что для комплексной экспоненты абсолютная и относительная погрешность совпадают. Точный результат можно найти, используя встроенную функцию Exp[I u]. Например

Для нахождения функции (11b) можно использовать следующий код

Эту функцию можно вызвать следующим образом

Первый пункт задания. С помощью функции для трех значений v = 0.1, 5.0, 10.0 нужно построить следующие графики.

Рис. 1, v = 0.1 Рис. 2, v = 5.0

Рис. 3, v = 10.0

Если для построения этих графиков использовать функцию Plot[], то это может занять слишком много времени. Удобней сначала создать двойной массив точек с помощью функции Table[]. Затем построить график с помощью функции ListPlot[]. Далее приведен пример получения графика синусоиды.

Рис. 4,

Второй пункт задания. Из рисунков 1, 2, 3 видно, что для фиксированного значения параметра v , функция Lr меняется в достаточно узком интервале, причем самое большое значение (наихудший вариант) принимает на границах интервала -v ≤ u ≤ +v. Поэтому мы упростим нашу задачу, будем искать зависимость Lr от u и v , предполагая, что u = v.

(13a)

Теперь построим график функции (13a).

Рис. 5,

Из рисунка 5 видно, что зависимость имеет в среднем линейный характер. Это имеет место и для больших значений параметра v. Например, на следующем рисунке область изменения v увеличена в 6 раз.

Рис. 6,

Поэтому возникает идея получить приближенную линейную формулу вычисления функции Lr, которую потом можно будет использовать для задания числа L членов ряда для заданного значения параметра v .

(14)

Задание будет состоять в нахождении чисел α, β, при которых прямая линия (14) наилучшим образом аппроксимирует кривую (13b), график которой приведен на рисунках 5 и 6.

Будем использовать минимизацию среднеквадратичного отклонения прямой линии (14) от заданной функции (13b). Выбираем N значений переменной v.

Вычисляем значения функции (13b) в этих точках

(14a)

Затем создаем функцию двух переменных α, β.

(14b)

Ищем экстремум функции

(14c)

Отсюда получаем линейную систему для нахождения коэффициентов α, β.

(14d)

Нужно создать программу для решения системы (14d) и найти коэффициенты α, β. Результат расчетов показан на следующем рисунке.

Рис. 7

Задание 2.

Графическое представление вычисленной плоской волны.

Научившись вычислять плоскую волну (10) с точностью δ = 0.01, построим графическое изображение этой волны. Так как плоская волна (10) обладает цилиндрической симметрией, то явно выделим две цилиндрические координаты z, ρ, а также волновое число электрона k.

(15)

Здесь коэффициенты α и β найдены в предыдущем задании, и будут использоваться во всех дальнейших заданиях в лабораторной работе 3 и лабораторной работе 4.

Обратим внимание на устранимую особенность в формуле (15) при kr→0. Из справочника специальных функций находим формулу для бесселевых функций при малых значения аргумента.

(15a)

Здесь Γ(p) - гамма функция. Используя формулу (15a) выясним поведение членов суммы (15) при kr→0.

(15b)

Из формулы (15b) видно, что все члены ряда (15) стремятся к нулю при kr→0, кроме члена с l=0. Этот член стремиться к единице.

(15c)

Отсюда следует, что ряд (15) равен единице при kr=0.

Напишем код функции (15) в пакете Mathematica, с учетом условия (15c), и с найденными значениями коэффициентов α и β.

Далее мы собираемся полученную функцию (15) представить в графическом виде. Для этого заметим, что функция зависит от двух пространственных переменных z и ρ. В цилиндрической системе координат радиус ρ лежит в пределах

Формально в силу симметрии в формуле (15) этот интервал можно расширить.

Это будет означать, что плоская волна рассматривается в точках некоторой плоскости проходящей через ось Oz , как это показано на рисунке 8.

Рис. 8

Построим зависимость действительной части функции (15) от переменных z и ρ с помощью функции Plot3D[].

Результат показан на рисунке 9.

Рис. 8

Изображение плоской волны получилось не слишком гладким. Это связано с тем, что функция Plot3D[]. строит 3D поверхность по умолчанию по 15 точкам вдоль каждой оси. Чтобы получить гладкую поверхность увеличим число точек вдоль оси z до 51 точки. Вдоль оси ρ число точек можно уменьшить до 5 точек, так как в этом направлении функция (15) почти не изменяется. Изменение числа точек можно выполнить с помощью опции PlotPoints→{51,5}.

Улучшенный результат показан на рисунке 9.

Рис. 9

Для анализа полученной 3D поверхности полезно построить линии уровня с помощью функции ContourPlot[]. Для примера построим поверхность для следующей функции

Применяем функцию Plot3D[].

Результат показан на рисунке 10.

Рис. 10

Применяем функцию ContourPlot[] для построения линий уровня.

Результат показан на рисунке 11.

Рис. 11

Задание состоит в следующем. Берем функцию Myexp[k, z, ρ], вычисляющую плоскую волну. Используя графические функции Plot3D[], и ContourPlot[].строим 3D поверхности и линии уровня для действительной части плоской волны в области -4 ≤ z ≤ 4 , -2 ≤ ρ ≤ 2 при различных волновых числах k = 0.1, 0.5, 1.0, 2.0, 5.0.

Кроме того, 3D поверхность на рисунке 9 надо рассечь плоскостью ρ = z/2, как это показано на рисунке 12.

Рис. 12

Затем необходимо изобразить полученную в сечении кривую с помощью функции Plot[].

Далее приведен пример расчета для k = 1.5.

Рис. 13

Лабораторная работа № 4

Приближение Борна в задаче рассеяния электрона

Если потенциал U₀ и радиус a рассеивающего центра имеют не слишком большие значения, то можно найти приближенное решение интегрального уравнения (7). Как показывается в теории квантового рассеяния, условия для приближенного решения имеют вид.

Если выполняются условия, приведенные выше, то для получения приближенного решения в интеграле (7) вместо искомой функции χ_l надо подставить j_l.

(16)

Такое приближение называют приближением Борна. Так как выбранная потенциальная энергия имеет простой вид

(17)

то интеграл в выражении (16) можно вычислить аналитически.

Подставляем (17) в формулу (16). В результате интеграл в (16) принимает вид.

(18)

Здесь для разных значений радиуса r возникают две формулы:

Для r ≥ a получаем следующую формулу

Используя явный вид сферической функции Ханкеля (8), получаем

(19)

Для r ≤ a получаем другую формулу

(20)

Здесь интегралы I₁, I₂ вычисляются по известным формулам, которые можно найти в пакете Mathematica. Для интеграла I₁ имеется простая формула

(21)

Чтобы вычислить интеграл I₂ надо немного потрудиться

(22)

Интеграл A вычисляется аналогично интегралу (21)

(23)

Для интеграла B получается следующее выражение

(24)

Здесь G(z) -гамма функция, в пакете она обозначается Gamma[], функция F(α, β, χ, z) называется гипергеометрическая функция и в пакете обозначается как

HypergeometricPFQ[{α}, {β, χ}, z]

В результате нам удалось получить аналитическое выражение для интеграла (18). Теперь, используя формулы (6), (7) можно найти волновую функцию электрона в задаче рассеяния на сферическом потенциальном барьере

(24)

Задание 1.

Ряд для волновой функции (24) очень похож на ряд для плоской волны. Например, при U₀ = 0 они совпадают.

Задание будет состоять в следующем. С помощью формулы (24) надо создать функцию Mypsi[a, U₀, k, z, ρ]. Надо использовать аналогичные выражения для плоской волны Myexp[k, z, ρ].

Первое, в ряде (24) вместо j_l(kr) надо поставить j_l(kr) - I(a, U₀, k, r).

Второе, число членов ряда L для случая r < a , надо выбирать по формул L = α (kr) + β, а для случая r > a, по формуле L = α (ka) + β.

Для полученной функции Mypsi[a, U₀, k, z, ρ] построить 3D поверхность в области

-a ≤ z ≤ 3a, -2a ≤ ρ ≤ 2a. Так как функция Mypsi[] комплексная, то интересно построить распределение плотности вероятности, т.е. надо строить поверхность для функции

Abs[Mypsi[]]^2. Далее надо построить линии уровня функции Abs[Mypsi[]]^2.

Параметры задачи a, U₀, k должны выбиться из условия применимости приближения Борна. Типичные значения параметров следующие

0 < U₀ < 0.1 (эВ)

0.1 < k < 5 (1/нм)

На следующем рисунке приведен пример для a = 2, U₀ = 0.02, k = 1.5.

Рис. 14

Красным кругом обозначено расположение центра рассеивания.

Ниже приводятся варианты выполнения лабораторной работы 4.