WWW.DISSERS.RU

БЕСПЛАТНАЯ ЭЛЕКТРОННАЯ БИБЛИОТЕКА

   Добро пожаловать!


Pages:     || 2 | 3 |
Саратовский Государственный Университет им. Н.Г. Чернышевского Кафедра теоретической и ядерной физики В. В. Серов Численные методы решения квантовомеханических задач Учебное пособие к курсу “Моделирование нелинейных процессов” для студентов специализации 010401 “Теоретическая физика” Саратов 2008 Оглавление Введение 3 1 Численное решение одномерной граничной задачи 6 1.1 Разностные схемы........................ 7 1.2 Прогонка............................. 9 Задание 1: Уравнение Пуассона................... 11 2 Решение временного уравнения Шредингера 12 2.1 Метод Кранка-Николсона.................... 14 Задание 2: Гармонический осциллятор во внешнем поле..... 15 2.2 Проблема нефизического отражения............. 17 Задание 3: Атомоподобная система во внешнем поле...... 20 а) Координатная калибровка.................. 20 б) Ускорительная калибровка................. 21 3 Эволюционные задачи, нелинейные по функции 22 3.1 Метод расщепления....................... 24 Задание 4: Солитоны......................... 26 4 Решение нелинейных стационарных задач 28 4.1 Непрерывный аналог метода Ньютона............ 29 Задание 5: Расчет стационарных состояний........... 32 а) Топологические моды одномерного Бозе-конденсата... 32 б) Возбужденные состояния атома гелия.......... 33 Приложения 35 А. Способы представления переменного внешнего поля..... 35 Б. Ортогональные полиномы..................... 37 Литература 38 2 Введение Содержанием данного курса является изучение методов моделирования нелинейных процессов в нерелятивистской квантовой механике. Под нелинейными процессами подразумеваются процессы, которые невозможно описать с помощью теории возмущений, т.е. нелинейные по внешнему полю. Актуальность таких задач связана с быстрым прогрессом лазерной технологии, позволяющей в настоящее время создавать сверхкороткие импульсы когерентного излучения, поле которых по интенсивности соизмеримо с полями внутри атомов. Кроме того, при приближенном описании больших квантовомеханических систем возникают уравнения, нелинейные по волновой функции, такие как уравнение Хартри–Фока или уравнение Гросса–Питаевского, описывающее динамику конденсата Бозе-Эйнштейна атомов. В отличии от стационарных задач и задач молекулярной динамики, для численного расчета нестационарных процессов нет стандартных программных комплексов, и программы обычно или создаются под конкретную задачу, или их приходится существенно модифицировать при применении к новой проблеме. Для этого необходимо знание базовых численных методов решения эволюционных уравнений.

Курс состоит из лабораторных работ, содержанием которых является написание программ, реализующих соответствующие алгоритмы решения, и сравнение результатов с тестовыми примерами. По причине практического характера курса, он ограничивается рассмотрением только тех методов, которые реализуемы от начала и до конца в течении разумного времени. Курс опирается на полученные ранее знания по математике (математический анализ, методы математической физики, численные методы) и физике (нерелятивистская квантовая механика [1], электродинамика [2]). В качестве основного языка программирования предлагается FreePascal версия языка Паскаль, имеющая свободно распространяемые реализации как под операционной системой Linux, так и под Windows. Помимо бесплатности, положительной стороной FreePascal является наличие модуля ucomplex для операций с комплексными числами, реализованного с помощью перегрузки арифметических операторов.

Так же допускается выполнение заданий на языках FORTRAN и С++.

Графическое представление результатов вычислений предполагается выполнять с помощью программы gnuplot.

Всего курс состоит из пяти лабораторных работ. Их порядок выбран таким образом, чтобы для выполнения каждой последующей работы использовались подпрограммы, написанные для предыдущих, что существенно ускоряет выполнение работ.

Поскольку в общефакультетских курсах численного моделирования на младших курсах методам решения граничных задач отводится слишком малое место, первая лабораторная посвящена изучению численных методов решения граничной задачи для обыкновенного дифференциального уравнения второго порядка. Далее рассматриваются методы решения временного уравнения Шредингера для осциллятора и для атомоподобной системы в сильном внешнем электрическом поле, а также временного и стационарного уравнения Гросса-Питаевского для бозе– конденсата как примера системы, нелинейной по волновой функции. Помимо изучения стандартных и широко распространенных численных методов, типа метода Кранка-Николсона и метода расщепления, в рамках курса рассматриваются методы, плохо отраженные до сего времени в учебной литературе, такие как непрерывный аналог метода Ньютона, или вообще малоизвестные в литературе на русском языке, типа метода внешнего комплексного скейлинга.

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

