Математическое моделирование и численный анализ нелинейных систем реакционно-диффузионного типа с запаздыванием тема диссертации и автореферата по ВАК РФ 05.13.18, кандидат наук Сорокин Всеволод Григорьевич

• Сорокин Всеволод Григорьевич
• кандидат наук
• 2018

Оглавление диссертации кандидат наук Сорокин Всеволод Григорьевич

Введение

Глава 1. Нелинейные модели реакционно-диффузионного типа с запаздыванием. Метод функциональных связей. Построение точных тестовых решений

1.1. Исследуемые модели. Метод функциональных связей

1.1.1. Нелинейные модели реакционно-диффузионного типа

с запаздыванием

1.1.2. Термин «точное тестовое решение»

1.1.3. РбХХТбНИЯ т ипа бегущей ВОЛНЫ

1.1.4. Метод функциональных связей

1.2. Точные решения нелинейных реакционно-диффузионных уравнений с постоянным запаздыванием

1.2.1. Кинетическая функция зависит от отношения ^/м

1.2.2. Кинетическая функция зависит от разности и —

1.2.3. Кинетическая функция зависит от линейной комбинации

и ж w

1.2.4. Кинетическая функция зависит от суммы квадратов и и w

1.3. Точные решения нелинейных реакционно-диффузионных уравнений с переменным запаздыванием общего вида т — т (г)

1.4. Точные решения нелинейных реакционно-диффузионных уравнений с постоянным запаздыванием и переменным коэффициентом переноса

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

1.6. Леммы о размножении решений для двух классов трехмерных реакционно-диффузионных уравнений с запаздыванием

1.7. Точные решения одномерных уравнений типа Фишера

с запаздыванием

1.8. Выводы к первой главе

Глава 2. Качественные особенности дифференциальных

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

2.1. Некоторые К^ЧбСТВбННЫб особенности дифференциальных

уравнений с запаздыванием

2.2. Подавление сингулярностей в решениях задач с обострением

путем введения запаздывания

2.3. Линейная устойчивость (неустойчивость) решений

реакционно-диффузионных задач с запаздыванием

2.3.1. Общая схема исследования линейной устойчивости (неустойчивости) решений

2.3.2. Случай стационарных решений. Анализ дисперсионного уравнения

2.3.3. Некоторые обобщения на случай нестационарных решений

2.3.4. Линейная неустойчивость решений для дифференциально-разностных моделей типа Каттанео — Вернотте

2.4. Нелинейная неустойчивость решений некоторых систем

реакционно-диффузионных уравнений с запаздыванием

2.4.1. Один подход для определения условий неустойчивости

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

2.4.3. Некорректность по Адамару некоторых задач с начальными данными

2.4.4. Некорректность по Адамару некоторых начально-краевых задач

2.4.5. Влияние запаздывания на неустойчивость решений

2.5. Принципы построения, выбора и использования тестовых задач

для уравнений с запаздыванием. Примеры тестовых задач

2.5.1. Предварительные замечания

2.5.2. Основные принципы построения, выбора и использования тестовых задач

2.5.3. Примеры построения тестовых задач для уравнений с квадратичными нелинейностями

2.5.4. Тестовые задачи со степенными нелинейностями. Прямое тестирование численных методов

2.6. Выводы ко второй главе

Глава 3. Численный анализ нелинейных

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

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

3.1.1. Краткое описание метода прямых

3.1.2. Интегрирование систем ОДУ с запаздыванием в среде МаЛета^са

3.2. Тестовые 3 сьд^сьч^ для ОДУ с запаздыванием. Сопоставление численных и точных решении

3.2.1. Предварительные замечания

3.2.2. Формулировки тестовых ЗсЬд^сЬч^. Точные и численные решения

3.3. Сопоставление численных и точных решении реакционно-диффузионных уравнений с запаздыванием параболического типа

3.3.1. Предварительные замечания

3.3.2. Формулировки тестовых задач для уравнений типа Фишера с запаздыванием

3.3.3. Сопоставление точных и численных решении тестовых

задач для уравнении типа Фишера с запаздыванием

3.4. Выводы к третьей главе

Глава 4. Численный анализ нелинейных

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

4.1. Метод прямых численного интегрирования обобщенных нелинейных реакционно-диффузионных уравнений с запаздыванием........ .....................ИЗ

4.2. Формулировки тестовых задач для уравнений типа Клейна — Гордона с запаздыванием

4.3. Сопоставление числ6нных и точных рбтттснии тестовых зсьд^сьч^

для уравнений типа Клейна — Гордона с запаздыванием

4.4. Общая схема численного интегрирования задач с запаздыванием методом прямых

4.5. Выводы к четвертой главе

Общие выводы и заключение

Список литературы

Введение диссертации (часть автореферата) на тему «Математическое моделирование и численный анализ нелинейных систем реакционно-диффузионного типа с запаздыванием»

Введение

Предварительные замечания

Математические модели с запаздыванием используются в динамике популяций, биологии, теории массо- и теплопереноса, биохимии, биомедицине, механике, гидродинамике, физике, химии, теории управления, математической теории искусственных нейронных сетей, экологии, экономике и других областях (Р. Беллман, С.А. Кащенко, В.Б. Колмановский, А.Д. Мышкис, А.Д. Полянин, Л.Э. Эльсгольц, J. Hale, J. Huang, G.E. Hutchinson, Y. Kuang, M.С. Mackey, A.J. Nicholson, H.L. Smith, K. Wang, J. Wu и др.) и позволяют учитывать предыдущую эволюцию процесса, либо его отдельные состояния в конкретные моменты в прошлом. Наиболее простые пространственно однородные модели с запаздыванием описываются обыкновенными дифференциальными уравнениями (ОДУ). Анализ и решение ОДУ с запаздыванием по сложности сопоставимы с анализом и решением уравнений в частных производных без запаздывания. Добавление диффузионного члена в модели с ОДУ дает возможность учесть РcLCi 1 рб^ДбЛбНИб ЧcLCТPIТобъектов или субстанции в пространстве и позволяет описывать более сложные явления или процессы.

Запаздывание может возникать ввиду различных причин. Например, в биологии и биомеханике запаздывание связано с ограниченной скоростью пере-д^ачи нервных и мышечных реакции в живых тканях, в з^дед^ицине в зад^ачах о распространении инфекционных заболеваний — время запаздывания определяется инкубационным периодом (промежутком времени от момента заражения до первых признаков проявления болезни); в динамике популяций запаздывание возникает из-за того, что особи вступают в репродуктивный период не сразу после рождения, а лишь по достижении определенного возраста; в теории управления запаздывание обычно связано с ограниченной скоростью распространения сигнала и ограниченной скоростью технологических процессов.

Многие методы численного интегрирования уравнений в частных производных без запаздывания допускают обобщения и модификации, после чего уже могут применяться для решения более сложных начально-краевых задач с запаздыванием (Д.А. Брацун, В.Г. Пименов, С.T. Baker, A. Bellen, Z. Jackiewicz, C.V. Pao, Q. Zhang, В. Zubik-Kowal и др.). В литературе наибольшее распространение получили два класса методов: конечно-разностные методы и метод

прямых (the method of lines). Метод прямых основан на аппроксимации рассматриваемого уравнения в частных производных с запаздыванием системой ОДУ с запаздыванием, которая может быть жесткой. Полученную систему ОДУ с запаздыванием можно решить с помощью различных специализированных численных методов, однако, по-видимому, в настоящее время наиболее перспективным и экономичным путем здесь является привлечение численных методов, встроенных в последние версии пакетов вычислительных программ Mathematica и Maple (или MATLAB), которые пока не справляются с решением уравнений в частных производных с запаздыванием, но достаточно хорошо интегрируют такие системы ОДУ.

Численный анализ нелинейных задач реакционно-диффузионного типа с запаздыванием (и более сложных моделей с запаздыванием) сопряжен с рядом специфических трудностей, к которым можно отнести возможную негладкость и неустойчивость решений, необходимость хранения большого массива данных и др. Важно отметить, что до сих не проводилось сопоставление результатов применения численных методов решения подобных задач с точными тестовыми решениями. Указанное обстоятельство связано с тем, что до середины 2012 года было известно всего лишь два уравнения реакционно-диффузионного типа с запаздыванием, допускающих невырожденные инвариантные решения (только одно из них выражалось через элементарные функции), отличные от решений типа бегущей волны (S.V. Meleshko, 2008). Поэтому большой теоретический и практический интерес представляет собой разработка и развитие методов построения точных тестовых решений и тестовых задач нелинейных уравнений в честных производных с зепезд^ь1ВсЪнис1\/1 j сшэлиз качественных особенностей таких уравнений и разработка эффективных численных методов их интегрирования, а также использование тестовых задач для оценки точности и области применимости численных решений уравнений с запаздыванием.

Обзор существующих математических моделей с запаздыванием

