close

Вход

Забыли?

вход по аккаунту

?

9894.Математическое моделирование

код для вставкиСкачать
Copyright ОАО «ЦКБ «БИБКОМ» & ООО «Aгентство Kнига-Cервис»
Министерство образования и науки Российской Федерации
Ярославский государственный университет им. П. Г. Демидова
С. Е. Биркган
Математическое
моделирование
Учебное пособие
Рекомендовано
Научно-методическим советом университета
для студентов, обучающихся по направлению
Электроника и наноэлектроника
Ярославль 2012
Copyright ОАО «ЦКБ «БИБКОМ» & ООО «Aгентство Kнига-Cервис»
УДК 519.8 (075.8)
ББК В183.5я73
Б 64
Рекомендовано
Редакционно-издательским советом университета
в качестве учебного издания. План 2012 года
Рецензенты:
кафедра математического анализа Ярославского государственного
педагогического университета им. К. Д. Ушинского;
Проказников А. В., д-р физ.-мат. наук, ведущий научный сотрудник ИМИ РАН
Б 64
Биркган, С. Е. Математическое моделирование: учебное
пособие / С. Е. Биркган; Яросл. гос. ун-т им. П. Г. Демидова. –
Ярославль: ЯрГУ, 2012. – 92 с.
ISBN 978-5-8397-0906-5
Учебное пособие посвящено моделированию взаимодействия
ионов с поверхностью – одной из фундаментальных проблем физики, микро- и наноэлектроники, и моделям цифровой обработки
сигналов, используемым в современной электронике. Первая из
этих проблем рассматривается с точки зрения изменения морфологии поверхности при ионном облучении и приводит к нелинейным дифференциальным уравнениям в частных производных.
Вторая – с точки зрения теории цифровой фильтрации – основной
составляющей цифровой обработки сигналов, приводящей к дискретным разностным уравнениям. Рассмотрены классические и
относительно новые результаты.
Пособие предназначено для студентов, обучающихся по направлению 210100.68 Электроника и наноэлектроника (дисциплины «Методы математического моделирования» и «Математическое моделирование в наноэлектронике», блок М1), очной формы
обучения.
УДК 519.8 (075.8)
ББК В183.5я73
ISBN 978-5-8397-0906-5
© Ярославский государственный университет
им. П. Г. Демидова, 2012
Copyright ОАО «ЦКБ «БИБКОМ» & ООО «Aгентство Kнига-Cервис»
Оглавление
Введение . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 5
1. Моделирование в задачах микро- и наноэлектроники . . . . . . . . . . . . . . . 7
1.1. Моделирование взаимодействия ионов с поверхностью . . . . . . . . . . . . 7
1.1.1. Уравнение Больцмана . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 10
1.1.2. Модели динамики рельефа поверхности, облучаемой ионами . . . 12
1.1.2.1. Общие закономерности выхода распыленных атомов . . . . . . 13
1.1.2.2. Зависимость коэффициента распыления от энергии . . . . . . . 16
1.1.2.3. Угловая зависимость выхода распыленных атомов . . . . . . . . 19
1.1.2.4. Распылительная модель и ее дифференциальные формы . . . 21
1.1.2.5. Уравнение Бредли-Харпера . . . . . . . . . . . . . . . . . . . . . . . . . . . 26
1.2. Исследование моделей динамики рельефа поверхностей . . . . . . . . . . . 30
1.2.1. Исследование распылительных моделей . . . . . . . . . . . . . . . . . . . . 31
1.2.1.1. Вывод основных расчетных формул . . . . . . . . . . . . . . . . . . . . 31
1.2.1.2. Компьютерное моделирование и визуализация результатов 33
1.2.1.3. 3D визуализация результатов моделирования . . . . . . . . . . . . 38
1.2.2. Возникновение “ripple” структур в модели Бредли-Харпера . . . . 39
Литература . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 43
2. Моделирование в задачах цифровой обработки сигналов . . . . . . . . . .
2.1. Моделирование цифровых устройств разностными уравнениями . . .
2.2. Начальные понятия теории дискретных разностных уравнений . . . . .
2.3. Разностные уравнения первого порядка . . . . . . . . . . . . . . . . . . . . . . . . .
2.4. Разностные уравнения второго порядка . . . . . . . . . . . . . . . . . . . . . . . . .
2.5. Методы решения разностных уравнений . . . . . . . . . . . . . . . . . . . . . . . .
2.5.1. Оценка скорости роста решений . . . . . . . . . . . . . . . . . . . . . . . . . . .
2.5.2. Применение Z-преобразования . . . . . . . . . . . .. . . . . . . . . . . . . . . . .
2.5.3. Алгебраический метод решения однородных уравнений . . . . . . .
2.5.4. Метод вариации произвольной постоянной . . . . . . . . . . . . . . . . .
2.5.5. Разностные уравнения с правой частью в виде квазиполинома
2.5.6. Применение функции Грина . . . . . . . . . . . . . . . . . . . . . . . . . . . . .
2.6. Устойчивость решений разностных уравнений . . . . . . . . . . . . . . . . . . .
2.6.1. Устойчивость уравнений с постоянными коэффициентами . . . .
2.6.2. Алгебраические методы анализа устойчивости . . . . . . . . . . . . . .
2.6.3. Частотные методы анализа устойчивости . . . . . . . . . . . . . . . . . . . .
2.6.4. Теорема Флоке-Ляпунова и устойчивость решений . . . . . . . . . . .
2.6.5. Устойчивость решений уравнений с малым параметром . . . . . . .
3
45
45
47
50
52
55
55
56
60
62
65
67
69
71
71
73
74
76
Copyright ОАО «ЦКБ «БИБКОМ» & ООО «Aгентство Kнига-Cервис»
2.6.6. Теорема об устойчивости по первому приближению . . . . . . . . . .
2.6.7. Устойчивость уравнений с гармоническими коэффициентами . .
2.6.7.1. Периодический случай . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .
2.6.7.2. Почти периодический случай . . . . . . . . . . . . . . . . . . . . . . . . .
2.7. Нелинейные разностные уравнения и хаос . . . . . . . . . . . . . . . . . . . . . . .
2.7.1. Логистическое уравнение . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .
2.7.2. Состояния равновесия логистического уравнения . . . . . . . . . . . .
2.7.3. Бифуркации удвоения периода и переход к хаосу . . . . . . . . . . . .
Литература . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .
4
78
79
79
84
85
86
87
87
90
Copyright ОАО «ЦКБ «БИБКОМ» & ООО «Aгентство Kнига-Cервис»
Введение
Идея замены реального физического объекта (явления, процесса) его математической, а затем, при необходимости, и компьютерной моделью лежит в основе методологии математического моделирования. Этот другой метод познания сочетает в себе многие достоинства теоретического и экспериментального
исследования. Современную науку невозможно представить себе без широкого
применения методов математического моделирования, которое дает в руки специалиста мощный, гибкий, удобный и универсальный инструмент для изучения
самых различных явлений, использующий новейшие физические представления, математические теории и достижения вычислительной техники. При этом
использование современных вычислительных методов и высочайшей производительности современных ЭВМ позволяет глубоко и полно изучать исследуемые объекты и получать результаты, недостижимые при чисто теоретическом
подходе, применяемом в недавнем прошлом. Особое значение математическое
и компьютерное моделирование приобретает при отработке новых технологий
в производстве, когда необходимо проведение массовых экспериментов для оптимизации того или иного технологического процесса при выборе наилучших
параметров и режимов процессов и достижения нужных характеристик выпускаемых изделий. Такие физические эксперименты могут быть весьма затратными и ведущими к серьезному удорожанию конечного продукта, снижающему
его конкурентоспособность. Разумеется, следует учитывать, что замена физического эксперимента компьютерным возможна лишь при наличии адекватных
физической и математической моделей.
Задача математического моделирования объекта или явления условно может быть разбита на три части.
• Математическая формулировка свойств и законов, которым подчиняется
объект или явление, и их описание на языке математики с помощью уравнений.
Этот этап называется выбором математической модели.
• Выбор адекватных аналитических и численных методов и алгоритмов
решения построенной модели, не искажающих ее основные свойства. Эти алгоритмы должны быть реализуемыми, экономичными и адаптированными к используемым ЭВМ.
• Разработка программы – компьютерной модели изучаемого объекта, которая, помимо адекватности и экономичности, должна иметь удобный графический интерфейс, позволяющий работать с ней конечному пользователю.
После проверки адекватности всех частей компьютерной модели исследователь получает в свое распоряжение универсальный, гибкий, мощный и недорогой инструмент, который может быть использован вместо трудоемких натурных экспериментов. При этом в процессе выполнения компьютерных экспериментов сама модель может изменяться и уточняться во всех своих частях.
Из сказанного выше ясно, что математическое моделирование не подменяет
собой математику, физику, биологию и другие научные дисциплины, а, наобо5
Copyright ОАО «ЦКБ «БИБКОМ» & ООО «Aгентство Kнига-Cервис»
рот, синтезирует и дает новые методы и подходы для их дальнейшего развития.
Следует отметить, что моделирование в определенном смысле присутствует
практически во всех видах творческой деятельности человека, однако использование математических методов позволяет подкрепить или опровергнуть поверхностные умозрительные заключения точными расчетами. Сфера применения математического моделирования чрезвычайно широка и существует много
книг и учебников, посвященных моделированию различных конкретных задач
физики, биологии, социологии, техники и т. д. При этом в силу огромного числа решаемых задач и используемых методов невозможно (да и не имеет смысла) в одном издании с достаточной полнотой отразить все известные модели.
Первая часть данного учебного пособия посвящена моделированию взаимодействия ионов с поверхностью – проблеме, являющейся одной из фундаментальных в физике и наноэлектронике. Во второй части пособия изучаются модели цифровой обработки сигналов, используемые в современной электронике.
Приведены классические и относительно новые результаты.
6
Copyright ОАО «ЦКБ «БИБКОМ» & ООО «Aгентство Kнига-Cервис»
1. Моделирование в задачах микро- и наноэлектроники
Одной из важнейших задач современной наноэлектроники является разработка новых технологий для создания элементов с характерными размерами
порядка 10 нм и менее. Такие приборы будут востребованы уже в самом ближайшем будущем для систем обработки информации, различного оборудования
электронного зондирования, источников СВЧ излучения и многих других целей. С недавнего времени ориентация на квантовые принципы работы рассматривается как основной путь для промышленных разработок сверхкомпактных
электронных устройств, элементами которых станут как нанотранзисторы и
элементы памяти, так и новые элементы, работа которых основана на этих
принципах. Результатом применения новейших наноматериалов и нанотехнологий станут сверхбыстрые и суперкомпактные компьютерные и коммуникационные системы терагерцового диапазона с минимальным потреблением
энергии, переводящие современную информатизацию на новый уровень.
Создание электронных приборов нового поколения немыслимо без всестороннего изучения свойств наноматериалов и структур и разработки способов и
технологий их изготовления. Данная задача распределяется на экспериментальные и теоретические исследования, причем роль последних постоянно возрастает. Это обусловлено рядом обстоятельств и в первую очередь значительно
выросшей мощью компьютерных систем, способных брать на себя проведение
огромной массы вычислительной работы, которая во многих случаях может заменить сложные, дорогие и долгие физические эксперименты. Следует отметить, что сами математические модели, используемые в наноэлектронике, становятся также все более сложными. Это связано еще и с тем, что наряду с традиционными факторами, такими как трехмерность, нестационарность, нелинейность, пространственно-временная неустойчивость и некорректность задач
появились новые особенности связанные со сложной геометрией, иерархией
пространственно-временных размеров, нелокальностью процессов и наличия
большого числа компонент. В итоге во многих новых математических моделях
наноэлектроники используются физические описания, для формализации которых требуется разный математический аппарат, и стыковка разных частей моделей вызывает дополнительные трудности. Еще одной весьма важной особенностью современных моделей наноэлектроники является сложность, а во многих случаях и невозможность анализа полученных расчетных данных без их
эффективной визуализации. Действительно, трудно представить себе, например, анализ изменения поверхности тела, распыляемого ионным лучом, с помощью таблиц из координат точек поверхности для разных значений времени.
Ясно, что для решения такого сорта задач требуется применение современных
методов визуализации, включая трехмерную и стерео визуализации, которые
позволят исследователю “погрузиться” в исследуемый процесс или явление.
При этом часть компьютерной модели, обеспечивающая визуализацию, может
быть весьма эффективно использована как для выбора параметров путем сравнения результатов расчетов с экспериментальными, так и для существенного
изменения самой модели и сравнения конкурирующих моделей.
7
Copyright ОАО «ЦКБ «БИБКОМ» & ООО «Aгентство Kнига-Cервис»
Основные задачи современной наноэлектроники, решаемые математическими и компьютерными методами, можно разделить на несколько классов. К
первому классу нужно отнести проблемы, связанные с расчетом состава и свойствами используемых материалов. Объектом исследований здесь являются
квантово-механические задачи, сформулированные для систем большой размерности, моделирующие физико-химические свойства используемых веществ
и материалов. Второй класс составляют задачи, связанные с созданием структур на поверхности материалов, предназначенных для будущих электронных
компонент. Здесь моделируются различные процессы создания поверхностной
топологии будущих электронных устройств. К третьему, традиционному, классу для электроники в целом относятся задачи, связанные с разработкой различных рисунков и схем, реализующих дизайн, отдельные компоненты и узлы устройств. Четвертый класс составляют задачи моделирования электрических процессов, происходящих в отдельных компонентах, узлах и в устройстве в целом,
обеспечивающие нужные свойства электронного изделия за счет оптимизации
режимов входящих в него элементов.
Трудно оценить количество научных публикаций, посвященных моделированию различных задач наноэлектроники, и этот список постоянно растет. В
данной части работы будут рассмотрены лишь некоторые модели развития
структур на поверхности материалов под воздействием ионной бомбардировки.
Эта проблема, относящаяся ко второму классу задач, описанных выше, на самом деле представляет собой часть одной из фундаментальных физических
проблем взаимодействия ионов с поверхностью.
1.1. Моделирование взаимодействия ионов с поверхностью
Исследование взаимодействия ионов с поверхностью материалов представляет большой научный и прикладной интерес. Впервые, по-видимому, на эту
проблему обратили внимание после открытия явления распыления в 1853 году
при исследовании электрических разрядов в газах. Проводя экспериментальные
исследования, Гроув [1] заметил вредный паразитный эффект загрязнения стенок разрядной трубки частицами материала катода. Этот эффект посчитали
вначале следствием испарения катода. Однако после проведения остроумного
опыта, когда в катоде было проделано отверстие, а за ним установлена металлическая мишень и после электрического разряда на стенках стеклянной трубки
за катодом были обнаружены частицы материала мишени, стало ясно, что дело
в другом. В качестве нового объяснения было выдвинуто предположение распыления материала мишени “каналовыми лучами”. В дальнейшем было установлено, что на самом деле эти “каналовые лучи” представляют собой поток
положительно заряженных ионов, которые образуются в газовом разряде и
притягиваются электрическим полем к катоду. Налетая на катод или мишень,
эти ускоренные ионы выбивают частицы материала, которые и оседают затем
на стенках разрядной трубки. Следует отметить, что, несмотря на отдельные
прорывы в исследовании эффекта распыления, решающие успехи в понимании
физики процесса были достигнуты лишь примерно через сто лет после открытия данного явления. Это объясняется несколькими причинами. Во-первых,
8
Copyright ОАО «ЦКБ «БИБКОМ» & ООО «Aгентство Kнига-Cервис»
только к этому времени серьезное развитие экспериментальной техники позволило получить обширный и достоверный экспериментальный материал для
кристаллических, поликристаллических и аморфных мишеней, что дало возможность составить достаточно полную картину изучаемого явления. Вовторых, распыление перестало быть узкой технической проблемой со своими
порой примитивными теоретическими концепциями и естественным образом
вписалось в круг задач радиационной физики. Теории распыления, созданные в
этот период, уже рассматривали процесс распыления как один из аспектов радиационного разрушения вещества при его облучении атомными частицами и
использовали математический аппарат, разработанный в радиационной физике.
В-третьих, в 1960 г. в Брукхейвенской лаборатории (США) было проведено успешное моделирование процессов, происходящих в твердых телах, подвергнутых облучению с помощью ЭВМ. При этом совпадение расчетных результатов
с экспериментальными было настолько хорошим, что это дало основание руководителю этих работ Дж. Виньярду [2] заявить о появлении в дополнение к
теоретической и экспериментальной “третьего вида (estate) физики”. Компьютерные эксперименты, которые в дальнейшем стали неотъемлемой частью исследований, позволили как бы проникнуть внутрь изучаемого явления. Вчетвертых, в шестидесятых годах были опубликованы фундаментальные теоретические исследования Линдхарда и его сотрудников [3-6], посвященные прохождению атомных частиц через конденсированные среды, и работы Зигмунда
[7,8], заложившие основу современной теории распыления.
Наиболее общей математической моделью процесса ионного распыления
поверхности в настоящее время является уравнение переноса, построенное на
теории линейных каскадов развитой в работах Зигмунда [7,8]. Однако ввиду
сложности исследование этого интегро-дифференциального уравнения проводят при дополнительных упрощающих предположениях, позволяющих свести
его к более простым дифференциальным уравнениям. При этом для описания
эрозии поверхности при ионном облучении в простых случаях успешно используются распылительные модели [9,10]. Появление сложной морфологии на распыляемой поверхности часто связывают со сложными математическими моделями, такими как уравнения Бредли-Харпера [11], Курамото-Сивашинского [12]
и Кадара-Паризи-Чанга [13]. Так, уравнение Бредли-Харпера было использовано
в [11] для обоснования возникновения волнообразных (“ripple”) структур. Однако в дальнейшем стало понятно, что многие эффекты появления волнообразного
рельефа (не говоря уже о рельефах других типов) не могут быть описаны в рамках этого уравнения (см., например, работу [12]). С целью устранения недостатков была предложена модификация модели Бредли-Харпера – уравнение КадараПаризи-Чанга [13], значительно усложняющее модель и отличающееся дополнительными нелинейными слагаемыми, введенными для ограничения роста периодических решений. Это уравнение, а также близкое к нему по духу уравнение
Курамото-Сивашинского в определенной мере сняли некоторые противоречия
модели Бредли-Харпера, но не смогли объяснить вновь обнаруженные эффекты,
в частности экспериментально наблюдаемое возникновение волнообразного
рельефа практически при любых наклонах направления ионного пучка к поверх9
Copyright ОАО «ЦКБ «БИБКОМ» & ООО «Aгентство Kнига-Cервис»
ности мишени [14]. Из других существенных недостатков упоминаемых уравнений следует отметить большое количество трудно определяемых физических параметров, вариации значений которых могут достигать десятичных порядков.
Кроме этого, следует отметить, что экспериментально наблюдаемый периодический рельеф в большинстве случаев имеет гребнеобразную, негладкую форму,
содержащую углы (кромки). Во избежание дополнительных проблем для описания таких режимов более естественным является использование дифференциальных уравнений, не содержащих высших производных, тем более что при
компьютерном моделировании, использующем численные методы решения
дифференциальных уравнений, наличие старших производных может также существенно усложнить задачу.
Таким образом, несмотря на существование достаточно большого количества
различных уравнений, описывающих процесс развития поверхности под воздействием ионного распыления, построение относительно простой универсальной и
адекватной модели, допускающей применение методов численного и компьютерного моделирования, по-прежнему является весьма актуальной задачей теории взаимодействия ионов с поверхностью.
1.1.1. Уравнение Больцмана. Как уже было сказано, для описания процесса распыления поверхности под воздействием ионного облучения в настоящее
время используется несколько моделей, которые можно рассматривать как
следствия уравнения переноса. Это уравнение, полученное Больцманом, возникает во многих разделах физики и его часто называют еще уравнением
транспорта. В работах Зигмунда [7,8] уравнение переноса было использовано
для описания взаимодействия атомов мишени с падающими на нее ионами.
Важными результатами, полученными на основе анализа этого уравнения, явились угловые и энергетические зависимости выхода распыленных атомов, их
дифференциальные угловые распределения и ряд других параметров. Главной
целью работ [7,8] являлось построение достаточно общей теории ионного распыления, объясняющей известные экспериментальные результаты и позволяющей прогнозировать новые эффекты и явления, возникающие при ионной бомбардировке поверхности. Вывод основных уравнений в [7] базируется на теории линейных каскадов, суть которой состоит в том, что падающий ион выбивает атомы мишени из состояния равновесия и приводит их в движение. Выбитые атомы движутся в веществе, сталкиваются с другими атомами и при определенных условиях некоторые из них, преодолевая поверхностный барьер, покидают мишень. Такой механизм распыления хотя и не является единственным,
но считается наиболее универсальным. При взаимодействии движущихся ионов
с атомами мишени обычно выделяют три разных случая: режим первичного
выбивания атомов, режим линейных каскадов и режим тепловых пиков. В режиме первичного выбивания бомбардирующий ион передает энергию атомам
мишени, которые могут выйти через поверхность после небольшого числа последующих столкновений при условии, что их скорость достаточна для преодоления поверхностных сил связи. Во втором и третьем случаях выбиваемые атомы приобретают большую энергию, которой достаточно для выбивания вторичных атомов. При этом определенная часть этих атомов достигает поверхно10
Copyright ОАО «ЦКБ «БИБКОМ» & ООО «Aгентство Kнига-Cервис»
сти и, преодолевая поверхностный барьер, покидает мишень. Режим линейных
каскадов характеризуется тем, что число движущихся атомов относительно мало, а для режима тепловых пиков это число велико. Отметим, что подавляющее
большинство атомов, покинувших мишень, имеют сравнительно небольшую
энергию. Атомы мишени, находящиеся вблизи траекторий первичных ионов и
выбитых атомов с высокой энергией, приобретая часть энергии быстро движущихся частиц, имеют небольшую длину свободного пробега. Следствием этого
является то, что они могут покинуть мишень только в случае, если изначально
располагались в пределах нескольких атомных слоев от поверхности. Но, с другой стороны, таких низкоэнергетичных атомов много, в связи с чем они дают
наибольший вклад в число распыленных частиц. Определение числа эмитированных частиц в работе [7] проводится в четыре этапа:
• вычисление энергии, вносимой ионами и высокоэнергетичными атомами
в приповерхностный слой;
• определение общего количества низкоэнергетичных частиц;
• вычисление количества низкоэнергетичных частиц, находящихся достаточно близко от поверхности мишени;
• определение количества частиц, находящихся близко к границе, способных преодолеть поверхностный барьер и покинуть мишень.
Решение всех четырех задач основано на одном уравнении переноса. Перейдем
к краткому обоснованию этого уравнения, следуя работе [7]. Обозначим через
среднее количество атомов, движущихся в момент времени в слое
со
скоростью из интервала
. Тогда число атомов с такими скоростями,
равно
проходящими через плоскость , лежащую внутри мишени, за время
– проекция скорости на ось , перпендикулярную плоской поверхности
где
мишени. При этом общее число атомов, распыленных с поверхности мишени
, может быть найдено по формуле
где интегрирование по скоростям производится по всем
с отрицательными
проекциями на ось (внутрь мишени) и таким, что этой скорости достаточно
для преодоления поверхностного потенциального барьера. Если выполнена так
называемая стохастическая гипотеза (гипотеза молекулярного хаоса) о том, что
после каждого столкновения устанавливается случайное распределение атомов
в пространстве координат и скоростей, т. е. начальные условия для последующих столкновений определяются только функцией распределения
11
Copyright ОАО «ЦКБ «БИБКОМ» & ООО «Aгентство Kнига-Cервис»
и не зависят от предыдущих столкновений, то согласно общим положениям
теории Больцмана эта функция удовлетворяет уравнению
где
,
– атомная плотность мишени, а функция
представляет собой дифференциальное сечение рассеяния ([4]),
- скорость
- скорость выбитого атома (атома отдачи),
.
отраженной частицы,
Уравнение (1.2) является “обращением” обычной формы уравнения переноса и
используется в теории переноса нейтронов. Оно может быть выведено стандартным способом, описанным, например, в работе [4]. Рассмотрим движущуюся частицу, которая в момент времени расположена в точке . В течение
времени она может испытать или не испытать столкновение. Следовательно,
Первый член в правой части представляет собой вероятность столкновения
двух частиц с таким исходом, при котором частицы приобретут скорости и , умноженную на сумму вкладов обеих частиц в функцию , проинтегрированную по всем возможным столкновениям. Второй член дает вероятность того, что столкновения не произойдет, умноженную на вклад в функцию
одной частицы с неизменившейся скоростью, но новыми начальным положением и начальным моментом времени. Раскладывая функции в этом равенстве
по степеням
и отбрасывая слагаемые выше первого порядка, приходим к
уравнению (1.2).
1.1.2. Модели динамики рельефа поверхности, облучаемой ионами. Рассмотренное выше уравнение переноса является общепризнанной моделью, описывающей процессы, происходящие в мишени при взаимодействии ионов с поверхностью облучаемого тела. Однако при математическом и компьютерном
моделировании изменения рельефа поверхности, подвергаемой ионной бомбардировке, обычно применяют дифференциальные уравнения в частных производных, которые получают из транспортной модели при различного рода упрощающих предположениях. Наиболее простые уравнения, которые могут быть
получены и непосредственно из экспериментальной угловой зависимости коэффициента распыления, будут рассмотрены ниже первыми. Эти уравнения, яв12
Copyright ОАО «ЦКБ «БИБКОМ» & ООО «Aгентство Kнига-Cервис»
ляющиеся по существу дифференциальными формами экспериментальной угловой зависимости коэффициента распыления, называют распылительными моделями. Исторически они являются первыми достаточно хорошо изученными
математическими описаниями процесса распыления, допускающими сравнительно простое компьютерное исследование и адекватно отражающими этот
процесс в простых случаях. Определенную сложность при использовании этих
распылительных моделей вносит появление особенностей на изучаемых поверхностях, в окрестностях которых соответствующие уравнения теряют адекватность. Сопутствующее этому экспериментально наблюдаемое появление
кромок на первоначально гладкой поверхности затрудняет моделирование и
требует уточнения распылительного закона для малых окрестностей особых точек. Одним из способов в какой-то мере исправить ситуацию может введение в
уравнения малой диффузии и переход к новой модели – уравнению Бюргерса.
Ценой такого перехода является существенное усложнение и смена типа рассматриваемого уравнения в частных производных. Еще одним вариантом
обобщения модели процесса распыления является переход к уравнению БредлиХарпера [11], которое в некоторых случаях позволяет обосновать возникновение волнообразного рельефа на облучаемой ионами поверхности тела. Это
уравнение в частных производных четвертого порядка будет рассмотрено ниже.
1.1.2.1. Общие закономерности выхода распыленных атомов. Экспериментальная проверка подтверждает, что при ионной бардировке поверхности мишени в большинстве случаев количество распыленных частиц пропорционально числу упавших на мишень ионов. Соответствующий коэффициент пропорциональности называют коэффициентом распыления (а иногда выходом распыления, от английского термина sputtering yield). Обычно его обозначают или латинской буквой Y (от слова yield), или буквой S (от слова sputtering). При этом
где в число распыленных атомов включаются как частицы, вылетевшие из мишени в виде ионов, так и в составе кластеров. Если распыляются атомы разных
видов, то для каждого из них вводят свой парциальный коэффициент распыления. Коэффициент распыления зависит от параметров мишени (состава, структуры, ориентации поверхности по отношению к кристаллографической решетке, шероховатости поверхности) и характеристик падающих ионов (типа, энергии, направления падения). Если энергия связи атомов мишени велика, то коэффициент распыления мал, а с увеличением атомного номера падающего иона
коэффициент распыления растет. При возрастании энергии падающего иона коэффициент распыления сначала растет, затем падает. При этом значение максимума приходится на энергию от нескольких кэВ до нескольких десятков кэВ.
Рост коэффициента распыления в начале этой зависимости до некоторого значения является следствием увеличения энергии, передаваемой атомам мишени,
а последующее убывание при больших энергиях объясняется тем, что слишком
“быстрый” ион пролетает мимо поверхностных атомов мишени, не успевая пе13
Copyright ОАО «ЦКБ «БИБКОМ» & ООО «Aгентство Kнига-Cервис»
редать им энергию. При этом, двигаясь внутрь мишени, первичный ион все
равно теряет свою кинетическую энергию, передавая ее атомам, но происходит
это уже на большей глубине, откуда атомы, выведенные из состояния равновесия, не могут покинуть мишень. Для тяжелых ионов, имеющих большее эффективное сечение, спад энергетической зависимости коэффициента распыления
наблюдается при более высоких энергиях, чем для легких. Приведенные в работе [16] экспериментальные и теоретические зависимости коэффициента распыления от энергии первичных ионов аргона для мишеней из алюминия и
кремния представлены на рис. 1.1.
Важнейшей характеристикой распыления, определяющей развитие рельефа
на поверхности, является ее зависимость от угла, измеряемого между нормалью
к поверхности и направлением ионного потока. При этом с возрастанием угла
от нуля до значений в 60-70 градусов коэффициент распыления растет, превосходя в максимуме до 10, а иногда и более раз значения при нуле градусов, а затем быстро падает, достигая нулевого уровня при значениях угла в 90 градусов.
Первоначальный рост выхода объясняется тем, что, когда движущийся атом
или ион взаимодействует с атомом мишени, происходит передача энергии тем
большая, чем ближе столкновение к лобовому. Но при лобовом столкновении
не происходит изменения направления движения. При касательном взаимодействии неподвижная частица получает меньше энергии. При этом потери энергии в цепочке атомных столкновений, начавшейся на падающем ионе и заканчивающейся на распыленном атоме, меньше в том случае, если меньше разница
в направлениях движения падающего иона и распыленного атома. Для первичного иона, падающего по нормали, разница в направлении его движения и направлении движения частицы, покинувшей мишень, должна быть больше 90
градусов. При увеличении угла падения ограничение на разницу в направлениях движения падающих и выбитых частиц ослабевает, и коэффициент распыления растет. При скользящих углах падения, превышающих 60-70 градусов, начинает преобладать другой эффект – отражение падающих ионов от мишени. В
14
Copyright ОАО «ЦКБ «БИБКОМ» & ООО «Aгентство Kнига-Cервис»
этом случае движущийся ион, постепенно приближаясь к поверхности распыляемого тела, успевает передать атомам, находящимся вблизи его траектории,
лишь небольшую часть своего импульса. Но как только компонента импульса
самого иона, перпендикулярная границе мишени, поменяет знак, ион начнет
удаляться от нее. В результате ион отражается, так и не передав сколько-нибудь
значительную часть своей энергии хотя бы одному атому мишени, что и вызывает уменьшение коэффициента распыления при скользящих углах. На рис. 1.2
показаны экспериментальные и теоретические зависимости коэффициента распыления от угла падения ионов аргона на поверхности мишеней из алюминия и
кремния, приведенные в работах [17,18].
При распылении монокристаллов для направлений ионного пучка близких к
какому-либо направлению плотной упаковки атомов, коэффициент распыления
падает. Этот эффект, называемый каналированием, объясняется тем, что
если ион попадает в промежуток между плотными рядами атомов – канал, и
при этом его импульс направлен вдоль
канала, то он может пройти достаточно большое расстояние в кристалле,
мало передавая свою энергию каждому из атомов, расположенных вблизи
его траектории. Поэтому для монокристаллов к зависимости от угла между вектором нормали к поверхности
и направлением движения первичных
ионов – полярного угла добавляется
еще зависимость коэффициента распыления от угла между проекцией направления ионного пучка на плоскую поверхность и каким-либо выбранным направлением в этой поверхности, называемого азимутальным. На рис. 1.3 показана экспериментальная угловая зависимость коэффициента распыления нитрида бора ионами Ar+ c энергией 1кэВ, полученная в работе [17], демонстрирующая эффект каналирования.
15
Copyright ОАО «ЦКБ «БИБКОМ» & ООО «Aгентство Kнига-Cервис»
Отметим, что важнейшими результатами, полученными в работах [7,8], являются теоретические зависимости коэффициента распыления от энергии и угла падения первичных ионов, основанные на формуле (1.1). При этом без преувеличения можно сказать, что именно эти зависимости в основном определяют характер рельефа, формируемого под воздействием ионной бомбардировки.
Задача исследования условий формирования того или иного типа рельефа на
поверхности, подвергаемой ионному облучению, чрезвычайно важна для многих приложений, в том числе и связанных с технологиями микро- и наноэлектроники. Следует отметить, что полученные в [7,8] соотношения во многих
случаях значительно расходятся с экспериментальными зависимостями коэффициента распыления и поэтому в дальнейшем предпринимались неоднократные попытки по уточнению теоретических формул из этих работ с целью их
лучшего согласования с постоянно накапливающимися экспериментальными
данными. Рассмотрим некоторые так называемые полуэмпирические модели
для энергетических зависимостей коэффициента распыления, полученные в работе [18] и использованные для расчета этих зависимостей для конкретных пар
мишень – ион в лаборатории NPL (National Physical Laboratory, UK [19]).
1.1.2.2. Зависимость коэффициента распыления от энергии. Анализируя равенство (1.1.) и используя результаты работ [3-6], Зигмунд обосновал формулу
для энергетической зависимости коэффициента распыления
, которая была
им представлена в виде
– поверхностная энергия связи атомов
где – энергия первичных ионов,
– коэффициент ядерного торможения атома с энергией , –
мишени,
безразмерный коэффициент, учитывающий возможность обратного отражения
падающего иона. Заметим, что при расчетах по формуле (1.3) величину
обычно заменяют теплотой сублимации атомов мишени. В результате всевозможных уточнений в дальнейшем выражение для коэффициента ядерного торможения
было приведено к виду
атомные номера первичного иона и атома мишени,
,
– их масгде ,
– диэлектрическая проницаемость вакуума. При
сы, – заряд электрона,
этом параметр
определяется равенством
16
Copyright ОАО «ЦКБ «БИБКОМ» & ООО «Aгентство Kнига-Cервис»
где
– боровский радиус. Величина
, как показано в статье Мацунами и
др. [20], может быть рассчитана с помощью приближений к уравнениям из
фундаментальной работы [4] по формуле
где приведенная энергия определена равенством
из соотношения (1.5) в формулы (1.4) и
Подставляя значение параметра
(1.7) и считая, что величина энергии измеряется в электрон-вольтах, получим
Формулы (1.3), (1.6), (1.8), (1.9) могут быть использованы для расчета энергетических зависимостей коэффициента распыления. В дальнейшем, ввиду того
что в этих соотношениях не учитывался вклад неупругого электронного торможения, роль которого возрастает с ростом энергии первичных ионов, в ряде
последующих теоретических работ были ведены дополнительные параметры
для учета этого фактора. С этой целью вместо формулы (1.3) было предложено
соотношение
где
– коэффициент неупругого электронного торможения, вклад которого
– пороговая энергия, снижающая распыление
растет для высоких энергий,
– коэффициент, определяемый подгонкой к эксперипри низких энергиях,
ментальным данным. При этом величины , в разных работах трактуются поразному. Согласно сложившимся в настоящее время представлениям, коэффициент неупругого электронного торможения
может быть найден по формуле
где коэффициент пропорциональности
определяется равенством
17
Copyright ОАО «ЦКБ «БИБКОМ» & ООО «Aгентство Kнига-Cервис»
в
При этом для определения значений параметра и пороговой энергии
разных работах предлагаются разные варианты. Так, Мацунами и др. в [20]
предлагают вычислять их с помощью равенств
В то время как в работе [21] для вычисления этих величин используются соотношения
где значение
выбирается из приводимой в [21] справочной таблицы, а параметр определяется формулой
в равенстве
В исследовании Мацунами и др. [20] для значения степени
(1.10) предлагается выбор величины , учитывающего возможность обратного
отражения падающего иона, в соответствии с формулой
В работе [21] та же степень для разных элементов берется из справочной табопределена
лицы, и ее значение равно либо 2.5, либо 2.8. При этом величина
соотношением
Последняя оставшаяся проблема определения подгоночного параметра в работах [20,21] натолкнулась на серьезные препятствия, связанные с большим
разбросом вычисляемых величин, превышающим для некоторых элементов несколько десятков процентов. В конечном итоге авторы этих работ вынуждены
были предложить определять этот параметр из справочных таблиц. Позднее в
работе [18] было получены формулы для величины , устранившие этот недос18
Copyright ОАО «ЦКБ «БИБКОМ» & ООО «Aгентство Kнига-Cервис»
таток для ряда элементов. С этой целью в [18] был введен параметр
атомное расстояние, для которого справедлива формула
– меж-
где – плотность мишени, а – число Авогадро. Тогда для параметра
следования [20] можно использовать представление
из ис-
Для расчетов по формулам работы [21] величину
нужно брать в виде
Соотношение (1.10), параметры которого определялись по формулам из исследования Мацунами и др. [20], дополненным уточненными параметрами из работы [18], были использованы для вычисления энергетических зависимостей
коэффициента распыления мишеней ионами инертных газов и опубликованы на
сайте NPL [19]. В настоящее время эти результаты дают хорошее приближение
к экспериментальным значениям и являются, по-видимому, наиболее достоверными.
1.1.2.3. Угловая зависимость выхода распыленных атомов. Важнейшей характеристикой для моделирования динамики морфологии поверхности, облучаемой ионами, является угловая зависимость коэффициента распыления, представляющая собой функцию выхода выбитых с поверхности атомов, приходящихся на один падающий ион от угла падения, отсчитываемого от локальной
нормали в данной точке поверхности. Одним из первых соотношений, описывающих эту зависимость, было предложенное Зигмундом в работе [8] равенство
где
– коэффициент распыления при нормальном падении ионов, θ – угол
между направлением ионного потока и локальной нормалью к поверхности –
постоянная, определяемая характеристиками ионов и материала облучаемой
поверхности. Эта теоретическая формула, справедливая лишь для достаточно
малых углов и практически неприменимая для компьютерных расчетов, послужила отправной точкой для последующих полуэмпирических формул, распространенных уже на более широкий диапазон углов наклона ионного потока.
Одним из вариантов простых формул такого вида была предложенная в работе
Оечснера [22] зависимость
19
Copyright ОАО «ЦКБ «БИБКОМ» & ООО «Aгентство Kнига-Cервис»
Данная квадратичная зависимость от угла падения первичных ионов лишь в самых общих чертах описывала поведения коэффициента распыления, и поэтому
даже в первых работах по численному моделированию эрозии поверхности
мишени под воздействием ионного потока исследователи предпочитали использовать более точные формулы. Значительно лучшую аппроксимацию угловой зависимости коэффициента распыления во многих случаях дает известная
полуэмпирическая формула Ямамуры, полученная в [23],
Простой анализ равенства (1.11) показывает, что параметры и однозначно
, при котором достигается максимальное
определяются значениями угла
, где
значение выхода распыленных атомов поверхности, и отношением
,
. Эта связь устанавливается соотношениями
Следует отметить, что приведенные формулы во многих случаях дают лучшее
приближение угловой зависимости коэффициента распыления к экспериментальной, чем значения, предсказываемые теорией Зигмунда [7,8]. Необходимо
отметить также, что сам вид функции, определяемой равенством (1.11), во многих случаях не годится для достаточно точной аппроксимации экспериментальных значений угловой зависимости выхода. Наглядным подтверждением этому
являются результаты, демонстрируемые на рис. 1.2, где приведены графики коэффициента распыления, полученные аппроксимацией угловой зависимости с
помощью формулы (1.11) экспериментальных значений для мишеней из алюминия и кремния при бомбардировке ионами аргона при оптимально выбранных значениях параметров и (аппроксимация по формуле (1.11) на правом
рисунке показана тонкой линией). В ряде последующих теоретических работ
(см., например, [15,24]) предпринимались усилия по уточнению формулы
(1.11). Так в [24] вместо этой формулы была предложена зависимость
где
20
Copyright ОАО «ЦКБ «БИБКОМ» & ООО «Aгентство Kнига-Cервис»
При этом поскольку в равенстве (1.12) по сравнению с исходной формулой
(1.11) изменения коснулись только значений входящих в эту формулу параметров, то проблемы точности аппроксимации экспериментальных угловых зависимостей с помощью соотношения (1.12) остались прежними. Следует отметить, что уже в ранних работах для численного моделирования процесса ионного распыления предпринимались попытки более точной аппроксимации экспериментальных угловых зависимостей выхода тригонометрическими многочленами. Так, в работе [16] для моделирования распыления поверхности кремния
ионами аргона с энергией 1 кэВ была использована формула
График этой зависимости показан на рис. 1.2 жирной линией. При этом благодаря наличию шести слагаемых в формуле (1.13) и весьма точно подобранным
коэффициентам, отклонение эмпирической кривой, определяемой (1.13), от
экспериментальных значений в целом оказалось достаточно малым. Исключение составили аппроксимированные значения из диапазона углов от 60 до 70
градусов, где наблюдается аномальная зависимость экспериментальных данных
от угла. Как уже было сказано выше, при распылении мишеней, имеющих кристаллическую структуру, возможно появление эффекта каналирования, проявлением которого является сильно выраженная аномальная зависимость коэффициента распыления от угла падения, показанная на рис. 1.3. Для численного
моделирования процесса ионного распыления таких мишеней, по-видимому,
вместо полуэмпирических формул типа Ямамуры (1.12) или аналогичных им
более оправданным было бы применение либо эмпирических формул, использующих разложения в ряды, либо использование гладких локальных аппроксимаций [25]. Отметим, что аномальное поведение экспериментальных угловых
зависимостей коэффициента распыления вблизи максимальных значений характерно и для некоторых других пар мишень – ион. Подобная ситуация рассматривалась в работе [25] для пары
, где моделирование динамики
рельефа проводилось с использованием гладких локальных аппроксимаций угловой зависимости коэффициента распыления.
1.1.2.4. Распылительная модель и ее дифференциальные формы. Наиболее
простыми, но во многих случаях достаточно эффективными моделями, описывающими динамику процесса распыления, являются нелинейные дифференциальные уравнения в частных производных, получаемые непосредственно из уг21
Copyright ОАО «ЦКБ «БИБКОМ» & ООО «Aгентство Kнига-Cервис»
ловой зависимости коэффициента распыления и ее связи со скоростью распыления точки поверхности, подвергаемой ионной бомбардировке. Такие уравнения называются дифференциальными формами распылительной модели. Для
вывода этих уравнений введем обозначения:
•
– управляемая интенсивность ионного потока,
•
– управляемый угол между направлением бомбардировки и
положительным направлением оси ,
•
– угол между локальной нормалью к поверхности и осью
в
момент времени ,
•
– коэффициент распыления (среднее количество атомов, выбиваемых одним ионом),
• – атомная плотность распыляемого вещества (количество атомов в единице объема).
Перейдем к выводу наиболее простой из дифференциальных форм распылительной модели. Следуя работе [26], предположим, что распыляемая цилиндрическая поверхность с образующими параллельными оси
определена уравне. Пусть ионный поток направлен вертикально вниз параллельно
нием
. Обозначим через
скорость распыления точки поверхности, наоси
правленную вдоль нормали к ней, тогда
Отсюда и из формулы
получим
Дифференцируя это равенство по , придем к соотношению
Меняя порядок дифференцирования в левой части этого равенства и используя
формулу
после простых преобразований придем к соотношению
22
Copyright ОАО «ЦКБ «БИБКОМ» & ООО «Aгентство Kнига-Cервис»
представляющему собой квазилинейное уравнение в частных производных
первого порядка. Уравнения такого типа возникают, например, в газовой динамике и описывают поведение ударных волн. Отметим, что конечной целью поставленной задачи является нахождение функции
, график которой
описывает динамику распыляемой поверхности. В связи с этим после решения
уравнения (1.16) найденную функцию нужно подставить в формулу (1.15) и интегрированием найти искомое уравнение распыляемой поверхности. Таким образом, постановка исходной задачи включает в себя систему уравнений (1.15),
(1.16), дополненную начальным условием
, определяющим распыляемую поверхность в начальный момент времени.
Более общей дифференциальной формой распылительной модели является
система уравнений, основанная на неявном представлении исследуемой поверхности в виде
, где
– координаты точки распыляемой
поверхности. Такая поверхность может иметь более общий вид, чем рассмотренная в предыдущем случае, и не обязательно допускать однозначное проектирование в целом на одну из координатных плоскостей. Вычисляя полную
производную функции
по , получим уравнение
которое может быть записано в виде
– вектор скорости точки поверхности с координатами
, направгде
ленный вдоль локальной нормали. Возводя в квадрат последнее равенство, получим
где скорость точки поверхности
определена формулой (1.14). Равенство (1.17) представляет собой нелинейное уравнение в частных производных первого порядка и называется уравнением эйконала. Это уравнение, описывающее волновой фронт, возникает, например, в геометрической оптике. Если исследуемая поверхность допускает представление в виде
, то
формула (1.17) может быть записана в виде
23
Copyright ОАО «ЦКБ «БИБКОМ» & ООО «Aгентство Kнига-Cервис»
Как и в предыдущем случае, для построения замкнутой системы уравнений,
описывающих процесс распыления, необходимо уравнение эйконала (1.17’)
(или (1.17), в более общей постановке) дополнить соотношением (1.14) и уравнением, связывающим угол с координатами
рассматриваемой поверхности и начальными условиями. Отметим, что такой подход, основанный на использовании уравнении эйконала, был использован при компьютерном моделировании процесса распыления в работах [27-29].
Еще один общий подход для изучения процесса распыления, основанный на
представлении исследуемой поверхности в параметрическом виде, приводит к
третьей дифференциальной форме распылительной модели. Как и в случае первой дифференциальной формы, при выводе уравнений ограничимся случаем
двух пространственных переменных. При этом можно считать, что имеет место
либо плоско параллельный случай, либо рассматриваемая поверхность обладает
осевой симметрией относительно оси
(в последнем случае нужно дополнимежду направлением бомбардировки и вертельно предполагать, что угол
тикальной осью равен нулю). Приведем вывод дифференциальной формы распылительной модели, приводящий к системе уравнений эрозии в переменных
цилиндрической системы координат (для декартовой системы достаточно заменить переменную на ), следуя результатам работы [30].
Рассмотрим процесс распыления поверхности аморфной мишени потоком
ионов, направление которого составляет угол
по отношению к оси
,
схематично показанный на рис. 1.4
(как уже упоминалось, для случая
осевой симметрии, приведенного
на этом рисунке,
), где
изображены последовательные положения исследуемой поверхности
в моменты времени t и
.
обозначено смещеЗдесь через
ние рассматриваемой точки поверхности вдоль локальной нормали за время , а через
на коордипроекции смещения
и
соответственно. Поскольку плотность ионного потока, панатные оси
к локальной нормали, равна
дающего на поверхность под углом
, то изменение
границы вдоль нормали может быть рассчитано
по формуле
Полагая
24
Copyright ОАО «ЦКБ «БИБКОМ» & ООО «Aгентство Kнига-Cервис»
запишем равенство (1.18) в виде
Из рис. 1.4 ясно, что проекции
и
гут быть найдены из соотношений
величины
на координатные оси мо-
Пусть исследуемая поверхность задана системой уравнений
Считая функции
приращений
,
дифференцируемыми, выпишем формулы для их
Отсюда с учетом равенств (1.19), (1.20) получим соотношения
зависит от
таким образом, что расОтметим, что в последней системе
сматриваемая точка границы поверхности, координаты которой получают приращения
и , перемещается вдоль своей локальной нормали. Предполагая
существование предела
разделим каждое из уравнений (5) на
зультате получим систему
и перейдем к пределу при
. В ре-
Поскольку для любого
соотношения (1.21) определяют поверхность, заданную параметрически, то имеет место равенство
25
Copyright ОАО «ЦКБ «БИБКОМ» & ООО «Aгентство Kнига-Cервис»
Умножая первое уравнение системы (1.22) на
, второе – на
и вычитая
из первого полученного равенства второе, с учетом последней формулы получим
Объединяя два последних соотношения, придем к системе уравнений эрозии
Основным достоинством полученной системы является то, что так же, как и
уравнение эйконала, она описывает процесс распыления поверхности общего
вида. В то же время, в отличие от этого уравнения, система (1.23) определяет
поверхность в явной форме (1.21), причем, как будет показано ниже, функции
,
в некоторых случаях могут быть найдены аналитически.
1.1.2.5. Уравнение Бредли-Харпера. Рассмотренные выше распылительные
модели весьма привлекательны с точки зрения численного исследования динамики рельефа поверхности, подвергаемой ионной бомбардировке. Однако они
не в состоянии объяснить появление сложных структур, экспериментально наблюдаемых при равномерном распылении мишеней ионным потоком. Одной из
таких структур, возникающих при эрозии облучаемой поверхности, является
волнообразный рельеф (“ripple” – в англоязычной литературе). Исследованию
условий возникновения морфологии поверхности такого вида посвящено множество публикаций, поскольку данные структуры чрезвычайно интересны с
прикладной точки зрения и особенно для новых технологий микро- и наноэлектроники. Объяснению появления периодически модулированного рельефа посвящено большое число теоретических работ. Первым исследованием, в котором данное явление получило удовлетворительное объяснение, была работа
Брэдли и Харпера [11], основанная на теории Зигмунда [7,8,31]. Основой результатов, полученных в [11], явилась формула для коэффициента распыления
поверхности ионным потоком, учитывающая морфологию самой поверхности.
Точнее, была установлена зависимость выхода количества распыленных атомов
мишени от локальной кривизны ее поверхности. Это позволило получить уточненную модель динамики распылительного процесса в виде дифференциального уравнения в частных производных четвертого порядка. Построенная модель
не только объяснила само появление волнообразного рельефа, но и позволила
получить формулу, удовлетворительно описывающую длину его волны. Рассмотрим кратко вывод основных результатов работы [11].
26
Copyright ОАО «ЦКБ «БИБКОМ» & ООО «Aгентство Kнига-Cервис»
В теоретических исследованиях поверхностной эрозии почти всегда предполагается, что коэффициент распыления не зависит от кривизны поверхности
(см. например, [8]). Это предположение, обычно приводящее к распылительным моделям, оказалось продуктивным при изучении временной эволюции
бомбардируемой ионами поверхности и является разумным приближением, если радиус кривизны в произвольной точке поверхности существенно больше
длины пробега иона. Тем не менее, если выход зависит только от угла падения
и другие эффекты игнорируются, малые возмущения плоской поверхности не
могут расти и поверхностные волны не образуются. Используя приближения
теории Зигмунда [31], найдем зависимость коэффициента распыления от кривизны. Согласно [31], скорость, с которой материал распыляется в произвольной точке поверхности, определяемой уравнением
, пропорциональна энергии, внесенной туда ионами при их столкновениях с атомами мишени,
результатом которых являются каскады взаимодействий, описанные выше.
Средняя энергия, переданная в точку
тела
ионом, кото, согласно [31]
рый до столкновения с поверхностью двигался вдоль оси
можно взять в виде
Здесь – полная внесенная энергия, – средняя глубина вклада энергии, а и
– размеры области вклада энергии, параллельные и перпендикулярные пучку
соответственно. При этом значения а, и соизмеримы по величине. Отметим,
что гауссова аппроксимация (1.24) часто используется во многих приложениях,
а доводы в поддержку ее применения в
данном случае приведены в работе [8].
Для определения зависимости коэффициента распыления от кривизны рассмотрим скорость эрозии в произвольной точке поверхности . Поместим
начало координатной системы в точку
, а ось
направим вдоль локальной
нормали распыляемой поверхности.
Направление ионного потока возьмем
лежащим в плоскости
, как показано на рис. 1.5. Для простоты высоту
поверхности выберем не зависящей от
, т. е.
(т. е. рассмотрим плоско параллельный случай). В результате
задача станет двумерной (случай произвольной поверхности будет рассмотрен
позже). Предположим также, что
изменяется достаточно медленно, так что
радиус кривизны в точке существенно больше . В этом случае членами
можно пренебвторого или более высокого порядка малости относительно
в точречь. Рассчитаем нормальную компоненту скорости поверхности
27
Copyright ОАО «ЦКБ «БИБКОМ» & ООО «Aгентство Kнига-Cервис»
ке когда равномерный поток с плотностью бомбардирует ее под углом
относительно оси .
, передаИон, падающий на поверхность в точку с координатами
имеет
ет заметное количество энергии в точку только если величина
, то высота поверхности в такой точке
может
порядок 1. Поскольку
быть приближенно найдена по формуле
. В первом приближении
относительно
плотность ионного потока в точке может быть представле, а длина дуги, лежащая между точками
на в виде
и
, равна . Учитывая эти формулы, из равенства
(1.24) получим соотношение
Здесь
, а – коэффициент пропорциональности между энергией, переданной в точку , и скоростью эрозии в этой точке.
Для приближенного вычисления интеграла в правой части (1.25) разложим
подынтегральное выражение по степеням
, удерживая члены нулевого и
первого порядка малости. В результате такого преобразования правую часть
равенства (1.25) можно будет представить в виде суммы гауссовых интегралов.
Выполним замену переменной интегрирования, положив
,
тогда соотношение (1.25) преобразуется к виду
где коэффициенты , , ,
найдены с помощью формул
,
не зависят от
и их значения могут быть
Заметим, что эти коэффициенты положительны при
. При этом
доминирующий вклад в значение интеграла получается, когда меняется вблизи единицы. Раскладывая подынтегральное выражение в (1.26) по степеням
и оставляя главные слагаемые, получим
28
Copyright ОАО «ЦКБ «БИБКОМ» & ООО «Aгентство Kнига-Cервис»
Вычисление интеграла в последней формуле приводит к равенству
где
Пусть теперь
– угловая зависимость коэффициента распыления плоской поверхности. Тогда
где – число атомов в единице объема тела. Отсюда и из формулы (1.28) получим
Из определения параметров
и
следует, что
и поэтому величина
является возрастающей функцией , так же как и
. Следовательно,
представляет собой возрастающую функцию угла падения . В реальных системах коэффициент распыления начинает убывать, когда угол падения превышает некоторое критическое значение и влияние отражения ионов становится
доминирующим. Поскольку в рассматриваемом приближении эффект отражения ионов от мишени не учитывается, то для скользящих углов падения полученные формулы не применимы.
При переходе к трехмерной модели объем вычислений значительно возрастает. В связи с этим приведем лишь окончательные результаты, полученные в
исследовании [11]. Используя прежнюю координатную систему, изображенную
на рис. 1.5, и вычисляя скорость распыления поверхности в точке с помощью
рассуждений, аналогичных приведенным выше, в [11] получена формула
29
Copyright ОАО «ЦКБ «БИБКОМ» & ООО «Aгентство Kнига-Cервис»
где
При этом, как и в двумерном случае, коэффициенты и управляют зависимостью скорости эрозии от кривизны поверхности. Заметим, что при
формула (1.29), справедливая для общей поверхности, переходит в соотношение (1.27), выведенное для плоско параллельного случая.
Пусть теперь невозмущенная поверхность лежит в плоскости
и ее
возмущение
достаточно мало. Предположим, что высота поверхности меняется достаточно медленно, так что для малых времен можно учитывать только производные от первых порядков. Пусть направление ионного
, а его направление составляет угол с осью .
потока лежит в плоскости
Тогда из формулы (1.29) получим
– скорость эрозии невозмущенной плоской погде
верхности. Заменяя последнее приближенное соотношение точным равенством,
получим уравнение
описывающее динамику распыляемой поверхности при сделанных выше допущениях. Отметим, что основное уравнение распылительной модели (1.14) может рассматриваться теперь как более грубое приближение процесса эрозии, не
учитывающее зависимость коэффициента распыления от локальной кривизны
поверхности в плоско параллельном случае, и является формальным следствием соотношения (1.29) при
.
1.2. Исследование моделей динамики рельефа
В этой части работы будет проведен анализ построенных моделей процесса
ионного распыления поверхности мишеней. Результаты, полученные с помощью более простых распылительных моделей, рассмотрим первыми. Как уже
30
Copyright ОАО «ЦКБ «БИБКОМ» & ООО «Aгентство Kнига-Cервис»
отмечалось выше, система уравнений эрозии в некоторых случаях допускает
аналитическое решение, что позволяет достаточно легко построить соответствующую компьютерную модель динамики процесса распыления поверхности.
Переход к значительно более сложному уравнению Брэдли-Харпера серьезно
усложняет численный анализ, а ограничения, накладываемые при его выводе,
не позволяют провести адекватное исследование динамики высокоамплитудного рельефа и изучение процесса при больших углах бомбардировки. Тем не менее приближенные аналитические расчеты, проведенные в рамках этой модели,
дают важную информацию об условиях возникновения периодических возмущений и связи периода с параметрами распылительного процесса.
1.2.1. Исследование распылительных моделей. Несмотря на сравнительную простоту распылительных моделей, в их рамках удается изучать динамику
достаточно сложных структур, реально возникающих в экспериментальных исследованиях. При этом возможно применение адекватного численного исследования и построение компьютерных моделей. Следует отметить, что ввиду
сложности возникающих в процессе моделирования поверхностных структур
весьма важной задачей является визуализация результатов проводимых расчетов. Ниже будет проведено аналитическое исследование системы уравнений
эрозии, для случая равномерного ионного потока получено его общее решение
и продемонстрированы некоторые результаты компьютерного моделирования.
1.2.1.1. Вывод основных расчетных формул. Рассмотрим поверхность, распыляемую ионным потоком, определяемую параметрически системой уравнений (1.21). Как было показано выше, в этом случае функции
и
удовлетворяют системе уравнений эрозии (1.22). Пусть график, определяемый этой системой, однозначно проектируется на ось
, тогда рас. Посматриваемая поверхность может быть определена уравнением
, приведем систему (1.22) к виду
лагая в этом случае
Считая для простоты ионный поток однородным, продифференцируем первое
равенство в (1.31) по . После упрощений получим
Выполняя дифференцирование в левой части, придем к формуле
31
Copyright ОАО «ЦКБ «БИБКОМ» & ООО «Aгентство Kнига-Cервис»
Это соотношение при
совпадает с уравнением (1.16) – первой дифференциальной формой распылительной модели, и является частным случаем
системы уравнений эрозии (1.22). Преобразуя правую часть равенства (1.32),
получим
Пусть теперь исследуемая поверхность, определяемая системой (1.21), задана
уравнением
. Полагая в этом случае в системе уравнений эрозии
, как и выше, придем к формуле
(1.22)
Считая, что при t=0 исследуемая поверхность определяется равенствами
и решая уравнения (1.33), (1.34) методом характеристик, получим искомую поверхность в параметрической форме
Если считать, что параметры и
уравнений эрозии (7) примет вид
не зависят от времени, то решение системы
Отметим, что в частном случае формулы, аналогичные (1.36), (1.37), были получены в работе [10]. Для неоднородного ионного потока решения системы
(1.22) в аналитической форме получить не удается.
Формулы (1.36), (1.37) позволяют проследить эволюцию распыляемой посистемой уравневерхности, определяемой в начальный момент времени
ний (1.35). При этом основными параметрами, определяющими динамику решений системы уравнений эрозии, являются атомная плотность вещества распыляемой поверхности – , интенсивность ионного потока –
и коэффициент распыления –
. Всюду ниже для расчетов будут использоваться сле32
Copyright ОАО «ЦКБ «БИБКОМ» & ООО «Aгентство Kнига-Cервис»
дующие значения этих физических постоянных:
,
, отвечающие распылению поверхности кремния ионами азота.
Следует отметить, что наиболее существенное влияние на результаты оказывает зависимость коэффициента распыления от угла бомбардирови
. Ниже
используется аппроксимация коэффициента распыления по формуле Ямамуры
(1.11). При этом, согласно результатам работы [10] значения
параметров формулы Ямамуры,
оптимально
соответствующие
экспериментальным данным распыления кремния ионами азота,
определяются
равенствами:
,
,
. График функции
, отвечающий этому набору
параметров, и экспериментально
полученные точки из [10] приведены на рис. 1.6. В дальнейшем
все расчеты в данной работе будут проводиться при выбранных
выше значениях атомной плотности мишени и интенсивности ионного потока и аппроксимации коэффициента распыления
формулой (1.11) при указанных значениях параметров.
1.2.1.2. Компьютерное моделирование и визуализация результатов. В этой
части работы будут представлены некоторые результаты визуализации численного исследования решений системы уравнений эрозии (1.22) для поверхностей
с осевой симметрией и однородного ионного потока, направленного вдоль оси
симметрии, а также в плоско параллельном случае ([30, 32]).
С целью тестирования используемых численных методов и алгоритмов и
для проверки адекватности самой системы уравнений эрозии (1.22) в трехмерном случае рассмотрим вначале стандартный пример распыления шара однородным вертикальным ионным потоком. Соответствующие теоретические и
экспериментальные результаты приведены, например, в работе [33]. Для обеспечения лучшей возможности сравнения на рис. 1.7 показаны одновременно
экспериментальные результаты распыления оловянного шарика, полученные в
работе [34] методами электронной микроскопии (рис. 1.7 a-c), теоретические
профили распыления сферы, рассчитанные с помощью кривой замедления эрозии приведенные в монографии [33] (рис.1.7 j), и результаты компьютерного
моделирования (рис. 1.7 a′-f ′) из работы [30]. Приведенная иллюстрация позволяет на тестовом примере распыления поверхности шара оценить хорошее качественное совпадение численных, теоретических и экспериментальных результатов. Достаточно интересной особенностью показанных на рис. 1.7 компьютерных графиков является появление кромок на распыляемой поверхности,
которые часто наблюдаются при экспериментальных исследованиях процесса
33
Copyright ОАО «ЦКБ «БИБКОМ» & ООО «Aгентство Kнига-Cервис»
ионной бомбардировки мишени. Следует отметить, что появление кромок на
первоначально гладкой поверхности при использовании распылительной модели является одной из серьезных проблем, решение которой существенно влияет
на дальнейшую динамику процесса. При этом в компьютерный эксперимент
появление кромок также вносит дополнительные сложности. Имея в виду
большое значение наличия указанных особенностей, на рис. 1.8 продемонстрированы детали компьютерного моделирования процесса распыления сферы, позволяющие понять и объяснить появление и взаимное слияние кромок в процессе ионной бомбардировки мишени.
34
Copyright ОАО «ЦКБ «БИБКОМ» & ООО «Aгентство Kнига-Cервис»
Дальнейшие возможности применения системы уравнений эрозии для исследования процесса распыления рассмотрим на примере поверхности с начальным рельефом синусоидальной формы, описываемой в цилиндрической
системе координат уравнением вида
с амплитудой A и длиной
волны . Результаты компьютерного моделирования такой поверхности, прове-
денные с помощью численного решения системы уравнений эрозии (1.22), приведены на рис. 1.9. Как видно из графиков, показанных на этом рисунке, в процессе распыления первоначально гладкая поверхность постепенно переходит в
структуру, состоящую из четырех конусов. Затем в результате борьбы четырех
мод вначале происходит поглощение верхнего наименьшего конуса следующим
за ним. В процессе дальнейшего распыления конус, ставший верхним, снова
поглощается расположенным ниже. В итоге формируется структура, состоящая
из двух конических поверхностей, которая существует достаточно продолжительное время. Как и в предыдущем случае, такая динамика морфологии объясняется тем, что распыление поглощаемых поверхностей идет быстрей ввиду
больших углов наклона нормалей к вертикальной оси, которой параллельно направление ионной бомбардировки.
Весьма интересным видом рельефа, как с точки зрения динамики, так и с
точки зрения практического применения, является так называемая структура
“конус в лунке”, часто возникающая в экспериментальных исследованиях при
ионной бомбардировке мишеней полупроводниковых и некоторых других ма35
Copyright ОАО «ЦКБ «БИБКОМ» & ООО «Aгентство Kнига-Cервис»
териалов ионами азота и аргона. Один из таких экспериментальных результатов, показанный на рис. 1.10, приведен в работе [35]. Для моделирования и
изучения условий его возникновения на исследуемой поверхности с помощью
компьютерной модели был “создан” кратер-затравка с достаточно крутыми
склонами. Для его “формирования” использовался сфокусированный ионный
луч с плотностью, распределенной по закону Гаусса. На левом верхнем снимке
рис. 1.11 показан завершающий момент создания этого кратера. После этого
компьютерная модель переключалась на “распыление” равномерным ионным
потоком. Несколько промежуточных снимков и окончательный результат мо-
делирования показаны на рис. 1.11 ([30]). Следует отметить близость экспериментального и рассчитанного результатов на качественном уровне. При этом
достаточно большую наблюдаемую разницу в соотношениях амплитуды и высоты экспериментальной структуры и ее компьютерной модели можно объяснить отличием угловых зависимостей коэффициента распыления, поскольку
36
Copyright ОАО «ЦКБ «БИБКОМ» & ООО «Aгентство Kнига-Cервис»
экспериментальная структура была получена при бомбардировке ионами
поверхности GaP, в то время, как уже упоминалось выше, моделирование распыления проводилось для кремниевой мишени.
Как показано выше, система уравнений эрозии (1.22) может быть использована для моделирования распыления достаточно сложного поверхностного
рельефа однородным ионным лучом. Следующий результат демонстрирует
адекватность этой системы и для случая неоднородного и нестационарного
ионного потока, вызывающего возникновения и развития и более сложной
морфологии. На рис. 1.12 показан снимок окна программы, моделирующей
распыление поверхности ионным лучом с плотностью, распределенной по за-
кону Гаусса. При этом скорость движения луча изменялась по гармоническому
закону таким образом, что в течение периода луч проходил расстояние, равное
его гауссовой ширине. Из рисунка видно, что результатом моделирования является волнообразный рельеф, период которого равен гауссовой ширине луча.
При этом другие аналогичные результаты компьютерного моделирования показали, что попытки уменьшения периода волнообразного рельефа путем уменьшения амплитуды синусоидально меняющейся скорости луча при сохранении
его ширины вызывают быстрое затухание аспектного числа этого рельефа и
даже при наличии произвольного начального рельефа приводят к его фактической полировке. Следует отметить, что распределение плотности ионного потока может оказывать существенное влияние на динамику морфологии распыляемой поверхности. При этом измерение истинного распределения плотности
ионного луча в реальном эксперименте сопряжено со значительными сложностями. Поэтому при интерпретации результатов экспериментальных работ
обычно считают, что плотность ионного потока распределена по закону Гаусса.
Проведенное исследование демонстрирует невозможность объяснения феномена волнообразного рельефа с малым периодом в рамках распылительной модели при использовании широкого ионного пучка с гауссовым распределением
37
Copyright ОАО «ЦКБ «БИБКОМ» & ООО «Aгентство Kнига-Cервис»
плотности. Отметим, что анонсированный в работе [36] результат показывает,
что если плотность ионного потока распределена по специальному закону, то
при определенных условиях все-таки есть возможность получения волнообраз
ного рельефа с периодом много меньшим ширины ионного луча.
1.2.1.3. 3D визуализация результатов моделирования. Трехмерные изображения результатов компьютерного моделирования процесса ионного распыления поверхности мишени, приведенные выше, демонстрируют высокую информативность, которую не могут дать обычные двумерные графики. Еще
большую перспективу в этом направлении открывает применение стереоизображений. При изучении динамики морфологии поверхности, подвергаемой
ионному облучению, большую роль могут играть детали, первоначально
имеющие очень малые размеры, которые, развиваясь, в дальнейшем становятся
определяющими глобальную структуру всего исследуемого объекта. В связи с
этим при компьютерном моделировании самых различных процессов все большее внимание начинают уделять не только точности математических расчетов,
но и реалистичности визуализации изучаемых объектов. При хорошем разрешении и качестве снимков в стереоизображении проявляются мельчайшие детали объектов, что позволяет судить об их физических размерах и положении в
пространстве. С точки зрения информативности эти методы наиболее полно
раскрывают и наглядно демонстрируют все нюансы моделируемого объекта. Во
многих случаях методы стереовизуализации просто незаменимы. В настоящее
время стереоизображения постепенно начинают использоваться во многих приложениях, использующих компьютерное моделирование как его важную, заключительную часть. При этом применение динамической стереовизуализации
привносит в процесс моделирования, кроме удобства, еще и эффекты “присутствия” и “погружения” исследователя “внутрь” изучаемого явления, создавая
иллюзию работы с реальными физическими объектами. При наличии адекватной компьютерной модели, снабженной развитым графическим интерфейсом,
могут создаваться виртуальные лаборатории и комплексы, аналогичные реальным физическим. При этом стоимость таких моделей может быть несравнимо
меньше физических установок, а скорость получения результатов намного выше.
Среди известных способов получения стереоизображений наиболее простым и весьма популярным является анаглифный метод (от греч. anagliphos –
рельефный). Суть этого метода состоит в наложении соответствующим образом
смещенных друг относительно друга изображений для каждого глаза. Каждое
из этих двух изображений, составляющих стереопару, окрашивается в дополнительные цвета, например в малиновый и бирюзовый. То, что данные цвета находятся на разных концах спектра, упрощает их сепарацию и при использовании специальных очков с малиновой линзой для левого глаза и бирюзовой для
правого у пользователя возникает иллюзия объемного изображения. При этом
указанный выбор цветов линз обеспечивает проявление стереоэффекта наилучшим образом. Кроме того, красный, синий и зеленый цвета являются основными, и при аддитивной модели синтеза из этих цветов могут быть получены
все остальные цвета.
38
Copyright ОАО «ЦКБ «БИБКОМ» & ООО «Aгентство Kнига-Cервис»
Применение стереовизуализации, и в частности анаглифного метода как
наиболее простого и дешевого способа получения в высокой степени информативных стереоизображений объектов и процессов для моделирования в микрои наноэлектронике, становится актуальной задачей для научных исследований
и разработки новых технологий электронных приборов следующего поколения.
1.2.2. Возникновение “ripple” структур в модели Бредли-Харпера. Как
уже было отмечено выше, одной из важнейших задач современной наноэлектроники является разработка новых технологий для создания элементов нанометрового размера. Хорошую перспективу в этом направлении открывает экспериментально наблюдаемый эффект спонтанного формирования волнообразного рельефа при ионном облучении поверхности однородным ионным потоком. При этом при разных условиях можно наблюдать волны различной формы,
амплитуды и периода. Задача выяснения условий возникновения таких структур
и управление их параметрами является весьма актуальной для разработки новых
технологий, основанных на данном явлении. Согласно сделанному выше замечанию, данное явление не может быть объяснено в рамках распылительной модели, и для его обоснования в разное время предлагались различные теории.
Первое достаточно удовлетворительное объяснение спонтанного формирования
волнообразного рельефа было дано в работе Брэдли и Харпера [11] на основе
анализа полученного ими уравнения (1.30).
Рассмотрим зависимость скорости эрозии
от кривизны поверхности
в рамках предположений пункта 1.1.2.5. Коэффициент
при
в (1.27)
. Это значит, что при нормальном падении пучка выстуотрицателен при
пы на поверхности распыляются медленнее, чем впадины (радиус кривизны
впадин отрицательный). Таким образом, распыление увеличивает амплитуду
возмущений на поверхности и приводит к неустойчивости. Чтобы наглядно
продемонстрировать, почему это происходит, рассмотрим процесс распыления
пика (рис. 1.13(a)) и впадины (рис. 1.13(b)). Энергия, переданная в точку ионом, упавшим в точке , равна
энергии, передаваемой в точку
ионом, попавшим на поверхность в этой точке. Но средняя
энергия, переданная в точку
ионом, который попал в точку ,
больше, чем выделенная в точке
ионом, упавшим в . Аналогичная ситуация наблюдается и
для точек
и . Поэтому скобольше,
рость эрозии в точке
чем в точке .
Поскольку
непрерывна, поверхность неустойчива по отношению к
развитию периодических возмущений, волновой вектор которых при малых углах падения параллелен проекции ионного пучка на поверхность мишени. Тем
не менее при углах, близких к скользящему падению ионов, нестабильность от39
Copyright ОАО «ЦКБ «БИБКОМ» & ООО «Aгентство Kнига-Cервис»
сутствует, поскольку
. Этот факт имеет простое физическое объяснение.
Рассмотрим эрозию выступа (рис. 1.14.a) и впадины (рис. 1.14.b) при ионном пучке с углом от нормального
падения. В среднем наибольшее количество энергии в точку
( )
приходит от ионов, упавших на поверхность, “вверх по течению” от
точки ( ). Хотя ионы, упавшие в
точке , передают в точку больше
энергии, чем те, которые ударили
поверхность в точке , передают в
число ионов, попавших на единицу площади за единицу времени
меньше в , чем в . Этот эффект
преобладает при углах , близких к
, поэтому скорость эрозии в
точке больше, чем в точке , когда направление пучка близко к скользящему
падению.
Найденная в пункте 1.1.2.5 зависимость коэффициента распыления от локальной кривизны поверхности позволяет изучить временную эволюцию малоамплитудного рельефа. Поскольку любое возмущение на плоской поверхности
в момент времени
может быть записано в виде суперпозиции гармонических модуляций высоты, то для малых времен каждая из гармоник будет эволюционировать независимо. Поэтому достаточно изучить временную эволюцию периодического возмущения
. Подставляя
функцию
где
,
,
и – вещественные числа, в уравнение (1.30), получим
всегда отрицательна, то амплитуда периодического
Поскольку величина
возмущение с волновым вектором
растет экспоненциально. Для углов
, меньших некоторого критического значения, величина
также отрицательна. Отсюда следует, что поверхность будет неустойчива к периодическим
возмущениям с произвольными волновыми векторами
. Разумеется, если пренебречь зависимостью коэффициента распыления от кривизны
поверхности
, то будет равен нулю. В этом случае амплитуда
малых возмущений не будет меняться с ростом времени.
40
Copyright ОАО «ЦКБ «БИБКОМ» & ООО «Aгентство Kнига-Cервис»
Следует отметить, что в данной модели быстрее всего растут периодические
возмущения с наименьшей длиной волны. Таким образом, следовало бы ожидать, что наблюдаемая длина волны периодической структуры будет сопоставима с микроскопической длиной отсечки по теории, т. е. в диапазоне вклада
энергии (т. е. порядка величины ). Тем не менее, как известно из экспериментальных работ, значение может на два порядка превосходить величину .
Для объяснения этого эффекта включим в теорию поверхностную самодиффузию. В результате вместо соотношения (1.30) получим уравнение
При этом если поверхностная диффузия активирована температурой, то коэффициент дается равенством
– коэффициент самодиффузии, – поверхностная плотность свободной
где
энергии, а – поверхностная плотность диффундирующих атомов. Отметим,
что это выражение дает заниженное значение для величины при низких температурах и высоких плотностях ионного потока в случае, когда важную роль
играет диффузия, вызванная ионной бомбардировкой. Учет эффекта поверхностной самодиффузии в рассматриваемую модель не приводит к изменению равенства (1.38), но формула (1.39) преобразуется к виду
Следует отметить, что экспериментально наблюдается тот волновой вектор ,
который соответствует максимальному значению . Если
, то вол, где
новой вектор определяется выражением
Выражение (1.38) показывает, что волны перемещаются по поверхности со скоростью
. Если теперь
, то наблюдаются структуры с волновым вектором
, где
Причем эти волны стационарны, поскольку
и
. Определим теперь
и
. Из соображений
условия, при которых
. Раскладывая функции
и
в
симметрии следует, что
ряды по и пренебрегая слагаемыми более второго порядка малости, получим
формулы
41
Copyright ОАО «ЦКБ «БИБКОМ» & ООО «Aгентство Kнига-Cервис»
справедливые при малых значениях . Отсюда следует, что
малых , если
для
Условие (1.45) выполняется, если, например,
, т. е. когда продольное
,
рассеяние выражено не очень сильно. С другой стороны, для , близких к
положителен, а
– отрицателен. Условие
выполнено,
, где
– первый критический угол, а неравенство
например, при
справедливо при
, где
– второй критический угол
(
). Если выполняется (1.45), то
. Заметим, что значения и
не
должны быть равны, но если это не выполнено, то существует третий критический угол
, лежащий между и , для которого
.
Объединяя сказанное, можно заключить, что если выполняется (1.45), то
волновой вектор гармонической структуры параллелен проекции ионного пучка на поверхность мишени, если угол не слишком велик. Если же этот угол
близок к скользящему, то волновой вектор перпендикулярен ионному потоку.
При нормальном падении значения и равны, и направление волнового вектора произвольно. В этом случае могут возникнуть несколько волн, ориентация
которых определяется дефектами поверхности и влиянием граничных эффектов, что в конечном итоге может привести к формированию на поверхности
сетки пиков и впадин. Такое поведение часто наблюдается в экспериментах. Развитая теория предсказывает,
что длина волны периодического
рельефа определяется формулами
и зависит от
угла падения. Графики этих зависимостей для
,и
приведены на рис. 1.15. При этом наблюдаемая длина волны
определяется равенством
42
Copyright ОАО «ЦКБ «БИБКОМ» & ООО «Aгентство Kнига-Cервис»
где
. Из последней формулы и равенства (1.41) следует,
что при низких плотностях ионного потока и высоких температурах справедлива формула
где
– энергия активации процесса поверхностной самодиффузии.
Отметим, что при высоких температурах и низких плотностях ионного потока полученная формула дает хорошее совпадение с экспериментальными результатами работы [37]. Таким образом, уравнение (1.40), полученное в исследовании [11], может служить теоретическим обоснованием возникновения волнообразного рельефа при ионном распылении поверхности равномерным ионным потоком.
Литература
1. Grove W. R. Phil. Mag. 1853. № 5. Р. 203.
2. Vineyard G. H. Computer Experiments with Lattice Models. Interatomic Potentials and Simulation of Lattice Defects. N. Y.: Plenum Press, 1972.
3. Lindhard J., Nielsen V., Scharff M., Thomsen P. V.Vidensk. Selsk Mat. Fys.
Medd. 1963.V. 33, № 10.
4. Lindhard J., Scharff M., Schiott H. E. K. Dan. Vidensk. Selsk. Mat. Fys. Medd.
1963.V. 33, № 14.
5. Lindhard J., Nielsen V., Scharff M. K. Dan. Vidensk. Selsk. Mat. Fys. Medd.
1968.V. 36, № 10.
6. Lindhard J. K. Dan. Vidensk. Selsk. Mat. Fys. Medd. 1965.V. 34, № 7.
7. Sigmund P., Radiat. Eff. 1969.V. 1, № 15.
8. Sigmund P., Phys. Rev. 1969.V. 184, № 2.
9. Kim H.-B., Hobler G., Steiger A., Lugstein A., Bertagnolli E. Nanotechnology.
2007. № 18.
10. Birkgan S. E., Bachurin V. I., Rudy A. S., Smirnov V. K. Rad. Eff. 2004.
V. 159, № 3.
11. Bradley R. M., Harper J. M. E. Vac. Sci. Technol. A. 1988.V. 6, № 4.
12. Makeev M. A., Barabasi A.-L. Nucl. Instr. B. 2004. P. 222.
43
Copyright ОАО «ЦКБ «БИБКОМ» & ООО «Aгентство Kнига-Cервис»
13. Kardar M., Parisi G., Zhang Y. C. Phys. Rev. Lett. 1986. № 56.
14. Ziberi B., Frost F., Hoche T., Rauschenbach B. Phys. Rev. B. 2005. № 72.
15. Behrisch R., Eckstein W. Sputtering by Particle Bombardment. Springer. 2007.
16. Ducommun J.P., Cantagrel M., Moulin M. J. Mat. Sci. 1975. V. 10.
17. Yurasova V. E., Elovikov S. S., Zykova E. Y. J. Surf. Invest. 2007. V. 1, № 3.
18. Seah M. P., Clifford C. A., Green M. F., Gilmore I. S. Surf. Interface Anal. 2005.
№ 37.
19. URL:http://www.npl.co.uk/science/
20. Matsunami N., Yamamura Y., Itakawa Y., Itoh N., Kazumata Y., Miyagawa S.,
Morita K., Shimizu R. Radiat. Eff. Lett. 1980. № 57.
21. Yamamura Y., Tawara H. At. Data Nucl. Data Tables. 1996. V. 62. Р. 149.
22. Oechsner H. Appl. Phys. 1975. V. 8.
23. Yamamura Y. Radiat. Eff. 1984. № 80.
24. Yamamura Y., Itikawa Y., Itoh N. IPPJ-AM-26. Nagoya University, 1983.
25. Биркган C. Е. Вестник ЯрГУ. Сер. «Физика. Радиотехника. Связь».
2010. № 1.
26. Nobes M. J, Colligon J. S., Carter G. J. Mat. Sci. 1969. № 4.
27. Tagg M. A., Smith R., Walls J. M. J. Mat. Sci. 1986. № 21.
28. Nobes M. J., Katardjiev I. V., Carter G., Smith R. J. Phys. D. Appl. Phys.
1987. № 20.
29. Katardjiev I. V., Carter G., Nobes M. J. J. Phys. D. Appl. Phys. 1989. № 22.
30. Биркган С. Е. 3D моделирование развития структурированной поверхности
с радиальной симметрией под воздействием ионного облучения //Трехмерная
визуализация научной, технической и социальной реальности: тр. Первой межд.
конф. Т. 2. Ижевск, 2009.
31. Sigmund P. J. Mat. Sci. 1973.№ 8.
32. URL:http://www.nd.uniyar.ac.ru/index.php/birkgan.
33. Behrisch R. Sputtering by Particle Bombardment II. Springer, 1983.
34. Volkert C. A., Minor A. M. Focused Ion Beam Microscopy and Micromachining.
MRS bulletin. 2007. V. 32.
35. Сошников И. П., Берт Н. А. ЖТФ. 2000. V. 70, № 9.
36. Birkgan S. E. SIMS Europe. Munster, 2006.
37. Vasiliu F., Teodorescu I. A., Glodeanu F. J. Mat. Sci. 1975. № 10.
44
Copyright ОАО «ЦКБ «БИБКОМ» & ООО «Aгентство Kнига-Cервис»
2. Моделирование в задачах цифровой обработки сигналов
Линейные и нелинейные уравнения в частных производных во многих случаях являются классическими моделями физических процессов и используются
в науке со времен появления дифференциального и интегрального исчисления.
Значительное число практически значимых задач, описываемых этими уравнениями, было решено аналитическими методами. Однако в своем большинстве
наиболее интересные в настоящее время проблемы, как правило, уже не могут
быть решены без применения численных методов. Чаще всего при численном
исследовании задач для дифференциальных уравнений используется метод конечных разностей или метод сеток, реализующий переход от задач математической физики к дискретным разностным уравнениям. Последние уже сравнительно просто могут быть решены на ЭВМ. Следует отметить, что, помимо задач математической физики, разностные уравнения самостоятельно и с успехом
используются для моделирования процессов и явлений во многих отраслях
науки и техники. Весьма плодотворно применение дискретных разностных
уравнений при моделировании процессов в устройствах цифровой обработки
сигналов. Этой теме и посвящена данная глава.
2.1. Моделирование цифровых устройств разностными уравнениями
В настоящее время все большее число электронных устройств переходит на
цифровые методы работы с электрическими сигналами, которые во многих
случаях несут огромные преимущества по сравнению с еще недавно безальтернативно использовавшимися методами для аналоговых сигналов. Одним из
наиболее мощных средств цифровой обработки сигналов является цифровая
фильтрация. Цифровые фильтры широко используются в телекоммуникации,
для адаптивной фильтрации, например при подавлении эха в модемах, шума и
для распознавания речи. Цифровые фильтры способны удовлетворить самым
высоким техническим требованиям, причем характеристики цифрового фильтра могут быть легко изменены программно. Согласно принятой в настоящее
время классификации, различают два основных типа цифровых фильтров:
фильтры с конечной импульсной характеристикой (КИХ) и фильтры с бесконечной импульсной характеристикой (БИХ). Как следует из названия, данная
классификация относится к импульсным характеристикам фильтров. С точки
зрения эффективности и сложности математических моделей этих устройств
более интересными являются БИХ-фильтры, которые называют еще рекурсивными. Математические модели таких фильтров, свойства решений соответствующих им разностных уравнений и будут обсуждаться в этой главе. БИХфильтры получили такое название, потому что их импульсные характеристики
растянуты на бесконечном временном интервале. Это объясняется тем, что они
используют обратную связь. БИХ-фильтры обычно реализуются с помощью
звеньев второго порядка, которые называются биквадратными фильтрами потому, что описываются биквадратными уравнениями в -области. Фильтры вы45
Copyright ОАО «ЦКБ «БИБКОМ» & ООО «Aгентство Kнига-Cервис»
сокого порядка проектируют, используя каскадирование биквадратных звеньев.
Структура биквадратного БИХ-фильтра представлена на рис. 2.1.
Рис.2.1.
Математической моделью процессов, происходящих в этом устройстве, является дискретное разностное уравнение
В общем случае фильтр высшего порядка может быть описан разностным уравнением
Еще одна эквивалентная схема биквадратного фильтра, приведенная на рис. 2.2,
называется второй прямой формой реализации и требует использования только
двух регистров.
Рис.2.2.
46
Copyright ОАО «ЦКБ «БИБКОМ» & ООО «Aгентство Kнига-Cервис»
Можно показать, что уравнение, описывающее биквадратный БИХ-фильтр второй прямой формы реализации, может быть сведено к уравнению первой прямой формы реализации.
2.2. Начальные понятия теории дискретных разностных уравнений
Во многих случаях моделирование процессов в дискретных и цифровых динамических системах приводит к уравнениям вида
– заданная
Такие уравнения называются разностными. Здесь
переменных, – аргумент, принимающий целые значения,
функция
– неизвестная числовая последовательность. Число
называется поряд– линейная отноком разностного уравнения. Если функция
, то уравнение (2.1) может быть записано в
сительно аргументов
виде
где
– известные последовательности. Уравнение (2.2) называется линейным разностным уравнением. Исследованию таких уравнений
посвящено значительное число работ. Основные результаты теории разностных
уравнений содержатся в монографиях [1–3]. Введем в рассмотрение разности
формулами
тогда равенство (2.2) может быть записано в форме
линейно выражаются через
где элементы последовательности
коэффициенты уравнения (2.2). Соотношение вида (2.3) называется уравнением
в конечных разностях. Форма записи и свойства решений этого уравнения во
многом аналогичны свойствам линейного дифференциального уравнения -го
порядка [1].
Разностное уравнение -го порядка может быть приведено к системе
уравнений первого порядка. Покажем это на примере линейного уравнения
(2.2). Выполним в этом уравнении замены
47
Copyright ОАО «ЦКБ «БИБКОМ» & ООО «Aгентство Kнига-Cервис»
С учетом формул (2.4) разностное соотношение (2.2) может быть заменено эквивалентной системой
В векторно-матричных обозначениях последняя система может быть записана в
виде
где матрица
и векторы
и
определяются равенствами
Линейное разностное уравнение (2.2) называется однородным, если
для всех , и неоднородным – в противном случае. Аналогичная терминология
используется и для системы (2.5). Решением разностного уравнения (2.2) (системы уравнений (2.5)) назовем числовую последовательность
(векторную
), которая обращает в тождество это уравнение (эту
последовательность
систему).
В дальнейшем в этой главе в основном будут рассматриваться лишь линейные разностные уравнения (2.2) и системы уравнений типа (2.5) с матрицей
и вектором
общего вида. Нелинейные уравнения вида (2.1) и системы таких уравнений будут рассмотрены в конце настоящей главы в связи с задачей об устойчивости движения.
Важными свойствами решений уравнения (2.2) и системы (2.5), вытекающими из линейности этих уравнений, являются следующие утверждения, носящие общий характер.
48
Copyright ОАО «ЦКБ «БИБКОМ» & ООО «Aгентство Kнига-Cервис»
• Любая линейная комбинация решений линейного однородного уравнения
(системы уравнений) также является решением этого уравнения (этой системы).
• Разность любых двух решений линейного неоднородного разностного
уравнения (системы уравнений) является решением соответствующего однородного уравнения (однородной системы уравнений).
Назовем общим решением линейного разностного уравнения (системы
уравнений) решение, зависящее от некоторых произвольных постоянных, такое, что любое решение этого уравнения (системы, уравнений) может быть получено из этого решения при некоторых значениях постоянных. Следствием
приведенных выше свойств являете следующее утверждение.
• Общее решение линейного неоднородного разностного уравнения (системы уравнений) может быть представлено в виде суммы некоторого фиксированного (частного) решения этого уравнения (системы уравнений) и общего
решения соответствующего однородного уравнения (однородной системы).
Разностное уравнение (2.2) или система уравнений (2.5) может рассматриваться как рекуррентная формула, позволяющая определять решение для любых номеров n, если только заданы некоторые начальные значения. Так, для того чтобы определить решение уравнения (2.2) для
, достаточно задать
. Для системы (2.5) нужно задать козначения
. Для уравнения в конечных разностях (2.3) естестординаты вектора
венной формой задания начальных условий является задание величин
. Таким образом, дополняя уравнение (2.2) условиями
где
– произвольные постоянные, приходим к начальной задаче для
этого уравнения. Аналогично система разностных уравнений (2.5) вместе с условиями
дает начальную задачу для системы разностных уравнений. Для уравнения (2.3)
в конечных разностях начальная задача состоит в отыскании решения этого
уравнения, удовлетворяющего дополнительно условиям
– произвольные начальные данные.
Пространством состояний, или фазовым пространством уравнений (2.1),
(2.2) является
– мерное пространство, элементами которого являются сово. Любая такая совокупность
купности из действительных чисел
чисел может рассматриваться в качестве начального условия вида (2.7) или
(2.8) для определения единственного решения (2.1), (2.2) или системы (2.5).
Другое определение фазового пространства можно получить, если в качестве элементов -мерного пространства выбирать величины
в
где
49
Copyright ОАО «ЦКБ «БИБКОМ» & ООО «Aгентство Kнига-Cервис»
формулах (2.9). Выбором начального значения и произвольного элемента фазового пространства однозначно определяется все решение уравнений (2.1) –
(2.3) или системы (2.5).
Для автономных уравнений и систем (т. е. в случае, когда функция в правой
части (2.1) и коэффициенты и правые части уравнений (2.2), (2.3) и системы
(2.5) не зависят от ) каждому решению отвечает некоторая траектория в фазовом пространстве. Траектории разных решений либо не пересекаются, либо
совпадают. Так же, как и в теории дифференциальных уравнений, исследование
фазового пространства автономных разностных уравнений используется для
качественного анализа динамики решений этих уравнений.
Изучение разностных уравнений целесообразно начать с линейного разностного уравнения вида (2.2), в котором
, и аналогичной системы вида (2.5) с постоянной матрицей. Вначале будет подробно
рассмотрено простейшее разностное уравнение вида
Оно описывает процессы, происходящие в линейной цифровой системе, называемой рекурсивным цифровым фильтром первого порядка [4]. После этого будет исследовано разностное уравнение второго порядка
Это уравнение описывает процессы в линейной цифровой системе, называемой
рекурсивным цифровым фильтром второго порядка. Затем будут рассмотрены
методы решения разностных уравнений высшего порядка вида (2.2) и системы
уравнений вида (2.5). Решения таких уравнений и систем могут быть получены
разными способами. Для объяснения других методов решения воспользуемся
методом интегральных преобразований, идея которого состоит в том, что разностное уравнение путем применения прямого интегрального преобразования
приводится к алгебраическому для функции-образа. Разрешая полученное
уравнение относительно этой функции и применяя обратное преобразование,
получают решение разностного уравнения. При этом начальные условия, которым должно удовлетворять решение, легко учитываются, поскольку входят в
уравнение для функции-образа. Следует отметить, что лишь для сравнительно
небольшого класса уравнений, включающего линейные разностные уравнения с
постоянными коэффициентами, описанным методом удается получить решение
в явном виде. Однако во многих случаях анализ функции-образа помогает исследовать свойства решений исходного уравнения. Суть метода и результаты
будут изложены на примере однородного уравнения второго порядка. Заметим,
что полученные результаты могут быть легко перенесены на общее линейное
разностное уравнение вида (2.2).
2.3. Разностные уравнения первого порядка
Сначала рассмотрим однородное разностное уравнение первого порядка:
50
Copyright ОАО «ЦКБ «БИБКОМ» & ООО «Aгентство Kнига-Cервис»
где постоянные
отличны от нуля. Переписывая его в виде
и рассматривая как рекуррентное соотношение для определения
, получим
Таким образом, общее решение уравнения (2.10) может быть записано в виде
где C – произвольная постоянная. Заметим, что число
уравнения
является решением
которое называется характеристическим для разностного уравнения (2.10).
Рассмотрим теперь линейное неоднородное разностное уравнение
– заданная последовательность. Это уравнение описывает процессы в
где
динамических устройствах под воздействием внешней силы. Для решения
(2.11) воспользуемся методом вариации произвольной постоянной, аналогичным известному для линейных дифференциальных уравнений первого порядка.
Согласно этому методу, частное решение неоднородного уравнения (2.11) будем искать в виде
– новая неизвестная последовательность. Подставляя (2.12) в (2.11),
где
получим
Преобразуя последнее равенство, придем к соотношению
Поскольку выражение в квадратных скобках равно нулю ( – решение однородного разностного уравнения (2.10)), то последнее равенство может быть записано в виде
51
Copyright ОАО «ЦКБ «БИБКОМ» & ООО «Aгентство Kнига-Cервис»
Суммируя в последнем соотношении от 1 до n, получим
Выражая отсюда
и подставляя в (2.12), придем к формуле для общего решения неоднородного разностного уравнения (2.11)
Здесь
– начальное значение последовательности
.
2.4. Разностные уравнения второго порядка
Рассмотрим линейное однородное разностное уравнение с постоянными коэффициентами второго порядка
отличны от нуля. Результаты предыдущего пункта навогде коэффициенты
. Действительдят на мысль искать решение уравнения (2.13) в виде
но, подставляя
в (2.13), получим
Предполагая, что
, и сокращая на
, придем к уравнению
которое называется характеристическим для разностного уравнения (2.13).
Корни этого уравнения определяются формулой
Рассмотрим различные случаи расположения корней (2.14) в зависимости от
величины дискриминанта
.
Пусть
, тогда
– различные действительные корни (2.14). При
– решения разностного уравэтом последовательности
нения (2.13), и общее решение этого уравнения может быть записано в виде
52
Copyright ОАО «ЦКБ «БИБКОМ» & ООО «Aгентство Kнига-Cервис»
где
– произвольные постоянные. Решение уравнения (2.13) с начальными
имеет вид
условиями
(строгое обоснование последних двух формул вытекает из результатов следующего пункта). Если
, то
– кратный корень характеристического
уравнения (2.14) и разностное уравнение (2.13) имеет лишь одно линейно независимое решение вида
. Выполним в (2.13) замену
.
в (2.13), получим
Подставляя
или
Поскольку
– кратный корень уравнения (2.14), то
Подставляя эти значения в последнее соотношение, придем к равенству
.
Это уравнение можно представить в виде
Отсюда следует, что
– арифметическая прогрессия. Таким образом, реше.
нием последнего уравнения является, например, последовательность
разностное уравнение (2.13)
Следовательно, наряду с решением
. Отсюда и из принципа суперпозиции вытекает, что
имеет решение
общее решение (2.13) в случае, когда характеристическое уравнение (2.14) имеет кратный корень, может быть записано в виде
– произвольные постоянные. При этом решение уравнения (2.13) с
где
начальными условиями
имеет вид
, тогда характеристическое уравнение (2.14) имеет комПуть теперь
плексно сопряженные корни
53
Copyright ОАО «ЦКБ «БИБКОМ» & ООО «Aгентство Kнига-Cервис»
Представим их в виде
где для модуля и аргумента комплексного числа
верны формулы
. В рассматриваемом случае общее решение
При этом ясно, что
(2.13) в комплексной форме по-прежнему может быть записано в виде (2.15), но
теперь постоянные
нужно считать комплексными. Формула (2.15) даст
действительное решение (2.15), если эти постоянные выбрать комплексно сопряженными. Чтобы записать это решение, рассмотрим метод выделения действительных элементарных решений из комплексных. Последние имеют вид
Справедливо следующее утверждение, имеющее общий характер.
Лемма 2.1. Пусть уравнение (2.13) с действительными коэффициентами
имеет комплексное решение
, тогда последовательности
и
также являются решениями (2.13).
в уравнение
Доказательство. Подставим функцию
(2.13), придем к тождеству
Разделяя действительные и мнимые части, получим
Последнее соотношение означает, что каждая из последовательностей
и
является решением уравнения (2.13). Лемма доказана.
Доказанное утверждение позволяет выделить из комплексных решений
(2.16) действительные, которые имеют вид
и записать общее решение разностного уравнения (2.13) для случая комплексных корней формулой
54
Copyright ОАО «ЦКБ «БИБКОМ» & ООО «Aгентство Kнига-Cервис»
где
– произвольные, уже действительные постоянные. При этом решение
имеет вид
с начальными условиями
Так же, как и в предыдущем пункте, соответствующее неоднородное разностное уравнение второго порядка может быть решено методом вариации произвольной постоянной, который в дальнейшем будет изложен для более общего
случая.
2.5. Методы решения разностных уравнений
Ниже рассматриваются основные методы решения линейных разностных
уравнений, при обосновании которых будет использоваться оценка скорости
роста решений, полученная в следующем пункте.
2.5.1. Оценка скорости роста решений. Здесь будет получена оценка скорости роста решений уравнения (2.2) и системы (2.5), необходимая для обоснования применения -преобразования к этим уравнениям. Будем рассматривать
систему разностных уравнений (2.5) с матрицей
и вектором
общего
вида. Таким образом, система (2.5) считается не связанной с уравнением (2.2).
элементы матрицы
, а через
– координаты векОбозначим через
тора
. Введем в рассмотрение векторную и матричную нормы равенствами
Лемма 2.2. Пусть коэффициенты матрицы
и координаты вектора
, тогда для любого решения разностной системы (2.5)
ограничены при
справедливо неравенство
где
55
Copyright ОАО «ЦКБ «БИБКОМ» & ООО «Aгентство Kнига-Cервис»
Доказательство. Рассматривая систему разностных уравнений (2.5) как реи выражая
последокуррентную формулу для определения вектора
вательно через
, получим
Оценивая последнее соотношение по норме, получим неравенство
Откуда с учетом обозначений (2.18) для
получаем оценку
Последнее неравенство легко преобразуется в (2.17). Лемма доказана.
Аналогичное утверждение справедливо и для решений уравнения (2.2). Для
этого достаточно заменить это уравнение эквивалентной системой вида (2.5)
так, как это сделано в начале главы, а затем воспользоваться результатом леммы 2.2. С учетом формул (2.4)-(2.6) получим для этого случая следующее утверждение.
Лемма 2.3. Пусть последовательности
и
ограничены
при
, тогда для любого решения разностного уравнения (2.2) справедливо
неравенство
где
2.5.2. Применение Z-преобразования. Одним из наиболее мощных (но не
самых простых) методов решения дифференциальных, разностных, дифференциально-разностных и некоторых других уравнений является операционный
метод. С помощью этого метода могут быть обоснованы и другие (например,
изложенные в этой главе) способы решения разностных уравнений. Рассмотрим
один из вариантов операционного метода для разностных уравнений, называемый методом Z-преобразования (соответствующим аналогом для систем с непрерывным временем является преобразование Лапласа).
56
Copyright ОАО «ЦКБ «БИБКОМ» & ООО «Aгентство Kнига-Cервис»
Предположим, что последовательность
определена при
и растет
не быстрее некоторой геометрической прогрессии. Тогда для достаточно больших по модулю комплексных можно определить функцию
равенством
Функция
аналитична вне некоторого круга с центром в начале координат
и может быть аналитически продолжена внутрь круга всюду, за исключением
некоторого числа точек, которые являются полюсами
. Оказывается, что
может быть восстановлена по функции
с помопоследовательность
щью формулы
где – контур, содержащий все особые точки
. Это утверждение называется теоремой об обращении [2]. Заметим, что интеграл в правой части последней
формулы может быть вычислен, например, с помощью вычетов.
Перейдем к решению разностных уравнений с помощью Z-преобразования.
Чтобы не загромождать изложение большими выкладками, ограничимся разностным уравнением второго порядка. Поскольку уравнение (2.15) является частным случаем уравнения (2.2), то из леммы 2.2 следует, что решения (2.2) растут
не быстрее некоторой геометрической прогрессии. Поэтому для отыскания решений уравнения (2.15) можно воспользоваться Z-преобразованием. Умножим
равенство (2.15) на
и просуммируем от 0 до , получим
Обозначим через
Z-преобразование последовательности
. Используя
формулы задержки (запаздывания), преобразуем суммы, стоящие в равенстве
(2.19)
57
Copyright ОАО «ЦКБ «БИБКОМ» & ООО «Aгентство Kнига-Cервис»
Отсюда и из формулы (2.19) получим для определения функции
ческое уравнение
алгебраи-
Из последнего соотношения получается формула для функции
Применим к обеим частям последнего равенства формулу обращения Zпреобразования, в результате получим
где в качестве контура интегрирования может быть выбрана окружность с центром в начале координат, содержащая все особенности подынтегральной функции в правой части (2.20) (отметим, что в качестве радиуса такой окружности
может быть выбрано число из леммы 2.3). Контурный интеграл в (2.20) может
быть вычислен с помощью вычетов [5]. При этом подынтегральная функция
имеет два простых полюса в точках
и , которые являются корнями характеристического многочлена, если эти корни различны, и полюс второго порядка
при
. Если и – различные корни, то, используя формулу
из равенства (2.20) получим
58
Copyright ОАО «ЦКБ «БИБКОМ» & ООО «Aгентство Kнига-Cервис»
Если корни и совпадают, то для вычисления интеграла в (2.20) используем
формулу для вычета функции в полюсе 2-го порядка
Это применительно к равенству (2.20) дает
Из приведенного анализа вытекает следующее утверждение.
Теорема 2.1. Пусть уравнение (2.14) имеет действительные корни
тогда общее действительное решение (2.13) имеет вид
,
где
– произвольные действительные постоянные. Если же
– комплексные числа, то общее решение (2.23) может быть приведено к форме
,
. Если характеристическое уравнение (2.16) имеет кратгде
ный корень , то общее решение разностного уравнения (2.15) может быть записано в виде
Для доказательства теоремы достаточно заметить, что формула (2.23) полу,
, если положить
чается из (2.21) при
Отсюда вытекает, что, ввиду произвольности начальных условий, постоянные
можно считать произвольными. Аналогично из равенства (2.22) получается формула (2.25). Если считать, что в формуле (2.23) комплексные числа, то,
используя тот факт, что (2.16) – уравнение с действительными коэффициента59
Copyright ОАО «ЦКБ «БИБКОМ» & ООО «Aгентство Kнига-Cервис»
ми, получаем, что
, – комплексно сопряженные числа. Полагая
, где
,
, из (2.23) получаем при
и
, что уравнение (2.15) имеет следующие комплексные элементарные решения:
Поскольку линейное уравнение (2.15) с действительными коэффициентами
имеет комплексные решения, то действительные и мнимые части также являются решениями (2.15). Поэтому уравнение (2.15) имеет линейно независимые
действительные решения вида
Отсюда следует, что общее решение уравнения (2.15) для случая комплексных
чисел имеет форму (2.24). Теорема 2.1 доказана.
2.5.3. Алгебраический метод решения однородных уравнений. Приведенный в предыдущем пункте метод отыскания решений линейных однородных разностных уравнений с постоянными коэффициентами второго порядка
нетрудно обобщить на соответствующие уравнения порядка m и системы разностных уравнений. Рассмотрим вначале уравнение m -го порядка
– постоянные, причем
где
(2.26) характеристическое уравнение
Пусть
– корни этого уравнения, a
ней. Рассмотрим систему последовательностей
. Введем соответствующее
– кратности этих кор-
Обобщением теоремы 2.1 является следующее утверждение.
Теорема 2.2. Пусть характеристическое уравнение (2.27) имеет корни
кратности
соответственно. Тогда каждая из последовательностей системы (2.28) является элементарным решением линейного разностного однородного уравнения (2.26), и общее решение этого уравнения может
быть записано в виде
60
Copyright ОАО «ЦКБ «БИБКОМ» & ООО «Aгентство Kнига-Cервис»
где
– элементы системы (2.28),
– произвольные постоянные.
Сформулируем теперь аналогичный общий результат для системы линейных однородных разностных уравнений вида
,
– m-мерный векгде постоянная квадратная матрица размеров
тор-столбец. Для системы (2.29) введем характеристическое уравнение
где – m-мерная единичная матрица. Соотношение (2.30) представляет собой
алгебраическое уравнение порядка m относительно переменной . Пусть известна полная система собственных и присоединенных векторов матрицы ,
которую запишем в виде
– собственные векторы, отвечающие собственным значениям
матрицы , а все остальные векторы системы (2.31) являются присоединенными к предыдущим векторам этой системы. Для простоты формулиотвечает нуровки предположим, что лишь первой серии векторов
левое собственное значение матрицы . Обозначим через
символ Кронекера
и рассмотрим систему векторных последовательностей
где
Справедливо следующее утверждение.
61
Copyright ОАО «ЦКБ «БИБКОМ» & ООО «Aгентство Kнига-Cервис»
Теорема 2.3. Пусть полная система собственных и присоединенных векто, опредеров матрицы , соответствующих собственным значениям
лена равенствами (2.31), причем
, а все остальные собственные значения
ненулевые. Тогда формулы (2.32) дают полную систему элементарных решений
разностной системы (2.30). Причем решения, отвечающие нулевому собственному значению, обращаются в нулевой вектор не более чем за шагов вправо и
продолжимы в сторону убывания времени не более чем на шагов. Общее ре, может быть запишение этой системы, задаваемое в начальный момент
сано в виде
определены равенстгде векторные последовательности
– произвольные постоянные.
вами (2.32), а
В заключение этого пункта отметим, что в случае, когда коэффициенты
уравнения (2.26) или системы (2.29) действительны, то общее решение также
может быть представлено в действительной форме. Процесс выделения действительных решений аналогичен рассмотренному в пункте 3.2.
2.5.4. Метод вариации произвольной постоянной. Перейдем к решению
неоднородных уравнений и систем. Как уже было отмечено выше, такие уравнения и системы возникают при моделировании процессов в цифровых устройствах под действием внешней силы. Из результатов пункта 1.1 вытекает, что
поскольку общее решение неоднородного уравнения может быть представлено
в виде суммы частного решения этого уравнения и общего решения соответствующего однородного, то, по крайней мере, для уравнения с постоянными коэффициентами достаточно научиться находить хотя бы одно решение неоднородного уравнения. Рассмотрим один из способов построения такого решения,
называемый методом вариации произвольной постоянной. Вначале рассмотрим
случай неоднородной системы (2.5) с постоянной матрицей общего вида и произвольной векторной последовательностью
. Пусть известна система элементарных решений (2.32) однородной системы (2.29). Определим фундаментальную матрицу решений (2.29) равенством
– элементарные решения вида (2.32)
Здесь векторы
соответствующей (2.5) однородной системы (2.29). Общее решение этой системы, определенное формулой (2.33), может быть записано в виде
. Считая, что все собстгде – вектор-столбец с координатами
венные значения матрицы отличны от нуля, получим, что система элементарных решений (2.32) линейно независима, и, следовательно, матрица – неособая. Будем искать решение неоднородной системы (2.6) в виде
62
Copyright ОАО «ЦКБ «БИБКОМ» & ООО «Aгентство Kнига-Cервис»
– новая неизвестная векторная последовательность. Подставляя (2.34)
где
в неоднородное уравнение (2.5), получим
Преобразуем последнее соотношение к виду
Учитывая, что столбцами матрицы
являются решения однородного
уравнения (2.29), приходим к уравнению для определения векторной последовательности
:
где
полагая
– матрица, обратная к
, получим
. Суммируя в равенстве (2.35) от
до
и
Таким образом, общее решение неоднородного уравнения (2.5) может быть записано в виде
Применим теперь метод вариации произвольной постоянной к неоднородному
разностному уравнению второго порядка
Для определенности будем считать, что соответствующее характеристическое
. В этом случае однородное уравнение
уравнение имеет различные корни
имеет общее решение вида (2.23). Метод вариации произвольной постоянной
состоит в том, что решение неоднородного уравнения (2.36) ищут в том же самом виде, что и решение однородного, но при этом произвольные постоянные
заменяют новыми неизвестными функциями дискретного аргумента . Другими словами, в уравнении (2.36) делают замену
где
,
– новые неизвестные последовательности. Вместо уравнения
(2.36) рассмотрим эквивалентную систему
63
Copyright ОАО «ЦКБ «БИБКОМ» & ООО «Aгентство Kнига-Cервис»
где
Решение системы (2.38) будем искать в виде
где
Заметим, что столбцами матрицы
являются решения однородной разностной системы. Подставляя (2.39) в уравнение (2.38), получим
Используя матричное равенство
преобразуем соотношение (2.40) к виду
слева, придем к формуле
Умножая последнее соотношение на
При этом матрица
имеет вид
Записывая формулу (2.41) в координатах, получим для неизвестных последовательностей
,
два разностных уравнения первого порядка
Решая эти уравнения и считая, что
,
64
заданы, получим
Copyright ОАО «ЦКБ «БИБКОМ» & ООО «Aгентство Kнига-Cервис»
Отсюда следует, что общее решение неоднородного уравнения (2.36) может
быть представлено в форме
В случае, когда характеристическое уравнение (2.16) имеет кратный корень,
решение неоднородного уравнения (2.36) ищется в виде
При этом, проводя вычисления, аналогичные рассмотренным выше, для определения
,
получим формулы
Общее решение уравнения (2.36) в этом случае имеет вид
2.5.5. Разностные уравнения с правой частью в виде квазиполинома.
Хотя рассмотренный в предыдущем пункте метод вариации произвольной постоянной в принципе решает вопрос об отыскании частного решения неоднородного разностного уравнения и системы таких уравнений, применение этого
метода во многих случаях сопряжено с громоздкими вычислениями. В связи с
этим рассмотрим еще один метод построения частного решения для неоднородных разностных уравнений с правыми частями специального вида
где
– произвольные числа,
– многочлены целого аргумента . Очевидно, что неоднородное уравнение с правой частью такого вида
можно решать "по частям", на каждом этапе отыскивая частное решение с неоднородностью вида
. Отметим, что неоднородности (2.42) встречаются
в приложениях особенно часто. Рассмотрим вначале случай неоднородного
разностного уравнения порядка вида
65
Copyright ОАО «ЦКБ «БИБКОМ» & ООО «Aгентство Kнига-Cервис»
Следующее утверждение позволяет отыскивать частное решение этого уравнения.
Теорема 2.4. Линейное неоднородное разностное уравнение (2.43) имеет частное решение вида
где
, если число не является корнем характеристического уравнения
– многочлен той же
(2.27) и s – кратность этого корня в противном случае,
с неизвестными коэффициентами, которые могут быть найстепени, что и
денными методом неопределенных коэффициентов.
Обоснование этого утверждения в общем случае достаточно громоздко. В
связи с этим продемонстрируем данный метод на примере линейного неоднородного разностного уравнения второго порядка
Соответствующее характеристическое уравнение
и
. Таким образом, частное решение (2.44), согласимеет корни
но теореме 2.4, нужно искать в виде
(число
Подставляя
является корнем характеристического уравнения кратности 1).
в соотношение (2.44), получим
Выделяя коэффициенты при одинаковых степенях n и решая полученную систему линейных уравнении, найдем
,
. Таким образом, общее
решение уравнения (2.44) имеет вид
.
Приведем формулировку утверждения, аналогичного теореме 2.4, применительно к системе разностных уравнений
где
– векторный многочлен степени .
Теорема 2.5. Система линейных неоднородных разностных уравнений (2.45)
имеет частное решение вида
66
Copyright ОАО «ЦКБ «БИБКОМ» & ООО «Aгентство Kнига-Cервис»
– векторный многочлен степени
, число имеет тот же смысл,
где
могут быть
что и в теореме 2.4, а векторные коэффициенты многочлена
найдены из уравнения (2.45) методом неопределенных коэффициентов. Заметим, что сформулированные в этом пункте утверждения переносятся на неоднородности более общего вида, дополнительно содержащие синусы и косинусы
аргумента . Сформулируем соответствующее утверждение для уравнения
Теорема 2.6. Линейное неоднородное уравнение (2.46) имеет частное решение
вида
, если число не является корнем характеристического уравнения
где
(2.27) и – кратность этого корня в противном случае,
– многочлен той
с неизвестными коэффициентами, которые могут быть
же степени, что и
определены из уравнения (2.46) методом неопределенных коэффициентов.
Следствием данной теоремы является условие разрешимости неоднородного разностного уравнения
в классе тригонометрических многочленов.
Теорема 2.7. Для того чтобы уравнение (2.47), где
имело
решение в виде тригонометрического многочлена, необходимо и достаточно,
чтобы среди корней характеристического уравнения (2.27) не было чисел вида
.
2.5.6. Применение функции Грина. Рассмотрим еще один метод отыскания частного решения неоднородных разностных уравнений, позволяющий
строить решения определенные и ограниченные при всех целых . Этот метод
опишем применительно к разностному уравнению второго порядка вида (2.36).
Пусть
– ограниченная последовательность. Для этого случая будет построено ограниченное решение уравнения (2.36) при условии, что характеристическое уравнение (2.16) не имеет корней, равных по модулю единице. Введем в рассмотрение функцию Грина
как решение уравнения (2.36) с правой частью
, где
– символ Кронекера. Найдем явный вид
для
различных случаев расположения корней характеристического уравнения
(2.16). Заметим, что определенная таким образом функция Грина удовлетворяет
системе уравнений
67
Copyright ОАО «ЦКБ «БИБКОМ» & ООО «Aгентство Kнига-Cервис»
,
Теорема 2.8. Пусть характеристическое уравнение (2.16) имеет корни
, тогда существует единственная функция Грина, отвечающая
причем
слеуравнению (2.36), которая определяется в зависимости от значений
дующими формулами
I.
,
,
II.
,
,
III.
,
,
IV.
,
,
V.
,
,
Для доказательства теоремы 2.8 достаточно непосредственно проверить, что
, определенная в этой теореме для различных случаев, удовлетвофункция
ряет условиям (2.48).
Следующее утверждение позволяет находить ограниченное частное решение уравнения (2.36).
68
Copyright ОАО «ЦКБ «БИБКОМ» & ООО «Aгентство Kнига-Cервис»
Теорема 2.9. Пусть
, тогда для любой ограниченной последовательности
существует единственное ограниченное решение уравнения
(2.36), которое определяется формулой
Доказательство. Прежде всего, заметим, что из формул (2.49) – (2.53) для
функции
вытекает, что ряд в правой части (2.54) сходится. Тот факт, что
формулой (2.54) действительно определяется решение, проверяется непосредственной подстановкой с учетом формул (2.48).
2.6. Устойчивость решений разностных уравнений
При исследовании разностных уравнений важным является вопрос о поведении решений при неограниченном возрастании аргумента, и в частности вопрос об устойчивости решений, тесно связанный с физической реализуемостью.
Особенно важным этот вопрос является применительно к стационарным и периодическим решениям. Сформулируем определение устойчивости (см. также
[3]) для разностных уравнений вида
Определение 2.1. Решение
разностного уравнения (2.55), определенное при
, называется устойчивым по Ляпунову, если для любого
, что для любого другого решения
уравнения
существует такое
(2.55), удовлетворяющего условию
для всех
выполнено неравенство
Определение 2.2. Решение
уравнения (2.55) называется асимптотически устойчивым по Ляпунову, если оно устойчиво, по Ляпунову и существует
такое
, что для любого решения
, удовлетворяющего условию
справедливо равенство
69
Copyright ОАО «ЦКБ «БИБКОМ» & ООО «Aгентство Kнига-Cервис»
Определение 2.3. Решение
уравнения (2.55) называется экспоненциально устойчивым по Ляпунову, если оно устойчиво по Ляпунову и существуют такие положительные постоянные и , причем
, что для некоторого
для любого решения
, удовлетворяющего условию
при
следует неравенство
Приведенные определения легко переносятся на случай систем разностных
уравнений. Сформулируем, например, определение устойчивости решений линейной разностной системы
где
) – квадратная
матрица,
и
– -мерные векторы.
системы разностных уравнений (2.57) наОпределение 2.4. Решение
существует такое
, что для
зывается устойчивым, если для любого
системы (2.57), удовлетворяющего условию
любого решения
при всех
выполнено неравенство
обозначена некоторая векторная норма в -мерном простран(Здесь через
стве.) Следующие утверждения аналогичны хорошо известным из теории устойчивости линейных систем дифференциальных уравнений.
Теорема 2.10. Для того чтобы любое решение линейной разностной системы
(2.57) было устойчиво (асимптотически устойчиво), необходимо и достаточно,
чтобы нулевое решение соответствующей однородной системы было устойчиво
(асимптотически устойчиво).
Сформулированное утверждение позволяет сводить вопрос об устойчивости
произвольного решения линейной разностной системы уравнений (2.57) (как
однородной, так и неоднородной) к исследованию устойчивости нулевого решения соответствующей однородной системы. В связи с этим дадим еще одно
определение.
Определение 2.5. Система линейных разностных уравнений (2.57) называется устойчивой (асимптотически устойчивой), если все решения этой системы
устойчивы (асимптотически устойчивы).
Из сказанного выше вытекает, что система (2.57) будет устойчивой (асимптотически устойчивой), если нулевое решение соответствующей однородной
системы будет устойчиво.
70
Copyright ОАО «ЦКБ «БИБКОМ» & ООО «Aгентство Kнига-Cервис»
2.6.1. Устойчивость уравнений с постоянными коэффициентами. Перейдем к вопросу об устойчивости решений линейных разностных уравнений с
постоянными коэффициентами вида
и системы разностных уравнений
с постоянной матрицей . При этом из утверждений предыдущего пункта ясно,
что достаточно ограничиться случаем однородных разностных уравнений и
систем (2.58), (2.59). Кроме этого, результаты пункта 2.5 данной главы позволяют в явном виде получить общие решения уравнения (2.58) и системы (2.59).
Вид этих решений устанавливается теоремами 2.2, 2.3. Сказанное позволяет
сформулировать следующие утверждения, вытекающие из самого вида общих
решений для разностного уравнения (2.58) и системы (2.59).
Теорема 2.11. Для того чтобы разностное уравнение (2.58) было устойчиво,
необходимо и достаточно, чтобы все корни соответствующего характеристического уравнения по модулю не превосходили единицу, а корни с модулем, равным единице, были простые.
Теорема 2.12. Для того чтобы разностное уравнение (2.58) было асимптотически устойчиво, необходимо и достаточно, чтобы все корни соответствующего
характеристического уравнения были по модулю меньше единицы.
Теорема 2.13. Для того чтобы система уравнений (2.59) была устойчива, необходимо и достаточно, чтобы все корни соответствующего характеристического уравнения по модулю не превосходили единицу, а геометрическая кратность корней с модулями, равными единице была бы равна единице.
Теорема 2.14. Для того чтобы система разностных уравнений (2.59) была
асимптотически устойчивой, необходимо и достаточно, чтобы все корни характеристического уравнения были по модулю меньше единицы.
Таким образом, для линейных разностных уравнений с постоянными коэффициентами приведенные утверждения сводят вопрос об устойчивости решений к исследованию расположения корней соответствующего характеристического уравнения. Задача исследования расположения корней таких уравнений
для разностных уравнений и систем порядка не выше третьего решается просто.
Для уравнений и систем более высокого порядка исследование может быть
проведено с помощью алгебраических критериев типа Рауса-Гурвица и частотных типа Найквиста-Михайлова.
2.6.2. Алгебраические методы анализа устойчивости. В этом пункте будет рассмотрен способ определения расположения корней многочлена
71
Copyright ОАО «ЦКБ «БИБКОМ» & ООО «Aгентство Kнига-Cервис»
относительно единичной окружности в комплексной плоскости переменной .
Чтобы свести эту задачу к известной проблеме определения расположения корней многочлена относительно мнимой оси, выполним в уравнении (2.60) замену
переводящую внутренность единичного круга комплексной плоскости в левую полуплоскость комплексной плоскости . В результате вместо (2.60) получим уравнение
После раскрытия скобок и приведения подобных членов придем к соотношению
где коэффициенты
линейно выражаются через
. Из
вида замены (2.61) следует, что корни уравнения (2.60), лежащие внутри единичной окружности, переходят в корни (2.62), лежащие в левой полуплоскости,
а корни (2.60), лежащие вне единичной окружности, переходят в корни (2.62),
лежащие в правой полуплоскости. Сформулируем теперь критерий РаусаГурвица применительно к многочлену (2.62). Для этого введем матрицу Гурвица
и определим главные угловые миноры этой матрицы равенствами
Теорема 2.15. Для того чтобы все корни многочлена (2.62) имели отрицательные действительные части, необходимо и достаточно, чтобы все главные
миноры матрицы Гурвица были положительными.
Следует отметить, что условия положительности всех главных миноров
матрицы Гурвица при условии положительности коэффициентов многочлена
(2.62) не являются независимыми, поскольку из положительности миноров четного порядка в этом случае следует положительность миноров нечетного порядка и наоборот. В связи с этим можно сформулировать критерий Льенарда-
72
Copyright ОАО «ЦКБ «БИБКОМ» & ООО «Aгентство Kнига-Cервис»
Шипара, использующий меньшее число проверяемых детерминантных условий.
Теорема 2.16. Для того чтобы все корни многочлена (2.62) имели отрицательные действительные части, необходимо и достаточно, чтобы все коэффициенты этого уравнения были положительны и выполнялись неравенства
.
Используя сформулированные утверждения, можно получить коэффициентные условия устойчивости решений уравнения (2.58) и системы (2.59). Для
характеристического уравнения (2.60) эти условия имеют следующий вид ([6]);
при
:
при
при
:
:
Таким образом, теоремы 2.6, 2.7. полностью решают вопрос об устойчивости
решений линейных разностных уравнений и систем, причем в случае уравнения
или системы порядка меньше чем четыре, можно использовать коэффициентные условия (2.63)–(2.65).
2.6.3. Частотные методы анализа устойчивости. Наряду с алгебраическими методами исследования устойчивости на практике широко используют и
частотные методы. Преимущества последних обычно проявляются при исследовании уравнений и систем высокого порядка. Рассмотрим один из вариантов
частотного метода, являющегося аналогом критерия Михайлова. Выполним в
уравнении (2.60) замену
. В результате этого уравнение преобразуется к
виду
Непосредственно из принципа аргумента [5] вытекает следующее утверждение.
Теорема 2.17. Для того чтобы все корни уравнения (2.60) были по модулю
меньше единицы, необходимо и достаточно, чтобы угол поворота против хода
73
Copyright ОАО «ЦКБ «БИБКОМ» & ООО «Aгентство Kнига-Cервис»
часовой стрелки вектора
при изменяющемся от до был бы равен
, а сам вектор
не обращался бы в нулевой.
Приведем другую формулировку частотного критерия, основанную на анализе поведения действительной и мнимой частей функции
. Положим
и
– действительная и мнимая части функции
. Из теоремы
где
2.8 вытекает следующее утверждение.
Теорема 2.18. Для того чтобы все корни уравнения (2.60) были по модулю
меньше единицы, необходимо и достаточно, чтобы функции
и
имеперемежающихся нулей.
ли на отрезке
При практическом применении теоремы 2.17 следует проследить за изменением вектора
при
и определить угол, на который повернется
этот вектор. Число корней характеристического уравнения (2.60), находящихся
внутри единичного круга, можно определить, разделив найденный угол на .
Для использования результатов теоремы 2.18 нужно построить графики действительной и мнимой частей функции
и проследить за чередованием нулей этих функций.
2.6.4. Теорема Флоке-Ляпунова и устойчивость решений. Рассмотрим
линейное однородное уравнение с периодическими коэффициентами
будем считать определенными для всех
где последовательности
для некоторого натуральцелых и такими, что а
для всех целых . Наряду с уравнением (2.67)
ного числа , причем
рассмотрим систему разностных уравнений
матрицы
определены для всех целых и
где элементы квадратной
. Кроме этого, будем предпериодичны с периодом так, что
для всех . В предыдущей главе было показано, что
полагать, что
уравнение (2.67) может быть сведено к некоторой системе вида (2.68) с матрицей специального вида. В связи с этим получим сначала общий результат для
системы (2.68). Обозначим через
линейно независимую
. Таким
систему решений (2.68) и положим
являются линейно независиобразом, столбцами построенной матрицы
мые решения системы (2.68). Отсюда следует, что для любого целого
. Матрицу
будем называть фундаментальной матрицей ре– неособая, то
шений системы разностных уравнений (2.68). Поскольку
. Сформулируем тедля любого существует обратная к ней матрица
перь утверждение, устанавливающее общий вид решения системы (2.68).
74
Copyright ОАО «ЦКБ «БИБКОМ» & ООО «Aгентство Kнига-Cервис»
Теорема 2.19 (Флоке-Ляпунова). Общее решение однородной системы липериода
нейных разностных уравнений (2.68) с периодической матрицей
может быть представлено в виде
– периодическая матрица периода
где – некоторая постоянная матрица,
, – произвольный постоянный вектор.
Доказательство. Рассмотрим соответствующее (2.68) матричное уравнение
Из самого определения фундаментальной матрицы
следует, что
также удовлетворяет
удовлетворяет (2.70). Покажем, что матрица
в систему (2.70), приходим к тождест(2.70). Действительно, подставляя
, придем
ву. Выбирая в нем значение независимого переменного равным
к формуле
Используя периодичность матрицы
, из последнего равенства получим
также удовлетворяет системе (2.70) и, следоваОтсюда вытекает, что
является фундаментальной матрицей. Таким образом, для нетельно,
которой постоянной неособой матрицы выполнено соотношение
Полагая в последнем равенстве
рицы
, получим формулу для определения мат-
нормированная матрица (
), получаем
В частном случае, когда
. Введем теперь в рассмотрение матрицу
равенством
Докажем, что матрица
– периодическая матрица. Действительно
75
Copyright ОАО «ЦКБ «БИБКОМ» & ООО «Aгентство Kнига-Cервис»
Вводя обозначение
и выражая фундаментальную матрицу
отношения (2.72), получаем
из со-
Такое представление фундаментальной матрицы решений эквивалентно формуле (2.69). Теорема доказана.
Матрица , определенная соотношением (2.71), называется матрицей монодромии системы (2.68). Из теоремы 2.19 вытекает, что свойства устойчивости
решений этой системы определяются собственными значениями матрицы .
Для отыскания этой матрицы можно построить фундаментальную систему решений для значений
от до
и затем воспользоваться формулой (2.71).
Применительно к уравнению (2.67) полученный результат может быть сформулирован в виде следующего утверждения.
Теорема 2.20. Любое решение линейного однородного разностного уравнения (2.67) с периодическими коэффициентами периода может быть записано
в виде
– некоторые постоянные,
– периодичегде
ские последовательности периода .
Сформулируем теперь утверждение об устойчивости решений разностных
уравнений и систем с периодическими коэффициентами. Будем называть собственные значения матрицы монодромии системы разностных уравнений (2.68)
мультипликаторами этой системы.
Теорема 2.21. Для того чтобы решения линейного разностного уравнения
(2.67) (системы (2.68)) с периодическими коэффициентами были устойчивыми,
необходимо и достаточно, чтобы все мультипликаторы этого уравнения (системы) не превосходили по модулю единицу, а мультипликаторам с модулями,
равными единице, отвечали бы лишь собственные векторы матрицы монодромии.
Теорема 2.22. Для того чтобы решения линейного разностного уравнения
(2.67) (системы (2.68)) с периодическими коэффициентами были асимптотически устойчивыми, необходимо и достаточно, чтобы все мультипликаторы этого
уравнения (системы) по модулю были меньше единицы.
2.6.5. Устойчивость решений уравнений с малым параметром. В этом
пункте будет рассмотрен алгоритм исследования устойчивости решений разностных уравнений с малым параметром в критическом случае, аналогичный известному для дифференциально-разностных уравнений [7,8]. Рассмотрим разностное уравнение
76
Copyright ОАО «ЦКБ «БИБКОМ» & ООО «Aгентство Kнига-Cервис»
где коэффициенты
,
уравнения имеют вид
),
причем элементами рядов
,
,
являются тригонометрические многочлены и эти ряды абсолютно сходятся равномерно относительно при достасоотношение (2.73) переходит в уравнеточно малых по модулю . При
ние с постоянными коэффициентами
характеристический многочлен (2.74).
Обозначим через
Если корни
по модулю меньше единицы, то решения уравнения (2.74)
есть хотя бы один корень, по модуасимптотически устойчивы. Если же у
лю превосходящий единицу, то решения уравнения (2.74) неустойчивы. При
этом оказываются устойчивыми и неустойчивыми, соответственно, и решения
уравнения (2.73) для малых по модулю . Таким образом, в случае, когда у характеристического многочлена корни по модулю меньше единицы или есть хотя бы один корень по модулю больше единицы, вопрос об устойчивости решений уравнения (2.73) при малых по модулю решается просто. Фактически в
этом случае он эквивалентен вопросу об устойчивости решений соответствующего уравнения с постоянными коэффициентами (2.74). Из сказанного ясно,
что остается открытым вопрос об устойчивости решений уравнения (2.73) в
случае, когда оба корня характеристического многочлена не превосходят по
модулю единицу и хотя бы один из корней имеет модуль равный единице. Такие случаи для уравнения (2.73) называются критическими. Для исследования
устойчивости решений уравнения (2.73) в критических случаях воспользуемся
алгоритмом, предложенным в [7,8] для дифференциально-разностных уравнений. Рассмотрим этот алгоритм в критическом случае простого единичного
корня. Введем в рассмотрение формальные ряды
где
– некоторые числа,
– тригонометрические полиномы. Будем искать решение уравнения (2.73) в виде
Подставляя
в соотношение (2.73), получим
77
Copyright ОАО «ЦКБ «БИБКОМ» & ООО «Aгентство Kнига-Cервис»
Нетрудно видеть, что, выделяя в уравнении (2.77) коэффициенты при различных степенях начиная с первой и используя условия разрешимости получающихся уравнений в классе тригонометрических полиномов, можно последовательно определять элементы рядов (2.75), (2.76) любого порядка. При этом вопрос об устойчивости может быть решен с помощью следующего утверждения.
Теорема 2.23. Пусть согласно описанному алгоритму вычислены коэффициенты рядов (2.75), (2.76) до номера
включительно, причем
,
, тогда решения уравнения (2.73) устойчивы при достаточно
малых по модулю
, если
, и неустойчивы, если
.
Обоснование этого утверждения может быть проведено в идейном плане
работ [7,8]. Ввиду сложности доказательства, приводить его здесь не будем.
Аналогично может быть исследована устойчивость решений уравнения (2.73) в
других критических случаях.
2.6.6. Исследование устойчивости по первому приближению. Рассмотрим нелинейное разностное уравнение второго порядка
– дифференцируемая функция переменных
. Пусть это уравнегде
(состояние равновесия).
ние имеет стационарное решение
. В результате вместо (2.78) полуВыполним в (2.78) замену
чим следующее нелинейное уравнение
a
– бесконечно малая функция, имеюгде
щая в нуле порядок малости выше первого. Последнее означает, что выполнено
следующее предельное соотношение
Наряду с нелинейным разностным уравнением (2.79) рассмотрим линейное
уравнение вида
Соотношение (2.80) будем называть уравнением первого приближения или линеаризацией уравнения (2.78) на состоянии равновесия . Можно показать, что
свойства асимптотической устойчивости нулевого решения уравнений (2.78) и
(2.80) одинаковы. Это означает справедливость следующего утверждения.
Теорема 2.24. Пусть
– состояние равновесия уравнения (2.78). Это состояние равновесия будет асимптотически устойчиво в том и только том случае, если асимптотически устойчиво нулевое решение линеаризованного урав78
Copyright ОАО «ЦКБ «БИБКОМ» & ООО «Aгентство Kнига-Cервис»
нения (2.80). Другими словами, состояние равновесия устойчиво, если оба корня характеристического уравнения, соответствующего (2.80), по модулю меньше единицы. Если же хотя бы один из корней этого уравнения по модулю
больше единицы, то состояние равновесия уравнения (2.78) неустойчиво.
Отметим, что последнее утверждение без труда переносится на случай нелинейных разностных уравнений и систем произвольного порядка.
2.6.7. Устойчивость уравнений с почти периодическими коэффициентами. В этом пункте будет исследована устойчивость решений разностного
уравнения
где
– постоянные, и построены области устойчивости в плоскости переменных
в зависимости от значений параметров и . Следует отметить,
что в случае, когда частота рационально соизмерима с числом
(периодический случай), область устойчивости оказывается зависящей от фазы . В
связи с этим для каждого фиксированного полезно выделить общую для всех
область устойчивости решений уравнения (2.81). В случае, когда и
рационально несоизмеримы (почти периодический случай), оказывается, что общая для всех
область устойчивости решений (2.81) является в некотором
смысле предельной для областей устойчивости в периодическом случае.
Для исследования устойчивости решений уравнения (2.81) перейдем от этого уравнения к эквивалентному в смысле устойчивости соотношению
Положим
тогда уравнение (2.82) может быть записано в виде
Ниже будет подробно исследовано поведение многочлена
, полностью определяющего свойства устойчивости решений уравнения (2.81), отдельно в периодическом и почти периодическом случаях.
2.6.7.1. Периодический случай. Пусть , – натуральные взаимно простые
числа. Тогда из формулы (2.83) вытекает, что многочлен
не
зависит от и что многочлены
и
отличаются лишь порядком сомножителей и поэтому равны. В связи с этим всюду
ниже в данном пункте будем считать
и многочлен (2.83) будем обо79
Copyright ОАО «ЦКБ «БИБКОМ» & ООО «Aгентство Kнига-Cервис»
значать
пеням
через
. Наряду с представлением многочлена
в виде произведения (2.83) рассмотрим его разложение по сте-
При этом для коэффициента
согласно [9], справедлива формула
Теорема 2.25. Коэффициенты
зависят от и определяются равенствами
где через
обозначена целая часть числа.
Из формул (2.81), (2.84) следует, что при
уравнения (2.81) справедливо представление
из разложения (2.85) не
для мультипликатора
и решения (2.81) устойчивы при
. Поэтому область устойчивости решений уравнения (2.81) –
, отвечающая
в плоскости параметров
, определяется соотношениями
При этом общая для всех значений область устойчивости решений уравнения (2.81) –
определяется условиями
Заметим, что в силу формул (2.85) – (2.88) при каждом фиксированном максимум по
в левой части (30) достигается либо при
, либо при
, что позволяет упростить построение областей устойчивости. Отметим
еще, что в силу соотношений (2.89), (2.90) области устойчивости симметричны относительно начала координат и поэтому достаточно ограничиться по80
Copyright ОАО «ЦКБ «БИБКОМ» & ООО «Aгентство Kнига-Cервис»
строением этих областей лишь в полуплоскости
. Вычисленные по формулам (2.86) – (2.88) многочлены
для
имеют вид:
а построенные по ним области устойчивости
значений в полуплоскости
приведены на рис. 2.3.
для некоторых
Рис. 2.3.
При этом непосредственно из формулы (2.81) вытекает, что на прямых
решения уравнения (2.81) суперустойчивы. Области устойчивости
для
показаны на рис. 2.4.
81
Copyright ОАО «ЦКБ «БИБКОМ» & ООО «Aгентство Kнига-Cервис»
Рис. 2.4.
Следующее утверждение устанавливает существование предельной области , являющейся объединением областей
. Причем из последовательности
можно выделить монотонно "возрастающие" подпоследовательности областей
и т.д.), "схо(например,
дящиеся" к области .
Теорема 2.26. Для любого действительного числа
для многочлена
справедливо предельное равенство
Доказательство. Пусть сначала
Положим
где коэффициенты
гда многочлен
– фиксированное число из отрезка
определены равенствами (2.87), (2.88). Томожет быть представлен в виде
причем для коэффициента
справедлива формула (2.86), а
висит от . Из равенства (2.83) следует, что
для некоторого
.
. Отсюда и из формулы (2.92) получим
82
не за-
Copyright ОАО «ЦКБ «БИБКОМ» & ООО «Aгентство Kнига-Cервис»
Поэтому многочлен
может быть представлен в виде
Из данной формулы и равенства (2.86) при
вытекает оценка
Отсюда следует справедливость формулы (2.91) при
, тогда из равенства (3) при
получим
. Пусть теперь
Данное выражение представляет собой интегральную сумму для интеграла
Отсюда и из периодичности подынтегральной функции получим (см. [9]) при
предельное равенство
которое эквивалентно (2.91) при
. Теорема доказана.
Из данного утверждения и соотношения (2.90) следует, что предельная область устойчивости G определяется условиями
Исключая из них параметр , получим неравенства
связывающие координаты
точек области
(для
) показана на рис. 2.5.
83
. Cама предельная область
Copyright ОАО «ЦКБ «БИБКОМ» & ООО «Aгентство Kнига-Cервис»
Рис. 2.5
2.6.7.2. Почти периодический случай. Пусть в уравнении (2.81) частота
несоизмерима рационально с . Для любого натурального представим число
в виде
где – натуральное число,
. Из [10] следует, что числовая последовательность
всюду плотна на отрезке
и, более того, является
равномерно распределенной на этом отрезке. Последнее означает, что если через
обозначить число точек
, лежащих на некотором отрезке
длины , принадлежащем отрезку
, то справедливо предельное равенство
Построение области устойчивости решений уравнения (2.81) в рассматриваемом случае основывается на следующем утверждении.
Теорема 2.27. Пусть функция
интегрируема на отрезке
, хотя
бы и в несобственном смысле, а последовательность
равномерно распределена на
и ее элементы не совпадают ни с одной из особых точек
, тогда
Пусть – произвольное действительное число, тогда следствием теоремы
2.27 является предельное равенство
84
Copyright ОАО «ЦКБ «БИБКОМ» & ООО «Aгентство Kнига-Cервис»
где
– многочлен (2.83), причем если
ется таким образом, что
, то число
выбира-
Из равенства (2.96) и формулы
(см. [9]) вытекает соотношение
аналогичное (2.91). При этом область асимптотической устойчивости решений
уравнения (2.81)
определяется условиями
а общая для всех
область асимптотической устойчивости этого уравнения
, являющаяся пересечением областей
, удовлетворяет соотношениям
Отсюда и из формулы (2.97) получим, что область
, как и предельная область , построенная в предыдущем пункте, определяется формулами (2.95).
Отметим еще, что для любых фиксированных и существует счетное, всюду
плотное на отрезке
множество чисел
такое, что на прямых
в плоскости параметров
решения уравнения
(2.81) суперустойчивы. Таким образом, для фиксированных и область устойчивости решений уравнения (2.81) состоит из объединения области
,
совпадающей с областью , представленной на рис. 2.5, и области суперустойчивости , показанной на этом же рисунке и состоящей из прямых
, которые всюду плотно заполняют конус
в плоскости параметров
.
2.7. Нелинейные разностные уравнения и хаос
Рассмотренные выше линейные разностные уравнения демонстрировали
сравнительно простой вид решений, которые определяли достаточно простую в
целом динамику моделируемых систем. В то же время известно, что многие явления, наблюдаемые в природе и в физических экспериментах, показывают
очень сложное поведение. В огромном большинстве случаев это связано с тем,
85
Copyright ОАО «ЦКБ «БИБКОМ» & ООО «Aгентство Kнига-Cервис»
что данные явления уже не могут быть адекватно описаны линейными уравнениями. Ниже будет рассмотрена простая нелинейная детерминированная модель динамической системы, обнаруживающая сложное поведение решений.
Большинство явлений природы нелинейны по своей сути, однако во многих
случаях они адекватно могут быть описаны хорошо выбранными линейными
математическими уравнениями. С другой стороны, при определенных условиях
те же явления могут вдруг начать проявлять свойства непредсказуемого поведения, отсутствия порядка или, как еще говорят, хаоса. Естественно было бы
предположить, что в этом случае в рассматриваемый процесс вмешиваются дополнительные факторы, оказывающие случайное воздействия на какие-либо
величины, определяющие поведение системы. Однако в действительности это
не так. Рассматриваемое ниже детерминированное нелинейное разностное
уравнение при определенных условиях демонстрирует сложное, хаотическое
поведение решений, называемое динамическим хаосом. При этом решения системы проявляют существенную зависимость от начальных условий. Дополнительной особенностью хаотических систем является появление странных аттракторов – областей притяжения с фрактальной структурой.
2.7.1. Логистическое уравнение. Ниже будет рассматриваться простейшее
нелинейное разностное уравнение
где
называемое логистическим отображением. Хотя уравнение (2.98) может быть
использовано в качестве математической модели самых разных по своей природе явлений (нелинейный цифровой фильтр, рост банковских сбережений и
т. д.) обычно при интерпретации (2.98) используется введенная еще в 1845 году
П. Ф. Ферхюльстом модель для описания динамики популяции в замкнутой
среде. При этом
трактуется как нормированная численность особей в
-й год, которая пропорциональна численности в предыдущий год, а так.
же пропорциональна свободной части жизненного пространства
Значение параметра зависит от плодовитости, свойств ареала обитания и т. д.
Для того чтобы численность популяции никогда не была отрицательной, параметр должен принимать значения от до . Биологическая трактовка возможного поведения решений уравнения (2.98) в зависимости от величины параметра сводится к следующим вариантам:
1. Популяция может вымереть:
.
2. Численность популяции может стабилизироваться на некотором стационарном значении:
.
3. Численность может периодически меняться через определенное число лет,
т. е.
. Здесь – период изменения численности.
4. Численность популяции может меняться хаотически.
86
Copyright ОАО «ЦКБ «БИБКОМ» & ООО «Aгентство Kнига-Cервис»
Отметим, что само уравнение (2.98) при этом может иметь одновременно стационарные, периодические и хаотические решения.
2.7.2. Состояния равновесия логистического уравнения. Найдем состояния равновесия уравнения (2.98). Подставляя в него
и используя формулу (2.99), получим
Это уравнение имеет два решения
Следует отметить, что, в силу требования неотрицательности , это стацио. Из теоремы об устойчивонарное решение существует лишь при
сти по первому приближению (см. теорему 2.24) следует, что для выяснения
вопроса об асимптотической устойчивости решений
и
достаточно найти
условия, при которых производная функции
в точках
и
по модулю
меньше единицы. Эти производные имеют вид:
уравнение (2.98) имеет лишь нулевое соОтсюда следует, что при
, которое асимптотически устойчиво. При
стояние равновесия
увеличении параметра это состояние равновесия становится неустойчивым и
у уравнения (2.98) появляется второе состояние равновесия
, которое определено для
и устойчиво при
.
изменения параметра принято называть окном устойчивоИнтервал
. Аналогично окном устойчивости стациости состояния равновесия
нарного решения
является интервал
. Отметим
производная
положительна и решения, отеще, что при
монотонно. Если же
личные от состояний равновесия, приближаются к
, то
отрицательна, и решения, отличные от состояния равновесия, приближаясь к , принимают то большие, то меньшие, чем , значерешение
теряет устойчивость колебания. При
тельным образом, т. е. у уравнения (2.98) появляется периодическое решение
(цикл) периода
.
2.7.3. Бифуркации удвоения периода и переход к хаосу. Формально потерю устойчивости состояния равновесия
, которое может рассматри, и появление периодичеваться как периодическое решение с периодом
, описанное в предыдущем пункте, при увелиского решения с периодом
чении параметра можно назвать бифуркацией удвоения периода. При этом
может быть найден из уравнения
появляющийся цикл периода
87
Copyright ОАО «ЦКБ «БИБКОМ» & ООО «Aгентство Kнига-Cервис»
в котором нужно считать
ние с формулой (2.99), придем к равенству
. Сопоставляя последнее уравне-
и, следуя [11], преобразуем его к виду
Откуда получим формулу
Решениями последнего уравнения, помимо уже найденных раньше стационаров
и , являются также числа и , определенные формулой
Отсюда следует, что периодическая последовательность, составленная из чередующихся чисел
и , является решением уравнения (2.98). Это периодиче. Ясно, что вместе с решением
ское решение удобно обозначить
. Этим двум
уравнение (2.98) имеет еще одно периодическое решение
, состоящий из двух точек
последовательностям отвечает цикл периода
, который будем обозначать
. Для исследования устойчивости этого
цикла нужно воспользоваться условием
которое применительно к данному случаю приводит к неравенству
Отсюда следует, что окно устойчивости для цикла
лено неравенствами.
88
периода
опреде-
Copyright ОАО «ЦКБ «БИБКОМ» & ООО «Aгентство Kнига-Cервис»
При этом цикл
притягивающийся и все траектории уравнения (2.98), за исключением конечного числа слипающихся со стационаром
притягиваются
циклом
. Последующие бифуркации удвоения периода происходят по тому
же сценарию.
Пусть интервал
– окно устойчивости цикла периода , тогда, как
из этого интервала, при котором этот
установлено в [11], найдется число
цикл суперустойчив. Более того, в [11] показано, что при возрастании параметра до значения
циклически повторяются бифуркации удвоения периода,
сопровождающиеся следующими событиями:
1. При
cтановится неустойчивым предыдущий цикл периода
и
рождается цикл периода .
2. При
цикл периода становится суперустойчив.
3. При
цикл периода cтановится неустойчивым и рождается цикл
.
периода
При этом для границ интервалов, определяющих окна сходимости, справедливы формулы
В работе [11] установлена общая формула связывающая границы окон устойчивости
из которой нетрудно получить ограниченность последовательности чисел .
Отсюда и из монотонного возрастания этой последовательности вытекает, что
она имеет предел , который может быть найден из последней формулы, если
в ней число устремить к бесконечности. Осуществляя этот предельный переход, получим
Отсюда получим приближенное значение
При этом существует предел
89
Copyright ОАО «ЦКБ «БИБКОМ» & ООО «Aгентство Kнига-Cервис»
Постоянная называется универсальной константой М. Фейгенбаума [12], который открыл ряд универсальных закономерностей для систем вида (2.98).
Проведенные в [12] расчеты для постоянных
и дали следующие значения
определяет для уравнения (2.98) верхнюю границу параПри этом величина
метра , при превышении которой все циклы, описанные выше, становятся неустойчивыми, причем существует множество решений, которые слипаются с
этими циклами. В то же время появляются решения, поведение которых становится хаотическим.
Литература
1. Гельфонд А. О. Исчисление конечных разностей. М.: Наука, 1967. 376 с.
2. Ципкин Н. З. Теория линейных импульсных систем. М.: Физматгиз, 1963.
968 с.
3. Мартынюк Д. И. Лекции по качественной теории разностных уравнений.
Киев: Наукова думка, 1972. 248 с.
4. Рабинер Л., Гоулд Б. Теория и применение цифровой обработки сигналов.
М.: Мир, 1976. 848 с.
5. Лаврентьев М. А., Шабат Б. В. Методы теории функции комплексного переменного. М.: Наука, 1973. 736 с.
6. Цыпкин Я. З. Основы теории автоматических систем. М.: Наука, 1977.
559 с.
7. Кащенко С. А. Об одном алгоритме исследования устойчивости решений
линейных дифференциальных уравнений n-го порядка с близкими к постоянными почти периодическими коэффициентами // Исследования по устойчивости и теории колебаний: сб. Ярославль, 1977. С. 60–70.
90
Copyright ОАО «ЦКБ «БИБКОМ» & ООО «Aгентство Kнига-Cервис»
8. Биркган С. Е. Экспоненциальная дихотомия и устойчивость решений линейных дифференциально-разностных уравнений нейтрального типа // Исследования по устойчивости и теории колебаний: сб. Ярославль, 1981. С. 90–112.
9. Прудников А. П., Брычков Ю. А., Маричев О. И. Интегралы и ряды. М.,
1981. 800 с.
10. Федорюк М. В. Обыкновенные дифференциальные уравнения. М., 1980.
11. Бурд В. Ш. Введение в динамику одномерных отображений. Ярославль:
ЯрГУ, 2006. 104 с.
12. Фейгенбаум М. Универсальное поведение нелинейных систем // УФН. 1983.
Т. 141, вып. 3. С. 343–374.
91
Copyright ОАО «ЦКБ «БИБКОМ» & ООО «Aгентство Kнига-Cервис»
Учебное издание
Биркган Сергей Ефимович
Математическое
моделирование
Учебное пособие
Редактор, корректор М. В. Никулина
Компьютерный набор, верстка С. Е. Биркган
Подписано в печать 25.01.13. Формат 60х84/16.
Гарнитура «TimesNewRoman». Бумага офсетная.
Усл. печ. л. 5,35. Уч.-изд. л. 5,0. Тираж 35 экз. Заказ
Оригинал-макет подготовлен
в редакционно-издательском отделе
Ярославского государственного университета
им. П. Г. Демидова
150000 Ярославль, ул. Советская, 14.
Отпечатано на ризографе
Ярославский государственный университет
им. П. Г. Демидова
Документ
Категория
Без категории
Просмотров
11
Размер файла
4 537 Кб
Теги
моделирование, математические, 9894
1/--страниц
Пожаловаться на содержимое документа