Данное учебное пособие создано при финансовой поддержке из средств совместного гранта Минобразования РНП2.2.2.3.415 и CRDF BRHE REC006 SR-006-X1/B75M06 Y3-P-06-08.

Глава Численное решение одномерной граничной задачи Цифровой компьютер не может решать дифференциальные уравнения напрямую, поэтому необходимо превратить их в алгебраические, точнее, свести задачу к решению системы линейных уравнений. Простейшим и наиболее “прозрачным” методом такого сведения является использование разностных схем[3, 4], которые мы рассмотрим ниже. Главным их плюсом является то, что они приводят к системам линейных уравнений с сильно разряженной (ленточной) матрицей, которые могут быть решены на компьютере весьма быстро, а минусом сложность корректного описания дифференциальных операторов нетривиального вида и граничных условий.



Альтернативой разностным схемам при решении дифференциальных уравнений является метод Галеркина[5], или попросту разложение функции по базису. В квантовой механике при расчетах часто используется разложение по глобальному базису(т.е. набору функций, неравных нулю в всей области определения искомой функции). Такой метод обычно приводит к системам уравнений с полностью заполненной матрицей, что компенсируется тем, что при удачном выборе базисных функций их количество, требующееся для достаточного приближения к точному решению, невелико. Промежуточным вариантом между сеточными(разностными) и “базисными” методами является использование локального базиса, строящегося на сетке, и приводящего к разряженной матрице. Наиболее часто используемые типы локальных базисных функций B-сплайны [6] и конечные элементы[7]. Оба этих типа функций строятся из кусков полиномов. Преимуществом B-сплайнов является то, что они дают максимальный порядок аппроксимации при минимальной ширине ленты матрицы. Одним из преимуществ конечных элементов является то, что коэффициенты разложения функции по ним совпадают со её значениями в узлах сетки, но главное их достоинство, приведшее к практической монополии метода конечных элементов при решении инженерных задач, является возможность постановки граничных условий на границах сложной формы в многомерном случае.

1.1 Разностные схемы Пусть у нас есть функция y(x), x [a..b]. Введем однородную сетку по x b - a xi = h(i - 1) + a; i = 1..N, h = (1.1) N - Обозначим yi = y(xi) и попытаемся выразить приближенное значение n-й производной в i-м узле сетки через значения функции в ближайших 2Nb + 1 узлах, j=i+Nb (n) yi = aijyj. (1.2) j=i-Nb т.е. заменить дифференцирование умножением на ленточную матрицу с полушириной ленты Nb. Разложим функцию в ряд Тейлора вблизи точки xi 1 1 y(x) = y(xi)+y(xi)(x-xi)+ y(xi)(x-xi)2+ y(xi)(x-xi)3+ y(4)(xi)(x-xi)4+...

2 6 (1.3) и подставим это разложение в (1.2). Чтобы получилось тождество, должны обнулится все члены разложения (1.3), кроме n-го порядка, что эквивалентно j=i+Nb aij(xj - xi)k = k!kn, k = 0,..., 2Nb;

j=i-Nb Это система 2Nb + 1 линейных уравнений относительно коэффициентов aij. Очевидно, что чтобы она имела однозначное решение(и вообще смысл), должно быть 2Nb n. Решая её для случая Nb = 1 (трехдиагональная матрица), получаем для первой и второй производных yi+1 - yi- yi = ; (1.4) 2h yi+1 - 2yi + yi- yi = ; (1.5) hЭти трехточечные формулы дают значения производных с точностью до членов порядка h3, что при решении уравнений дает точность порядка h2, поскольку число узлов сетки обратно пропорционально h. Заметим, что матрица для второй производной - симметрична, а для первой - антисимметрична. Это особенность сетки с равным шагом. Для сетки с непостоянным шагом полученные разностные операторы не будут обладать этим полезным свойством, и, если исходный дифференциальный оператор эрмитов, разностный оператор эрмитовым не будет, что является одним из главных недостатков разностных схем. Можно пользоваться симметричной разностной формулой [4] 2 yi+1 - yi yi - yi- yi = - ; (1.6) xi+1 - xi-1 xi+1 - xi xi - xi-но она для сеток с неравномерным шагом в общем случае даёт только первый порядок аппроксимации.

