автореферат диссертации по информатике, вычислительной технике и управлению, 05.13.18, диссертация на тему:Разработка двухстадийных схем Розенброка с комплексными коэффициентами и их применение в задачах моделирования образования периодических наноструктур

На правах рукописи

ЛИМОНОВ Александр Георгиевич

П * П 004616950

Разработка двухстадииных схем Розенорока с комплексными коэффициентами и их применение в задачах моделирования образования периодических

наноструктур

05.13.18 - Математическое моделирование, численные методы и комплексы программ

АВТОРЕФЕРАТ

диссертации на соискание учёной степени кандидата физико-математических наук

- 9 ДЕК 2010

Екатеринбург — 2010

004616950

Работа выполнена на кафедре Мультимедиа технологий, ФГАОУ ВГЮ "Уральский федеральный университет — имени первого Президента России Б.Н. Ельцина"

Научный руководитель: доктор физико-математических наук,

профессор А.Н. Красовский

Официальные оппоненты: доктор физико-математических наук,

профессор, член-корреспондент РАН H.H. Калиткин

доктор физико-математических наук, А.Л. Агеев

Ведущая организация: Институт вычислительного моделирования

Сибирского отделения РАН, г. Красноярск

Защита состоится « £■£» QSXGLbpia 2010 г. в часов на заседании диссертационного совета1 Д 212.286.10 при ГОУ ВПО "Уральский государственный университет им. A.M. Горького'по адресу: 620000, г. Екатеринбург, пр. Ленина 51, комн. 248.

С диссертацией можно ознакомиться в научной библиотеке ГОУ ВПО "Уральский государственный университет им. A.M. Горького".

Автореферат разослан « » Upobpg

2010 г.

Учёный секретарь диссертационного совета, доктор физико-математических наук, профессор / ' " В.Г. Пименов

Общая характеристика работы

Актуальность

Жёсткие системы обыкновенных дифференциальных уравнений (ОДУ) неизбежно возникают при математическом моделировании задач со многими диссипативными процессами, характерные скорости которых сильно различаются. Примером таких задач являются: моделирование работы двигателей внутреннего сгорания и ракетных двигателей, расчёт широкодиапазопной радиоаппаратуры, моделирование процессов в химических и нейтронных реакторах, а так же многие задачи микроэлектроники.

В практической части данной работы смоделирован процесс образования периодических наноструктур на поверхности оксида алюминия. Подбирать режимы формирования пор экспериментальным путём - длительное и дорогостоящее занятие. Моделирование процесса на компьютере позволяет это сделать быстрее и гораздо дешевле.

Процесс образования нанопор описывается уравнением Курамото - Сива-шинского. Одномерное уравнение Курамото - Сивашинекого в классической монографии Ханера-Ваннера отмечено как один из самых сложных тестов для численных методов. В дайной диссертации рассмотрен ещё более сложный -двумерный случай.

Жёсткие системы ОДУ требуют применения специальных численных методов, обладающих повышенной устойчивостью. Разработка численных методов, сочетающих требования высокой точности аппроксимации, устойчивости и экономичности, которой посвящена данная диссертация, по сей день является актуальной задачей.

Цель работы

Целью данной работы является построение семейства двухстадийных схем Розенброка с комплексными коэффициентами для численного решения жёстких систем ОДУ. Данная задача включает:

1. Разработку алгоритма автоматического построения системы уравнений -условий порядка для двухстадийных схем Розенброка с комплексными коэффициентами;

2. Решение полученной системы уравнений и нахождение коэффициентов схем;

3. Исследование свойств построенных схем.

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

Научная новизна

До последнего времени автоматическое построение условий порядка применялось только для схем типа Рунге-Кутта. В данной работе впервые разработан и применён алгоритм автоматического построения системы уравнений-условий порядка для двухстадийных схем Розенброка.

Получено семейство двухстадийных схем Розенброка, 4 из которых впервые получены в данной работе, реализующих четвёртый порядок точности. Ранее четвёртый порядок точности в схемах Розенброка мог быть получен при числе стадий не менее четырёх*.

Впервые доказана теорема о максимальном порядке точности схем данного класса.

Положения, выносимые на защиту

• Алгоритм автоматического построения условий порядка двухстадийных схем типа Розенброка с комплексными коэффициентами, позволяющий получать условия аппроксимации до любого заданного порядка точности.

• Теорема о сходимости и теорема о максимально возможном порядке точности для двухстадийных методов Розенброка с комплексными коэффициентами.

• Семейство построенных двухстадийных схем Розенброка с комплексными коэффициентами, удовлетворяющих условиям аппроксимации 4-го порядка на произвольных задачах и сочетающих требования устойчивости и повышенной степени аппроксимации на линейных задачах.

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

• Комплекс программ для моделирования образования периодических наноструктур на поверхности оксида алюминия.

'Здесь имеется в виду известный пакет 1Ю54, реализующий чегырёхстадийний метод Розенброка с действительными коэффициентами

Достоверность полученных результатов

Сходимость предложенных методов доказана теоретически. Кроме того, результаты подтверждены тестированием на нескольких задачах различной жёсткости с апостериорной оценкой погрешности.

Результаты математического моделирования процесса образования периодических наноструктур согласуются с результатами физических экспериментов.

Теоретическая и практическая значимость

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

По результатам тестирования, полученные в данной работе схемы обеспечивают точность в 10-1000 раз лучше чем известные схемы входящие в пакет 1Ю54 при одинаковом числе узлов сетки. Более того, чем больше жёсткость задачи, тем больше преимущество построенных схем. Это существенно сокращает трудоёмкость вычислений и позволяет рекомендовать их к практическому применению для жёстких и сверхжёстких задач.

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

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

Публикации

Основные результаты диссертации отражены в работах [1-Ю]. Работы [I-5] опубликованы в журналах, входящих в перечень ВАК. В совместных публикациях [1-4] Е.А. Альшиной и А.Б. Альшину принадлежат постановки задач, общие методики исследований и идеи доказательств основных утверждений, а диссертанту - доказательства основных утверждений, реализация численных методов, разработка программных средств и проведение численных экспериментов.

Апробация работы

Результаты работы были представлены на конференции 1CNAAM 2008, посвященной 75-ти летию Д.Бутчера в г.Кос, Греция [6], на международной конференции NASConi'08 в Ростове-на-Дону [7], дважды на конкурсе молодых ученых в области наноэлектроники в рамках международного форума по нанотехнологиям РОСНАНО '08 [8] и РОСНАНО'09 [9], на конференции Микроэлектроника и наноинженерия-2008 в Московском институте электронной техники [10].

В 2008-2010 годах сделаны доклады на семинаре Института вычислительного моделирования Сибирского отделения РАН, на семинаре члена-корреспондента РАН H.H. Калиткина в Институте математического моделирования РАН, на семинаре кафедры ВМ1 Московского государственного института электронной техники и на семинаре кафедры Вычислительной математики Уральского государственного университета им. A.M. Горького.

Структура диссертации

Диссертация состоит из введения, четырёх глав и одного приложения, посвященных изложению оригинальных результатов автора. Главы пронумерованы арабскими цифрами. Каждая глава разбита на разделы. В номере раздела первая цифра обозначает номер главы, а вторая цифра, отделённая точкой, номер раздела внутри главы. Диссертация состоит из 92-х страниц, содержит в общей сложности 19 таблиц и 22 рисунка. Список цитируемой литературы включает 33 наименования.

Благодарности

Благодарю научного руководителя, д.ф.-м.н. А.Н. Красовского за помощь и поддержку на протяжении долгого времени.

Особенную благодарность выражаю к.ф.-м.н. Е.А. Алыииной и к.ф.-м.н. А.Б. Альшину за ценные советы, без которых данная работы бы не состоялась.

Благодарю сотрудников кафедры материаловедения и физической химии МИЭТ: д.т.н. С. А. Гаврилова и к.т.н. А.Н. Белова за предоставление результатов физических экспериментов.

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

Содержание работы