Динамика популяций. В динамике популяций время запаздывания т обычно отвечает за период взросления, во время которого особи не способны к размножению и могут быть неподвижны. Например, у многих насекомых это стадии яйца, личинки и куколки. Рост популяции в момент времени t может зависеть как от текущей плотности популяции u = u(x,t) (из-за проявления

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

Рассмотрим, например, задачу для диффузионного логистического уравнения с запаздыванием

П = апхх + Ьп(1 — и) + ](х, г), 0 < х < I, г > 0, (1)

с граничными и начальными условиями

п(0,г) = д1(г), п(1,г) = д2(г), г> 0; п(х,г) = ф(х,г), 0 ^ х ^ I, —т ^ г ^ 0.

Здесь моделируется зависимость собственной скорости роста популяции п от

п(х, г)

насыщения, зависящими, в свою очередь, от плотности популяции в прошлом и(х,г) = п(х,г — т). Коэффициент диффузии а во многих случаях является достаточно малым и отвечает за эффективность перехода из состояния с высокой плотностью особей в состояние с низкой плотностью. Отметим, что под словом «особь» понимается не только некое живое существо, но и, например, клетка крови.

ат

фузионное логистическое уравнение с запаздыванием имеет ограниченные и неограниченные решения. В [3] показано, что при достаточно большом врет

запаздыванием становятся неустойчивыми и что существует последователь-

т

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

п = Пхх + п/(и).

Для численного решения задач, описываемых диффузионным логистическим уравнением с запаздыванием, в [5] использовалось псевдоспектральное пространственное разложение Чебышева, в [6] — псевдоспектральное пространственное разложение Лежандра.

Замечание. В литературе встречаются различные названия уравнения (1): диффузионное логистическое уравнение с запаздыванием, уравнение Фишера с запаздыванием, диффузионное уравнение Хатчинсона. Достаточно часто

диффузионное уравнение (1) называют просто уравнением Хатчинсона (см, например [1, 2, 5]). Однако, исторически термин «уравнение Хатчинсона» закрепился за обыкновенным дифференциальным уравнением с запаздыванием с нелинейным слагаемым вида и(1 — т), которое рассматривалось Хатчинсоном в [7] (см. также [8-10]). Поэтому в данной работе, чтобы не вносить путаницу, используется термин «диффузионное логистическое уравнение с запаздыванием».

В [11-13] изучаются решения типа бегущей волны для уравнений типа диффузионного логистического уравнения с запаздыванием. Там также рассматривается уравнение Кобаяши

иг = ихх + Ц1 — и),

которое моделирует зависимость скорости роста популяции от плотности популяции в момент £ — т, а эффектов насыщения от тбкущбй плотности популяции. Впервые это уравнение исследовалось в работе [14]. Явный и неявный методы Эйлера для численного решения уравнения Кобаяши, в том числе условия устойчивости, описаны в [15].

Диффузионное логистическое уравнение с запаздыванием и уравнение Кобаяши являются обобщениями, с одной стороны, классического уравнения Фишера, или уравнения Колмогорова-Петровского-Пискунова (КПП), без запаздывания

иг = ихх + и(1 — и) и, с другой стороны, уравнения Хатчинсона

и'г = и(1 — т).

Существование решений типа бегущей волны для кооперативной модели Лотки-Вольтерры с четырьмя запаздываниями изучается в работе [16]. Эта модель описывается уравнениями

ди1(х,£) д2и1(х,1) / чг-, т / ч / м

-д£-= а1—дХ2--+ Г1и1(х,г)[1 — Ь1Щ (х,Ь — п) + С1щ(х,£ — Т2)],

ди2(х,£) д2и2(х,1) л,/ м

-д£- = а2-дХ2--+ Т2и2(х,г)[1 + С2Щ(х,г — Тз) — Ь2и2(х,г — ТА)],

где гг, Ьг, сг, т^ > 0, % = 1, 2, ] = 1, 2,3, 4, и при с1 = с2 = 0 распадается на два независимых диффузионных логистических уравнения с запаздыванием.

Модель обобщает уравнение Фишера без запаздывания и уравнение Хатчинсона на случай кооперативного взаимодействия нескольких видов особей.

Диффузионная модель хищник-жертва Лотки-Вольтерры с двумя запаздываниями имеет вид

дп\(х,г) д2п\(х,г)

—дг— = а^—дх2--+ пп1(х,г)[1 — Ь\п1(х,г) — с\щ(х,г — п)\,

дп2(х,г) д2п2(х,г) / м- -< / ч 7 / м -дг- = а2-дх2--+ Г2П2(х,г)[ — 1 + 02Щ(х,г — т2) — 02Щ(х,г)\,

где Ь\2,т\ ^ 0 г1,2,с1,2,т2 > 0 п\ и п2 — плотности популяций жертвы и хищника, соответственно. В отсутствие хищника, при п2(х, г) = 0, плотность популяции жертвы подчиняется уравнению Фишера-КПП без константа Ь\ является строго положительной и отвечает за самоограничение роста численности популяции. В отсутствие жертвы плотность популяции хищника также подчиняется уравнению Фишера-КПП, причем численность хищников уменьшается. Вопросам устойчивости решений этой системы посвящены работы [17, 18]. Уравнение вида

1 — Ьи ,

п = апхх + гп—-, г,Ь,с> 0,

1 + си

описывает рост популяции большой дафнии. Существование решений типа бегущей волны этого уравнения изучается в [19]. Это уравнение является следствием обыкновенного дифференциального уравнения без

3 йп 9)зд ы в йн и я

1 — Ьп

п = гп-,

г 1 + сп

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

сь НО) стадии насыщения толькхз но жизнеобеспечение. Удельная скорость роста популяции, пропорциональная отношению п[/п., лип

пп

скорость роста становится практически постоянной. Указанное качественное поведение удельной скорости роста популяции хорошо согласуется с результатами экспериментов [20].

Еще одной распространенной моделью динамики популяций является диффузионное уравнение Николсона

с запаздыванием для описания поведения

и

популяции мясных мух

и = ихх — Ьи + рт)е—ч™,

где Ь > 0 смертность взрослых особей за один день в пересчете на одну особь, р > 0 максимальное количество яиц, производимое за один день, в пересчете на одну особь, 1/д > 0 размер популяции, при котором она воспроизводится с максимальной скоростью. Это уравнение изучается, например, в [21-27]. Соответствующее обыкновенное дифференциальное уравнение с запаздыванием было впервые ^в^ы^в^зд-^^зно .в работе [28] с опорой на эксперименты Николсона [29,30].

Биология и биомедицина. Процесс кроветворения, гомопоэз, описывается реакционно-диффузионным уравнением с запаздыванием [31-33]

ввпт

щ = ихх — 6и + —-, т = и(хЛ — т),

вп + к ;

где и = и(х, £) — плотность зрелых стволовых клеток в потоке крови, т — время запаздывания между началом производства стволовых клеток в костном мозге и выходом зрелых клеток в кровь. Более сложная модель

ввпи 2ввпт

щ = ихх — ои —---+ --е 1

вп + ип вп + тп

содержит экспоненциальное слагаемое, которое отвечает за поток клеток из костного мозга в кровь и зависит от плотности клеток в момент времени £ — т.

Эти уравнения являются обобщениями предложенной в [34] модели кроветворения Мэкки-Гласса на основе обыкновенного дифференциального уравнения с запаздыванием. Введение в уравнение диффузионного слагаемого позволяет учесть перемещение клеток в пространстве из области с высокой концентрацией в область с низкой. Физический смысл и диапазоны изменения констант 6, в, в, 7 > 0 и п > 1 обсуждаются в [35]. В работе [31] исследуются осцилляционные свойства решений и устойчивость равновесных состояний п одоб н ых уравнений, в [32] изучаются решения типа бегущей ВОЛНЫ, в зз рассматривается локальный метод Галеркина для численного решения соответствующих ЗсЬД^сЬЧ^ •

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

щ = X — Ьи(х, £) — /(и(х, £),т(х, £),у(х, £))у(х, £),

Щ = /(п(х, г — т\),и(х, г — т\),у(х, г — т\))у(х, г — т\)е—11Т1 — си(х, г), уг = аАу + ки(х, г — т2)е—12Т2 — ту(х, г).

Здесь п = п(х,г), и = и(х,г), У = у(х,г) плотности неинфицирован-ных клеток, инфицированных клеток и свободных вирусов, соответственно^ Л — скорость мобилизации неинфицированных клеток, к — скорость воспроизводства вирусов инфицированными клетками, Ь, с, т ~ скорости гибели неинфицированных клеток, инфицированных а

— коэффициент диффузии, А — оператор Лапласа. Время запаздывания т\ это время с момента заражения клетки до начала производства вирионов (вирусных частиц). Множитель е—11Т1 отвечает за вероятность выживания клетки за время от г — т\ до г, где 7\ — скорость гибели инфицированных, но еще не

т2

та образования вириона до момента, когда он станет инфекционной частицей. Вероятность выживания вириона определяется множителем е—12Т2., а средняя продолжительность его жизни равна 1/^2- Функция заболевав мости / (п,и,у) является непрерывно-дифференцируемой по своим аргументам и удовлетворяет трем гипотезам:

