Содержание к диссертации
Введение
1. Аналитические методы решения прямых и обратных задач для уравнения диффузии дробного порядка по времени с постоянным коэффициентом 20
1.1. Определения и основные свойства интегродифференциальных операторов дробного порядка 20
1.2. Уравнение диффузии дробного порядка по времени 28
1.2.1. Случайное блуждание при непрерывном времени 28
1.2.2. Вывод уравнения диффузии дробного порядка по времени 31
1.3. Аналитический подход к решению прямых и обратных задач для уравнения диффузии дробного порядка по времени 35
1.3.1. Решение краевой задачи без начальных условий с периодическим источником для уравнения диффузии дробного порядка по времени 35
1.3.2. Фундаментальное решение уравнения диффузии дробного порядка по времени 40
1.3.3. Решение первой краевой задачи для уравнения диффузии дробного порядка по времени 43
1.3.4. Связь между фундаментальным решением уравнения диффузии дробного порядка по времени и функцией Грина первой краевой задачи 45
1.3.5. Решение краевой задачи с однородными граничными условиями для уравнения диффузии дробного
порядка по времени методом Фурье 47
1.3.6. Постановка и решение обратных задач для уравнения диффузии дробного порядка по времени 49
1.4. Выводы 52
2. Численные методы решения прямых и обратных задач для уравнения диффузии дробного порядка по времени с постоянным коэффициентом 54
2.1. Метод конечных разностей численного решения краевых задач для уравнения диффузии дробного порядка по времени с постоянным коэффициентом 55
2.1.1. Разностная схема для уравнения диффузии дробного порядка по времени с постоянным коэффициентом 56
2.1.2. Обобщенный метод прогонки решения первой краевой задачи для уравнения диффузии дробного порядка
по времени с постоянным коэффициентом 58
2.1.3. Получение оценок решения разностной краевой задачи 60
2.1.4. Исследование устойчивости разностной схемы
для уравнения диффузии дробного порядка по времени 63
2.1.5. Определение порядка аппроксимации разностной схемы для уравнения диффузии дробного порядка по времени 67
2.2. Метод Монте-Карло численного решения краевых задач для уравнения диффузии дробного порядка по времени с постоянным коэффициентом 71
2.2.1. Связь дробно-устойчивых распределений с фундаментальным решением уравнения с дробной производной по времени 73
2.2.2. Классическая модель случайного блуждания 76
2.2.3. Построение дискретной модели случайного блуждания для диффузии дробного порядка по времени 78
2.2.4. Разложение случайного блуждания на диффузионную и дисперсионную составляющие 81
2.3. Результаты вычислительных экспериментов 85
2.3.1. Численное решение краевой задачи с однородными граничными условиями 85
2.3.2. Решение уравнения диффузии дробного порядка по времени методом Монте-Карло 86
2.4. Выводы 89
3. Численные методы решения прямых и обратных задач для уравнения диффузии дробного порядка по времени с переменным коэффициентом 90
3.1. Метод конечных разностей численного решения уравнения диффузии дробного порядка
по времени с переменным коэффициентом 90
3.1.1. Интегро-интерполяционный метод построения разностной схемы для уравнения диффузии дробного порядка с переменным коэффициентом 90
3.1.2. Обобщенный метод прогонки численного решения краевой задачи для разностного уравнения диффузии дробного порядка по времени с переменным коэффициентом 95
3.1.3. Получение оценок решения разностной краевой задачи 98
3.1.4. Исследование устойчивости разностной схемы для уравнения диффузии дробного порядка
по времени с переменным коэффициентом 99
3.2. Метод минимизации функционала невязки решения обратных задач для уравнения диффузии дробного порядка по времени с переменным коэффициентом 102
3.2.1. Постановка прямых и обратных задач 102
3.2.2. Свойства решения разностной прямой задачи 103
3.2.3. Построение и основные свойства функционала невязки 107
3.2.4. Метод Левенберга—Марквардта минимизации функционала невязки 113
3.2.5. Метод секущих минимизации функционала невязки . 118
3.2.6. Метод Флстчера—Ривса минимизации функционала невязки 120
3.3. Результаты вычислительных экспериментов 124
3.3.1. Выбор критериев останова алгоритмов в методе минимизации функционала невязки 125
3.3.2. Исследование свойств функционала невязки 127
3.3.3. Численное решение обратных задач для уравнения диффузии дробного порядка по времени 132
3.3.4. Случай постоянного коэфициента 136
3.4. Выводы 137
4. Эволюционные методы решения обратных задач для уравнения диффузии дробного порядка по времени с переменным коэффициентом 139
4.1. Применение генетических алгоритмов для минимизации функционала невязки 139
4.2. Применение алгоритма поиска по шаблону для улучшения приближенного решения 145
4.3. Результаты вычислительных экспериментов 147
4.3.1. Численное решение обратных задач для уравнения диффузии дробного порядка по времени 147
4.3.2. Сравнительный анализ эффективности классических и эволюционных алгоритмов 148
4.3.3. Пример последовательного использования ГА и метода Левенберга—Марквардта 150
4.4. Выводы 151
Заключение 152
Список использованных источников
- Случайное блуждание при непрерывном времени
- Разностная схема для уравнения диффузии дробного порядка по времени с постоянным коэффициентом
- Обобщенный метод прогонки численного решения краевой задачи для разностного уравнения диффузии дробного порядка по времени с переменным коэффициентом
- Применение алгоритма поиска по шаблону для улучшения приближенного решения
Введение к работе
В последнее время значительный интерес представляют физические исследования диффузионных процессов аномальной природы, отклоняющихся от классической гауссовской диффузии, которые встречаются во множестве физических систем (см. [77, 78] и цитированную там литературу). Для аномального диффузионного процесса характерно в первую очередь то, что зависимость среднеквадратического смещения от времени
отклоняется от «нормального» линейного закона (ж2()) = 2Kit. Здесь К\ и Ка — «обычный» и обобщенный коэффициенты диффузии размерности см2 с-1 и см2с~а соответственно (при а = 1 имеем Г(2) = 1). Показатель аномальной диффузии аф\ определяет, будет ли процесс классифицирован как субдиффузионный (дисперсионный, медленный) при 0 < а < 1 или супердиффузионный (ускоренный, быстрый) при а > 1. Обычно рассматривается область 1 < а < 2, где а = 2 — баллистический предел, описываемый волновым уравнением.
Соотношение (0.1) описывает так называемые «странные» процессы переноса в нелинейных динамических системах [21], то есть негаус-совые процессы, допускающие корреляции на сколь угодно больпшх пространственно-временных масштабах. Так, например, броуновское движение, для которого функция плотности вероятности имеет вид
1 f х1 \
Рвм{х,) = —===ехр
допускает обобщение с помощью модели случайного блуждания при непрерывном времени (Continuous Time Random Walk — CTRW) на случай субдиффузии или переноса с дисперсией, где
PsD{x,t) ~ СіГ«/2-(1-а)/(2-а) ехр
при f = И/Г/2, 0 < а < 1.
_С22/(2-а)
Рассматриваемым моделям могут быть поставлены в соответствие дифференциальные уравнения дробного порядка, из которых видно, что данные процессы являются сильно нелокальными и характеризуются широкими корреляциями во времени и/или пространстве, представимыми в виде медленно убывающих но степенному закону ядер в соответствующих интегро-дифференциальных уравнениях.
В работе [78] подробно рассматриваются примеры систем из различных областей (геофизики, геологии, физики, химии, астрономии, экономики), в которых были обнаружены процессы аномальной природы (см. также цитированную там литературу). В частности, с точки зрения аномальной диффузии было предложено объяснение водородного эффекта в структуре силиконовых электродов при некоторых электрохимических условиях, так же, как и в контексте нелинейного электрофореза. Методы дробного исчислеїшя применялись при анализе процессов аномальной диффузии, обнаруженных в аморфных электроактивных материалах. Аномальная диффузия катионов была обнаружена при изучении механизмов роста образцов наплавленного оксида молибдена, кинетики переноса электронов в пленках из поли-3,4-этилендиокситиофена, а также атомных процессов переноса и химических реакций в диэлектрических пленках с высокой диэлектрической постоянной. Дробная динамика (fractional dynamics) может лежать в основе статистики объединенной функции плотности вероятности «скорость—координата» движения частицы в турбулентном поле. Обобщение закона Ричардсона с позиций дробного исчисления было предложено для описания переноса воды в ненасыщенном грунте.
Функция плотности вероятности для диффузии дробного порядка описывает самодиффузионные профили на перколирующих фрактальных структурах. С этих же позиций рассматривается 1//а-шум (т. е. сигнал или процесс со степенной спектральной плотностью, пропорциональной l/fa, где / — частота) и соответствующее промежуточное состояние при моделировании молекулярной динамики замерзания воды. Замечено также, что для процессов «старения» в неупорядоченных системах, равно как и в динамических системах, характерна «нелокальность» во времени, подчи-
няющаяся степенному закону вида (0.1).
Процессы вида (0.1) обнаружены, в частности, в электроактивных материалах [53], лабораторной плазме, турбулентных жидкостях выше некоторого порогового значения числа Рейнольдса, замагпиченных вихревых потоках (см. [21] и указанные там ссылки), в средах с временной дисперсией, таких как биологические ткани и материалы, имеющие самоподобную «архитектуру». Известны многочисленные примеры физических структур, обладающих самоподобием на уровне, промежуточном между микро- и макроуровнями (на мезоуровне), которые стали объектами пристального внимания со стороны физиков-экспериментаторов [23]. Оказалось, что материалы, мезоскопическая структура которых обладает свойством масштабной инвариантности всего 5—8 порядков, имеют уникальные физические свойства, являющиеся результатом их внутренней самоподобной «архитектуры» [21]. В связи с этим актуальной становится проблема построения адекватных математических моделей таких сред.
С макроскопической точки зрения, диффузионный процесс описывается уравнением диффузии
ut(x,t) ~ Kiuxx(x,t), (0.2)
где u(x,t) представляет собой плотность вероятности обнаружить частицу в точке х в момент времени t. С микроскопической точки зрения, диффузия представляет собой марковский процесс, в котором микроскопические частицы выполняют случайные «прыжки» конечной длины и конечной дисперсии. С другой стороны, если рассматривается немарковский процесс, в котором «прыжки» частиц, выбираются из распределения с длинным временным хвостом t~a~l, то диффузионный процесс является аномальным (см., например, работу G. Rangarajan и М. Ding [96]). Как показано в работах М. Giona и Е. Roman [57], R. Metzler и J. Klafter [77, 78], W. Schneider [98], W. Wyss [104], R. Gorenflo и F. Mainardi [61], А.И. Саиче-ва и С.Г. Уткина [38, 39], функция плотности вероятности u(x,t), которая описывает движение частиц в случае аномальной диффузии, удовлетворя-
10 ет уравнению с дробной производной вида
Уравнения вида (0.3) часто встречаются в современной литературе при описании различных аномальных диффузионных процессов [77, 78, 108], в частности процессов, протекающих в сильно пористых (фрактальных) средах [86, 102]. Показатель а может интерпретироваться как характеристика фрактальной размерности «активного» времени [21], в котором реальные блуждания частиц выглядят как случайный процесс; интервал активного времени пропорционален ta.
Важный класс объектов с самоподобной структурой образуют множества, описывающие геометрию перколяции (или протекания), то есть случайного распространения жидкости через среду, причем абстрактные понятия «жидкость» и «среда» могут быть интерпретированы в соответствии с физическим смыслом рассматриваемой задачи. Диффузионные процессы на перколирующих фрактальных структурах существенно не гауссовы и не согласуются с традиционными представлениями о процессах переноса как о случайном броуновском движении частиц в среде. Уравнения вида (0.3) учитывают эффекты памяти и нелокальное, выходящие далеко за пределы традиционной гауссовой статистики (0.2).
Примеры физических систем, которые описываются в терминах дробного исчисления, приведены в работах Л. М. Зеленого и А. В. Мило-ванова [21], A. Le Mehaute [69], F. Mainardi [73], P. P. Нигматулли-на [87, 88, 89, 90], А. И. Олемского [93], R. Hilfer [62], I. Sokolov [100]. Однако в соответствующей литературе часто отсутствует достоверная информация о результатах проверок адекватности получаемых моделей в ходе реальных экспериментов. В диссертации предлагается подход, использующий предположение, что процесс аномальной диффузии описывается уравнением с дробной производной по времени вида (0.3), однако порядок дифференцирования не известен, то есть ставится задача нахождения вида уравнения аномальной диффузии по реально измеряемым данным. Такие задачи не являются классическими коэффициентными обратными задачами для
уравнений в частных производных, вид которых заранее известен [36, 37]. В книгах А.Н. Тихонова и А.А. Самарского [45] и А. Зоммерфельда [22] рассматривалась известная задача Фурье о распространении температурных волн в почве и было показано, что получить информацию о строении среды можно, измерив ее реакцию на воздействие температурной волны. С другой стороны, в ряде оптических экспериментов по исследованию распространения света в биологических средах было замечено, что источник света с периодически модулированной интенсивностью порождает волну плотности энергии. Эта волна, распространяясь в оптически плотной среде, имеет сферический волновой фронт, форма которого зависит от функции плотности среды. Этому явлению было дано название «диффузионная волна» [74]. Таким образом, по реакции среды на воздействие диффузионной волны можно определить основные параметры среды. Применение данного подхода на мезоуровне позволяет по данным измерений провести проверку адекватности предлагаемой математической модели и определить вид соответствующего дифференциального уравнения.
Актуальность темы диссертации определяется возрастающим интересом к исследованиям в области мезоскопического моделирования процессов переноса. Трудности, связанные с применением аналитических методов решения прямых задач для неоднородных сред па мезоуровне, приводят к необходимости построения эффективных численных методов, чему способствуют и высокие темпы развития компьютерной техники. Разработанные в диссертации методы численного решения обратных задач позволяют на основе реальных данных построить математическую модель аномального диффузионного процесса.
Цель работы заключается в построении точных решений обратных задач для уравнения диффузии дробного порядка по времени с постоянным коэффициентом, а также разработке численных методов решения прямых и обратных задач для уравнения диффузии дробного порядка по времени с переменным коэффициентом.
В рамках указанной цели были поставлены следующие задачи:
1. Получить аналитическое решение краевой задачи без начальных усло
вий для дифференциального уравнения с дробной производной по времени,
на основании которого построить точные формулы нахождения обобщен
ного коэффициента диффузии и показателя аномальной диффузии.
Разработать метод статистического моделирования (Монте-Карло) для численного решения начально-краевых задач для уравнения диффузии дробного порядка по времени с постоянным коэффициентом.
На основании классической теории разностных схем построить разностную схему с весами и предложить обобщенный метод прогонки численного решения начально-краевых задач для дифференциального уравнения с дробной производной по времени с постоянным и переменным коэффициентом.
Получить условия устойчивости разностных схем и оценки решений разностных краевых задач для дифференциального уравнения с дробной производной по времени с постоянным и переменным коэффициентом.
Разработать модифицированные методы минимизации функционала невязки для численного решения обратных задач, состоящих в определении (переменного) обобщенного коэффициента диффузии и показателя аномальной диффузии, и провести их сравнительный анализ.
Методика исследования. При решении поставленных задач использовались методы математической физики, дифференциальных уравнений в частных производных, конечных разностей, статистического моделирования, теории вероятностей, дробного, операционного и вариационного исчислений, а также классические и эволюционные методы безусловной оптимизации.
Научная новизна работы состоит в следующем:
1. Впервые получено решение краевой задачи без начальных условий с периодическим источником методом разделения переменных и построено интегральное преобразование, связывающее решение краевой задачи
без начальных условий с периодическим источником для уравнения диффузии дробного порядка по времени с решением аналогичной задачи для уравнения параболического типа.
Предложены новые постановки обратных задач для уравнения диффузии дробного порядка по времени, заключающихся в восстановлении обобщенного коэффициента диффузии и дробного показателя диффузии. В некоторых простых случаях решения обратных задач представлены в виде точных аналитических формул, а в общем случае предложен алгоритм численного решения.
На основании существующей дискретной модели случайного блуждания для уравнения диффузии дробного порядка по времени разработан метод Монте-Карло численного решения первой краевой задачи. Впервые в явном виде получено разложение случайного блуждания на диффузионную и дисперсионную составляющие.
Впервые для уравнения диффузии дробного порядка по времени с переменным коэффициентом с помощью интегро-интерполяционного метода построена разностная схема с весами и разработан обобщенный метод прогонки численного решения первой краевой задачи и получены равномерные оценки решения. Исследована устойчивость разностной схемы с весами для уравнения диффузии дробного порядка по времени с переменным коэффициентом и условие устойчивости получено в явном виде.
Впервые предложены постановки обратных задач для уравнения диффузии дробного порядка по времени, заключающихся в восстановлении обобщенного коэффициента диффузии и дробного показателя диффузии, и разработан комплекс программ и реализованы оптимизационные методы ньютоновского типа, сопряженных градиентов, а также эволюционные алгоритмы для их решения.
Теоретическое значение работы заключается в том, что предложенные в ней постановки обратных задач являются новыми, а разработанные для их решения методы представляют собой обобщение существующих подходов на случай уравнений с дробной производной по времени.
Практическая ценность работы состоит в том, что результаты исследований, проведенных в диссертации, позволят на практике по данным измерений строить математические модели сред с неизвестными характеристиками.
Достоверность и обоснованность результатов подтверждена сравнением результатов решения прямых задач теории аномальной диффузии сеточными методами, методом Монте-Карло и аналитическими методами.
На защиту выносятся следующие положения:
Разработаны методы аналитического решения обратных задач для уравнения диффузии дробного порядка по времени с постоянным коэффициентом, заключающихся в восстановлении его параметров: коэффициента и порядка временной производной.
Разработан алгоритм статистического моделирования (метод Монте-Карло) для решения прямых задач теории аномальной диффузии в однородных средах.
Разработан метод конечных разностей численного решения краевых задач для уравнения диффузии дробного порядка по времени как с постоянным, так и с переменным коэффициентом.
С помощью классической техники получены условия устойчивости разностных схем с весами для уравнения диффузии дробного порядка по времени с постоянным и переменным коэффициентом, а также оценки решений краевых задач для данного уравнения.
Представлены новые постановки обратных задач для уравнения диффузии дробного порядка по времени с переменным коэффициентом, заключающихся в восстановлении его параметров: коэффициента и порядка временной производной. Для их решения разработаны и реализованы в рамках программного комплекса INPRANDI оптимизационные алгоритмы ньютоновского типа, представляющие собой модификации методов Левенберга—Марквардта, секущих и Флетчера—Ривса, а также гибридные эволюционные алгоритмы.
Апробация работы. Основные положения диссертации и отдельные ее результаты докладывались и обсуждались в рамках семинаров, проводимых на кафедре Высшей Математики НГТУ; на семинаре академика В.Н. Монахова в Институте гидродинамики СО РАН; на семинарах чл.-корр. В.Г. Романова и проф. A.M. Блохина в Институте математики им. С.Л. Соболева СО РАН, а также на следующих конференциях:
Региональная научная конференция «Наука. Техника. Инновации». Новосибирск, 2001, 2002 гг.
Всероссийская научная конференция молодых ученых «Наука. Технологии. Инновации». Новосибирск, 2003, 2004, 2005, 2006 гг.
Korea—Russia International Symposium on Science and Technology. Novosibirsk, 2002; Tomsk, 2004; Novosibirsk, 2005.
Международная конференция «Дифференциальные уравнения, теория функций и приложения», посвященная 100-летию со дня рождения академика Ильи Несторовича Векуа. Новосибирск, 2007 г.
Всероссийская конференция но вычислительной математике КВМ-2007. Новосибирск, 2007 г.
Международная конференция «Обратные и некорректные задачи математической физики», посвященной 75-летию академика М. М. Лаврентьева. Новосибирск, 2007 г.
Публикации. По теме диссертации опубликовано 15 работ, в том числе 1 работа в журнале из перечня ВАК.
Структура работы. Диссертация состоит из введения, четырех глав, заключения, приложения и списка литературы из 109 наименований. Общий объем диссертации составляет 187 страниц, в том числе основной текст 163 страницы.
Краткое содержание работы.
В первой главе разработан аналитический подход к решению прямых и обратных задач для уравнения диффузии дробного порядка по времени с постоянным коэффициентом. В разд. 1.1 представлены определения и основные свойства интегродифференциальных операторов дробного поряд-
ка, используемые в диссертации. В разд. 1.2 представлен вывод уравнения диффузии дробного порядка по времени на основе концепции случайного блуждания при непрерывном времени (CTRW). В разд. 1.2.1 приводится описание концепции CTRW [64, 81, 82, 83, 97, 103], а в разд. 1.2.2 дан обзор существующих методов вывода уравнения диффузии дробного порядка по времени [77]. В разд. 1.3 обсуждается аналитический подход к решению прямых и обратных задач для уравнения диффузии дробного порядка по времени с постоянным коэффициентом. В разд. 1.3.1 представлена реализация метода разделения переменных решения краевой задачи без начальных условий с периодическим источником [45] для рассматриваемого уравнения, а в разд. 1.3.2 на основании существующих методов [75] получено его фундаментальное решение. В разд. 1.3.3 рассматривается первая краевая задача для уравнения диффузии дробного порядка по времени, на основании существующих методов [75] построена функция Грина данной задачи, и решение первой краевой задачи представлено в виде свертки функции Грина с функцией, определяющей граничное условие. В разд. 1.3.4 показано, что решение краевой задачи без начальных условий с периодическим источником может быть представлено в виде свертки специального вида функции, определяющей граничное условие, с ядром, представляющим собой функцию Грина первой краевой задачи, а также установлена связь между функцией Грина и фундаментальным решением уравнения диффузии дробного порядка по времени. В разд. 1.3.5 рассматривается краевая задача с однородными граничными условиями для уравнения диффузии дробного порядка по времени и строится ее решение с использованием метода Фурье. В разд. 1.3.6 представлены постановки обратных задач для уравнения диффузии дробного порядка по времени с постоянным коэффициентом и решения представлены в виде точных аналитических формул.
Во второй главе рассматриваются численные методы решения прямых и обратных задач для уравнения диффузии дробного порядка по времени с постоянным коэффициентом. В разд. 2.1 на основе классической теории разностных схем [40] разработан метод конечных разностей численного решения краевых задач для рассматриваемого уравнения. В разд. 2.1.1 по-
строєна разностная схема с весами для уравнения диффузии дробного порядка по времени с постоянным коэффициентом, а в разд. 2.1.2 разработан метод прогонки для ее компьютерной реализации. В разд. 2.1.3 с помощью принципа максимума и формул прогонки получены оценки решений ряда разностных краевых задач. В разд. 2.1.4 аналитически проведено исследование устойчивости схемы с весами и условие устойчивости получено в явном виде. В разд. 2.1.5 на основании существующих результатов [40, 56] в случаях явной и неявной схем определена погрешность разностной аппроксимации как в произвольном внутреннем узле сеточной области (локально), так и на всей сетке. В разд. 2.2 рассматривается метод Монте-Карло численного решения краевых задач для уравнения диффузии дробного порядка по времени с постоянным коэффициентом. В разд. 2.2.3 на основании результатов [59] построена дискретная модель случайного блуждания для диффузии дробного порядка по времени, а в разд. 2.2.4 представлено разложение данного случайного блуждания на диффузионную и дисперсионную составляющие. В разд. 2.3 представлены результаты вычислительных экспериментов: в разд. 2.3.1 представлены результаты численного решения краевой задачи с однородными граничными условиями для уравнения диффузии дробного порядка по времени; в разд. 2.3.2 приводятся результаты решения уравнения диффузии дробного порядка по времени методом Монте-Карло.
В третьей главе рассматриваются численные методы решения прямых и обратных задач для уравнения диффузии дробного порядка по времени с переменным коэффициентом. В разд. 3.1 на основе классической теории разностных схем [40] разработан метод конечных разностей численного решения краевых задач для рассматриваемого уравнения. В разд. 3.1.1 представлен интегро-интерполяционный метод построения разностной схемы с весами для уравнения диффузии дробного порядка по времени с переменным коэффициентом, а в разд. 3.1.2 разработан метод прогонки для ее компьютерной реализации. В разд. 3.1.3 с помощью формул прогонки получены оценки решения разностной краевой задачи для рассматриваемого уравнения. В разд. 3.1.4 аналитически проведено исследование устойчи-
вости схемы с весами и условие устойчивости получено в явном виде. В разд. 3.2 рассматривается метод минимизации функционала невязки численного решения обратных задач для уравнения диффузии дробного порядка по времени с переменным коэффициентом. В разд. 3.2.2 исследуются свойства решения разностной прямой задачи и получена для него оценка. В разд. 3.2.3, опираясь на результаты [19], построен функционал невязки, исследованы его основные свойства и представлена общая схема минимизации. В разд. 3.2.4 и 3.2.5 рассматриваются метод Левепберга—Марквардта и метод секущих минимизации функционала невязки. В разд. 3.2.6 рассматривается метод Флетчера—Ривса минимизации функционала невязки, в рамках которого используются метод Ньютона и метод золотого сечения одномерной минимизации. В разд. 3.3 представлены результаты вычислительных экспериментов: в разд. 3.3.1 исследуется проблема выбора критерия останова алгоритма минимизации, в разд. 3.3.2 изучаются дальнейшие свойства функционала невязки, в частности, топографическая структура поверхности, им определяемой. В разд. 3.3.3 приведены результаты численного решения обратных задач, состоящих в восстановлении обобщенного коэффициента диффузии q(x) и показателя аномальной диффузии а, и проведен сравнительный анализ эффективности используемых алгоритмов; в разд. 3.3.4 — результаты численного решения обратных задач, состоящих в восстановлении обобщенного (постоянного) коэффициента диффузии А и показателя аномальной диффузии а, когда данные обратной задачи рассчитываются методом Монте-Карло.
В четвертой главе рассматриваются эволюционные методы решения обратных задач для уравнения диффузии дробного порядка по времени с переменным коэффициентом. В разд. 4.1 представлены общая схема и принцип работы эволюционного алгоритма, а в разд. 4.2 рассматривается алгоритм поиска по шаблону, который может быть использован для улучшения полученного приближенного решения. В разд. 4.3 представлены результаты вычислительных экспериментов: в разд. 4.3.1 показаны решения обратных задач восстановления обобщенного коэффициента диффузии, иллюстрирующие преимущества генетических алгоритмов; разд. 4.3.2 содер-
жит результаты сравнительного анализа эффективности эволюционных и классических алгоритмов, а в разд. 4.3.3 приведен пример их последовательного использования.
В заключении представлены основные результаты диссертационной работы.
Случайное блуждание при непрерывном времени
При классическом броуновском блуждании на прямой предполагается, что частица через равные промежутки времени At делает прыжок длины Лес в один из соседних узлов сетки, который выбирается случайным образом. Процесс подобного типа, локальный во времени и пространстве, можно смоделировать при помощи так называемого мастер-уравнения [77]: , . х и(х — АхЛ) и(х + АхЛ) и(х, t + At) = -А —-U. + -І _2-2. (1.39)
Уравнение (1.39) определяет плотность распределения вероятности того, что частица попадет в точку х в момент времени t + At, в зависимости от того, в каком из соседних узлов х ± Ах она находится в момент t. Коэффициент 1/2 означает, что выбор одного из соседних узлов равновероятен. При переходе к пределу при At — 0 и Ах — 0 разложения в ряд Тейлора функции и относительно At и Ах д и(х, t + At)= и(х, t) + At o7 (z, t) + О ([At]2) (1.40) и д (Ах)2 д2 и(х ± Ах, t) = и(х, t) ± Ах j-u(x, t) + 2 2и(х, t) + О ([Аж]3) (1.41) приводят к уравнению диффузии -u(x,t) = K1—u(x,t), (1.42) где Kl = lim &- (L43) AxAt- o 2 At — коэффициент диффузии. Для больших интервалов времени, соответствующих достаточно большому числу прыжков, плотность распределения вероятности пребывания в точке х в момент времени t дается уравнением диффузии (1.42) и имеет вид 1 / х2 \ и(хЛ) = , ехр( -——- . (1-44)
При использовании подхода на основе концепции случайного блуждания при непрерывном времени предполагается, что длина прыжка частицы, так же как и время ожидания между двумя последовательными прыжками, выводятся из плотности распределения вероятности прыжка if)(x,t). Таким образом, плотность распределения длины прыжка и плотность распределения времени ожидания имеют вид 00 оо М = /«х,«) и ,(«)=/ (,, ) (1.45) О -оо соответственно. Значит, і?(ж) dx дает вероятность того, что длина прыжка лежит в интервале (ж, ж + dx), a w(t)dt дает вероятность того, что время ожидания лежит в интервале (t,t + dt). Если длина прыжка и время ожидания являются независимыми случайными величинами, то ф{х,Ь) = d{x)w(t). При классификации различных типов процессов CTRW обращают внимание на то, являются ли характеристическая функция времени ожидания оо т = и дисперсия длины прыжка Е2 = w(t)tdt (1.46) / ${x)x2dx (1.47) -со зо конечными или бесконечными величинами. Процесс CTRW описывается следующим мастер-уравнением [67]: 00 оо r}{x,t) = [ dx J ф ,ї) ф(х - x ,t - t ) dt + S(x)5(t), (1.48) -сю О которое связывает величину rj(x:t) — плотность распределения вероятности попасть в точку х в момент времени t — с величиной r)(x ,t )} отвечающей за попадание в точку х в момент t . Второе слагаемое в уравнении (1.48) характеризует начальное условие случайного блуждания, в данном случае выбранное в виде 5(х). Следовательно, плотность распределения вероятности нахождения в точке х в момент времени t дается формулой и (ж, t) = f т/(ж, t )V{t - t ) dt . (1.49)
Здесь подразумевается, что частица попала в точку х в момент времени t и с тех пор не двигалась. Последнее определяется интегральной вероятностью t Ф( ) = 1- f w{t )dt , (1.50) о соответствующей вероятности того, что прыжок частицы не произойдет в течение временного интервала (0,). В работе [67] показано, что u(x,t) удовлетворяет следующему соотношению: з 1-V(f,s) где й(, s) — Фурье- и Лаплас-образ функции и(х, і) и щ{) — Фурье-образ начального условия и(х, 0); в случае (1.48) имеем йо() — 1 В качестве примера рассмотрим модель CTRW для броуновского движения в случае, когда длина прыжка и время ожидания — независимые случайные величины, т. е. i{j(x,t) — &(x)w(t), и величины Т и Е2 — конечны. Пусть время ожидания имеет пуассоновскую плотность распределения w(t) = -ехр (--) , (1-52) причем Т = т, а длина прыжка — гауссовскую плотность распределения )=тйрехр(-3 (1-5з) причем Е2 = 2сг2. Тогда соответствующие преобразования Лапласа и Фурье имеют вид С {w(t); s} = —i-j 1 - sr + 0(т2) , (1.54) W ); Д = ехР (-eV) - 1 - V + 0(a4) . (1.55)
Действительно, любая пара функций распределения, при которой величины ТиЕ2 конечны, приводит к тому же результату [64, 67]. Подставляя (1.54) и (1.55) в (1.51), при начальном условии и(х,0) 5(х) в пределе (С? s) (0;0) получим Фурье- и Лаплас-образ функции распространения: Заметим, что при обратном переходе к переменным (ж, t) получим гауссовскую функцию распространения (1-44), а преобразовав выражение (1.56) к виду -КгЄй( s) = su(, s) - «o(0 (1-57) и воспользовавшись свойствами преобразований Фурье (1-21) и Лапласа (1.36), получим уравнение диффузии (1.42).
Разностная схема для уравнения диффузии дробного порядка по времени с постоянным коэффициентом
Рассмотрим по аналогии с [59] дифференциальное уравнение с дробной производной по времени: д2 tV%+u(x,t) = \z7nu(x,t), 0 а 1, (2.7) где х Є R, t Є RQ , а оператор fX g+ — это оператор дробного дифференцирования Капуто [94, 59]. Будем рассматривать краевую задачу для уравнения (2.7) в прямоугольнике Q = {Q x 1, 0 t Т}. (2.8) Требуется найти непрерывное в прямоугольнике Q решение и = u(x,t) задачи д2 tV%+u{x,t) = Л2—-J и(х, t), 0 х 1, 0 t Т, (2.9) м(ж,0) = щ(х), 0 ж 1, (2.10) гх(о, ) = 4 i(t), и{1,г) = ф2{г), o t T. (2.11) Введем сетки uh = {ХІ = гТг, і = 0,1,..., iV, /г = 1/iV} , (2.12) "т = ft = З , 3 = 0,1,..., М, г = Т/М} (2.13) и сетку в Q,: - -} 0, = cjhx шт t{Xi,tj) ХІ Є wh, tj Є и) (2.14) Обозначим через y\ значение в узле (xi,tj) сеточной функции у, определенной на Г2. Заменяя производную t o+ ее разностным аналогом (2.6), а д2/дх2 — второй разностной производной и вводя вещественный параметр О а 1, получим разностную схему с весами VI (а) = ЛА (ayj+1 + (1- а)уї) , (2.15) где оператор А\ действует по правилу г Ллй - L УІ-г - 2yl + y3i+l (2.16) Заметим, что данная разностная схема является многослойной с переменным числом слоев. Начальные и краевые условия аппроксимируем точно: У і =У{ХІ,0) =Щ(ХІ), (2.17) УІ = ФІ Ум = ФІ (2-18) Рассмотрим схему с весами (2.15), где в случае а = 0 получим формулу, позволяющую явным образом выразить значения у на следующем слое через значения на текущем слое: Я-1 УІ = а а AV Ь2 УІ + Q h2 (й-і + !/Ін)+ 7, (2-19) k=0 V fc= j+l-k [?/: (2.20) При (7 0 приведем выражение (2.15) к каноническому виду Я-1 2/i — л V+1 _ r7V+1 4- 7? V+1 — — Fj — АгУі-1 гУі + ПгУм — -г,- , (2.21) а именно vj+1 -Уі-і h2 Л2Та + J+1 . -.Я-i Уг Уг+1 -F3 (2.22) ah/ h2 F? = 1 а\2та 1-(7 Y? + сг\2та a При «7=1 (чисто неявная схема) имеем h2 + 2 2/І-І j+i . j+i Л2т« УЇ + a /г5 rL Л2 УІі + УІ+і) - (2.23) (2.24) У І + «2// На рис. П. 1, стр. 165 представлены шаблоны разностной схемы с весами (2.15) для случаев а=1и0 а 1.
Обобщенный метод прогонки решения первой краевой задачи для уравнения диффузии дробного порядка по времени с постоянным коэффициентом
Рассмотрим для разностного уравнения (2.22) краевую задачу (2.17), (2.18). Для определения yj+ нужно решить линейную систему уравнений (2.22). Матрица данной системы трехдиагональная, поэтому для решения системы воспользуемся методом прогонки. Заметим сразу, что элементы матрицы, лежащие на главной диагонали, преобладают: \СІ\ \АІ\ + \ВІ\. Это условие обеспечивает единственность разностного решения и устойчивость прогонки [26].
Построим алгоритм, реализующий метод прогонки численного решения прямых задач для уравнений вида (2.7). Предположим, что J+1 .7 + 1 і (2.25) и, соответственно, i+i i+i , (2.26) Подставим (2.26) в (2.22): ,і+і т-іУі + т-і tf аХ2та + УІ+1+УІ:І = -ГІ (2.27) i+l _ Уі+1 Уі Тогда і УІЇІ + Vi-i + Fj p + 2- fii-i P h2 crAV (2.28) откуда следует, что № p + 2 - m-i (2.29) здесь г = 2,3,..., iV — 1. Возьмем г = 1: p + 2 - /ij_i yV (p + 2)yi+i+vi+i = -F{, (2.31) откуда уГ = + , (,32) и, следовательно, « = 2 щ = ТТ (2 33) Из (2.26) также следует, что УІЇ-i = Им-іФІ+1 + m-i- (2-34) Таким образом, для обобщенного метода прогонки выведены формулы «прямого хода» (2.33), (2.29), (2.30), где г = 1, 2,..., N — 1 и формулы «обратного хода» (2.34), (2.25), в которых і = N — 1, N — 2,..., 1. Заметим, что т. к. ( I = 1, I ] = 1 и ( ] = 0, к 1, то при а = 1 в правой части (2.22) останется F? = А2 h2 -2±ZJL ал т УІ + — - (УІ-І + Уі+і)
Алгоритм 2.1. Обобщенный метод прогонки численного решения прямой задачи для уравнения диффузии дробного порядка по времени с постоянным коэффициентом 1) Задаем краевые условия (2.17), (2.18) и полагаем j = 0; 2) Определяем начальные значения коэффициентов ii\ и щ согласно формулам (2.33); 3) Вычисляем коэффициенты щ и гц по формулам прямого хода (2.29), (2.30); 4) Вычисляем значения у3+ по формулам обратного хода (2.34), (2.25); 5) Увеличиваем индекс j на единицу и переходим к гиагу 2.
Обобщенный метод прогонки численного решения краевой задачи для разностного уравнения диффузии дробного порядка по времени с переменным коэффициентом
Отметим некоторые полезные свойства решений разностной прямой задачи (3.69)-(3.71). Обозначим через QY величину Y? из (3.30) при условии (3.70). Таким образом, о ! = (-WfH+1" = ІШ+1 к- (з-73) к=2 fc= Представим сеточную область Q, определяемую из (3.68), в виде суммы областей, содержащих только граничные узлы, Т , и только внутренние узлы, 9? , так что Q = Х + 9. Воспользовавшись техникой книги [40], сформулируем принцип максимума для уравнения (3.69) при условиях (3.70), (3.71). Пусть Р Q — точка сетки Q и у(Р) — функция, заданная на сетке Q. Если Р Є Q т0 У(Р) принимает значения [л{Р) согласно (3.70), (3.71). Пусть Р Є 9 — внутренний узел сетки, а узел Q Є ЇЇ8(Р), где 2її(Р) Є 9 — окрестность узла Р, т. е. множество внутренних узлов сетки, не содержащее узла Р. Покажем, что уравнение (3.69) представимо в канонической форме А(Р)у(Р)= Y, S(P,Q)2/(Q), РЄЩ, (3.74) QeW(P) где для любых Р Є 9 и Q Є 2В(Р) выполнены условия А(Р) 0, (P,Q) 0, (3.75) D(P) = A{P)- ]Г Q) 0. (3.76) Qe2U(P) Рассмотрим уравнение (3.69) в виде (3.32), взяв Y? в форме (3.73). Заметим, что Р = P(xi,tj+i), Qi = {xi,tj), Q2 = (ari_i, +i), Q3 = (xi+1,tj+1), QA = (rcf_i, ty), Q5 = (xi+1,tj), Q4+k = {xi,tj+i-k), к = 2,3,..., j. Граница 3 включает узлы (ж,-, 0), (0, tj), (TV, fy), г = 0,1,.:., N, j = 0,1, 2 ..., M. Далее, нетрудно видеть, что а при условиях (3.63) и 0 а 1 выполняется i?(P, Q) 0. Кроме этого получаем, что L (P) = — сгт 0, fc=l причем D(P) — 0, если а = 1, а в случае а 1 -О(-Р) - 0 при j — М, где М — достаточно велико (2.90). Таким образом оказывается, что разностная схема вида (3.69) допускает представление в канонической форме (3.74), а значит, для нее может быть сформулирован принцип максимума, применяемый для оценок в С решений разностных параболических уравнений. Определение 3.2. Сеточная область 1 называется связной, если для любых точек Р ,Р" Є 9 найдутся такие узлы Pi, Р2, , Рт из 9? ; что РіЄЩР ), Р2ЄЩРі), ..-, РтЄЩРт-і), Р"єЩРт), причем В(Р{, Рі+1) ф 0, і = 1, 2,..., т - 1, Р(Р , А) ф 0, В(Рт,Р") ф 0. (3.77) Суть принципа максимума раскрывает следующая Теорема 3.4. Пусть у(Р) ф const — сеточная функция, определенная на связной сетке Q, и пусть выполнены условия (3.75) и (3.77). Тогда из условия L[y(P)] = A(P)y(P)- J2 B(P,Q)y(QJ 0 ( 0) QeW{P) на 9 следует, что у{Р) не может принимать наибольшего полооїси-тельного (наименьшего отрицательного) значения во внутренних узлах Ретй.
Доказательство в техническом отношении опирается на результаты [40]. Пусть L [у{Р)] 0 во всех внутренних узлах Р Є 9. Предположим, что у(Р) принимает наибольшее положительное значение во внутреннем узле fe% у(Р ) = тахт/(Р) = М0 0. Покажем, что существует внутренняя точка Р, в которой Я у(Р\ 0, что противоречит условию &[у(Р)] 0. Так как у(Р ) y{Q) для всех Q Є ЩР ), то &ЫР )] = D{P )y{P )+ ]Г B(P\Q)[y(P )-y(Q)] ЯєЩР1) D(P )y(P ) 0, 106 так как D(P ) 0 и у(Р ) 0. Заметим также, что [у(Р )] = 0 только если D(P ) = 0 и y(Q) = у(Р ) для всех Q Є ЩР ). Теперь нужно повторить рассуждения для узла Рі Є 2її(Р ), в котором у{Рі) у{Р ) — 3Vto- Так как у(Р) ф const на 9 и сетка связная, то существует такая последовательность узлов Pi, Pi,..., Pm, Р, для которых выполнены условия (3.77), что у(Рт) = у(Р ) — Мо, а у(Р") Мо, Р"є2ЇЇ(Рт). Тогда Л[у(Рт)} D(Pm)y(Pm) + B(PmiP")[y(Pm)-y(P")] В(Рт,Р»)[у(Р )-у(Р")] 0, т. е. Р = Рт. Первое утверждение теоремы доказано, а второе сводится к первому, если заменить у(Р) на —г/(Р). Следующая теорема позволяет получить оценку решения разностной прямой задачи (3.69)-(3.71). Теорема 3.5. Решение разностной прямой задачи (3.69)-(3.71) допускает оценку 1М1с 1М1ст, (3-78) з где Ы\с = тах \уЛ с, IMIcv = тах (р \, j = 0,1, 2,..., М. з Доказательство. Воспользуемся представлением (3.74), (3.75) и (3.77) краевой задачи (3.69)-(3.71), обозначив Г , г = 0, j-0,l,2,...,M; Q(P)=l 0, i = N, j = 0,1,2,...,М- (3.79) [ 0, j = 0, г = 0,1,2,..., . Заметим, что у(Р) 0, Р Є Q в силу принадлежности г/(Р) классу ІЇГ. Рассмотрим функцию у(Р) Є #, удовлетворяющую (3.74), (3.75) и (3.77), такую, что у = \\g\\c = max \g(P)\ и
Применение алгоритма поиска по шаблону для улучшения приближенного решения
Другой подход к реализации формулы (3.91), основная идея которого подробно изложена в [19], состоит в аппроксимации гессиана по секущим. Как правило, на практике много усилий приходится направлять на достижение выигрыша в эффективности, численной устойчивости и глобальной сходимости, обусловленного симметричностью и положительной определенностью матрицы Н = У2Ф(х)- Бройденом, Флетчером, Голдфарбом и Шанно была предложена следующая формула пересчета аппроксимаций гессиана, позволяющая сохранить его необходимые характеристики: А(УФ) А(УФ)Т Н АХ Ахт Я A(W)TAx АхНАх1 где А (УФ) = УФ(х+) — УФ(хс) а ньютоновский шаг Ах может быть получен путем решения линейной системы уравнений ЯАХ = -УФ(х). (3.97) Обратимся к вопросу о сходимости метода секущих. Основной результат принадлежит Бройдену, Дэннису и Морэ.
Теорема 3.10. Пусть 7] : Rn —»- R дважды непрерывно дифференцируема в открытом выпуклом мнооюестве G С Rn, и пусть У2т? Є ІлрДС). Предположим, что найдется х Є G, такое, что У 77 (х ) — О и У2т? (х ) не вырождена и полооїсительно определена. Тогда существуют положительные постоянные є и 5, такие, что если \х Х 2 — є и \HQ — У2ї](х )о ії г е матрица Н симметрична, то полооїсительно определенный метод секущих (3.96) корректно определен, а последовательность {х } остается в G и сходится сверхлииейно к х , т- е. выполняется x fc+1 — X ck \х Х \ где {(} — некоторая последовательность, сходящаяся к .
Доказательство представлено в [19]. 119 Алгоритм 3.3. Метод секущих, модифицированный для минимизации функционала невязки Ф(х) 1) Задаем вектор данных обратной задачи ф, устанавливаем счетчик итераций п = 0 и задаем нулевое приближение х 2) Вычисляем начальное значение функционала невязки ф(); а) решаем прямую задачу относительно пулевого приближения, получаем вектор данных C(x )i б) формируем вектор невязки i](x ) — Ф (х ) в) вычисляем s\ = Ф(х ); 3) Вычисляем нулевое приближение для гессиана: а) вычисляем якобиан J = УФ(х ); б) полагаем Н = JTJ + jE, где 7 — константа, которую можно выбрать аналогично алгоритму 3.2; 4) Строим градиент УФ(х): для каоїсдой j-u компоненты вектора х а) вычисляем Xj + 5х г е &Х 0 некоторое малое приращение аргумента; б) решаем прямую задачу относительно Xj + $Х получаем вектор данных {XJ + Sx); в) формируем вектор невязки r](xj + b x) — Ф (Xj + $x)i г) вычисляем Ф(ХІ + х) л) Y7 Tff \ (Xj + 5х) - 51 а) заполняем вектор Vw(Xj) — —т ; 5) Вычисляем шаг Ах — —#_1УФ(х); 6) Вычисляем Д(УФ): а) решаем прямую задачу относительно х + Ах; б) формируем вектор невязки г)(х + Ах)/ в) вычисляем S2 = Ф(х 4- А%); 120 г) строим градиент УФ(х + Ах) аналогично шагу 4 , д) полагаем А (УФ) = УФ(х + Ах) - УФ(х)/ 7) Выполняем пересчет гессиана, по формуле (3.96); 8) Полагаем s\ = S2, УФ(х) = УФ(х + Ах), X — X + Ах и переходим к шагу 5.
Здесь нужно обратить внимание на то, что при использовании метода секущих значение функционала Ф(х) не всегда монотонно уменьшается от итерации к итерации, а в зависимости от выбора 7 время от времени скачкообразно возрастает, а затем вновь стабильно убывает. Принцип построения нулевого приближения для гессиана в методе секущих целесообразно использовать такой же, как в методе Лсвенбсрга—Марквардта, с 7(о) = 10з.
Метод Флетчера—Ривса минимизации функционала невязки Данный метод использует последовательность направлений поиска, каждое из которых является линейной комбинацией антиградиента в текущей точке и предыдущего направления спуска. Он является методом первого порядка, но скорость его сходимости квадратична [27, 31]. В работе используется следующая модификация метода Флетчера—Ривса: х(А+1)=х( )+7( у ); (3.98) = -УФ (х(/с)) + "V -1), к 1; (3.99) (0) = _w х(о) . (з.ЮО) „ -!) = (v (x(fc)),v U(fc))-v (x(fc-1))) 71 (v -1)), УФ(х ))) Текущий шаг метода 7 определяется путем решения задачи одномерной минимизации: ф (x(fc) + 7( VW) = nun Ф (х(к) + 7/ (А) , (3-102) 121 для чего в работе используются метод Ньютона и метод золотого сечения. Характер сходимости метода Флетчера—Ривса устанавливает следующая
Теорема 3.11. Пусть Ф(х) сильно выпуклая функция на множестве X, т. е. существует постоянная $ О, такая, что Ф (аХ + (1 - а)х) аФ (Х) + (1 - «)Ф Ш - Ф - а)$ \Х - X? при всех х-, X Є % и всех а, О а 1. Тогда последовательность {х } построенная по методу (3.98)-(3.102), сходится к минимуму х функции Ф(х); причем имеет мест,о оценка