Во введении обоснована актуальность данной работы. Сформулированы требования, предъявляемые к численным методам решения жёстких ОДУ. Этими требованиями традиционно являются условия Л-устойчивости и L-устойчивости. В данной работе также рассматривается более детальная классификация, предложенная H.H. Калиткиным — так называемая Lp-устойчивость.

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

Первая глава посвящена алгоритму символьных вычислений для построения системы уравнений-условий порядка для двустадийных схем Розенброка с комплексными коэффициентами. Схема Розенброка позволяет построить численное решение задачи Коши:

y'(t)=f(ym (1)

y{to) =Уо-

Здесь у - искомая вектор-функция, a f(y(t)) - заданная вектор-функция той же размерности.

Двухстадийная схема Розенброка имеет вид:

Уп+i = Уп + Re{(nki + b2k2), (2)

где к\ и к2 находятся из соответствующих систем линейных уравнений:

\Е - Taify(yn)\k! = т/(уп), [Е - Ta2fv(yn + тЯе{ик1))]к2 = rf(yn + r/?e(cfca)).

Здесь у(п} - численное решение задачи ( 1 ), т - шаг по времени, Е - единичная матрица, /у - матрица Якоби системы (1), a ai,a2>bi,b2,a и с - комплексные параметры, определяющие свойства схемы.

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

Задача построения системы уравнений-условий порядка является достаточно трудоемкой, и многие исследователи прибегали к некоторым упрощающим соображениям для построения условий порядка. В данной работе все выкладки выполняются с помощью специально разработанной программы, позволяющей производить символьные вычисления. Работа программы основана на представлении элементарных дифференциалов в разложении точного и численного решения в виде деревьев.

В разделе 1.1 определён метод преобразования символьной записи элементарных дифференциалов задачи (1) в запись в виде графов (деревьев). Для автоматического определения коммутирующих тензоров, введено понятие топологической эквивалентности деревьев, соответствующих коммутируемым тензорам в разложении точного и численного решения. Построен алгоритм нахождения топологически эквивалентных деревьев с помощью вычисления функции стоимости вершин.

Раздел 1.2 описывает процесс построения разложения точного решения с использованием представления элементарных дифференциалов в виде деревьев. Выведено правило дифференцирования в терминах деревьев и представлено разложение точного решения.

Разложение точного решения не зависит от выбора численного метода и может быть получено несколькими разными способами, включая методы комбинаторики. В данной работе используемый метод дифференцирования деревьев, приведённый в раздел 1.2 обобщён для построения разложения численного решения.

Раздел 1.3 посвящен разложению численного решения. Численное решение, задаваемое схемой Розенброка (2), на новом временном слое представляется через приращения первой и второй стадии. В подразделах 1.3.1 и 1.3.2 рассматриваются разложения первой и второй стадии отдельно.

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

В части 1.4 объединяются ранее полученные разложения точного и численного решения и строятся условия порядка для двухстадийного метода Розенброка (2). Полученные условия порядка приведены в таблице 1.

Порядок Деревья Условия порядка

1 _ 1 = Яе (¡>1 + Ьг)

^ = Яе (г>1«1 + Ь2и-2) + Щс)йе (62)