(!) ](0, и, у) = 0 при и,у ^ 0, (11) ¡и(п,и,у) > 0 при п> 0, и, у ^ 0, (Ш) (п,и,у) ^ 0, ¡у(п,и,у) ^ 0 при п,и,у ^ 0.

Начальные и граничные условия имеют вид

п(х,г) = ф\(х,г), и(х,г) = ф2(х,г), у(х,г) = ф3(х,г),

х Е П, —тт ^ г ^ 0, тт = т.ах{т\,т2}; ду

-7^ = 0, х е дп, г> 0, дп

где П — ограниченная область в К3 с гладкой границей дП, д/дп — производная по направлению внешней нормали к границе дП, ^¡(х,г) ^ 0 заданные функции (I = 1, 2,3).

Эта модель обобщает многие диффузионные модели динамики протекания вирусной инфекции гепатита В, которые можно встретить в литературе. При т\ = т2 = 0 /(п,и,у) = вп где в > 0 некоторая постоянная инфекционного процесса, получаем диффузионную модель без запаздывания, которая впервые

была рассмотрена в [37]. В этой модели предполагается, что рассматриваемые клетки печени являются неподвижными, а вирусы свободно передвигаются по печени согласно закону диффузии Фика. При ^^ = 0,т2 = 0 I(и,т,у) = ви получаем предложенную в [38] диффузионную модель с запаздыванием, которая является обобщением модели [37]. Запаздывание отвечает за учет времени между заражением клетки и началом производства вирионов. При т2 = 0, I(и,т,у) = ви/(и + т) получаем диффузионную модель из статьи [39]. Выбор функции I(и,т,у) = ви/(1 + 61и + 62у + 63иу)} где 6122з3 ^ 0 некоторые константы, приводит к модели из работы [40], частными случаями которой являются модели из статей [41,42].

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

В [43] (см. также похожую модель в [44]) рассматривается реакционно-диффузионная модель эпидемии с запаздыванием

л п / л/-, / \\ Яоит

щ = а1Аи + иК, (и + т)(1 — (и + т)) — ии — —--,

ил \ ии

. я0ит , ч _ / ч

= а2Ат + --- — т, и = и(х,£ — т ), гл = — т ),

и + IV

где и = и(х,£), т = т(х,£) — плотности популяции восприимчивых и заразных особей, х = (х1,х2) € П, О некоторая ограниченная связная область в К2 с гладкой границей

дО А двумерный оператор Лапласа, а1, а2

т

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

я,V=Яо = в

6 + ш' 6 + 6 +

Ьо

ш

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

п(х, г) = ф1(х, г) ^ 0, и(х, г) = ф2(х, г) ^ 0, х е п, —т ^ г ^ 0;

дп ди

— = — = 0, х е дп, г> 0,

дп дп

д/дп дп

Рассматриваемая модель является следствием системы двух обыкновенных ьных уравнений без запаздывания (см. [43, 44]), одно из которых основано на логистическом уравнении, а значит, рост популяции восприимчивых особей подчиняется логистическому закону. При Яа > 1 происходит рост популяции, при Яа ^ 1 популяция погибает. Динамическое распространение болезни осущвст вляет благодаря взапмодепствпю восприимчивых и заразных особей, которые передвигаются случайным образом. Широкое распространение инфекции происходит при Я0 > 1 и не происходит при Я0 < 1. С луча й Я = 1 является граничным, при котором процесс может пойти или в ту, или в другую сторону.

Экология и биология (модель пищевой цепи). В работе [45] предложена система уравнений с запаздыванием параболического типа для описания простой пищевой цепи, состоящей из п + 1 вида живых организмов — зоопланктона г7 фитопланктона пп и мпкроорганизмов щ, г = 1,... ,п — 1а также из п видов растворенных органических и неорганических веществ , ] = 1,... ,п и детрита d (останков всех рассматриваемых живых организмов).

Пусть микроорганизм щ, г = 1,... ,п — 1, питается веществом фитопланктон пп питается вещеетвом ип, зоопланктон 2 питается фитопланктоном пп. Вещество — это концентрация растворенных органических веществ, образующихся в результате частичного распада мертвых организмов d, а также в результате жизнедеятельности фитопланктона пп и зоопланктона Вещество , ] = 2,..., п, является продуктом метаболизма микроорганизма п^Моделирование п + 1 уровня живых организмов, детрита и п веществ ведется в терминах содержания в них азота. Будем предполагать, что поток вещества от уровня к уровню изменяется согласно гипотезе Лотки-Вольтерры как произведение взаимодействующих компонент, что приводит к нелинейности системы дифференциальных уравнений. Таким образом, поведение простой пищевой

цепи можно описать системой реакционно-диффузионных уравнений с запаздыванием [45]:

dt aui + Ui(x,t Tui)Ui(wi(x,t Tui))

- Ui(x,t - Tei)Ei (wi(x,t - Tei)) - Ui(x,t)Mi(wi(x,t)), i = 1,...,n - 1;

dt aun + Un(xi t Tun)Un(wn(x, t Tun))

- Un(x, t - Ten)En(Wn(x, t - Ten)) - Un(x, t)Mn(Wn(x, t))-Z (x, t Tuz )Uz {Un(x, t Tuz ));

3w\ d 2w\ ,п, ,

= + Z(x,t - Tez )Ez (Un(x,t - Tez )) +

+ Un(x,t - Ten)En(wn(x,t - Ten)) + Kid(x,t - Tkl)-

- Ui(x,t - Tui)Ui(wi(x,t - Tui));

dwj д 2wj ,п , ,

= + U-i(x,t - Tej-i)EJ-i{wj-i{x,t - Tej-i))-

- Uj(x,t - Tuj)Uj(wj(x,t - Tuj)), j = 2,... ,n; dz д 2 z

dt = aidx2 + z(x,t - Tuz)Uz(Un(x,t - Tuz))-

- Z (x, t - Tez )Ez (Un(x,t - Tez )) - z(x,t)Mz (Un(x,t));

dd d2d n dt a2 dx2

= a2+ UiMi(wi) + zMz(Un) - Kid(x,t - Tki);

i=i

0 < x < 1, t> 0, с граничными условиями второго рода

dz dd

= 0

dUi dwi dz dd

dx x=0,i dx x=oi dx x=0,i dx

x=0i

и начальными условиями

Ui(x,t) = фi(x) ^ 0, i = 1,...,n, wj(x,t) = <n+j(x) ^ 0, j = 1,...,n, z(x,t) = <P2n+i(x) ^ 0,

d(x,t) = <2n+2(x) ^ 0, 0 < x < 1, -Tm < t < 0, Tm = m^{r„},

где Ui(x,t) wj(x,t) z(x,t) d(x,t) — концентрации перерабатываемой материн в микроорганизмах, доступных питательных веществах, зоопланктоне и детрите, соответственно; Tui и Tuz — времена запаздывания при потреблении веществ

г-м организмом (г = 1,..., п) и зоопланктоном, те^ и тег — времена запаздывания при выделении веществ г-м организмом (г = 1,..., п) и зоопланктоном, ть\ время запаздывания при распаде детрита, К — некоторый коэффициент. Функции и и иг определяют скорости потребления веществ г-м организмом (г = 1,... , п) и зоопланктоном, функции Е¡и Ег — скорости выделения веществ г-м организмом (г = 1,..., п) и зоопланктоном, функции Ы^ и Ыг — скорости смертности г-го организма (г = 1,... ,п) и зоопланктона. В работе [45] исследуется устойчивость решений данной задачи, приводятся примеры численного моделирования.

Физика (теплопроводность, диффузия и гидродинамика). Дифференциальные уравнения в частных производных параболического типа, применяемые для описания многих нестационарных процессов, обладают существенным недостатком — бесконечной скоростью распространения возмущений, что не наблюдается на практике. Это приводит к необходимости разработки более адекватных моделей, которые дают конечную скорость распространения возмущений. Ниже приведены примеры подобных дифференциально-разностными моделей для описания процессов переноса тепла и массы, а также движения вязкой несжимаемой жидкости.

В одномерном случае дифференциально-разностная модель теплопроводности имеет вид (см. [46,47] и ссылки в них)

Я\г+т = —Лпх, (2)

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

т

ально-разностное уравнение теплопроводности

уг = апхх, у = п(х,г + т),