1.2 Прогонка Пусть нам нужно решить дифференциальное уравнение типа y(x) + (x)y(x) + (x)y(x) = f(x) (1.7) В электродинамике уравнение такого типа - это уравнение Пуассона в неоднородной среде, аналогичный вид имеют уравнения стационарной теплопередачи и диффузии. Оно должно быть дополнено граничными условиями y(a) + ay(a) = 0; y(b) + by(b) = 0. (1.8) Частными случаями являются граничные условия Дирихле y(a) = 0 и/или y(b) = 0.

и Неймана y(a) = 0 и/или y(b) = 0.

Пользуясь трехточечными разностными формулами, уравнение (1.7) превращается в систему aiyi+1 + biyi + ciyi-1 = fi (1.9) где 1 i 2 1 i ai = + ; bi = - + i; ci = -. (1.10) h2 2h h2 h2 2h Представим решение (1.9) в виде yi = iyi+1 + i (1.11) Переписав это выражение для yi-1, подставив в (1.9) и сравнив то, что получилось с выражением для yi (1.11), получаем реккурентные формулы ai fi - cii-i = - ; i =. (1.12) bi + cii-1 bi + cii-Граничные условия, в свою очередь, дают y1 - y0 yN+1 - yN + ay1 = 0; + byN = 0.

h h где мы ввели фиктивные узлы y0 и yN+1. Сравнивая с (1.11), получаем 0 = ah + 1; 0 = 0; (1.13) N yN =. (1.14) 1 + (bh - 1)N Таким образом, алгоритм прогонки следующий: мы берем 0, 0 из (1.13) и по формуле (1.12) вычисляем коэффициенты i, i, i = 1,..., N, затем берем yN из (1.14) и вычисляем решение yi, i = N - 1,..., 1 по формуле (1.11). В этой процедуре граничные условия как бы “прогоняются” сначала слева направо, а потом справа налево, откуда и произошло название метода в отечественной литературе. Метод прогонки полностью эквивалентен методу LU-декомпозиции для трехдиагональных матриц.

Задание 1: Уравнение Пуассона Написать программу, решающую уравнение (1.7) с граничными условиями (1.8) с помощью разностной схемы 2го порядка (1.4,1.5) и метода прогонки (1.13,1.12,1.14,1.11).

В качестве примера решить уравнение Пуассона для потенциала электронного облака в атоме водорода 1 e r2 = -4(r); (r) = exp(-2r/aB), r2 r r aB где aB боровский радиус, e заряд электрона, в атомной системе единиц aB = 1 и e = -1, с помощью замены (r) = y(r)/r приведя уравнение к виду y(r) = f(r); f(r) = -4r(r), с граничными условиями y(0) = 0; y(rmax aB) = 0.





Полученное решение сравнить с аналитическим e r + aB 2r (r) = 1 - exp -.

r aB aB Вывести ri, (ri), i = 1... N в файл и построить график с помощью программы gnuplot.

Зафиксировав верхнюю границу сетки b = 20, вычислить отклонение N = max |(ri) - calc| i i=1,...,N вычисленного решения calc от точного (ri) для разного числа узлов i сетки N = 1001, N = 2001, N = 4001 и определить порядок точности, т.е. показатель степени p в выражении N hp.

Глава Решение временного уравнения Шредингера Эволюция квантовой системы, состоящей из N частиц движущихся со скоростями много меньше скорости света, без учета спина описывается временным уравнением Шредингера (..., t) r1, rN, i = (..., t)(..., t). (2.1) r1, rN, r1, rN, t где гамильтониан системы(оператор полной энергии), постоянная Планка, волновая функция системы, |(..., t)|2 есть плотr1, rN, ность вероятности того, что в момент времени t первая частица будет находится в точке вторая частица в..., N-я частица в точке r1, r2, rN.

Точное численное решение такого уравнения даже на самых мощных современных компьютерах возможно только для квантовых систем, состоящих из небольшого числа частиц. Мы для простоты будем рассматривать одномерную квантовую систему. Одномерные задачи получаются при разделении переменных в многомерных уравнениях, возможном для некоторых задач, например о движении в центрально-симметричном потенциале. Одномерное временное уравнение Шредингера в атомных единицах измерения имеет вид (x, t) i = (x, t)(x, t). (2.2) t Атомные единицы определяются таким образом, что в них постоянная Планка, заряд электрона e и его масса me равны единице. Помимо уравнения движения, нужно еще задать начальное состояние системы (начальное условие) (x, 0) = 0(x).