\ =Яе(т,)Яр(Ы + Яе(с№ (Ь2п2) +

и

+ Яе (Ь2а^ + «561)

V <1= + Яе ь

^ =Ле(ш?)Яс(г>2) + Яс(са,)Ле (Ь2«2) + + Не {Ь2аг2) Яе(с) + Яс (а^Ь, + г>2о2)

4 = Яе (Ь,о2) ^Яс(е)2 + Яе (Ь2г^) Яе(а) ¿4 2

- =Яе (<;«!) Яе(е)Яе (Ь2) + Яе (Ь2<»2) Я« («<»,) +

О

+ Яе (Ь2а2) Яе (а) Яе (с) + Яе (^а2) Яе (а) \|/ ~ = ^ (с)3 Яе (62) + Яе (62о2) Яс (а)2

\ 40

1

30

\ Ш =Ле + Яе ^ Яс ^ Нс +

+ Яе (Ь2а2) Яе (си?) + Яе (Ь2о2) Яе (с) + + Яе(М2) Яе^пО

\ ¿о = \Пе ^ Яе Сс)2 + Яс (а)

^ =Яе (Ь2«2) Яе (паО + Яе (Ь2о*) Яе (а) + + Яе (Ь2аг) Яе (с) Яе (сгц) + + Яе (62а2) Яе (а) Яе (с) =Яе (Ьгаг) Яе (аа?) + Яе (Ь2»|) Яе (а) +

+ Яе (с) Яе (Ь2) Яе («»?) + Яе (с) Яе (Ь2а^) Яе (а) + + Яе (Ь2о2) Яе (а) Яе (етц)

X Й = 1Яе (г,2а2'Яе +Не Не ^

\/ 4 = (б2°2) Яе (а) Кс (с)2 + Лс Де (а)2

V ~ (Ы Яе («л)2 + Яе (62а2) Яе (аа,) + + Яе (Ь2а2) Яе (с) Яе (ааг])

\/ 5) (&2) ЙС (ССП) Ле (")2 + (''2"2) Лб (С) ЙК (а)2 +

+ Яе (г»2«2) Яе (а) Яе (001) + ^Яе (Ь2с$ Яе (а)2

Ч/^' й = \Ле (Ьг) Ле (с)4+Ле (Ь2<1а) йе (а)3

Таблица 1: Условия 5-го порядка для двухстадийной схемы Розенброка

Ранее условия 4-го порядка для схемы (2) были получены П.Д.Ширковым в 1992 году.

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

В разделе 1.5 доказана теорема о сходимости численного метода (2) при

условии выполнения уравнения-условий порядка.

Теорема 1. Пусть правая часть уравнения (1) четырежды непрерывно дифференцируема, а комплексные числа а^аг,^!, Ь2,аис таковы, что выполнены первые 8 уравнений из таблицы I, тогда двухстадийный метод Розенброка (2) сходится с четвёртым порядком, то есть ||уп — у (£п)||с = О (г4) для £„ — пт < Т.

Условия 5-го порядка содержат 17 уравнений, а схема (2) имеет всего б комплексных (12 действительных) параметров. Таким образом, без введения дополнительных параметров в схему (2), эту систему решить нельзя. В разделе 1.6 рассматривается естественное обобщение двухстадийного метода Розенброка с введением дополнительных параметров, которое позволит в будущем построить схемы более высокого порядка всего на двух стадиях:

у =у + Яе Ьхкг + Ке к2, [Е - та^у] кг =т/,

7 (3)

\Е - та2/у (у + т Ие (а**))] к2 =т ]Г Ь28/(у + т Яе (с^)).

4 = 1

Схема (2) является частным случаем схемы (3) при J — 1. Условия порядка для схемы (3) получаются из условий порядка для схемы (2) путем замены

./ .7

&2 и с на соответствующие суммы Ь^ и В разделе 1.6 основного

5=1 5=1

текста работы приводятся уравнения-условия порядка для обобщённого метода Розенброка.

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

Теорема 2. Невозможно построить двухстадийную схему Розенброка (3), имеющую точность выше, чем 0(т6).

Вторая глава посвящена решению системы уравнений-условий порядка до 4-го порядка включительно, полученной в Главе 1, и нахождению коэффициентов схемы (2). Эти уравнения связывают 6 комплексных (или 12 действительных) параметров. Таким образом, получившаяся система 8-ми уравнений недоопределена. Систему можно дополнить уравнениями обеспечивающими затухание функции устойчивости порядка 0 < р < 4 и обеспечивающими повышенный порядок точности 4 < тп < 8 на линейных задачах. При этом

порядок р затухания функции устойчивости на бесконечности и максимальный порядок точности т на линейных задачах связаны соотношением (р + т) = 8.

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

Обозначение р тп

СЦ052_4 4 4

саоэ2 3 3 5

СЯОБ2 2 2 6

1 7

С1ЮБ2_0 0 8

Таблица 2: Обозначения построенных схем и их свойства: порядок затухания функции устойчивости р и порядок точности на линейных задачах га.

Для двухстадийной схемы (2) функция устойчивости имеет вид:

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

где А, В, С и И так же выражаются через коэффициенты «1 и »2-Для получение нужного порядка затухания функции устойчивости, коэффициенты при степенях £ в числителе функции устойчивости(4) необходимо приравнять к нулю. Тогда, при А — В — С~ И — 0 получим схему с порядком затухания р = 4, при В ~ С = И — 0 — с порядком р = 3, при С = .О = 0 — с порядком р = 2, и при И — 0 — с порядком р = 1.

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

Схема СЯОБ2_4

«1 =0.4573733434972976 + г • 0.2351004879985425, «2 =0.04262665650270241 + i • 0.3946329531721134, =0.7893434641361923 + г ■ 0.9821367946107931, Ь2 =0.2106565358638077 - г • 0.5705215732509971, с =0.6444138212147357 - г • 1.143956305335963, а =0.5250462591428808 + i ■ 1.453646467184172.

Схема СЯОЭ2 3

«1 =0.3074021043872249 + г • 0.1292532396046484, о>2 =0.09259789561277514 + г • 0.2576121583025594, =0.8644582665498726 + г • 0.9366952975243449, Ь2 =0.1355417334501275 - г • 1.154171181438793, с =0.3353594637740966 - г • 0.4983420242149068, а =0.5132472378039463 + г • 0.2267734198731172.

Схема СНОв2_2

а! =0.2334763488700170 + г • 0.08527040833242157, а2 =0.09985698446331641 + г • 0.1870544254177949, Ьг =0.9248875101862942 + г ■ 0.7077449395923038, Ь2 =0.07511248981370576 - { ■ 1.698741848884691, с =0.2549862725007512 - г • 0.3381738431416763, а =0.5049068514817424 - г • 0.4325579331793709.

Схема СИ082_1

оц =0.09705048233513194 + г • 0.1441824711215367, а2 =0.1886638033791538 + г • 0.06177441689689114, Ъх =0.04833419895509594 - г • 0.3205959705202483, 62 =0.9516658010449041 - г • 1.696774337833587, с =0.1730887968652113 --г • 0.1694095699539014, а =0.5359744564304916 - г • 0.9665922748484184.

Схема СЦ082_0

щ =0.09156624026571748 + г ■ 0.1156626130131271 а2 =0.1584337597342825 + г ■ 0.04744101257110980 Ъх =0.2803648780046792 - I ■ 0.1985114563571360 Ь2 =0.7196351219953208 - г ■ 2.479090148032283 ( '

с =0.3053528612690534 - г • 0.2319031769536116 а =0.5747096314647993 - г • 1.092069699677101

График сравнения степени затухания функций устойчивости для схем СГ{082_4, С 1^.0Б2_3, СЯОБ2_2, СЯОЭ2_1 и СЯОБ2_0 приведён на

рисунке 1.

Рис. 1: Сравнение затухания функций устойчивости с различной степенью р

Из графика, приведённого на рисунке 1, видно, что модуль функции устойчивости схем СЯОЭ2_4 и С7? ОБ2_3 превышает 1. Таким образом, для этих схем в небольшой области нарушается условие А-устойчивости. И наоборот, для схем СЯОБ2_2, С1{032__1 и СКОБ2_0 модуль функции устойчивости на мнимой оси (а в силу принципа максимума и на всей левой полуплоскости) не превышает 1. Следовательно, схемы С1^0Б2_2, С11,082_1 и СЯОБ2_0 являются Ь2,Ы и Л-устойчивыми схемами соответственно.

В третьей главе производится оценка погрешности для полученных ранее схем и выполняются тестовые расчёты трёх задачах с апостериорной оценкой погрешности.

В разделе 3.1 сравниваются невязки в системе уравнений-условий 5-го порядка. В условиях 5-го порядка наибольший модуль невязок имеет схема

С1Ю82_4, а наименьший Тем самым, в общем случае на нежёст-

ком участке точность схемы СИ.082_0 будет выше всех остальных схем.

В разделе 3.2 рассматривается метод апостериорной оценки погрешности Эйткена. Для такой оценки погрешности выполняется серия расчётов на сгущающихся равномерных сетках с шагом 2т и 4т. В совпадающих узлах t локальная погрешность по методу Эйткена выражается следующей формулой:

Л2 = «<*>(*)-«<*>(*); (10)

Здесь и - решение на текущем слое с шагом т, 2т

и 4т соответственно, а е - локальная погрешность. Оценка погрешности £ асимптотически точна при т —> 0.

Для проверки эффективного порядка точности полученных схем в подразделе 3.2.1 проводится расчёт следующей задачи:

у[ = - 1лвт(у2), у\ (0) = 10;

у'г =зт(У1), у2( 0)=0; (11)

ц =20, 0 < « < 3.

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

Для задачи (11) был выполнен расчёт на десяти сгущающихся сетках с увеличением шага в два раза. На графике 2 в двойном логарифмическом масштабе изображена зависимость максимума модуля локальной погрешности е, вычисленной по формуле Эйткена (10), от шага сетки т для схем С1Ю82_4 -

СЯОв2_0.

Начиная с достаточно большого шага по времени т < 0.5 все точки графиков хорошо ложатся на прямые с углом наклона ф, где ¿Д (</->) = 4. Это является экспериментальным подтверждением того, что все полученные схемы реализуют теоретический четвёртый порядок точности.

В подразделе 3.2.2 рассматривается стандартная задача Протеро-Робин-сона, которая в монографии Хайер-Ваннера включена в список двенадцати жёстких задач-тестов для разностных схем:

У\ = -РЫ ~

^2 = 1» (12) 2/1 (0) = 0, у2 (0) = 0.

Здесь /х - коэффициент жесткости.

Рис. 2: Погрешность численного решения задачи (11)

На рисунке 3 приведены графики зависимости максимально го модуля локальной погрешности от размера шага т. Погрешность оценена по методу Эйткена (10).

а) б)

Рис. 3: Погрешность численного решения задачи Протеро-Робинсона: а) ¿t = 20, б) ß = 200.

Для сравнения на рисунке 3 приведены результаты чстырёхстадийной схемы Розенброка с действительными коэффициентами, входящей в пакет ROS4. Среди всех схем этого пакеты была выбрана схема с оптимальными параметрами, предложенная Д,С. Гужевым и H.H. Калиткиным. Только схема CROS2_4

уступает по точности лучшей схеме пакета Я,084. Остальные схемы, построенные в данной работе, на данной тестовой задаче обеспечивают в 10-1000 раз лучшую точность.

Для демонстрации преимущества Ьр-устойчивых методов с большими р была проведена серия расчётов задачи (12) с постоянным шагом т — 0.075, но с увеличением жесткости ц. Зависимость максимального модуля локальной погрешности от увеличения жёсткости ¡л приведена на рисунке 4.

100-. 100.1, О.01. 1Е-3.

1Е-6 -1Е-7,

1Е-8 ] ^

1Е-Ю4-100

Рис. 4: Погрешность численного решения задачи Протеро-Робинсона при постоянном шаге и увеличении жёсткости ц.

При маленьких значениях (I точность схемы ()Б2__0 выше, поскольк>' эта схема имеет меньшие невязки в условиях аппроксимации 5-го порядка и, следовательно, меньшую погрешность. Но с ростом жёсткости [I погрешность схемы СЯ082_0 быстро растёт. Для 1Л-устойчивой схемы С/?052_/ погрешность нарастает гораздо медленнее, а для схем с р > 2 погрешность вообще практически не увеличивается. При ¡1 > 3000 видно, что, например, Ь2~ устойчивая схема дает погрешность на порядок лучше 2/1-устойчивой схемы, а решению Л-устойчивой при таком же значении ¡А. и шаге т вообще не стоит доверять.

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

В подразделе 3.2.3 рассматривается осциллятор Ван-дер-Поля:

У\ =3/2, У1(0) = 0;

у'2 (1 - У\) У2 ~ УЪ Ш (0) = 2. К }

Фазовая траектория уравнения (13) при /х = —100 изображена на рисунке 5. Участки с плавным изменением фазовых координат сменяются резкими поворотами, что является проявлением жёсткости.

Рис. 5: Фазовая траектория задачи Ван-дер-Поля при ц — —100

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

В качестве механизма для автоматического выбора шага был использован следующий алгоритм. На каждом шаге производится оценка локальной погрешности по формуле Эйткена (10). Для этого выполнялся один шаг с заданным г, 2 шага с шагом т/2 и 4 шага с шагом т/4.

Задав желаемую точность расчёта toi, на каждом шаге будем сравнивать априорно оценённую глобальную погрешность со значением toi. Если оценка глобальной погрешности сильно отличается от желаемой точности toi будем проводить корректировку шага т.

Приведённый выше метод автоматического выбора шага позволил существенно сократить трудоёмкость численного решения. Например точность 10~4 в расчёте с автоматикой выбора шага требует в 100 раз меньшего числа узлов сетки. Заданное значение toi использовалось для априорной оценки глобальной погрешности. Глобальную погрешность можно оценить, проведя серию расчетов на сгущающихся сетках по методу Эйткена, и такая оценка, как известно, является асимптотически точной. Сравним априорную й апостериорную оценки глобальной погрешности.

На рисунке 6а по оси абсцисс отложено заданное значение toi, а по оси ординат — апостериорная оценка глобальной погрешности. Хорошо видно, что эти величины отличаются на 2-3 порядка. Априорная оценка может и должна быть использована в стратегии выбора шага, хотя она и не является надёжной.

Также была оценена погрешность от числа шагов метода выполненных при заданном toi. Эта зависимость показана на рисунке 66. Угол наклона графика

■1S0'

•2

О

tz±)

Рис, 6: Апостериорная оценка глобальной погрешности численного решения задачи Ван-Дер-Поля при /л = 100 и времени расчёта Т = 200. а) Зависимость погрешности от заданного значения toi. б) Зависимость погрешности от числа шагов N.

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