где а = Л/(рср) — коэффициент температуропроводности. В [47,48] показано ^ что рбшбния ЗОДОЧИ 1а( шш для дифференциально-разностного уравнения теплопроводности, как правило, неустойчивы (за исключением некоторых спет

Дифференциальная модель теплопроводности Каттанео — Вернотте (см. [46-48] и ссылки в них)

Я + тяг = -Хпх (3)

приводит к уравнению теплопроводности гиперболического типа

тпи + п = апхх,

которое также дает конечную скорость распространения тепловой волны. Формально разложив левую часть дифференциально-разностной модели (2) в ряд Тейлора и оставив два первых члена, можно получить модель Каттанео — Вернотте (3). Однако, этот часто встречающийся в литературе прием не является математически корректным: член тqг при т ^ 0 играет значительную роль лишь при малых временах £ ~ т, однако при £ ~ т разложение зависимости (2) в ряд Тейлора невозможно. В [48] также показано, что свойства решений дифференциально-разностного и гиперболического уравнений теплопроводности различаются, так как различаются наборы собственных значений уравнений.

В [46] описан ряд^ точных решении нелинейных дифференциально-разностных уравнений теплопроводности вида

= сИу [/(п)Уп\ + д(у), V = п(х,£ + т).

Аналогичные дифференциально-разностные и гиперболические модели и

п

а

Дифференциально-разностный аналог уравнений Навье — Стокса для описания движения вязкой несжимаемой жидкости имеет вид (см. [49] и ссылки в ней)

+ (V • = -Ур + VДи, V- и = 0,

где и = и(х, £) — вектор скорости жидкости, V = и(х,£ + т), т — время релаксации, р — отношение давления к плотности жидкости, V — кинематическая вязкость, V — оператор градиента, Д — оператор Лапаса. Отметим, что левая

часть уравнения движения вычисляется со сдвигом по времени, при г + т. Граничные условия ставятся так же, как и для классических уравнений Навье — Стокса, а начальные условия должны быть заданы на отрезке 0 ^ г ^ т.

Химия. Химическую реакцию Белоусова — Жаботинского (см. [11,50]) можно описать системой

п = пхх + п(1 — п — гий), ий = и(х,г — т), = ихх — Ьпи,

V <Лу <Лу }

где г,Ь > 0 — константы, функции п = п(х,г) и и = и(х,г) отвечают за

т

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

п(—ж,г)=0, и(—ж,г) = 1, п(+ж,г) = 1, и(+ж,г) = 0.

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

Теория управления. Процесс термической обработки металлических листов можно описать уравнением [51-54]

п = апхх + д (иг)пх + с(/(и2) — п), и>г(х,г) = п(х,г — т¡), 0 <х< 1, г > 0,

с граничными условиями

п(0,г) = п(1,г) = 0, г ^ 0,

и начальными условиями

п(х,г) = ф(х,г), 0 <х< 1, —тт ^ г ^ 0, тт = шах{т1,т2}.

Здесь т\2 > 0 — времена запаздывания; п(х, г) — распределение температуры в металлическом листе; д (и\) — скорость ли ста, / — функция распределенного источника.

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

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

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

Многие модели искусственных нейронных сетей основаны на обыкновенных дифференциальных уравнениях с запаздыванием (см., например, [55-63] и ссылки в них

). Однако в некоторых случаях требуется диффузион-

ного члена, который позволяет учесть движение электронов в асимметричном магнитном поле. В [64-69] изучается класс моделей искусственных нейронных сетей, который

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

о / о \ п

> д ( ощ \ ,

Сг4 Т 3 ( и3

дщ з д ( ди \ ' '

~ = кС зХк\агкдх~к) — ЬгЩ + ^ Сгзи) + ^ ^9з+(4)

дг

изг = из(х,г — Гц (г)), х е и, г > о,

где щ = щ(х, г) — функция состояния г-го нейрона сети (г = 1,...,и)7 а1к = aik(х,г, и) ^ 0 гладкие функции, моделирующие диффузионные операторы трансмиссии вдоль г-го нейрона, ^ > 0 скорость, с которой потенциал г

яния покоя, Сз — константы, характеризующие взаимодействие нейронов, Тз (из) и дз (Щ^ — функции активации '-ого нейрона, 1() — функции внешнего воздействия на г-й нейрон, (г) — времена запаздывания, и — замкнутая

Список литературы диссертационного исследования кандидат наук Сорокин Всеволод Григорьевич, 2018 год

/ / ^ -

^^ ^ » "

1000

100

10

1

0 1 2 3 г 0 0.5 1 1.5 2 г

Рис. 2.1. Точное решение (2.2) задачи без запаздывания (2.1) (сплошная линия

и численные решения (а) задачи с запаздыванием (2.3) и (б) задачи с запаздыванием (2.4) при т = 0.1 (штриховая линия) и т = 0.5 (штрих-пунктирная

линия).

Из Рис. 2.1,а и 2.1,6 видно, что в рассмотренных тестовых задачах введение запаздывания полностью подавляет в решении сингулярную особенность с обострением.

Нетрудно доказать следующее утверждение. Пусть задача Коши

и[ = $(и), г > 0; и(0) = а,

где а > 0 и $ (и) > 0 при и > а, имеет решение с обострением. Тогда решен И^З модифицированной задачи с запаздыванием

и'г = $(ш), ш = и(г — т), г > 0; и(г) = а, —т ^ г ^ 0,

где т > 0, не имеет сингулярностей при конечном t.

2.3. Линейная устойчивость (неустойчивость) решений реакционно-диффузионных задач с запаздыванием

2.3.1. Общая схема исследования линейной устойчивости

(неустойчивости) решений

Опишем общую схему исследования линейной устойчивости (или неустойчивости) решений начально-краевых задач на отрезке 0 ^ x ^ 1 для следующего класса нелинейных реакционно-диффузионных уравнений с запаздыванием !

£utt + aut = auxx + f (u,w), w = u(x,t — т) (2.5)

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

Пусть u0 = uo(x,t) — нетривиальное решение рассматриваемой задачи,

т.е.

s(uo)tt + a(uo)t = a(uo)xx + f (uo,wo), wo = uo(x,t — т) (2.6)

Для анализа линейной устойчивости/неустойчивости решения uo рассмотрим возмущенные решения вида

u = uo(x,t) + ôu(x,t), (2.7)

СЧИТйЯ^ что 6 — малый параметр и u/uo = O(1) для —т ^ t < œ. Подставим (2.7) в уравнение (2.5). Учитывая равенство (2.6) и отбрасывая члены порядка S2 и выше, получим линейное уравнение для функции u = u(x,t):

eût + ащ = auxx + fu(uo,wo )u + fw (uo,wo)w, w = u(x,t — т ). (2.8)

На границах рассматриваемой области искомая функция u должна удовлетворять однородным краевым условиям Дирихле

u(0,t) = u(1,t) = 0. (2.9)

Если для любых функций u, удовлетворяющих уравнению (2.8) и граничным условиям (2.9), при выполнении соотношения

\и\ + \щ\ < Cl при — т < t < 0, (2.10)

имеет место неравенство

\и\ + \Ut\ < C2 при t> 0, (2.11)

то решение и0 будет устойчивым (Ci и C2 — некоторые постоянные, не зависящие от выбора U). Другими словами, для устойчивоети решения и0 решения линеаризованной задачи (2.8)-(2.9) с ограниченными начальными данными при —т ^ t ^ 0 должны быть ограничены и при всех t > 0.

Замечание. Для задач, описываемых уравнениями параболического типа (2.5) при £ = 0и а = 1 в неравенствах (2.10) и (2.11), следует опустить второе слагаемое \Ut\.

Если сущвст вубт функция и, удовлетворяющая уравнению (2.8) и условиям (2.9) и (2.10), которая в некоторой точке x* (0 < x* < 1) имеет свойство lim (\и\ + \Ut\) = œ (для уравнений параболического типа следует опустить

t—>со

слагаемое то решение ио является неустойчивым.

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

2.3.2. Случай стационарных решений. Анализ дисперсионного

уравнения

Пусть константа и = ио является решением уравнения (2.5), т.е. имеет место равенство

$ (ио,ио) = 0. (2.12)

Константа и0 является также решением задачи на отрезке 0 ^ х ^ 1 для уравнения (2.5) с граничными условиями первого рода и начальными данными

и(0, г) = и(1 , г) = ио г > 0;

и(х,г) = и0, иДх,г) = 0 при — т ^ г ^ 0.

Поиск возмущенных решений вида (2.7) в данном случае приводит к линейному уравнению в частных производных с запаздыванием с постоянными коэффициентами

еии + ащ = аихх + $и(щ,ио)и + (ио,ио)ш, ш = и(х,г — т), (2.14)

точные решения которого можно искать методом разделения переменных в виде произведения функций разных аргументов и = X(x)T(t). Удовлетворяя однородным краевым условиям Дирихле (2.9), находим X(x) = sin(nnx)7 где n = 1, 2, ... Для функции T(t) из (2.14), получим обыкновенное дифференциальное уравнение с запаздыванием

eT''(t) + аГ(t) = [-an2n2 + ¡и(щ,щ)]Т(t) + fwЩщ)Т(t - т), (2.15)

где штрих обозначает производную по t решения уравнения (2.15)

ищем в виде Т(t) = e~xt. В результате получим дисперсионное (характеристическое) уравнение для определения спектрального параметра Л:

-еЛ2 + аЛ - a(nn)2 + fu(uo,uo) + fw(ио,щ)еХт = 0, (2.16)

где n = 1, 2, ... (левая часть уравнения (2.16) является квазиполиномом специального вида [76, 164]).

Таким образом показано, что дисперсионное уравнение (2.16) можно было получить путем линеаризации исходного уравнения в частных производных с запаздыванием (2.5), если искать возмущенные решения в виде

и = щ0 + öe~xt sin(nnx), n = 1, 2, ... , (2.17)

где 6 — малый параметр, Л — спектральный параметр, подлежащий определению, и0 = const.

Если действительная часть хотя бы одного корня трансцендентного уравнения (2.16) отрицательна, то решение и0 = const рассматриваемой начально-краевой задачи для уравнения (2.5) является неустойчивым (поскольку lim Щ = ж в любой точке x*, удовлетворяющей условию sin(nnx*) = 0). Для линейной устойчивости этой задачи действительные части всех корней трансцендентного уравнения (2.16) должны быть положительны.

Отметим, что некоторые результаты, касающиеся вопроса о распределении действительных частей корней трансцендентных уравнений вида (2.16), можно найти, например, в [48, 7б, 164].

Рассмотрим подробнее частный случай е = 0 а = 1,

СЧИТО)Я^ что

кинетическая функция f (и, w) не зависит явно от т. Тогда дисперсионное уравнение (2.16) можно представить в виде

Л + p + qeXT = 0, 2

p = fu(uo,uo) - a(nn) , q = fw(щ,щ),

где действительные константы р и д не зависят от т. Зафиксируем целочисленное значение п и положим Л = а + г/3, где а и в ~ действительная и мнимая части спектрального параметра, г = —1. Уравнение (2.18) сводится к троне цендентной системе уравнений для определения а и в'-

а + р + деат ^(вт) = 0, в + деат вт(вт) = 0. (2.19)

Система (2.19) не меняется при замене в на —в■ Поэтому, без ограничения общности, будем считать, что в ^ 0.

а=0

а=0

1 ±аГСсов,-Р/«. + 2™ ,

(2.20)

в = Р, т = 1 ±arccos (-p/q) + 2nm

P

\p/q\ ^ 1, m = 0, 1, 2, ... ,

где надо отбросить отрицательный корень т, соответствующий значению m = 0. Упорядочив и переобозначив значения времени запаздывания, которые определяются формулой (2.20), имеем

Т1 <Т2 < ••• <TS < Ts+1 < . . .

pq

а = а(т) и в = в(т). Поэтому при изменении т в области Ts < т < Ts+1, знак действительной части спектрального параметра не меняется. Указанное обстоятельство позволяет делать выводы относительно линейной устойчивости или неустойчивости стационарных решений рассматриваемой задачи при разных т

При \p/q \ > 1 действительная часть спектрального параметра Л не меняет знак при любых т. В этом случае, взяв 0 < т ^ 1 из (2.19) получим а = —p — q, откуда следует sign а = —sign (p + q).

Посмотрим поведение корней системы (2.19) при высокочастотных возмущениях. При n ^ то имеем p ^ —то, а из первого уравнения (2.19) получим

а > 0 (n ^ то). (2.21)

Поэтому для рассматриваемого класса задач неустойчивость решений может появляться только при возмущениях с небольшой частотой (для практических целей целесообразно исследовать случай с минимальным значением n = 1)-

Укажем простой (но достаточно грубый) критерий неустойчивости стационарного решения u0 = const. Для этого рассмотрим левую часть дисперсионного уравнения Q(X) = X + p + qeXr. Имеем Q(X) ^ —ж при Л ^ —ж. Если Q(0) = p + q > 0 то существует действительное значение X = X* < 0, для которого Q(X*) = 0. Этому значению соответствует неустойчивое возмущенное решение (2.17). Возвращаясь к исходным обозначениям (2.18), описанныи достаточный признак линейной неустойчивости можно представить в виде неравен ст в cL

fu(uo,uo) + fw(uo,uo) > an2, (2.22)

при записи которого было взято минимальное значение n = 1. Отметим, что критерий неустойчивости стационарного решения (2.22) остается справедливым также для более общего дисперсионного уравнения (2.16) при £ ^ 0 и ^ 0

(И + И > 0).

Пример. Уравнение (2.5) с квадратичной кинетической функцией

f (u, w) = b\u + b2w — b3u2 — b4uw + (b3 + b4 — b\ — b2)w2

допускает стационарное решение u0 = 1. Используя (2.22), получаем достаточный критерий неустойчивости этого решения

С\

—b\ — b2 > an .

Отметим, что для численного определения а и в систему (2.19) полезно преобразовать к следующему виду:

в

exp

sin(er )

вт COs(eT) sin(er )

PT „gos(BT )

= —e—pT q, а = — P, 2.23

Бт(вт)

где первое уравнение является изолированным трансцендентным уравнением

вв а уравнений (2.23) ДсШ в разделе z. о. .

2.3.3. Некоторые обобщения на случай нестационарных решений

Дисперсионное уравнение (2.16) можно использовать также для анализа более широкого класса начально-краевых задач для уравнения (2.5), решения которых u0 = uo(x,t) отличны от константы, но удовлетворяют асимптотическому свойству lim u0 = u* = const (поскольку неустойчивость определяется

поведением возмущающей функции U при больших временах). В этом случае для анализа линейной неустойчивости (устойчивости) решение щ = щ(x,t) в уравнениях (2.6) и (2.16) заменяется его асимптотическим значением Щ при t ^ ж.

Отметим, что функции, входящие в граничные условия

u(0,t) = ^o(t), u(1,t)= ^i(t)

начально-краевых задач, которые можно исследовать указанным способом, должны удовлетворять условиям: lim ^0(t) = lim ^i(t) = и0-

Дисперсионное уравнение (2.16) при п = 1 используется далее для определения области неустойчивости точных решений некоторых нелинейных задач реакционно-диффузионного

типа с запаздыванием в главах 3 и 4.

2.3.4. Линейная неустойчивость решений для дифференциально-разностных моделей типа Каттанео — Вернотте

Рассмотрим теперь нелинейное дифференциально-разностное уравнение

щ = awxx + f (u,w), w = u(x,t — т) (2.24)

с граничными условиями Дирихле и некоторыми начальными данными (граничные условия и начальные данные могут быть произвольными и здесь опускаются). Уравнение (2.24) получено формальной заменой слагаемого ихх в уравнении (2.5) на слагаемое wxx. В частном случае f = 0 уравнение (2.24) собой линейную дифференциально-разностную модель тепло-проводноетИ Каттанео — Вернотте, в которой введено обозначение u(x, t) = 0(x,t + т), где в — температура (некорректность нэ/чэльно краевых зад^ач д^ля этого уравнения обсуждалась в [46-48]).

Пусть константа и = и0 (щ ф 0) является решением уравнения (2.24), т.е. имеет место равенство (2.12). Как и ранее, возмущенные решения ищем в виде (2.17). В результате линеаризации уравнения (2.24) приходим к дисперсионному уравнению

Л — а(пп)2вХт + fu(uo,uo) + fw (ио,ио)вХт = 0, (2.25) где п =1, 2, ... , которое удобно представить в виде

Л + p + двХт = 0,

2 (2-26) p = fu(uo,Uo), q = fw(ио,Uo) — a(nn) ,

где действительные константы р и д не зависят от т. Положив Л = а + гв и разделяя действительную и мнимую части в (2.26), для определения а и в получим трансцендентную систему уравнений (2.19).

Исследуем поведение корней системы (2.19), где р и д определяются из (2.26), при высокочастотных возмущениях. При п ^ ж имеем д ~ —а(пп)2 ^ —ж. Исключив из уравнений (2.19) экспоненциальный член, выразима через в'-

Г1СОБ(вт) „„ч

а = в .гй — Р. 2.27

Бт^вт)

а

в

в

exp

Бт(вт)

вт COS^T) Бт(вт)

= —e—pT q. (2.28)

При q ^ —ж из уравнения (2.28) следует, что вт ^ п — 0. Используя связь (2.27) имеем а ^ —ж (при n ^ ж, что эквивалентно q ^ —ж). Отсюда с учетом представления решения в виде (2.17) следует линейная неустойчивость стационарных решении для обобщенной модели Каттанео — Вернотте (2.24) для любых кинетических функций f (u,w) (глобальная неустойчивость).

Полученные u=

u(x,t) более общих начально-краевых задач для уравнения (2.24), которые при t ^ ж стремятся к постоянной величине (u ^ u0 = const), также будут неустойчивы для любых f(u, w)

Если кинетическая функция линейна по своим аргументам, т. е. f (u, w) = c\u + c2w7 то данный результат будет точным для любых решений u = u(x,t) (в этом случае неустойчивость имеет место без линеаризации уравнения), а соответствующие начально-краевые задачи будут некорректны по Адамару.

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

2.4.1. Один подход для определения условий неустойчивости

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

Пусть u0 = Uo(x,t) — вектор-функция, являющаяся частным решением некоторой нелинейной параболической системы уравнений в частных производных с запаздыванием. Пусть удалось найти другое точное решение этой системы

в виде суммы

и = По(х, £) + у(х,Ь,8),

где у(х,1,5) — некоторая достаточно гладкая по всем аргументам функция, ограниченная всюду при конечных £ и зависящая от параметра 5, который не входит в рассматриваемую систему уравнений.

Решение и0 = и0(х,£) является неустойчивым, если функция V удовлетворяет двум условиям:

^(х,£,5)\^ 0 при 0 < £ < т, 5 ^ 0, ^(х, £, 5)\ ^ то при £ ^ то, 5 = 0 — любое.

Действительно, из условия (2.29) следует, что при 0 ^ £ ^ т в силу непрерывности V по параметру 5 для любого достаточно малого А можно выбрать такое значение 5, что выполняется неравенство \и — и0\ ^ А, а при £ ^ то величина \и — и0\ становится неограниченной.

Для гиперболической системы уравнений в частных производных с запаздыванием к условиям (2.29)

добавить условие для производных по времени \^(х,£,5)\ ^ 0 при 5 ^ 0 (для 0 ^ £ ^ т).

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

Применим описанный метод для анализа нелинейной неустойчивости систем реакционно-диффузионных уравнений гиперболического типа с запаздыванием следующего вида [119]:

£х (п\)и + а1(п1)г = а1(и1)хх + Ьщ + Г (их — kwí,u2,w2), ^

£2(и2)й + а2(и2); = а2(и2)хх + 0(щ — kwí,u2, w2), к> 0,

где их,2 = их,2(х,£), Wí,2 = их,2(ж,£ — т); ах,2 > 0 £х,2 ^ 0 ах,2 ^ 0 (£х,2 + ах,2 = 0); Г(и,-и^), С(и, V, w) — произвольные функции трех аргументов (Ги ф т — время запаздывания.

В общем случае система (2.30) (при т = 0) допускает простейшие решения: стационарное (не зависящее от £), однородное (не зависящее от х) и типа бегущей волны их = их(г), и2 = и2(г), где г = ах + в£- Устойчивость таких и некоторых других решений различных реакционно-диффузионных уравнений с запаздыванием и систем таких уравнений рассматривается, например, в [24, 165-167].

Лемма 2.1 (обобщает лемму 1.2 на системы вида (2.30)). Пусть

Що = Що(х,г), и 20 = П2о(х,г) (2.31)

есть некоторое частное решение системы (2.30). Тогда система (2.30) имеет также решение

и1 = и10(х,1) + вс1У(х,г), и2 = и20(х,г), с = — 1пк, к> 0, (2.32)

т

где У = У(х, £) — любое т-периодическое решение линейного гиперболического уравнения

еУи + (( + 2б1 с)У = а\Ухх + (Ь — £гс2 — агс)У, У(х,1) = У(х,1 — т). (2.33)

Формулы (1.23)-(1.24) определяют общий вид функции У(х,1) = У1(х,1,а1 + 2е1с, Ь — £\с2 — а1с).

Воспользуемся леммой 2.1 для получения условий неустойчивости нелинейной реакционно-диффузионной системы (2.30). Для этого возьмем стационарное пространственно-периодическое решение уравнения (2.33):

У = 5бш^х + ц), 7 = \/(Ь — е1с2 — а1с)/а1, Ь — £1с2 — а1с ^ 0, (2.34)

где 5, ц — произвольные постоянные.

Анализ соотношений (2.32) и (2.34) показывает, что для системы (2.30) условие (2.29) будет выполняться при

к> — , т > 0, Ьт2 — ег(1п к)2 — ахт 1п к ^ 0, (2.35)

и тогда любое решение и10(х,1)щ20(х,1) системы (2.30) будет неустойчивым. Условия (2.35) можно представить в более н^глядном виде •

к> — , Ь > 0, т ^ то, то =Ь2Ьк(а1 + а2 + (2.36)

Соотношения (2.36) являются следствием решения квадратного неравенства т

т ^ т*, п =

1п к

*' '* ~ пи Г^ \ и1

((I —

2Ь

т > 0

Физический смысл условий (2.36) состоит в том, что в области параметров к > 1,Ь > 0 неустойчивость возникает за счет запаздывания, которое должно быть достаточно большим: т ^ т0.

Поскольку вид кинетических функций Г и О не влияет на условия неустойчивости (2.36) решений реакционно-диффузионной системы (2.30), бу-д^ез^/г их н аз ыв ат ь глобальными условиями неустойчивости.

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

Замечание 2.2. Система

£х(их)ы + ах(щ)г = ах(щ)хх + Ьхих + díWí + гх(щ — kwí,u2,w2),

у/.о (;

£2(и2)й + 02(и2^ = а2(щ)хх + О(щ — kWí,U2, и^), может быть записана в виде (2.30), где

dí dí

Ь = Ьí + , Г (г,и2 ,W2) = Fí(z,U2,W2) — — г. (2.38)

Поэтому условия неустойчивости решений системы (2.37) получаются из (2.36) заменой параметра Ь на ^ + ^/Н).

Замечание 2.3. Полученные результаты по неустойчивости распространяются на нелинейные многокомпонентные системы вида

£\_ии + аи = а]_ихх + ^и + Г(х, £,и — кп), и2, w2,..., ит, wm), £n(un)tt + а n(un)t ап(ип )хх + О п(х, £, и — kw,u2,w2,... ,um,wm), (2.39) п = 2,... ,т, к > 0,

где и = и(х, £), w = и(х,£ — т), ип = ип(х, £), wn = ип(х,£ — тп); Г(...), О(...) —

т тп

рые могут быть разными).

Замечание 2.4. Полученные результаты по неустойчивости распространяются также на одиночные реакционно-диффузионные уравнения с запаздыванием вида

еии + ащ = аихх + Ьи + Г(х, £,и — kw). (2.40)

2.4.3. Некорректность по Адамару некоторых задач с начальными

данными

Пусть (2.31) — решение задачи типа Коши с начальными условиями общего вида

U1 = фц(х,г), dtUi = ^12(x,t), U2 = ^2i(x,t), dtU2 = ^22(x,t) При 0 < t < T

(2.41)

для нелинейной системы уравнений с запаздыванием (2.30) на всей области изменения пространственной переменной —то < x < то. Здесь и далее dt обо-

t

Из леммы 2.1 следует, что при k > 0 система (2.30) имеет также решение, которое определяется формулами (2.32) и (2.34). Обозначив это решение Ui, U2, получим

U1 = U1o + Sect sin(Yx + ¡i), U2 = U2o, (2.42)

где S и ¡i — произвольные постоянные, c = 1 ln k, а коэффициент y опреде-

t

0 ^ t ^ t имеем

|Ui — Uio| < SeCT, \dtUi — dtUio| < SceCT, |U2 — U20I = 0, \dtU2 — OtU2o\ = 0.

(2.43)

При фиксированных t и k (при k > 1, что соответствует c > 0) разности (2.43)

t

S

будут мало различаться при 0 ^ t ^ т. С другой стороны, при выполнении

условий (2.36) их = -(П — l)

7v 2

имеем

\Uí — щ0\ = 5ес ^ то при £ ^ то,

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

Указанная неустойчивость решений нелинейной системы уравнений с запаздыванием (2.30) относительно н^ч^льных данных д^ел^сьет Зсьд^сьч^у^ х\.оттти^ для системы (2.30) некорректно поставленной по Адамару (в случае выполнения условий (2.36)). Важно отметить, что эта неустойчивость носит общий харак-

ГО

2.4.4. Некорректность по Адамару некоторых начально-краевых

задач

Покажем, что при выполнении условий (2.36) может иметь место глобальная неустойчивость решений некоторых нелинейных начально-краевых задач с граничными условиями первого, второго, третьего рода в ограниченной области 0 ^ х ^ Н [119].

Пусть (2.31) — решение начально-краевой задачи для системы уравнений с запаздыванием (2.30) с начальными условиями (2.41) и общими граничными условиями первого рода:

щ(0,г) = фи(х), и1 (н,£) = фП(г), щ(0,г) = щ(н,г) = Ф22(^), (2.44)

где Н = п/^, а коэффициент 7 определен в (2.34).

Формула (2.42) при ¡1 — 0 Д-^сКЗтг решение системны (2.30), которое точно удовлетворяет граничным условиям (2.44). Это решение ЗсЬ сч^ет 5 можно сделать сколь угодно близким к решению (2.31) в области начальных данных 0 ^ £ ^ т (для доказательства используются рассуждения полностью аналогичные рассуждениям в разделе 2.4.3). Однако при выполнении глобальных условий неустойчивости (2.36) первоначально близкие решения (2.31) и (2.42) рассматриваемых начально-краевых задач будут экспоненциально расходиться при £ ^ ж в центре рассматриваемой области (при х = Н/2). Такая неустойчивость решений системы (2.30) относительно начальных данных делает начально-краевую задачу для системы (2.30) некорректно поставленной по Адамару (в случае выполнения условий (2.36)).

Замечание 2.5. В случае краевых условий второго рода, когда на границах области задаются производные (и\)х\1 (и2)х7 решение (2.31) следует сравнивать с решением, полученным с помощью леммы 2.1 и формулы (2.34) при ¡1 = п/2.

Замечание 2.6. Вопросы устойчивости решений дифференциальных уравнений с запаздыванием рассматриваются в [3,4,18,24,47,48,161,165-169].

2.4.5. Влияние запаздывания на неустойчивость решений

Для выбора ад^екватных численных зугетод^ов необходимо хорошо понимать качественные особенности задач с запаздыванием. Введение 3 1 сЬЗгЛ^Ы В Э)Н и я в стандартные математические модели (описываемые дифференциальными уравнениями без запаздывания) может превратить корректную задачу в некорректно поставленную. Поясним сказанное на конкретных примерах.

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

, ,/n — kw\ £(U — kw\ /n

entt + out = auxx + bu — bl k ) + Я 1 _ k I, k = 1, (2.45)

где b = const, a /(z) — произвольная функция. Это уравнение является частным

T = 0 k = 0

переходит в уравнение

entt + out = aUxx + f (n). (2.46)

Любой корректной задаче Коши для уравнения (2.46) можно поставить в соответствие некорректную задачу с начальными данными для уравнения (2.45),

bkT

Проиллюстрируем вышесказанное наглядным примером. Рассмотрим некоторое семейство задач с начальными данными для уравнения (2.45) при т ^ 0 и фиксированных ПОСТ оя н н ых k > 1, b > 0. При т = 0 имеем обычную задачу Коши без (корректную по Адамару при достаточно общих предположениях), решения которой устойчивы. Постепенно увеличивая т, можно прийти к некорректной задаче по Адамару, которая при т ^ to, где to = ^nr (o + л/o2 + 4eb), характеризуется глобальной неустойчивостью решений. Таким образом, обобщение математической модели путем введения запаздывания может превратить соответствующую ей хорошую (корректную) задачу в некорректно поставленную.

Примеры некорректности по Адамару

начально-краевых задач для уравнений параболического и гиперболического типов с запаздыванием можно найти в [167, 170].

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

КйЧбСТВбННЫб

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

Используемые численные схемы целесообразно тестировать на модельных задачах, допускающих точные решения в замкнутом виде (см. [171]).

2.5. Принципы построения, выбора и использования тестовых задач

для уравнений с запаздыванием. Примеры тестовых задач

2.5.1. Предварительные замечания

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

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

2.5.2. Основные принципы построения, выбора и использования

тестовых задач

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

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

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

3. Рассматриваемые тестовые зи^ дол^ж^ны содбр^ж^^ть ряд свободных параметров, которые можно свободно варьировать в широких пределах.

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

5. Для оценки точности численных методов лучше использовать несколько различных тестовых зад^ач.

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

7. Полезно тестировать численные методы на быстро растущих решениях при достаточно больших временах.

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

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

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

2.5.3. Примеры построения тестовых задач для уравнений с квадратичными нелинейностями

Описанные в разделе 1.7 точные решения уравнений (1.115) и (1.124) содержат свободные параметры и могут быть использованы в качестве основы для формулировки тестовых задач ^ х 1 рна;3ненных^ для отт^енк^и^ точ^н^ости^ численных методов интегрирования реакционно-диффузионных уравнений с запаздыванием .

Рассмотрим уравнение (1.124). Выберем его параметры так, чтобы оно, как и диффузионное логистическое уравнение с запаздыванием (1) (при / = 0), имело стационарные решения и0 = 0 и и0 = 1. Первое тривиальное стационарное решение уже имеется. Чтобы совпало второе стационарное решение в (1.125) надо положить и0 = 1. В результате получим в = 1/(1 — к). Подставив

это значение в (1.124), приходим к уравнению

1 Л и — к,ш\ , . ,

щ = аихх + Ьи1 1--1—— ), = и(х,к — т). (2.47)

При к ^ 0 уравнение (2.47) переходит в уравнение Фишера, а при к ^ — в диффузионное логистическое уравнение с запаздыванием (уравнение Фишера с запаздыванием). Ниже приведены две группы наиболее простых точных решений уравнения (2.47).

1) Решения при к > 0 (к =1):

и = всЬ[Л еоБ^х) + В Бт^ж)], 7 = \/ (Ь — с)/а при Ь > с; и = есС,(Ле—тх + Ве{х), 7 = уДс—Ь/а при Ь < с;

и = есЬ(Лх + В), при Ь = с;

и = 1 + ес [Л еоБ^х) + В вт(7х)], 7 = \/—с/а при с< 0; и = 1 + ес(Ле—1х + Ветх), 7 = л/ф, при с > 0,