Если в начальный момент времени система была связанной, то волновая функция должна удовлетворять граничным условиям (-, t) = 0; (, t) = 0.

Это уравнение параболического типа, но с мнимой единицей, что приводит к тому, что в отличии от, например, уравнения нестационарной теплопроводности, уравнение Шредингера имеет осциллирующие решения.

Формально решение (2.2) можно записать как t (x, t) = exp -i (x, t)dt (x, 0) (2.3) где множитель перед (x, 0) в правой части называется оператором эволюции или пропагатором. Заметим, что оператор эволюции унитарный, т.е. комплексно сопряженный к нему равен обратному. Это приводит к тому, что норма решения сохраняется |(x, t)|2 dx = 1 (2.4) Это условие можно рассматривать как необходимое условие устойчивости любого приближенного метода решения (2.2). Под устойчивостью численного метода подразумевается, что ошибка (отклонение приближенного решения от точного) растет по степенному, а не по экспоненциальному закону. Чтобы оно гарантированно удовлетворялось, необходимо, чтобы приближенный оператор эволюции также был унитарным.

Необходимо подчеркнуть, что устойчивость и точность (порядок аппроксимации) независимые свойства численной схемы, и высокая точность не означает устойчивости. Примером точного, но неустойчивого метода является метод Рунге–Кутта высокого порядка. Но устойчивость схемы при решении (2.2) куда более важна, чем точность.

2.1 Метод Кранка-Николсона Начнем с того, что выберем однородную сетку по времени с шагом. Из (2.3) следует n+1 = exp -i(tn+1/2) n (2.5) Мы заменили интеграл его приближением с помощью формулы прямоугольников. Здесь tn = n, n = (x, tn). Теперь мы должны выбрать приближение для экспоненты, поскольку напрямую экспоненту от оператора, т.е. матрицы, мы подсчитать не можем. Обычное разложение Тейлора не годится, поскольку приближенный оператор будет неунитарным. На практике это будет оборачиваться экспоненциальным ростом нормы со временем, т.е. неустойчивостью.

В качестве приближения для экспоненты выберем дробно-рациональную аппроксимацию Паде i I - (tn+1/2) exp -i(tn+1/2) = + O(3) (2.6) i I + (tn+1/2) где I единичный оператор. То, что это верно, нетрудно проверить с помощью сравнения разложений Тейлора экспоненты и дроби. Унитарность дробно-рационального оператора очевидна. Подставляя дробнорационального оператор в (2.5) вместо экспоненты и домножая обе части i выражения на I + H(tn+1/2), получаем i i I + (tn+1/2) n+1 = I - (tn+1/2) n (2.7) 2 Если гамильтониан дискретизирован, т.е. превращен в матрицу, например с помощью разностной схемы, процедура вычисления сводится к решению линейной системы уравнений типа (1.9) относительно вектора n+1. Такой метод называется схемой Кранка-Николсона[4].

Задание 2: Гармонический осциллятор во внешнем поле Написать программу, решающую уравнение (2.2) с помощью метода КранкаНиколсона (2.7), разностной схемы (1.5) и прогонки, для гамильтониана 1 (x, t) = - + U(x, t) 2 xВ качестве примера рассмотреть гармонический осциллятор с единичной частотой под действием периодической вынуждающей силы с частотой и амплитудой ExU(x, t) = - E0x sin t а начальное состояние вида 1 x(x, 0) = exp - ; (2.8) 2aa Осцилляторным потенциалом приближается потенциал взаимодействия атомов в молекуле вблизи положения равновесия, так что физически сформулированная выше задача есть задача о двухатомной полярной молекуле в электрическом поле линейно-поляризованной электромагнитной волны, а x есть отклонение расстояния между атомами от равновесного.

Pages:     || 2 | 3 |










© 2011 www.dissers.ru - «Бесплатная электронная библиотека»

Материалы этого сайта размещены для ознакомления, все права принадлежат их авторам.
Если Вы не согласны с тем, что Ваш материал размещён на этом сайте, пожалуйста, напишите нам, мы в течении 1-2 рабочих дней удалим его.