В заключении главы 3 делается вывод о том, что на тестовых задачах лучшие результаты показали схемы CROS2__l и CROS2_2. И именно они рекомендуются к практическому применению.

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

В разделе 4.1 описана математическая модель процесса. Она построена по аналогии с работами С. Sample и A.A. Golovin, однако содержит некоторые уточнения, позволяющие получить количественное совпадение результата численного моделирования с физическим экспериментом.

Схема процесса анодного окисления алюминия изображена на рисунке 7.

В приближении слабых нелинейных возмущений эволюция толщины поверхности и(Х, Г) — — описывается начально-краевой задачей для уравнения Курамото-Сивашинского:

8\U + öi V2u + а^Ч^и — ~(Vu)2 — 0. (14)

¿i

u(x + dx,y,t)=u(x,y,t),

u(x,y + dy, t) =u (x, y, t),

u{x,y,t) =fo{x, y).

оксид

;электролит

СДХЛ) Цу.1)

Рис. 7: Схема процесса анодного окисления.

Здесь fo(x, у) - начальный профиль поверхности. Также выбраны периодические граничные условия по пространству.

Коэффициенты уравнения öi, а2 и vs зависят от параметров среды, напряжения на аноде V, кислотности pH и температуры Т. Точный вид зависимости громоздкий и здесь не приводится.

В разделе 4.2 описан численный метод решения задачи (14). Данная задача для уравнения в частных производных решалась методом прямых.

Как известно из эксперимента, поры имеют гексагональную структуру. Гексагональная симметрия ячеек возможна, если стороны прямоугольной области имеют размеры dx = d, dy — dy/3.

Была введена пространственная сетка:

х„ =hxn, hx = (1 < п < N),

dy/3

ут =hym, hy = , (1 <т< М).

Число узлов сетки по а; и по у выбрано таким образом, чтобы hx ~ hy ~ h. В практических расчётах были использованы сетки размеров N = 10 • 2', М = 17 • 2г, где I - целое число.

Пространственные производные в (14) были заменены разделёнными разностями, используя специально построенную для данной работы аппроксимацию точности 0(h4). Таким образом, уравнение в частных производных было сведено к системе ОДУ для неизвестной толщины поверхности в узлах сетки

^ (-£«) Ут) — Мпт-

Эта система ОДУ решалась с помощью схемы CROS2_2, построенной в данной работе и показавшей лучшие результаты при тестировании.

Выражения для коэффициентов а\ и vs были ранее получены в работах С. Sample и A.A. Golovin. Однако приведённая в этих работах формула для

коэффициента a<i даже по размерности не совпадает с '-^j-, как должно быть для соразмерности величин в уравнении (14).

В данной работе была проведена калибровка модели по результатам вычислительного эксперимента и установлена зависимость коэффициента а2 от напряжения на аноде V.

В статьях О. Jessensky, F. Muller и U. Gosele, а так же F.Y.Li, L.Zhang и R.M.Metzger описаны несколько данных физических экспериментов, где для заданных напряжения V, кислотности pH и температуры Т приведены снимки наноструктур на поверхности оксида алюминия, сделанные при помощи электронного микроскопа.

При тех же значениях напряжения V, кислотности pH и температуры Т были проведены серии расчётов с разными значениями параметра а2. Целью данного вычислительного эксперимента являлось достижение максимального количественного совпадения размера и формы периодических наноструктур, полученных математическим моделированием, с данными реальных экспериментов.

В результате этих экспериментов была получена следующая зависимость параметра 02 от напряжения V:

а-2 = 1,28395 + 0,87346 • V + 0,07272 • V2.

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

а) б)

Рис. 8: Наноструктуры, образовавшиеся при V = 25У,рН — 0.52, Т = 10°С. а) Микроснимок поверхности оксида алюминия, б) Результат численного моделирования.

(шт.;* ШФФФФФФФФФФФФФФФФФЛ1

(••••••••О*•*•«•«•«

• ••••••••О**»*»«**«

»»»»••«•«•«»»»»«а*** --шффффффффффффффффф

1ВС0 (пт)

б)

Рис. 9: Наноструктуры, образовавшиеся при и = 40V,рН = 0.52, Г = ГС. а) Микроснимок поверхности оксида алюминия, б) Результат численного моделирования.

б)

Рис. 10: Наноструктуры, образовавшиеся при V — 160У.рН = 0.52, Т = 3°. а) Микроснимок поверхности оксида алюминия, б) Результат численного моделирования.

Благодаря проведённой калибровке модели удалось добиться хорошего количественного совпадения результатов моделирования с результатами эксперимента. Так на рисунках 8, 9 и 10 период получившихся наноструктур равен б Опт, 95 пт и 420пт соответственно.

Основные результаты работы

1. Разработан алгоритм и написана программа, позволяющая получать систему уравнений-условий порядка для двухстадийных схем Розенброка с комплексными коэффициентами, обеспечивающих аппроксимацию до любого заданного порядка точности. Доказана теорема о сходимости.

2. Исследован максимально возможный порядок точности схем данного типа. Доказана теорема о барьере точности, аналогичная барьеру Бутчера для схем Рунгс-Кутта.

3. Построено семейство из 5-ти схем (4 из них новые) для решения жёстких систем ОДУ. Исследована их устойчивость и проведены серии сравнительных тестов, которые позволяют рекомендовать две из построенных схем к практическому применению.

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

Публикации по теме диссертации

Статьи в ведущих рецензируемых научных журналах, определённых ВАК:

[1] Л.Г. Лимонов Л.Б. Ллыиии, К-А. Альшипа. Автоматизированное символьное построение условий порядка для двухстадпйных комплексных схем типа Розенброка // Матем. моделирование, 2009. том 21, Л» 12 е. 70-88.

[2) Limonov A.G., Al'shina Е.А., Al'shin, А.В. Two-stage complex Rosenbrock schemes for stiff systems // Computational Mathematics and Mathematical Physics, MAIK Nauka/lnterperiodica, 2009, vol. 49, no. 2, pp. 261-278.