где с = (1п к)/т; Л, В произвольные постоянные.

к<0

и = Лпес1тКх еоБ(впг т 1иХ + Сп), /Зп = п(2п — 1),

т

вп Л N (Ь — с)2 + в! — Ь + с )1/2

_ гп л _

^п о л , Лп

С

2аЛп' п \ 2а )

и = 1 + ЛпвСтХпх ео$>(впк Т 1пХ + Сп), вп = П(2П — 1)

т

вп л {^с^ГЩ + с )1/2

_ гп л _

2аЛп \ 2а

,

где с = (1п |к|)/^ Лп7 Сп — произвольные постоянные; п = 1, 2, ...

Для тестирования численных методов интегрирования нелинейных реакционно-диффузионных уравнений можно выбирать указанные выше точные решения уравнения (2.47). Отметим к алеет вен^н^ ые особенности некоторые решений, которые полезно использовать для сопоставления с результатами численных расчетов.

Первое и четвертое решения из первой группы (1) являются периодических

качестве тестовых решений для начально-краевых задач с граничными условиями первого и второго рода на отрезке 0 ^ х ^ шп/^ (ш = 1, 2, ...).

Подходящим выбором свободных постоянных А и В можно сделать искомую функцию на границе равной нулю или единице (для краевых условий первого рода) или получить на границе нулевую производную по х (для краевых условий второго рода). В случае задач со смешанными граничными условиями удобно рассматривать эти решения на отрезках 0 ^ х ^ 2 тп/^ (т = 1, 2, ...). Начальные данные прИ /Т ^ £ ^ 0 (или 0 ^ £ ^ т) определяются из решений, используемых в качестве тестовых задач. Полезно сравнить численные и точные решения тестовых Зсьд^сьч^ п^эр! к, близких к единице (когда решения слабо

к

ются быстро).

