Karmanov_Trojanov

Проверка точности нахождения собственных векторов – расчет невязки системы уравнений:

Волновые функции первых состояний

Рассчитанные пробные функции (штриховые кривые), точные волновые функции (сплошные линии) и волновые функции основного и первого возбужденного состояний гармонического осциллятора (пунктир) показаны на рис. 2.

Точная волновая функция и ее нормировка:

ψ(n,x) := exp −y(x) y(x)β(n) mhyper(−n,2 β(n) + 1,y(x))

2

40

D(E) = H − E I = 0 ,

где I – единичная матрица 4-го порядка. Проведите графический анализ расположения корней и вычислите корни, используя встроенную функцию root().

6. Используя встроенную функцию lsolve(), найдите собственные векторы C из решения системы уравнений

j = 0

7. Считая известными собственные значения энергии, найдите собственные векторы матрицы H, используя метод квадратного корня (разложение Холецкого). Для выполнения разложения Холецкого можно воспользоваться встроенной функцией cholesky().

8.Рассчитайте собственные значения и собственные векторы для гамильтониана данной задачи итерационным методом, используя QR - разложение.

9.Проверьте нормировку и ортогональность рассчитанных волновых функций.

10.Используя точную волновую функцию Ψ(x)n для молекулы H2 в потенциале Морса, рассчитайте матричные элементы гамильтониана для первых четырех состояний. Дайте оценку обусловленности задачи вычисления собственных значений и собственных векторов гамильтониана H (см. [10]).

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

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

первых 4-х состояний молекулы H2 в потенциале вида (см. [9])

b) в потенциале Кратцера с параметрами, указанными выше, V(r) = −Vo ror 2 − ror .

42

Глава 4. Равновесное положение атомов ионного кристалла NaCl

(решение нелинейных уравнений)

Рассчитаем равновесное положение атомов и удельную энергию связи ионного кристалла NaCl.

Потенциал взаимодействия

Будем считать, что потенциал взаимодействия пары ионов решетки кристалла может быть представлен как сумма кулоновского потенциала и центрального потенциала отталкивания, записанного в форме БорнаМайера [13] :

где λ и ρ − эмпирические параметры потенциала Борна-Майера, ri,j – расстояние между ионами и qi , qj – их заряды.

Поскольку взаимное отталкивание является короткодействующим, то его вклад при суммировании достаточно учесть только для ближайших (соседних) атомов и тогда полная энергия кристалла представляется в виде

43

Используя режим Format-> Graph-> Trace получаем оценку равновесного положения и удельной энергии кристалла:

Уравнение, определяющее равновесное положение кристалла,

Метод половинного деления

Алгоритм метода половинного деления [11] реализован в виде подпрограммы Rmin(f,a,b), представленной ниже. Параметры этой подпрограммы: f – функция, задающая уравнение; a и b – начальные границы интервала, на котором ищется корень. Программа проверяет наличие одного действительного корня на заданном отрезке. Если корней нет или их четное число на данном отрезке, то программа выдает сообщение "No Root". Точность расчета корня eps

характеризуется шириной интервала локализации корня.

44

c ← a +2 b fc ← f(c)

b ← c if fa fc < 0

otherwise a ← c

fa ← fc rm ← a +2 b

Поскольку в данном случае корень находится на интервале [2,3], то, используя подпрограмму Rmin, можно вычислить корень урав-

нения f(R) = 0. Используя пункт меню Format-> Result-> Number Format-> Number of decimal places, установим режим представления результатов, например, с 5-ю десятичными знаками после запятой, для того, чтобы можно было следить за динамикой изменения погрешности разных методов при выбранном уровне точности eps = 0.001.

Rmin(f,2,3) = 2.81543

Оценка невязки для уравнения f(R) = 0 :

f(RmRm:=) =Rmin(f,2,3) f(Rm) = −1.21358 × 10− 6

а, например, для интервала [3,5] получаем ответ, свидетельствующий об отсутствии корня:

Rmin(f,3,5) = "No Root"

45

Метод Ньютона

Метод Ньютона [11] реализован в виде подпрограммы– функции newton (f,df,Rold,eps). Параметры этой подпрограммы:

f – функция, представляющая уравнение, df – функция, задающая производную функции f, Rold – начальное приближение к корню и eps – точность вычисления корня.

Интересно отметить, что, несмотря на одинаковые требования к точности вычисления корней (eps = 0.001) в методе половинного деления и методе Ньютона, невязка уравнения f(R) = 0 в данном случае много меньше.

f(newton(f,df ,2.0,0.001)) = 8.19403 × 10− 10

или f(Rn) = 8.19403 × 10− 10

Определение корней с помощью функции root(y,x)

Прежде всего, проверим используемое по умолчанию значение параметра TOL , определяющего "точность" вычисления корня,

46

Ясно, что требуемая точность вычисления корня не достигается и параметр TOL задает ограничение на точность, с которой выражение f(r) может считаться равным нулю!!!

Теперь результаты расчета корней тремя способами совпадают!

Итерационный процесс

Проследим за динамикой изменения положения корня и невязки в ходе итерационного процесса. Для этого модифицируем программу вычисления корня по методу Ньютона и построим графическую иллюстрацию к итерационному процессу.

47

Упражнения

1.Рассчитайте постоянную Маделунга кристалла NaCl (см. [13–15]).

2.Выполните суммирование энергий взаимодействия пар ионов по решетке NaCl и получите выражение Utot(R) для удельной энер-

гии кристалла (см. [13–15]).

3. Постройте графическую иллюстрацию итерационного процесса нахождения корня методом половинного деления. Оцените скорость сходимости.

4. На основе вычислительного эксперимента оцените границы области сходимости метода Ньютона в данной задаче. Сопоставьте Ваши результаты с достаточным условием

48

5. В данной задаче на основе вычислительного эксперимента оцените границы области монотоннойсходимости метода Ньютона и скорость его сходимости.

6.Напишите программу и найдите методом секущих ([10–12]) равновесные параметры кристалла NaCl. Постройте графическую иллюстрацию итерационного процесса.

7.На основе вычислительного эксперимента оцените границы области сходимости, скорость сходимости метода секущих в данной задаче. Сопоставьте Ваши результаты с результатами упражнений 4–5.

8.Предполагая, что атомы инертных газов взаимодействуют между собой посредством потенциала Ленарда-Джонса [13–15],

получите выражение для удельной энергии кристалла с гранецентрированной (ГЦК) решеткой

Здесь С6 , С12 – решеточные суммы для ГЦК – решетки, а ε и σ – параметры потенциала Ленарда-Джонса.

9. Рассчитайте решеточные суммы С6 , С12 для ГЦК– решетки кристаллов инертных газов (см. [13–15]).

10. Принимая следующие значения параметров потенциала ЛенардаДжонса для атомов Ne

рассчитайте любым из методов, обсуждавшихся в данном разделе, равновесное положение, удельную энергию и модуль всестороннего сжатия B для ГЦК решетки кристалла Ne, если

49