[3j Limonov A.G., Alshin А.В., Aishina F.A. Symbolic Derivation of Order Conditions for Two-Stage Complex Rosenbrock Scheme. //American Institute of Phisics, 2008, vol. 1018, pp. '17-51.

[■1] Л.Г. Лимонов A. 13. Альшин, E.A. Альшипа. Двухетаднйные комплексные схемы Розснброка для жестких систем. //Ж- вычисл. матем. и матем. физ., 2009, том 49, № 2. с. 270-287.

[5] А.Г. Лимонов. Моделирование образования гексагональных периодических наноструктур на поверхности оксида алюминия. //Матем. моделирование, 2010, том 22, № 8 с. 97-10§.

Другие публикации:

[6] A. G. Limonov, Е. A. Aishina, А.В. Alshin. Symbolic Derivation о! Order Conditions for Two-Stage Complex Rosenbrock Scheme. // International Conference on Numerical Analysis and Applied Mathematics, 2008.

[7] A.G. Limonov, E.A. Aishina, A.B. Alshin. Constructing Two-Stage Complex Rosenbrock Schemes By Using Automatic Symbolic Derivation Of Order Conditions.//international Workshop Numerical Analysis and Scientific Computing(NASCom'08), 2008.

[8] А.Г. Лимонов, E.A. Альшина, А.Б. Альшин, C.A. Гаврилов. Математическое моделирование процесса образования нано-пор на поверхности оксида алюминия. // Международный Конкурс Научных Работ Молодых Ученых в Области Нанотехнологий, Москва, Экспоцентр, 2008."

[9| А.Г. Лимонов, Е.А. Альшина, А.Б. Альшин, А.Н. Белов. Математическое моделирование процесса образования нано-пор па поверхности оксида алюминия. // Международный Конкурс Научных Работ Молодых Ученых в Области Нанотехнологий, Москва, Экспоцентр, 2009.

[10] А.Г. Лимонов, Е.А. Альшина, А.Б. Альшин, А.Н. Белов, Численное моделирование процесса образования нано-пор на поверхности оксида алюминия. // Микроэлектроника и наноинжцнерия-2008, Московский институт электронной техники.

Подписано в печать 29.10.2010. Формат 60x84 1/16 Гарнитура «Times». Усл. печ. л. 1,4. Тираж 100 экз. Заказ № 0 (Ч

Отпечатано в типографии ИПЦ «Издательство УрГ'У» 620000, Екатеринбург, ул. Тургенева, 4

Оглавление

Введение

Цель работы.

Жёсткие задачи.

Жёсткая устойчивость.

Численные методы для жёстких систем

Метод Розенброка.

Интерполяционные свойства схем.

Оценка погрешности

Символьные вычисления.

Пористый анодный оксид алюминия.

Уравнение Курамото-Сивашинского

Структура работы.

Публикации

Апробация работы.

1 Построение условий порядка

1.1 Запись в виде деревьев.

1.2 Разложение точного решения

1.3 Разложение численного решения.

1.3.1 Разложение первой стадии.

1.3.2 Разложение второй стадии.

1.4 Условия порядка.

1.5 Сходимость.

1.6 Обобщение двухстадийных схем Розенброка.

1.7 Барьер точности для двухстадийных схем Розенброка.

2 Вычисление коэффициентов схем и анализ устойчивости

2.1 Функция устойчивости.

2.2 Решение системы уравнений - условий порядка.

2.3 Коэффициенты схем и анализ функции устойчивости

3 Тестирование полученных схем и оценка погрешности

3.1 Априорная оценка погрешности.

3.2 Апостериорная оценка погрешности.

3.2.1 Пример 1. Проверка эффективного порядка точности.

3.2.2 Пример 2. Задача Протеро-Робинсона.

3.2.3 Пример 3. Задача Ван-Дер-Поля

3.3 Выводы.

4 Моделирование образования нанопор на поверхности оксида алюминия.

4.1 Математическая модель.

4.2 Численное решение.

4.2.1 Пространственная модель

4.2.2 Экспериментальное уточнение параметра аг.

4.3 Сравнение результатов моделирования с физическими экспериментами

4.3.1 Описание физического эксперимента.

4.3.2 Эксперимент 1.

4.3.3 Эксперимент 2.

4.3.4 Эксперимент 3.

Цель работы

Цель данной работы - построение семейства двухстадийных схем Розенброка с комплексными коэффициентами для численного решения жёстких систем обыкновенных дифференциальных уравнений (ОДУ). Данная задача включает: разработку алгоритма символьных вычислений для получение системы уравнений - условий порядка, нахождение решений этой системы, вычисление коэффициентов схем и исследование свойств построенных схем.

Практической целью работы является моделирование процесса образования нанопор на поверхности оксида алюминия - одной из актуальных задач современной микроэлектроники.

Жёсткие задачи

Жёсткие системы ОДУ неизбежно возникают при математическом моделировании задач со многими диссипативными процессами, характерные скорости которых сильно различаются.

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

Высокая точность при расчёте мощности ракетного двигателя является вполне обоснованным экономическим требованием. Так погрешность всего в 1% при расчёте тяги приводит к погрешности в 50% при вычислении массы полезного груза, который может данный ракетный двигатель вынести на орбиту.

Расчёт широкодиапазонной радиоаппаратуры приводит к жёсткой системе ОДУ. Любые нестационарные процессы в разветвлённых электрических цепях описываются дифференциально-алгебраической* системой законов Ома и Кирхгофа.

Моделирование процессов в химических и нейтронных реакторах приводит к задаче, где несколько реакций протекают одновременно с разной скоростью. Различие скорости протекания реакций хможет составлять 105 — Ю10 .

Дифференциально-алгебраическую систему можно трактовать, как задачу с бесконечной жёсткостью.

При решении начально-краевых задач для уравнений в частных производных часто применяют, так называемый, метод прямых. При этом все пространственные производные заменяют разностными аналогами. В результате каждое уравнение в частных производных заменятся на систему ОДУ для значений решения в узлах пространственной сетки. Получающиеся при такой замене система ОДУ большой размерности в большинстве случаев является очень жёсткой и требует применения специальных численных методов, пригодных для жёстких систем.

Исторически интерес к жёстким системам возник в середине XX в. при изучении уравнений химической кинетики с одновременным присутствием очень медленно и очень быстро протекающих химических реакций. Тогда неожиданно оказалось, что считавшиеся исключительно надёжными методы Рунге—Кутта [1] стали давать сбой при расчёте этих задач.

Строгого общепринятого математического определения жёстких ОДУ не существует. В 50-х годах Кёртиссом и Хиршфельдером[2] было предложено следующее определение жёстких систем: "Жёсткие уравнения - это уравнения, для которых определённые неявные методы, в частности ФДН, дают лучший результат, обычно несравненно более хороший, чем явные методы"

Жёсткость системы ОДУ принято характеризовать в терминах собственных значений матрицы Якоби её правой части. Традиционно при численном решении большую сложность представляют задачи, в которых спектр матрицы Якоби правой части разделяется на две части, где жёсткая часть спектра имеет хотя бы одно собственное значение с большой по модулю и отрицательной реальной частью, а мягкая часть спектра имеет собственные значения, существенно меньшие по модулю. Рассмотрим систему ОДУ: t = fi(t,yuV2>-,yn), dy 2 dt

1) г = fn{t,yuy2, •••,Уп),

Без ограничения общности можно считать задачу (1) автономной (правая часть не зависит от времени £ явно), так как любая неавтономная задача может быть сведена к автономной задаче путём добавления в систему одного уравнения и увеличения размерности вектора у на 1. Перепишем систему (1) в автономной форме: f = 1\{УиУ2,---,Уп,Уп+1), f = Vli V2t •■•> Uní Уп+l))

Г = ín{t, Уи г/2, Уп, Уп+l),

У'Ч i i dt ~ -*->

2)

Далее в тексте везде будет применяться более компактная запись автономной задачи Коши (2) в векторной форме: y'(t) =/(Ут

V(to) =Уо

Здесь у - искомая вектор функция, a f(y(t)) - заданная вектор функция той же раз мерности.

3)

Жёсткая устойчивость

Традиционно принято исследовать поведение разностных схем для жёстких задач на примере задачи Далквиста.

V'(t) = Ay(i) m t> 0,7/(0) = 1.

Она с одной стороны достаточно проста, чтобы выписать точное решение, а с другой хорошо иллюстрирует типичные трудности, возникающие при численном решении жёстких систем. Решением этой задачи является y{t) = eAt, которое стремится к 0 при t —» оо, когда Re(А) < 0.