Решения из второй группы (11) при к = —1 являются периодическими по времени и быстроосциллирующими по обеим переменным при т ^ 0. Такие решения полезно использовать для оценки точности численных методов в зад^сьчсьх с большими градиентами.

Тестовые задачи формулируются так: выбранное уравнение и его точное решение дополняют начальными данными при —т ^ £ ^ 0 и граничными условиями при х = 0 и х = 1, которые получают из используемого точного решения. Для иллюстрации сказанного в Таблице 1 приведены уравнения и точные аналитические решения трех из семи тестовых задач, которые будут более детально описаны далее в главах 3 и 4. Начальные данные п^зи^ т ^ £ ^ 0 и граничные условия при х = 0 и х = 1 для этих задач получаются из приведенных решений и здесь опускаются. Эти и другие задачи будут использованы в главах 3 и 4 для тестирования используемых численных методов.

Таблица 1.

Некоторые тестовые задачи с запаздыванием (0 ^ х ^ 1, £ > 0)

№ Уравнение Решение

1 щ = аихх + Ьи[1 — в (и — кт)] и = вс1[со$,(пх/2) + 2^\п(пх/2)}./ к > 0, Ь = (1п к)/т + ап2/4 с = (1п к)/т

2 щ = аихх + Ьи — в (и — кт)2 есг+\ и =1+ 2 (ех в-х), с =(1п к)/т 2 — 1 к > 0 к = 1 Ь = (1п к)/т — а в = Ь/(1 — к)2

