Содержание к диссертации
Введение
1 Основные положения методологии построения алгоритмов решения одного класса систем операторных уравнений 12
1.1 Основы общей методологии 13
1.2 Случай симметричного оператора А 17
2 Численное решение нестационарной системы Стокса, возмущенной кососимметрическим оператором 22
2.1 Постановка задачи для нестационарной системы Стокса и переформулировка ее как задачи оптимального управления . 22
2.2 Вариационные уравнения 26
2.3 Итерационные процессы 31
2.4 Методы уточнения приближенного решения 33
2.5 Численные эксперименты 35
3 Численное репіение системы уравнений динамики приливов в декартовых координатах 41
3.1 Постановка задачи 41
3.2 Схема расщепления 42
3.3 Решение стационарной системы 44
3.4 Задача оптимального управления 46
3.5 Итерационный процесс решения задачи 49
3.6 Результаты численных экспериментов 51
4 Численное репіение системы уравнений динамики приливов на сфере 58
4.1 Постановка задачи 58
4.2 Схема расщепления 61
4.3 Решение стационарной системы 63
4.4 Задача оптимального управления 66
4.5 Итерационные процессы 71
4.6 Некоторые свойства гладких решений 76
5 Численное исследование итерационных процессов на сфере 83
5.1 Стационарная линейная система уравнений динамики приливов с оператором До : сферический слои 84
5.2 Нестационарная линейная система уравнений динамики приливов с оператором До : сфера 90
5.3 Нестационарная линейная система уравнений динамики приливов с оператором Ді : сфера 92
5.4 Численное исследование ошибок алгоритма от замены оператора Ді на оператор До 96
5.5 Численное исследование влияния специальных условий в "полюсных точках" 98
5.6 Численные эксперименты для тестовых решений 101
5.7 Численные эксперименты с реальными данными 103
Заключение 113
Приложение. Алгоритмы решения уравнений Навье-Стокса. 115
Список литературы 124
- Случай симметричного оператора А
- Вариационные уравнения
- Решение стационарной системы
- Нестационарная линейная система уравнений динамики приливов с оператором До : сфера
Введение к работе
Большое количество физических процессов динамики океана описывается моделями, использующими различные модификации и упрощения классических задач гидродинамики, таких как система уравнений Навье-Стокса [16, 29, 31, 41, 15, 36, 30, 82].
Исследование и численное решение системы уравнений Навье-Стокса - одна из наиболее сложных задач вычислительной математики и гидродинамики, методы решения которой активно разрабатывались в течение последних 50 лет. Количество работ, опубликованных на эту тему, не поддается исчислению, отметим лишь малую часть монографий [24, 43, 37, 4, 45, 70, 71]. Большинство этих работ посвящено развитию численных методов решения уравнений Навье-Стокса и различных их модификаций. Трудности при решении этих задач связаны с недостаточной информацией о точных решениях (почти все найденные точные решения не несут в себе специфики нелинейной задачи). Дополнительные сложности возникают при учете реальных физических данных (геометрия области, разброс значений коэффициентов, специфика поведения решений, размерность задачи, расчет на долгий интервал по времени, ограниченные ресурсы ЭВМ и многое другое). Существующая, в то же время, большая востребованность решения данных задач при моделировании физических процессов гидродинамики оставляет актуальным вопрос о разработке эффективных методов решения в каждом конкретном случае.
Одним из подходов конструирования новых алгоритмов решения задач математической физики (в том числе и задач гидродинамики) является методология их построения, базирующаяся на методах теории оптимального управления. Вероятно, впервые эти подходы были предложены в работе [52] в применении к решению классической стационарной системы Стокса. Идея
построения таких методов при рассмотрении системы Стокса состоит в следующем: функция давления рассматривается в качестве "дополнительной" неизвестной по отношению к компонентам вектора скорости - "основным" компонентам решения, а уравнение неразрывности рассматривается в качестве "дополнительного уравнения" и относится к уравнению замыкания задачи. Затем задача рассматривается как обратная задача (или задача оптимального управления) н включается в семейство задач оптимального управления, зависящих от регуляризирующего члена. Далее исследуются задачи оптимального управления рі для их решения применяются классические численные методы. В последующем распространение данных подходов к построению численных алгоритмов и их формулировке для систем операторных уравнений было выполнено в [1].
В настоящей работе проведено исследование данной методологии в применении к нестационарной системе Стокса, возмущенной ограниченным ко-сосимметрическим оператором и системе уравнений динамики приливов, а также предложены направления развития этой методологии для решения нелинейных уравнений Навье-Стокса.
Приведем описание рассматриваемых в данной работе задач и обзор методов их решения.
Нестационарная система уравнений Стокса является линеаризацией уравнений Навье-Стокса. С физической стороны она описывает движения вязкой несжимаемой жидкости при малых числах Рейнольдса (т.е. при малых скоростях, большой вязкости или малых размеров рассматриваемой области), и в этом случае решение уравнений Стокса хорошо приближают решения соответствующих уравнений Навье-Стокса. Кроме того, она часто возникает в качестве одного из этапов решения уравнений Навье-Стокса, например, при дискретизации по времени конвективного слагаемого с предыдущего слоя.
Помимо упомянутой выше методологии для решения задачи Стокса существует большое количество методов. Приведем обзор этих методов по разным направлениям, цитируя работу [20]. Это алгоритмы типа Узавы-сопряженных градиентов [56, 73, 43, 65, 70, 77, 63, 33, 81], методы релаксации [56, 65, 43, 48, 21, 62], методы, основанные на идеях симметризации и предобусловливания [46, 47, 5, 61, 66, 84, 85], многосеточные и многоуровне-
вые алгоритмы [49, 60, 66, 87, 88, 89, 83], расщепление граничных условий для давления [34, 35], итерирование граничных условий для скорости [64], метод фиктивных областей [6, 58], использование гармонической составляющей скорости [86]. Рассматриваемый в диссертации алгоритм решения нестационарной задачи Стокса является новым по отношению к перечисленным. Из особенностей алгоритма можно отметить отсутствие жесткого требования удовлетворения условию divu = 0 последовательности приближенных решений, а также требования симметричности оператора задачи (в большинстве из перечисленных выше алгоритмов по существу используется свойство симметричности) .
Одной из возможных постановок для задачи динамической теории приливов является рассматриваемая в работе "гиперболо-параболическая" система уравнений [82, 30], которая относится к классу систем уравнений теории "мелкой воды" [13] и является упрощением уравнений Навье-Стокса. Отметим некоторые особенности этой системы. Она содержит "слабую" нелинейность, кососимметричный и "диффузионный" (диссипативный) операторы. При устремлении коэффициента диффузии к нулю и постановке соответствующих граничных условий система принимает "чисто" гиперболический вид. Специальная подобная система (без нелинейного слагаемого) исследовалась Соболевым С.Л. в работе [42], которая явилась основополагающей для нескольких направлений в исследованиях задач математической физики. Также системы подобного вида исследовались в работах [68, 76]. В отличие от [42], в настоящей работе рассматриваются услови5і "прилипания" на границе, а коэффициент диффузии считается положительным. Ряд результатов по существованию и единственности решения подобных задач приведен, например, в [30].
Для численного решения системы уравнений модели динамики приливов можно использовать HN-метод (сокращение от английского Hydrodynamic Numerical Method) [72], метод дробных шагов [79], метод конечных элементов [78] и др. Если функция рельефа дна Н является константой и отсутствует нелинейное слагаемое, то задача принимает вид классической задачи газовой динамики, описывающей движения вязкого слабосжимаемого газа. В этом случае можно предположить, что большое число способов исследова-
ния и решения такого рода задач (см. [40, 18]), разработанных в классической теории газодинамики, могут найти здесь применение. Однако нас будет интересовать случай сферической системы координат и реальной топографии дна. Здесь по существу возникают сложности, связанные с: а) аппроксимацией области, соответствующей акватории Мирового океана, и поведением решений возле границы этой области; б) реальной функцией рельефа дна, являющейся "сильно негладкой"; в) нелинейностью некоторых коэффициентов системы и специфическим их поведением в различных частях области, в которой решается задача, а также другими особенностями системы уравнений динамики приливов и её решений. Кроме того, можно предположить, что в России на сегодняшний день отсутствует национальная гидродинамическая глобальная модель приливов.
Это побудило автора к исследованию и разработке новых алгоритмов решения задачи динамики приливов с использованием исследуемой в данной работе методологии, уделяя особое внимание поведению гладкого решения в окрестности "полюсных точек" и проблеме аппроксимации исходных уравнений на равномерных сетках в сферической системе координат.
Цель работы - разработка и исследование новых алгоритмов решения класса задач геофизической гидродинамики (возмущенных уравнений Сток-са, задач динамической теории приливов в декартовых и сферических координатах) на основе методологии их построения, базирующейся на подходах теории оптимального управления и сопряженных уравнений.
Диссертация состоит из введения, пяти глав, заключения, приложения и списка литературы. Далее кратко рассмотрим содержание работы по главам.
В первой главе рассматривается методология конструирования новых методов решения следующего класса систем операторных уравнений:
Ьф = / + Ви, C(f> = h + Du, (1)
где /, h - заданные элементы, ф,и- неизвестные, L, B,C,D - линейные операторы. Для этого вводится соответствующая система пространств. Затем задача (1) рассматривается как обратная задача (или как задача управления, тогда и следует трактовать как "управление") и включается в семейство задач оптимального управления, зависящих от регуляризирующего члена.
Далее исследуются задачи оптимального управления и для их решения применяются классические итерационные методы. Доказывается сходимость последовательности приближений к решению задачи (1).
Во второй главе проведено исследование этой методологии в применении к нестационарной системе Стокса, возмущенной ограниченным косо-симметрическим оператором, в ее классической постановке: найти функции ф : П х [0,Т] -> R" и р : П х [О, Г] -> Ж, где П ограниченная область в Жп, такие, что
фг - аАф + \хф + Ьф = f -Vp в П х (О, Т) = QT,
divtf> = 0 в QT, ф\г = О V* Є (О, Г), 0Uo = Фо, (2)
Jpdn = 0\/t, а
где ф = (фі,..., фп), 1x0= (—1ф2,1фъ 0,..., 0), а, 6, I = const, а > 0, b ^ 0. Для этого мы переформулируем классическую постановку нестационарной системы Стокса как задачу оптимального управления, в которой давление р(х, t) выступает в роли "управления", а уравнение div(/> = 0 рассматривается в качестве "уравнения замыкания" и входит в функционал стоимости. Предложено несколько итерационных процессов решения поставленной задачи, которые имеют следующее свойство: в них отсутствует ограничение div = 0 для сеточных функций на каждой итерации. Реализация этих процессов основана на последовательном решении задач параболического вида. Доказана сходимость последовательности приближений к решению задачи (2). В конце главы приводятся результаты численных экспериментов.
В третьей главе рассматривается применение общей методологии совместно с классическими методами расщепления к построению численного метода решения системы уравнений динамики приливов в декартовых координатах и области прямоугольной формы. А именно, в Qt = Пх (0,Т), Г2 = [0, А] х [0, В] требуется найти функции U = (и, и) Є {L>2(Qt))2,C Є L2(QT) такие, что и, v Є L2(0,T; (1 (tt))2), ut,vt,(t Є L2(0,T; (Wj (Q))2y и удовлетворяющие следующей системе уравнений
U, - Div(z/VU) + KU +
G + divU = 0, (3)
UlcK2 = 0> U|t=0 = U(0), Clt=o = C(0),
где г/, Н > 0 в fi, Н,д = const, divU == g + g,
Div(z/VU) =
,K
*|U| -^ І ft|U|
/, г* = cons/ > О,
U = (її, її) - заданная вектор-функция, компоненты которой являются достаточно гладкими. Компоненты вектор-функции U имеют физический смысл составляющих полного потока, - отклонение поверхности океана относительно "поверхности равновесия", Н - функция рельефа дна, г* и v - коэффициенты придонного и горизонтального турбулентного трения, I - параметр Кориолиса.
Сначала для задачи (3) применяется схема расщепления, а затем для реализации одного из этапов схемы расщепления формулируется задача оптимального управления, доказывается её корректная разрешимость и приводится итерационный процесс минимизации функционала "стоимости". Получены оценки границ спектра оператора стационарной задачи и на их основе оценки скорости сходимости итерационных процессов. В заключении приводятся результаты численных экспериментов.
Четвертая глава посвящена применению исследуемой методологии совместно с классическими методами расщепления к построению численного метода решения системы уравнений динамики приливов на сфере.
Постановка задачи рассматривается в сферических координатах #, А, г, X Є [0, 27г], 9 Є [0,7г], где г полагается равным среднему радиусу Земли R, и выписана в терминах "средних" скоростей и симметризованном виде. Область fi соответствует поверхности Мирового океана на "сферической Земле" или некоторой её части. В Qt е fix (0,Г), Т < оо, требуется найти функции U == (щ и) Є (Z/2(Qt))2>C Є 1/2(Qr) такие, что u,vX являются достаточно гладкими на Qt функциями и удовлетворяют следующей системе уравнений
HVt - z/A/?U + KV + gHWC = f,
Ct + divJEnj = /3, (4)
uUn = o» uLo = u(o)3 CLo = C(o)-
Здесь
ApV = (Aw, Av) + /3DV, /З Є [0,1],
,. „ _, ( 1 дф \дф\ ,. тт 1 ди ldtvsine)
A = divV, V* = [j^^, Ш) , d1VU ее ^-^- + д^^р
1 f dv ди\
Ше , -^ + 2cosfl—,-?;-2cos6>— ,
i?2sin^ V oA ал/
1Ш= (r*|U|u-Zv,ZM + r*|U|v), -
її H Є W^Q), I = —2a; cos (9, cj - угловая скорость вращения Земли, остальные функции и коэффициенты были определены ранее. В приводимых выше обозначениях круглые скобки ( , ) определяют вектор-функцию.
Отдельной проблемой, которая обсуждается в данной главе, становится постановка задачи в сферической системе координат, в силу наличия сингулярности в полюсных точках. Как и в предыдущей главе, сначала для задачи (4) применяется схема расщепления, а затем для реализации одного из этапов схемы расщепления формулируется задача оптимального управления, доказывается её корректная разрешимость и приводится итерационный процесс минимизации функционала "стоимости". Получены оценки границ спектра "оператора давления" стационарной задачи и на их основе оценки скорости сходимости итерационных процессов. Исследована специфика поведения гладких решений задачи динамики приливов в окрестности "полюсных точек", получены специальные формулы уточнения численного решения в этих окрестностях.
Пятая глава посвящена численным аспектам рассмотренных в четвертой главе алгоритмов. Численно исследуются: а) скорость сходимости итерационных процессов применительно к линейной системе уравнений динамики приливов на сфере; б) ошибке при замене оператора А і в исходных уравнениях динамики приливов на сфере на оператор До, широко используемой в геофизической гидродинамике; в) аппроксимация уравнений в окрестности "полюсных точек"; г) решение нелинейной нестационарной системы уравнений динамики приливов на сфере. В конце главы приводятся расчёты в акватории Мирового океана с реальной функцией рельефа дна и приливным потенциалом.
В заключении кратко сформулированы основные результаты диссертационной работы.
В приложении изложены некоторые возможные направления для развития исследуемой методологии применительно к построению алгоритмов решения системы уравнений Навье-Стокса без детального изучения рассматриваемых алгоритмов и их численной проверки.
Апробация работы. Результаты диссертации обсуждались на семинарах Института вычислительной математики РАН, Вычислительного Центра РАН, кафедры Вычислительной математики Механико-математического факультета МГУ им. Ломоносова.
Результаты диссертации докладывались на Международной конференции "Тихонов и современная математика" (Москва, 2006), Международной конференции "Дифференциальные уравнения и смежные вопросы" (Москва, 2007), конференции "Тихоновские чтения" (Москва, 2007), Международной конференции "Ломоносов-2008" (Москва, 2008), Международной конференции EGU General Assembly (Вена, 2008).
Публикации. По теме диссертации опубликовано 7 печатных работ [2, 3, 7, 53, 54, 55, 59]. Из них 4 статьи в рецензируемых журналах. Во всех работах, подготовленных в соавторстве, диссертант совместно осуществлял теоретическое исследование задач, самостоятельно проводил обоснования итерационных процессов, планирование численных экспериментов, выполнил все расчеты и осуществлял анализ результатов.
Автор искренне благодарен профессору В.И. Агошкову за постановку задач, научное руководство и всестороннюю поддержку, профессору В.П. Шу-тяеву за обсуждение и полезные советы, а также Г.М. Кобелькову, В.Б. Залесному, Ю.В. Василевскому, Е.И. Пармузину, В.Б. Сухову, В.М. Ипатовоп и многим сотрудникам ИВМ РАН и преподавателям кафедры Вычислительной математики МГУ им. Ломоносова за постоянную возможность консультаций и готовность к обсуждениям.
Случай симметричного оператора А
При решении прикладных задач, относящихся к исследуемому нами классу, в ряде конкретных случаев удается переформулировать исходную постановку задачи таким образом, что оператор А становится симметричным (или самосопряженным при соответствующем выборе исходных пространств). В этом случае итерационный алгоритм (1.12) решения системы операторных уравнений (1.1) можно упростить за счет исключения этапа решения сопряженных задач.
Итак, пусть операторы L,B,C,D такие, что С = В , В = С , D — D 0, L — L 0 и существует L l (из чего следует, что Ker{L} = Ker{L } — 0 и (L-1) = (L )-1). Тогда оператор Л является положительным и симметричным. Кроме того, по-прежнему считаем, что оператор L l ограничен и справедлива оценка (1.2). Пусть также для упрощения изложения Хс = Не.
Здесь мы допускаем случай, когда задача (1.18) (или, что тоже самое, система (1.1)) либо некорректно поставлена (например, оператор А является вполне непрерывным), либо наличие свойства корректной разрешимости для этой задачи заранее неизвестно. Рассмотрим семейство регуляризовапных задач Лаи = аи + Аи = д, (1.19) a = const 0. Задача (1.18) включена в семейство задач (1.19) в качестве предельного случая при а = 0. Регуляризованное уравнение (1.19) представляет собой метод регуляризации М.М. Лаврентьева [23]. Оператор Аа является положительно определенным, сформулируем для решения задачи (1.19) некоторые классические итерационные методы. Они будут относится к двум классам алгоритмов: двухслойным и трехслойным.
Трехслойные градиентные методы. Свойства оператора А в уравнении (1.19) позволяют использовать трехслойные градиентные методы. Доказательство сходимости и(а), ф(сг) ки,ф- решению системы (1.1), при а — +0 можно провести по аналогии с соответствующими рассуждениями из [1]. Замечание 2. Если задача (1.18) является корректно поставленной, то во всех рассуждениях выше можно полагать а — 0. Однако если нижняя граница спектра Со близка к нулю, то, как легко заметить скорость сходимости итерационных процессов уменьшается и параметр регуляризации начинает играть более существенную роль в скорости сходимости итерационных процессов. Отметим, что в каждом случае применения изучаемой в данной главе методологии к конкретным задачам математической физики принципиальными являются следующие вопросы: исследование существование решения при а = 0и сходимости решений (Фа, Ua) ПРИ а — +0 к решению задачи (1.1); нахождение границ спектра оператора А, что позволяет получить условия сходимости итерационных методов и оценить их скорость сходимости; исследование плотной разрешимости задачи (1.1) для ускорения сходимости итерационных методов.
Изучение этих вопросов составляет большую часть обоснования всего алгоритма решения задач вида (1.1), а получаемые при этом утверждения представляют собой отдельные математические результаты. Именно таким исследованиям в применении к конкретным задачам математической физики посвящена данная диссертационная работа. Глава 2
Численное решение нестационарной системы Стокса, возмущенной кососимметрическим оператором В данной главе исследуется приложение общей методологии, изложенной выше, к построению итерационных алгоритмов для численного решения нестационарной системы Стокса, возмущенной кососимметрическим оператором. Для этого мы переформулируем классическую постановку нестационарной системы Стокса как задачу оптимального управления, в которой давление выступает в роли "управления", а уравнение div / = 0 рассматривается в качестве "данных наблюдений" и входит в функционал стоимости. Предложено несколько итерационных процессов решения поставленной задачи, которые имеют свойство: в них отсутствует ограничение div0 = 0 для сеточных функций на каждой итерации. Реализация этих процессов основана на последовательном решении задач параболического вида. Доказана сходимость последовательности приближений к решению исходной задачи. В конце главы приводятся результаты численных экспериментов.
Вариационные уравнения
На каждом шаге итерационного процесса (2.30)-(2.34) для решения систем параболических уравнений можно использовать, например, метод Фурье (или другие подходящие методы). Именно этот метод, возможно не являющийся экономичным, был выбран для проведения численных экспериментов, поскольку мы здесь не стремились использовать наиболее эффективные методы, так как нашей главной целью было исследование итерационного процесса (1.12) и его характеристик. При этом нередко возникает возможность достаточно точно аппроксимировать элемент ф весьма небольшим числом некоторых базисных функций {фк}- однако может оказаться, что в результате итерационного процесса в конечномерных подпространствах точности вычисления приближения к р будет недостаточно. Для получения более точного численного приближения к функции р при малом числе базисных функций {фк} можно предложить различные способы.
Численные результаты этого способа решения иллюстрируют рис. 2.3 и табл. 2.3, 2.4. Другой способ уточнения приближений к функции р можно осуществить путем решения дополнительного уравнения для р после окончания итерационного процесса. Изложим идею такого рода уточнения в применении к рассмотренному выше алгоритму решения двумерной нестационарной системы Стокса в прямоугольнике Q = (О, Л) х (0,5). Так, если решения u,v,p достаточно гладкие, то в каждый фиксированный момент времени t функция р удовлетворяет уравнению Ар — divf — ІтоЬф в О, др. д2и др. __ d2v у ЛА)
Таким образом, после окончания итерационного процесса с целью уточнения приближений к функции давления можно решить серию уравнений (2.42) на каждом временном шаге. Численные результаты этого способа решения иллюстрируют рис. 2.4 и табл. 2.5.
Заметим, что алгоритм уточнения путем решения вспомогательных задач (2.42) может быть выписан в терминах обобщенной постановки задачи (2.1). 2.5 Численные эксперименты. Главная цель численных экспериментов заключалась в подтверждении сходимости итерационных процессов (2.30)-(2.32) и в исследовании их скорости сходимости. Кроме того, нас интересовал вопрос об уточнении приближения к функции давления и применение подходов описанных в предыдущем разделе.
Все численные расчеты проводились для О, С М2, Q = (0, А) х (0, В). При численной реализации итерационных процессов для решения системы параболических уравнений на каждом шаге использовался метод собственных функций, который для сетки 50 х 50 узлов по пространству, 200 - по временному интервалу на гладких решениях давал результаты относительной ошибки « 0.001. На рис. 2.1, 2.2 и в табл. 2.1, 2.2 приведены результаты численных расчетов решения задачи (2.8) (пример 1, 2) без предварительного преобразования правой части и поправок функции давления. Эти результаты показывают, что итерационный процесс (2.30)-(2.32) сходится. На рис. 2.3 и в табл. 2.3, 2.4 приведено сравнение численного решения задачи (2.8) (пример 1, 3) с предварительным преобразованием правой части и без него. Из табл. 2.3 видно, что для функции давления относительная ошибка в норме L уменьшилась в два раза, точность определения функции скорости в целом не изменилась.
Результаты этих расчетов позволяют высказать следующее предположение: если мы знаем, что при численном решении задачи (2.1) функции скорости определяются более точно (как, например, в случае реализации итерационного процесса (2.30)-(2.32)), то становится целесообразным дополнительно решить уравнение (2.42) для уточнения приближения к функции давления. Как видно из таблиц, функция давления приближается менее точно, чем функции скоростей, поэтому для демонстрации различных зависимостей скорости сходимости итерационных процессов от входных параметров приведены графики значений относительной ошибки в норме 1/2 именно функции давления. На рис. 2.5 приведены графики численного расчета примера 1, демонстрирующие увеличение скорости сходимости итерационного процесса при выборе оптимальных параметров т . В табл. 2.6 приведены результаты расчета примера 2, демонстрирующие зависимость сходимости итерационного процесса от параметра I. Из табл. 2.6 видно, что при увлечении параметра I растет относительная ошибка, и для ее уменьшение требуется большее количества итераций. Все расчеты (за исключением расчета, результат которого приведен на рис. 2.5) проводились с выбором параметров 7 по формуле (2.35).
Сформулируем основные выводы, к которым приводят теоретические и экспериментальные результаты по исследованию итерационных алгоритмов (2.30)-(2.34). Результаты численных экспериментов подтверждают сходимость итерационных процессов (2.30)-(2.34). При этом выбор параметров Tk по формуле (2.35) существенно увеличивает скорость сходимости итерационных процессов. Кроме того, такой выбор снимает вопрос о величине выбираемых параметров т/-, которые зависят от спектра оператора в уравнении (2.17) ((2.25), (2.23)).
Решение стационарной системы
Перейдем к исследованию следующей стационарной системы, полученной после применения схемы расщепления -aAU + Ш + cVC = F, иш = 0 (3.5) cdivU + d( = G. Данная система уравнений при b — d — 0, с = 1, G = О представляет собой классическую постановку для задачи Стокса. Методы численного решения подобных задач и возникающие при этом трудности рассматривались многими исследователями (см. [43, 71, 48] и цитируемую там литературу). В данной главе исследуется новый итерационный алгоритм решения (3.5), построенный с помощью методологии, изложенной в главе 1.
Ясно, что если функции U, С являются решением (3.6), то они будут решениями и задачи (3.10). Согласно изложенной в главе 1 общей методологии, для того, чтобы гарантировать существование и единственность решение задачи (3.10) мы можем включить в функционал стоимости (3.10) регуляризи-рующий член — С2- Однако позже будет доказано, что задача (3.10) всегда имеет единственное решение (корректно поставлена). И это решение, в свою очередь, при а — 0 является решением (3.5). В этом случае можно полагать а = 0. Выпишем необходимые условия экстремума для функционала J(U, С) SJ = (cdivU + d - G, cdiv U + d5) = 0. (3.11) При этом зависимая вариация 6TJ связана с независимой вариацией 5С, посредством уравнения -aAJU + bSU = -cV8(. Условие (3.11) можно записать в терминах билинейных форм a[U, U] + 6(U, U)H - (С, cdivU) = (F, U)H VUG W, a[U , U ] + b(U , U )H = (cdivU + dt;-G, cdivU ) VU є W, (3.12) (cdivU , C) + d(cdivV + dC-G,O = 0 VC Є L2(Q). При наличии соответствующей гладкости решений (3.12) можно записать в виде системы вариационных уравнений -aAU + 6U + cVC = F, U\an = 0, -aAV + b\J = cV(cdivV + dC-G), U+ = 0, (3.13) -cdivU + d(cdivU + d( - G) = 0. Для исследования свойств разрешимости задачи (3.12), а, следовательно, и эквивалентности обобщенных постановок (3.10) и (3.6) для задачи (3.5), введем операторы (LU, U)w = a[U, U] + b(U, U)H VU, U Є W, L : W - W , L (L) = W, (?,U) = с (CdivU) VC Є L2(Q),U Є W, Б : L2(0) - W , () - L2(ft), (CU, C) = с (divU, C) VU W, С Є L2(fi), С : W -+ L2(fi), (C) = W, (U, L U)W = a[U, U] + b(U, U)H VU, U є W, L : W - W , D(L ) = W, (C, B U) = с (С, divU) VC Є L2(fl), U Є W, Б : W - L2(fi), (# ) = W, (U,C C) = c(divU,C) VU Є W,C Є L2(fi), C : L2(H) -. W , L (C ) - Ь2(П). Отмстим, что введенные операторы являются ограниченными. Кроме того L = L , (L-1) = (L )"1, В = С ,С = В .
Отметим также, что если в (3.16) принять: і cdivU(fc) + )-q2 1к 2 - cdivU ( ) + d(cdivU( + d k) - GQ2(0/ K } то получаем метод градиентного спуска процесса минимизации J(, U) с "оптимизированным" выбором релаксационных параметров {7А} [1] Реализация итерационного процесса (3.16) сводится к последовательному решению задачи Дирихле для "обычного" двумерного эллиптического уравнения второго порядка. Для решения данной задачи разработаны различные эффективные численные методы [39, 11, 9]. В приводимых ниже численных экспериментах мы использовали известный метод конечных разностей. Мы предполагаем, что шаги сеток при аппроксимации этих эллиптических задач достаточно малы, чтобы возникающими при этом численными ошибками можно было пренебречь. Таким образом, мы считаем, что основную роль играют ошибки итерационного процесса (3.16). По очевидным соображениям изложение метода конечных разностей в применении к данным задачам мы здесь не приводим (см, например, [39]).
Сделаем замечания о некоторых предельных случаях поведения коэффициентов в (3.5) и соответствующей им скорости сходимости итерационного процесса (3.16). Если вычислять {7А} согласно (3.17), то из приведенного выше выражения для q можно получить некоторые асимптотические значения этой величины, которые могут быть полезными при анализе сходимости (3.16). Так, если предположить, что At — +0 (т.е. Ъ —» +оо и d — +00), то приближенно можно принять: о = — и сделать предположение о аа-I- & быстрой сходимости рассматриваемых итерационных процессов типа (3.16). Этот асимптотический случай указывает на принципиальное отличие скорости сходимости изучаемых в данной работе алгоритмов в применении к задачам (3.1), (3.2)-(3.3) по сравнению с теми же алгоритмами, но в применении к задаче Стокса (см. [2, 52]).
Все численные эксперименты мы приведём для Q С М2, Г2 = (0,1) х (0,1) с использованием метода конечных разностей для реализации шагов схемы расщепления (3.2)-(3.3). Для этого введем сетку О, = {(ж , у3) ж,- = ihx, і — 0,...,NX, yj = jhy, j = 0,..., Ny}, hx = 1/NX, hy = 1/Ny. Неизвестными будем считать значения компонентов вектора скорости и(х{, yj), v(xi, yj) при і — 1,..., АГа; — 1, j = 1,..., Ny — 1 и значения уровня поверхности океана С(ЖІ, yj) при і — 0,..., Nx, j = 0,..., Ny. Все дифференциальные операторы в итерационном процессе (3.16) будем аппроксимировать конечными разио стями со вторым порядком точности. Для решения получаемых при этом симметричных систем алгебраических уравнений использовался метод сопряженных градиентов.
Главной целью численных экспериментов было исследование сходимости итерационного процесса (3.16). Предварительные численные эксперименты проводились при выборе итерационных параметров { Ук} п0 формулам (3.17) и (3.19). При вычислении параметров по формуле (3.17) были получены результаты, подтверждающие утверждение теоремы 3.1. Однако при этом было установлено, что выбор параметров {7fe} согласно (3.19) на порядок увеличивал скорость сходимости итерационного процесса (3.16) по сравнению с выбором этих параметров по формуле (3.17). Поэтому во всех приведенных ниже результатах численных экспериментов {jk} выбирались согласно (3.19). Критерием окончания итерационного процесса (3.16) являлось достижение функционала (3.10) значения, меньшего Ю-4. В численных экспериментах в качестве Uj-i в формуле для оператора Kj-i/2 = -K"(Uj_i) из (3.3) мы брали U , т. е. функцию U с предыдущего слоя. Тогда оператор Kj-i/2 описывает влияние сил придонного трения и имеет более практически значимый интерес.
Нестационарная линейная система уравнений динамики приливов с оператором До : сфера
Теперь мы приведём результаты численных экспериментов для линейной нестационарной системы (4.3) при [3 = О, Г2 = [0, 27г] х [0,7г] и отрезка по времени [0,Т], Т оо с использованием схемы Кранка-Николсон и метода конечных разностей для реализации итерационного процесса (4.17). Введем сетку по времени tn = nAt, п = 0,1,...,NT, At = T/NT, по простран ству будем использовать сетки Q,f\ Qh/2, определённые выше. Неизвестными будем считать значения вектора скорости U = (u(Ai,9j:tn), v(Aj,9j,tn)) в точках сетки Q,h и значения уровня поверхности океана С (ЛІ, 9J, tn) в точках сетки при п = 1,... ,NT- Все дифференциальные операторы в итера ционном процессе (4.17) будем аппроксимировать конечными разностями со вторым порядком точности. Для решения получаемых при этом систем алгебраических уравнений с разреженной матрицей, как и в предыдущем параграфе, использовалось LU разложение и реализующий его программный пакет SuperLU [10, 74].
Во всех приведенных ниже результатах численных экспериментов параметры {7/с} для итерационного процесса (4.17) выбирались согласно (4.21). На каждом шаге по времени tn = nAt начальное приближение для итерационного процесса (4.17) полагалось равным = (Л, #,n_i) и выполнялось фиксированное количество итераций - 10.
Более подробный анализ результатов показывает, что наибольшие ошибки присутствуют в окрестности точек 9 = 0 и 9 = 7Г. Поэтому, в следующем тестовом примере продемонстрирован один из простейших приемов сглаживания решения около "полюсных" точек, позволяющий уменьшить численные ошибки. Тест N7. Пусть Т = 86400, At = 30. Замыкание систем разностных уравнений осуществлялось с помощью соотношений (5.2)-(5.3), но теперь после реализации одной итерации алгоритма (4.17) решение при в = в\ принималось равным этим значениям в полюсах (т.е., фактически применяются одни из простейших способов "сглаживания" численного решения). В табл. 5.7 приводятся относительные ошибки на момент времени t — Т.
В данном параграфе мы приведём результаты численных экспериментов для линейной нестационарной системы (4.3) при (3 = 1, ГІ = [0, 27г] х [0,7г] и от резка по времени [0, Т], Т со с использованием полностью неявной схемы и метода конечных разностей для реализации итерационного процесса (4.22) (4.24). Как и прежде будем считать неизвестными значения вектора скорости U = (u(\j, 6j, tn), у(\г, 9j, tn)) в точках сетки QJ1 и значения уровня поверхно сти океана (Аг, 6j, tn} в точках сетки при п = 1,..., NT- Все дифферен циальные операторы в итерационном процессе (4.22)-(4.24) будем аппрокси мировать конечными разностями со вторым порядком точности. Для реше ния получаемых при этом систем алгебраических уравнений с разреженной матрицей, использовалось LU разложение и реализующий его программный пакет SuperLU [10, 74].
В качестве условий замыкания разностной схемы для оператора Ді значения сеточных функций uh, vh при 9 = 0 и 9 = 7Г полагалось равным известному точному решению. Такой способ замыкание непригоден для расчётов с реальными данными, однако позволяет выявить возникающее при численном счете проблемы. Тест N8. Расчёт примера 3 проводился при v = 1012. В этом случае быстрая скорость сходимости итерационного процесса наблюдалась при At = 3600 (1 час), что соответствует полученной оценке (4.26). Значение функционала J{u Vh, C/i) »на каждом временном шаге уменьшалось за 3 итерации с 6.1е-4 до 2.4е-10. В таблице 5.8 представлены относительные ошибки на различные моменты времени, на рис. 5.3 приведены их графики в различных нормах с течением времени.
С точки зрения численной реализации решения эллиптических уравнений в задачах геофизической гидродинамики, рассматриваемых на сфере в области Мирового океана, иногда весьма выгодной считается возможность использования конечномерного оператора (AQ)-1 вместо (А )-1. В таком случае уравнения для компонентов вектора скорости становятся независимыми, что существенно уменьшает как объемы памяти для хранения матриц, так и время, затрачиваемое на их обращение. В данном параграфе мы приведем результаты численных экспериментов, демонстрирующих разницу такой замены на тестовых примерах. Пусть Н = 4000, д = 9.80665, R = 6371000, г = 0, I = 0. При расчётах количество итераций итерационного процесса на каждом временном шаге было фиксированным - 3 итерации. Шаги сетки по пространству: h\ = he = 1.
Тест N10. В данном численном эксперименте правые части задавались подстановкой точного решения из примера 3 в исходное уравнение при (3=1. А при численной реализации решения эллиптических задач, возникающих на каждой итерации итерационного процесса (4.22)-(4.24), обращался сеточный оператор AQ . В качестве условий замыкания разностной схемы для оператора AQ значения сеточных функций uh, vh при 9 = 0 и 9 = тг полагались равными известному точному решению. Коэффициент v полагался равным 1012, шаг по времени At = 3600. Значение функционала J(uh, v , ) на каждом временном шаге уменьшалось за 3 итерации с 5.0е-4 до 1.2е-8. В таблице 5.10 представлены относительные ошибки на различные моменты времени. На рис. 5.5-5.7 приведены графики точных и вычисленных решений на конечный момент времени расчёта.