Для любой линейной схемы переход на следующий временной слой при решении задачи (4) имеет вид у = R(£)y, где /?,(<!;) называется функцией роста или функцией устойчивости зависящей от £ = А т.

Определение 1. Схема называется А-устойчивой, если |Д(0| < 1 при Re{0 < 0.

Это минимальное требование, предъявляемое к численным методам для решения жёстких задач. Если схема не обладает хотя бы А-устойчивостью, то она вообще не пригодна для жёстких задач. Но для большой жёсткости даже А-устойчивости оказывается недостаточно для надёжной работы численного метода. Желательно, чтобы при Яе(£) —> оо функция устойчивости также сильно затухала, иначе жёсткие компоненты хотя и будут демпфироваться, но достаточно медленно, что приведёт к длительному процессу стабилизации.

Определение 2. Метод называется L-устойчивым, если он является А-устойчивым и lim = 0. оо

Для оценки степени этого затухания вводят понятие Lp-устойчивости [3]. В этом случае функция устойчивости затухает пропорционально

Определение 3. Схема называется Lp-устпойчивой, если она является А-устойчивой и = при оо.

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

Определение А-устойчивости, с одной стороны, слишком слабое, но, с другой стороны, оно слишком сильное, что многие не так уж плохие методы не являются А-устойчивыми [2]. Следующее определение сужает область устойчивости метода до сектора Sa = {£; |arc/(-OI < а, ^ф 0}.

Определение 4. Схема называется Аа-устойчивой, если < 1 при

Ksr(-OI < а.

Аналогично определению Аа-устойчивости, вводится определение Lpa-устойчивости:

Определение 5. Схема называется Lpa-устойчивой, если |i?(£)| < 1 для < оси R{0 = 0(ГР) при оо.

Численные методы для жёстких систем

Начиная с 50-х годов прошлого века, для жёстких задач стали создавать специальные неявные методы. Наиболее полный обзор методов решения жёстких систем содержится в монографии Э.Хайера и Г. Ваннера [2].

Итерационные методы

Среди неявных схем широко распространены схемы Рунге-Кутта (IRK - implicit Runge-Kutta methods). He все схемы IRK подходят для решения жёстких задач. Огромное число работ посвящено схемам IRK, в [2] дан всеобъемлющий обзор IRK методов, пригодных для жестких задач, выделены методы, являющиеся жёстко точными (тем самым пригодными для ДАУ).

Любой метод IRK для перехода на новый временной слой требует неоднократного решения системы, вообще говоря, нелинейных уравнений при помощи итераций ньютоновского типа. Для s-стадийного IRK минимальное число возникающих нелинейных систем s соответствует диагонально неявным методам (DIRK - diagonal implicit Runge-Kutta methods). Именно они чаще всего и используются на практике.

В [2] исследованы диапазоны параметров методов DIRK, при которых схемы являются А и L-устойчивыми.

Безитерационные методы

Наличие итераций сильно усложняет использование схем IRK, так как к проблемам устойчивости добавляется проблема сходимости итерационного процесса. Альтернатива, которая обходит эту трудность - это методы типа Розенброка ROS и Розенброка-Ваннера ROW [2]. В том числе развитию таких методов посвящены работы Е.А.Новикова [4], Н.Н.Калиткина [5], [6] и других авторов.

Формально эти схемы неявные, но итераций в них не возникает и число арифметических действий для перехода на новый временной слой фиксировано и заранее известно (как в явных схемах). За это безусловное преимущество эти схемы получили название явно-неявных или полуявных.

Среди одностадийных схем Розенброка выделяется схема с комплексным коэффициентом (CROS). К сожалению, эта схема до сих пор малоизвестна. Например, в классической монографии [2] она даже не упомянута. Применение схем типа Розенброка для интегрирования ДАУ рассмотрено в [7], но там также исследованы только схемы с действительными коэффициентами.

Многошаговые методы

Многошаговые методы положены в основу популярных программ Гира [8]. Коэффициенты многошаговых методов подбираются так, что бы (/-шаговый метод имел точность 0{rq). Переход на новый временной слой для этих методов не требует итераций.

Но можно показать, что многошаговые методы в лучшем случае лишь Lpa-устойчивы с малым значением р = т.е. по надёжности они сильно уступают ROS, ROW и IRK методам.

Метод Розенброка

Среди схем для жёстких систем дифференциальных уравнений стоит выделить схемы Розенброка. Схемы эти по сути неявные, но для перехода на новый временной слой нужно решать линейную систему уравнений с хорошо обусловленной матрицей, что позволяет избежать итераций.

Рассмотрим задачу Коши (3) для системы дифференциальных уравнений . Будем полагать, что функция f(t,y(t)) не менее четырёх* раз непрерывно дифференцируема в окрестности решения задачи Коши .

Схемы Розенброка одношаговые - для вычисления значения численного решения у на новом временном слое £ + г используется лишь значение у с текущего слоя ¿. Для 5-стадийной схемы переход на новый временной слой происходит по формулам: у = у + Re ( ¿ bmk,

771 — 1

Е - тат/у ( у + rRe ( £ amih m=l т-1 кт = т/ [у + тЛе £ Cmik 1 1 (5)

1 < т < 6-, где Е - единичная матрица, /у - матрица Якоби системы (3), а ат,Ьт,ат1 и ст/ - комплексные параметры, определяющие свойства схемы.

Использование комплексной арифметики требует большего числа арифметических операций для решения линейных систем в (5) по сравнению с действительной схемой. Но при выборе коэффициентов комплексной схемы имеется в два раза большее число степеней свободны, что позволяет построить схему более высокого порядка точности и большей надёжности, по сравнению с аналогичной действительной схемой.

Одностадийная комплексная схема (5) была предложена Розенброком в 1963 году [9] Уп+1 =Уп + Яе{к), ^

Е - та/у(уп)]к = т/(уп).

Больший но сравнению с действительными схемами объем вычислений стал в те годы непреодолимым препятствием для ее широкого распространения. Для современного уровня вычислительной техники на первый план уже давно выходит не объем вычислений, а надёжность и точность схемы. Одностадийная схема (6) при а = ^ обладает уникальным сочетанием свойств: она имеет максимальный для одностадийной схемы -второй порядок точности и обладает Ьр-устойчивостыо с максимальным для одностадийных схем р = 2.

Целью данной работы является построение комплексной схемы (5) с я = 2, обладающей максимально возможной точностью и устойчивостью. Формула перехода на новый

Для построения схемы порядка р функция должна быть непрерывно дифференцируема р раз. В данной работе будет построена схема 4-го порядка точности. временной слой для 2-х стадийной схемы Розенброка имеет вид:

Уп+1 = Уп + Яе{Ь\к1 + Ь2к2), где кг и к2 находятся из соответствующих систем линейных уравнений:

Е - та1/у(уп)]к1 = т/{уп), [Е - та2/у(уп + тЯе{а2\к\))]к2 = т/(уп + тЯе^к^).

8)

9)

Параметры ац, «2! Ь2, а2х и с2\ будут находиться из условий порядка и дополнительных условий устойчивости схемы. В дальнейшем сокращённо будем записывать параметры а и с без индексов.

Интерполяционные свойства схем

В работе [10] приведены примеры, когда небольшие нарушения интерполяционно-сти приводят к неустойчивой работе численных методов, поскольку при экстраполяции погрешность нарастает быстрее.

Для схем типа Розенброка (5) с действительным коэффициентами условия интерполяционное™ формулируются достаточно просто:

При вычислении с использованием схемы с комплексными коэффициентами, коэффициенты определённым образом перемножаются и от произведения берётся реальная часть. Все реальные части полученных произведений должны лежать в пределах от 0 до 1. В этом случае схема будет интерполяционной.

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

Оценка погрешности

Льюис Фрай Ричардсон показал, как по двум расчетам с шагами /г/2 и к по схеме р-ого порядка точности можно получить хорошую оценку погрешности, а также повысить точность результата. Этот способ не требовал знания производных точного решения, и давал апостериорную оценку.Систематическое изложение этого метода для различных классов задач можно найти в монографии Г.И. Марчука и В.В. Шайдурова [11].

Так же асимптотически точную оценку погрешности можно получить с помощью метода Эйткена [12], проведя расчёт на трёх сетках с шагами и/а,к/2 и и. В этом случае, в отличие от метода Ричардсона, не нужно знать порядок точности метода.

Позднее такие процедуры были построены для произвольного числа расчётов на последовательно сгущающихся сетках. Они позволяют получать результат, по погрешности эквивалентный расчёту со схемой гораздо более высокого порядка точности, чем исходная схема Общие формулы такого расчёта сравнительно громоздки. Но Ромберг

0 < ат < 1,0 < Ьт < 1,0 < ат1 < 1,0 < с™, < 1.

10) в 1955 г. нашёл простой вариант этих формул, если при каждом очередном сгущении шаг сетки уменьшается ровно вдвое.

Эти методы применимы практически к любым задачам с сеточным представлением функции - к интерполяции, численному дифференцированию, численному интегрированию, решению задач Коши и краевых задач для обыкновенных дифференциальных уравнений, задач для уравнений в частных производных, для интегральных и интегро-дифференциальных уравнений. Более того, до сих пор это единственный способ получить для сеточных решений асимптотически точную оценку погрешности. Они служат надёжной основой для построения программ с контролем точности. Однако классическая формулировка этого метода рассчитана на равномерные сетки. Практические же расчёты зачастую требуют неравномерных сеток, адаптированных к конкретным задачам. В монографии [13] рассмотрен важный класс неравномерных сеток - квазиравномерные сетки, на который метод сгущения полностью распространяется.

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

В методе невязок погрешность приближённого решения оценивается через норму невязки в пространстве образов оператора соответствующей краевой задачи.

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

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

Указанные методы можно использовать для построения адаптивных алгоритмов, когда сетка последовательно улучшается на основе информации, полученной на предыдущей сетке. Эти методы можно использовать на неравномерных и неструктурированных сетках, а также для задач, точное решение которых не обладает высокой гладкостью. В этом главное преимущество указанных методов перед методами Ричардсона и Эйткена.

Недостатком рассмотренных методов является то, что получаемые оценки для погрешности решения нельзя использовать для уточнения решения. Хорошо построенная адаптивная сетка улучшает точность расчётов в несколько раз, а одно рекуррентное уточнение по Ричардсону на несколько порядков.

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

В главе 3 данной работы проведено тестирование построенных разностных схем на сгущающихся сетках для нескольких задач, оценён эффективный порядок точности и погрешность полученных схем. Эти тестовые расчёты подтвердили заявленные свойства построенных схем.

Символьные вычисления

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

7)

Получение условий порядка - довольно трудоёмкая задача. Многие исследователи прибегали к некоторым упрощениям, например, положив а = 0 в формуле (7) (см. схему Розенброка-Ваннера)[2], что существенно снижает сложность выкладок, но влияет на потенциальные возможности получаемой схемы.

При современном развитии компьютерной техники, имеет смысл использовать ЭВМ для выполнения рутинной работы. Пакеты символьных вычислений уже давно стали неотъемлемой частью программного обеспечения таких как Maple, Matlab, Mathematika и других [14].

Некоторые математические пакеты содержат модули для построения условий порядка для явных методов. Примером таких модулей являются модули построения условий порядка методов Рунге-Кутта. Например, для среды Mathematika [15] или для среды Matlab [16]. Но ни один из пакетов не содержит модулей построения условий порядка для неявных или полуявных схем и не оперируют комплексной арифметикой.

В данной работе будет создана специальная программа для автоматического символьного построения условий порядка двустадийных комплексных схем Розенброка.

Пористый анодный оксид алюминия

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

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

В настоящее время пористые плёнки из анодированного оксида алюминия широко используются в качестве матриц для получения различных наноматериалов. Исследователи из Японии таким способом научились изготавливать не что иное как нано-пробирки (carbon nano test tubes, CNTTs), т.е. трубки, один конец которых открыт, а другой закрыт. Таким образом, внутрь такой пробирки можно что-нибудь поместить, например, магнитные частицы. Кроме того, образующиеся пробирки хорошо диспергируются в воде, что открывает множество их возможных применений, например, в области доставки лекарств.

Анодные оксидные плёнки образуются на поверхности большинства металлов и полупроводников. В основе процесса лежит реакция взаимодействия ионов окисляемого материала с молекулами воды. В зависимости от характера взаимодействия образующегося оксида с компонентами электролита имеется возможность управления структурой оксида. Так при отсутствии растворения оксида на поверхности электрода образуются беспористые барьерные плёнки. При частичном растворении оксида электролитом имеет место рост пористых плёнок. Среди последних наибольшее практическое значение имеют пористые анодные оксидные плёнки алюминия.

Мембраны анодированного оксида алюминия (АОА) обладают однородной пористой структурой с гексагональной упаковкой цилиндрических каналов и узким распределением пор по размерам. Использование различных электролитов, напряжений и времени анодирования позволяет варьировать диаметр пор (Dp), расстояние между порами (Ant) н толщину пленки (Lf) в широких пределах (Dp = 15-200 нм; Dint = 50-500 нм; L/до нескольких сотен микрон).

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

В настоящее время наибольшей популярностью пользуются два способа синтеза оксидных плёнок. Анодирование в мягких условиях (Mild Anodization - MA), включающее две стадии, протекает при малых значениях напряжения (U=40 В для щавелевой кислоты и U=25 В для серной) и характеризуется малой скоростью роста (порядка 2 мкм/час). Таким образом, получение толстых плёнок требует существенных временных затрат. Жёсткие условия анодирования (Hard Anodization - НА) требуют больших напряжений (до 180 В в щавелевой кислоте и до 80 В в серной), что увеличивает скорость роста до 50 мкм/час и позволяет получать большие расстояния между центрами пор для аналогичного электролита по сравнению с первым способом.

Подбирать режимы формирования пор экспериментальным путём - длительное и дорогостоящее занятие. Моделирование процесса на компьютере позволяет это сделать быстрее и гораздо дешевле.

Математическая модель, описывающая процесс образования периодических наноструктур, сводится к уравнению Курамото-Сивашинского.

Уравнение Курамото-Сивашинского

Уравнение Курамото-Сивашинского Ut = —aAU — bA2U + c(VC/)2 было выведено независимо Курамото [17] для описания образования диссипативных структур в системах с диффузией и Сиваишнским [18] для изучения гидродинамической неустойчивости при ламинарном горении. Позднее оказалось, что это уравнение описывает неустойчивость во многих других нелинейных процессах. В данной работе с помощью этого уравнения моделируется процесс образования нанопор на поверхности оксида алюминия.

В классической монографии Хайера-Ваннера [2] уравнению Курамото-Сиващинского даётся следующее. описание: "Первое слагаемое в этом уравнении отвечает за нарастание неустойчивости, но для колебаний высоких гармоник она стабилизируется вторым слагаемым. Нелинейное слагаемое связывает различные моды колебаний и обеспечивает ограниченность решение. Все это порождает удивительный хаос".

Одномерное уравнение Курамото-Сивашинского входит в список двенадцати задач-тестов для жёстких систем. При этом она отмечена как один из самых сложных тестов.

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

Общая теория уравнения Курамото-Сивашинского далека от завершения. Так условия глобального существования, а также условия разрушения решения за конечное время были получены Похожаевым [19], Галактионовым и Митидиери [20], [21] относительно недавно (в 2008 г).

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

Структура работы

Диссертация состоит из введения, четырёх глав и одного приложения, посвящённых изложению оригинальных результатов автора. Каждая глава разбита на разделы.

4.4 Заключение

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

Кроме того, в данной работе, благодаря калибровке модели, удалось добиться не только качественного, но и количественного совпадения результатов моделирования с результатами физического эксперимента.

В настоящий момент установлена зависимость параметра а2 в уравнении Курамото-Сивашинского (4.26) от напряжения на аноде при постоянной кислотности и температуре. В будущем, при наличии достаточного количества экспериментальных данных, возможно установить зависимость параметра а2 от кислотности рН и температуры Т. Это уточнение модели позволит лучше моделировать процесс образования периодических наноструктур в широком диапазоне напряжений, при различной кислотности среды и температуре.

1. Ваннер Г. Хайрер Э. Решение обыкновенных дифференциальных уравнений. Жесткие и дифференциально-алгебраические задачи. Мир, Москва, 1999. 685с. 5, 6, 7, 11, 13, 52, 55.

2. Н.Н. Калиткин. Численные методы. М.Наука, 1978. 512с. 6]

3. Е.А. Новиков. Одношаговые безитерационные методы решения жестких систем. Докторская диссертация. Красноярск, ВЦ АН СССР, 1991. 7]

4. Н.Н. Калиткин Д.С. Гужев. Оптимальная схема для пакета ros4. Матем. моделирование, 6(11):128-138, 1994. 7, 54.

5. C.JI. Панченко Н.Н. Калиткин. Оптимальные схемы для жестких неавтономных систем. Матпсм. моделирование, 11(6):52-81, 1999. 7.

6. Michel Roche. Rosenbrock methods for differential algebraic equations. Numerische Mathematik, 52:45-63, 1987. 10.1007/BF01401021. 7.

7. C. William Gear. Numerical Initial Value Problems in Ordinary Differential Equations. Prentice Hall PTR, Upper Saddle River, NJ, USA, 1971. 7.

8. H. H. Rosenbrock. Some general implicit processes for the numerical solution of differential equations. The Computer Journal, 5(4):329-330, 1963. 8.

9. Н.Н. Калиткин. Численные методы решения жестких систем. Матем. моделирование, 7(5):8—11, 1995. 9]

10. Шайдуров В.В. Марчук Г.И. Повышение точности решений разностных схем. М., Наука, 1979. 319с. 9]