3 иы = аихх + и (и — т) и = $,т([3х/л/а) соъ(в£), в = 2п/т

Замечание 2.9. Уравнение (2.47) является обобщенным диффузионным логистическим уравнением с запаздыванием (обобщенным уравнением Фишера

с запаздыванием) и может быть записано в более удобном виде иг = аихх + Ьи[1 — (а^и + а2^)], а! + а2 = 1,

где а! = 1/(1 — к).

2.5.4. Тестовые задачи со степенными нелинейностями. Прямое

тестирование численных методов

Тестовые задачи со степенными нелинейностями. Вместо пятипа-раметрических семейств уравнений (1.115) и (1.124) для получения тестовых задач можно использовать точные решения более общих реакционно-диффузионных уравнении с запаздыванием вида

и = аихх + Ьи — в(и — Ы)п+!, (2.48)

иг = аихх + Ьи [1 — в(и — к-ш)п]т, (2.49)

к = 0 п = ш = 1 т = 0 п = ш = 1 Фишера. Уравнение (2.48) является частным случаем одномерного уравнения (1.96) при /(г) = —вгп+!, а уравнение (2.49) — частным случаем одномерного уравнения (1.109) при ](г) = к(г) = 0 9(г) = Ь(1 — вгп)т.

Точные решения уравнений (2.48) и (2.49), как и ранее, строятся с помощью лемм 1.4-1.6 и 1.8-1.10.

О прямом тестировании численных методов. Для получения тестовых задач помимо уравнений с запаздыванием вида (1.1), допускающих точные решения, можно использовать более широкий класс уравнений

еии + аиг = аихх + Г (и,,ш) + С(ж,к), = и(х,к — т). (2.50)