11. А.С. Aitken. On interpolation by iteration of proportional parts, without the use of differences. Proc. Edinburgh Math. Soc. Second ser., 1932. Vol. 3, p. 56-76. 9, 50]

12. Е.А. Альшина В.В. Рогов Н.Н. Калиткин, А.Б. Альшин. Вычисления на квазиравномерных сетках. Физматлит, 2005. 224с. 10, 57]

13. W. Gander and J. Hrebicek. Solving problems in scientific computing using Maple and MATLAB. Springer, 2004. 11]

14. I. Th. Famelis, S. N. Papakostas, and Ch. Tsitouras. Symbolic derivation of runge-kutta order conditions. J. Symb. Compute 37(3):311-327, 2004. 11, 17]

15. Frank Cameron. A matlab package for automatically generating runge-kutta trees, order conditions, and truncation error coefficients. ACM Trans. Math. Softw., 32(2):274-298, 2006. 11]

16. Y. Kuramoto and T. Tsuzuki. On the Formation of Dissipative Structures in Reaction-Diffusion Systems —Reductive Perturbation Approach—. Progress of Theoretical Physics, 54:687-699, September 1975. 12]

17. G.I. Sivashinsky. Nonlinear analysis of hydrodynamic instability in laminar flaines-i. derivation of basic equations. Acta Astronáutica, 4(11-12):1177 1206, 1977. 12]

18. Похожаев С.И. О разрушении решений уравнения Курамото-Сивашинского. Ма-тем. сб., 199:98, 2008. 97-106. 13]

19. Похожаев С.И. Галактионов В.А., Митидиэри Э. Существование и отсутсвие глобальных решений уравнения Курамото-Сивашинского. Доклады Академии наук, 2008. т.419 номер 4 с.439-442. 13] ;

20. Е. Mitidieri V.A. Galaktionov and S.I. Pohozaev. On global solutions and blow-up for Kuramoto-Sivashinsky-type models, and well-posed Burnett equations. 13]

21. А.Г. Лимонов А.Б. Алынин, E.A. Алынина. Автоматизированное символьное построение условий порядка для двухстадийных комплексных схем типа Розенброка. Матем. моделирование, 21(10):76-88, 2009. 14]

22. А.Г. Лимонов А.Б. Алынин, Е.А. Альшина. Двухстадийные комплексные схемы Розенброка для жестких систем. Ж. вычисл. матем. и матем. физ., 49(2):270-287, 2009. 14, 26]

23. A. Al'shin, Е. Al'shina, and A. Limonov. Two-stage complex rosenbrock schemes for stiff systems. Computational Mathematics and Mathematical Physics, 49:261-278, 2009. 14]

24. A. B. Alshin, E. A. Alshina, and A. G. Limonov. Symbolic derivation of order conditions for two-stage complex rosenbrock scheme. AIP Conference Proceedings, 1048(l):47-51, 2008. 14]

25. А.Г. Лимонов. Моделирование образования гексагональных периодических наноструктур на поверхности оксида алюминия. Матем. моделирование, 22(8):97-108, 2010. 14]

26. Е. A. Alshina А.В. Alshin and A. G. Limonov. Symbolic Derivation of Order Conditions for Two-Stage Complex Rosenbrock Scheme. International Conference on Numerical Analysis and Applied Mathematics, 2008. 14]

27. E.A. Alshina A.B. Alshin and A.G. Limonov. Constructing Two-Stage Complex Rosenbrock Schemes By Using Automatic Symbolic Derivation Of Order Conditions. International Workshop Numerical Analysis and Scientific Computing (NASCom'08),

28. E.A. Альшина А.Б. Альшин А.Г. Лимонов, С.А. Гаврилов. Математическое моделирование процесса образования нано-пор на поверхности оксида алюминия. Международный Конкурс Научных Работ Молодых Ученых в Области Нанотехноло-гий, Москва, Экспоцентр, 2008. 14]

29. Е.А. Альшина А.Б. Альшин А.Г. Лимонов, А.Н. Белов. Математическое моделирование процесса образования нано-пор на поверхности оксида алюминия. Международный Конкурс Научных Работ Молодых Ученых в Области Нанотехнологий, Москва, Экспоцентр, 2009. 14]

30. Е.А. Альшина А.Б. Альшин А.Г. Лимонов, А.Н. Белов. Численное моделирование процесса образования нано-пор на поверхности оксида алюминия. Микроэлектроника и наноинженерия-2008. 14]

31. Ваннер Г. Хайрер Э., Нёрсетт С. Решение обыкновенных дифференциальных уравнений. Нежесткие задачи. Мир, Москва, 1990. 512с. 15]

32. J. Riordan. Introduction to Combinatorial Analysis. Dover books on mathematics. Dover Publications, 2002. 22]

33. П.Д. Ширков. Оптимально затухающие схемы с комплексными коэффициентами для жестких систем ОДУ. Матем. моделирование, 4(8):47-57, 1992. 26, 41]

34. J. P. O'Sullivan and G. С. Wood. The morphology and mechanism of formation of porous anodic films on aluminium. Proc. R. Soc., 317(1531):511-543, 1970. 59]

35. C. Sample and A. A. Golovin. Formation of porous metal oxides in the anodization process. Phys. Rev. E, 74(4):041606, Oct 2006. 59]

36. G. K. Singh, A. A. Golovin, and I. S. Aranson. Formation of self-organized nanoscale porous structures in anodic aluminum oxide. Phys. Rev. B, 73(20):205422, May 2006. 59, 65]

37. J.O.M. Bockris and A.K.N. Reddy. Modern electrochemistry. Number v. 2 in Modern Electrochemistry. Plenum Press, 2001. 60]

38. O. Jessensky, F. Miiller, and U. Gosele. Self-organized formation of hexagonal pore structures in anodic alumina. Journal of The Electrochemical Society, 145(11) :3735-3740, 1998. 61, 63, 65]

39. Feiyue Li, Lan Zhang, and Robert M. Metzger. On the growth of highly ordered pores in anodized aluminum oxide. Chemistry of Materials, 10(9):2470-2480, 1998. 61, 63]

40. A. P. Li, F. Miiller, A. Birner, K. Nielsch, and U. Gosele. Hexagonal pore arrays with a 50-420 nm interpore distance formed by self-organization in anodic alumina. Journal of Applied Physics, 84(ll):6023-6026, ' ,67]2008. 14.