Выбирают (достаточно произвольно) некоторую функцию^ = п(х,к), удовлетворяющую заданным граничным условиям. Эта функция является точным решением уравнения (2.50) при

С(х,г) = пг — апхх — Г (ц,ц), п = п(х,к — т). (2.51)

Уравнение (2.50)—(2.51) вместе с соответствующими начальными и граничными условиями является тестовой задачей, которая имеет точное решение и = п(х^). Это решение сравнивают с численным решением этой задачи. Различные функции п = П(х,к) порождают различные уравнения (2.50) и раз-личны6 тестовые •

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

Замечание 2.10. При применении прямого метода в качестве функций п вполне допустимо использовать описанные в главе 1 точные решения реакционно-диффузионных уравнений с запаздыванием.

2.6. Выводы ко второй главе

1. Выявлены КйЧбСТВбННЫб особенности некоторых классов задач реакционно-диффузионного типа с запаздыванием, связанные с подавлением сингулярности и возможной неустойчивостью их решений.

2. Показано, что усложнение математических моделей путем введения за-11 йЗДЫ В 9iH И Я В 3 BjX С обострением может полностью подавить сингулярность

рвТТТвН li jiri •

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

4. Доказано, что любые стационарные решения нелинейных дифференциально-разностных моделей массопереноса типа Каттанео — Вернотте с произвольной кинетической функцией являются неустойчивыми (в линейном приближении).

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

6. Описаны принципы построения, выбора и использования тестовых задач с запаздыванием.

Приводятся примеры построения тестовых задач.

Глава 3. Численный анализ нелинейных реакционно-диффузионных систем с запаздыванием параболического типа

3.1. Метод прямых численного интегрирования реакционно-диффузионных уравнений с запаздыванием

параболического типа

3.1.1. Краткое описание метода прямых

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

щ = [р(х,г,и,и)их]х + я(х,г,и,и)их + /(х,г,и,и), г> 0, 0 ^ х ^ 1; и(х,г) = ф(х,г), -т ^ г ^ 0; (3.1)

и(0,г) = ф0(г), п(1,г) = ^(г), г> о,

где и = и(х,г — т). Второе слагаемое в правой части уравнения соответствует возможной конвективной (движущейся) составляющей модели; в частности, при р(х,г,и,и) = 1, д(х,г,и,и) = —и уравнение (3.1) является уравнением Бюргерса с нелинейным источником и запаздыванием.

Введем пространственную сетку хт = тН7 где т = 0,1,...,М, Н = 1/М — шаг сетки. Аппроксимируя производные по х разностными аналогами

хт

(итX = ^х[Рт^хит\ + Ят^т + /т, Ш =1,... ,М — 1, 0 <г < Т; ит(г) = ф(хт,г), т = 0,1,...,М, —т ^ г ^ 0; (3.2)

ио(г) = фо(г), им (г) = ^(г), 0 <г ^ т.

ЗдеСЬ ит ит(г) u(xm, г) рт р(хт, г, ит1 итдт 0_(хт1 г, ит1 /т

/(хт, г, ит, ит), ит = и(хт, г — т), Т — временной интервал вычислений, 5х — разностный оператор, который определяется так:

^т = НЩт.- ит,.

^х[рт^хит\ н2 [рт(ит+1 ит) рт—1(ит ит—1)]-

Система (3.2) содержит М — 1 неизвестных функций ит(г) и столько же уравнений, а также две известные функции и0(г), им (г). Отметим, что полученная

система ОДУ часто является жесткой и ее решение должно вестись соответствующими численными методами.

Замечание 3.1. Функция wm(t) является известной и обозначает функцию um(t — т), которая была вычислена на несколько временных слоев ранее. Численное интегрирование уравнений с запаздыванием, в отличие от уравнений без запаздывания, требует хранения данных со всех временных слоев из промежутка от tn — т до tn—i, где tn — расчетный временной слой, что приводит к существенным затратам оперативной памяти [162].

3.1.2. Интегрирование систем ОДУ с запаздыванием в среде

Mathematica

Основным методом численного интегрирования ОДУ и систем ОДУ, в том числе с запаздыванием, в среде Mathematica является команда (встроенная функция) NDSolve [172,173]. На текущий момент с помощью команды NDSolve можно решать только системы с постоянным запаздыванием (в том числе с несколькими запаздываниями) [94]. Без дополнительных опций команда NDSolve решает систему комплексным методом, при котором в процессе вычисления происходит автоматическая смена методов решения системы ОДУ и выбор значений параметров метода. В ранних версиях Mathematica такой подход не удавалось применять для решения жестких задач, поскольку программа требовала конкретизировать метод решения, однако в более поздних версиях (например, Mathematica 11.2.0) комплексный метод уже давал некоторые результаты. Тем не менее, в диссертации мы будем использовать команду NDSolve, самостоятельно выбирая определенный метод решения системы ОДУ с помощью опции Method, которая позволяет использовать один из встроенных методов решения жестких систем ОДУ: метод из класса неявных методов Рунге — Кутты [174, 175] или неявный многошаговый метод Гира, основанный на формуле дифференцирования назад (BDF — Backward differentiation formula [176]). Комплексный метод будет использоваться только в разделе 3.2 при тестировании задач Коши для ОДУ с запаздыванием.

Приведем краткое описание применяемых методов. Для начала запишем систему (3.2) в общем виде в векторной форме:

u't = F (t, u, w), w = u(x,t — т), 0 <t < T;

u(t) = p(t), —т < t < 0 ^

где u = u(t) = (u0(t), • • • ,uM(t)), p(t) = (^(xo,t),... , ф(хм,t)), компоненты вектора F определяются так:

Fo =(фо)'г + uo(t) - фо(t),

Fm &x\Pm^xum\ + qm$xum + fm i ^^ • • • i -M -1,

Fm = (Ф1); + um (t) - ^i(t).

Метод Гири встроен в Mathematica как часть пакета IDA, входящего в библиотеку методов SUNDIALS, которая разрабатывается Ливерморской национальной лабораторией им. Э. Лоуренса, США (IDA — Implicit Differential-Algebraic solver — неявный дифференциально-алгебраический решатель, SUNDIALS — SUite of Nonlinear and DIfferential/ALgebraic equation Solvers — набор нелинейных и дифференциальных/алгебраических решателей) [176]. Программный код методов IDA (см. руководство пользователя [177]) основан на DASPK [178,179] — программах на языке Фортран, позволяющих решать дифференциально-алгебраические системы больших размерностей. Опишем работу метода Гира, руководствуясь материалами [176,177].

Метод Гира основан на аппроксимации временной производной с помощью формулы дифференцирования назад порядка r с переменным шагом, которая имеет вид

(un)'t = an,Un-i, (3.4)

—1

¿=0

где Пп = и(^)

, 8п ^п ^п—1 шаг по времен и, / 8п т , tn текущий временной слой; коэффициенты однозначно определяются порядком г и значениями предыдущих шагов вп [177].

Применяя формулу (3.4) к системе (3.3), получаем систему нелинейных алгебраических уравнений:

r

G(un) = F(tn, Un, Wn) - s-1 ^ OLn,Wn-i = 0, (3.5)

n n n

i=0

которая на каждом слое 1п решается итерационным методом Ньютона

Ыиг] — и(п >) = —си >).

Здесь на каждой итерации ] решается линейная система, Лп — аппроксимация Якобиана системы (3.5), которая вычисляется по формуле

дС(ип) дЕ(ип, <Ъп) ап,о

J

n

dun dun

s

n

Функция wn _ u(tn — т) вычисляется с помощью интерполяции, если значение tn — т лежит вне точек сетки, и является известной функцией u значения которой были вычислены ранее на слое t*7 тел и tn — т _ t

На каждом временном слое метод Гира вычисляет оценку En локальной погрешности и автоматически выбирает размер шага sn и порядок r таким образом, чтобы выполнялось соотношение ||En/wn|| < 1, где m-я компонента ^n,m вектор a определяется по формуле

_ 1

Значения констант p и q определяются с помощью опций PrecisionGoal ^ p и AccuracyGoal ^ q команды NDSolve. Норма Ц^Ц по умолчанию автоматически выбирается командой NDSolve в зависимости от метода решения (но может быть задана вручную). Для метода Гира это норма вида ||x||2 _ \xi\2 [180].

Методы Рунге — Кутты основаны на использовании квадратурных формул. Интегрируя уравнения системы (3.3) по времени от tn—i до tn, получаем:

ptn

un _ un—i + F(t, u, w)dt. (3.6)

Jtn-1

Интеграл в формуле (3.6) аппроксимируем некой квадратурной формулой и

r

ты [92, 174, 175]:

r

un _ un—i + SnY, bKn, (3-7)

i=i

r

K'n _ F(tn + CiSn, un—i + Sn^2 a'iJKn, Wn), (3.8)

j=i

где sn — шаг сетки, bi — веса квадратурной формулы, ci — коэффициенты, оире-

aij

условию ci _ Y^r=i aij (i _ 1, •••, r). Функция Wn _ u(tn — т) вычисляется с