Содержание к диссертации
Введение
1 Сингулярное разложение в некорректных задачах 14
1.1 Понятие г-решения для систем линейных алгебраических уравнений 16
1.1.1 Регуляризация 17
1.2 Понятие s-чисел вполне непрерывного и ограниченного операторов 18
1.2.1 Понятие s-чисел вполне непрерывного оператора 18
1.2.2 Понятие s-чисел ограниченного оператора 22
1.3 Примеры некорректных задач и построение сингулярных чисел 26
2 Определение источника колебаний в волновом уравнении по данным о колебаниях на части границы области 28
2.1 Физическая и математическая модели 30
2.1.1 Физическая постановка 30
2.1.2 Математическая постановка прямой и обратньгх задач 31
2.1.3 Единственность решения обратной задачи 33
2.2 Численньге алгоритмьг решения прямой задачи 34
2.2.1 Численное интегрирование 35
2.2.2 Метод Фурье 35
2.2.3 Конечно-разностные схемы 37
2.3 Анализ некорректности обратной задачи термоакустики 39
2.3.1 Представление обратных задач в матричном виде 40
2.3.2 Анализ сингулярных чисел дискретных операторов обратных задач 41
2.3.3 Условная устойчивость обратной задачи 42
2.4 Оптимизационные методы численного решения обратных задач термоакустики 44
2.4.1 Выражения для градиентов функционалов, используемых в численных расчетах 46
2.4.2 Численные эксперименты решения обратных задач термоакустики градиентными методами 48
2.5 Сингулярное разложение и прием С.К. Годунова 55
2.5.1 Метод усеченного сингулярного разложения. Анализ численных расчетов 57
2.5.2 Прием С.К. Годунова. Анализ численных расчетов 58
3 Определение источника колебаний в волновом уравнении по данным о колебаниях в фиксированный момент времени 61
3.1 Постановка задачи 63
3.2 Неустойчивость и неединственность задачи Дирихле 64
3.3 Исследование обратной задачи в случае постоянных коэффициентов 66
3.3.1 Теорема единственности 67
3.3.2 Неустойчивость обратной задачи 68
3.4 Сингулярный анализ обратной задачи
3.4.1 Сингулярное разложение оператора обратной задачи 68
3.4.2 Сингулярный анализ дискретного аналога оператора обратной задачи 70
3.5 Вариационная постановка обратной задачи. Градиент целевого функционала 73
3.5.1 Формула градиента целевого функционала, используемая в численньгх
расчетах 74
3.6 Численньге экспериментьг 76
3.6.1 Метод сингулярного разложения 78
3.6.2 Метод простой итерации 80
3.6.3 Использование г-решения в качестве начального приближения для метода простой итерации 84
4 Определение источника колебаний в волновом уравнении по данным о колебаниях в конечном числе точек 86
4.1 Краткая история изучения задачи определения начального возмущения для линейных уравнений мелкой воды 87
4.2 Постановка задачи и ее разрешимость 88
4.3 Вариационная постановка обратной задачи. Градиент целевого функционала 91
4.4 Численное решение прямой и сопряженной задач: уровень вычислительной ошибки 93
4.5 Степень некорректности обратной задачи 97
4.6 Метод сопряженных градиентов: численные расчеты 100
4.7 Совмещенная постановка обратной задачи для уравнений мелкой воды в линейном приближении 105
4.7.1 Вариационная постановка совмещенной обратной задачи 105
4.7.2 Результаты численных расчетов 106
4.7.3 Преимущество совмещенных данных 108
5 Численный алгоритм определения амплитуды переднего фронта волны 109
5.1 Уравнение эйконала ПО
5.1.1 Уравнение эйконала в геометрической оптике. Принцип Ферма ПО
5.1.2 Схема С.К. Годунова 111
5.1.3 Метод бихарактеристик 113
5.1.4 Метод Рунге-Кутты решения уравнений Эйлера 114
5.1.5 Сравнение метода Годунова и метода бихарактеристик 115
5.2 Алгоритм определения амплитуды фронта волны в случае линейного источника 117
5.2.1 О существовании и единственности решения задачи определения амплитуды переднего фронта волны 119
5.2.2 Сравнение с акустическим фронтом в одномерном случае 120
5.3 Алгоритм определения амплитуды фронта волны в случае точечного источника 121
5.4 Численные эксперименты. Амплитуда переднего фронта волны 125
Заключение 127
Список рисунков 132
Список таблиц 133
Список работ, опубликованных по теме диссертации
- Понятие s-чисел вполне непрерывного и ограниченного операторов
- Математическая постановка прямой и обратньгх задач
- Теорема единственности
- Метод Рунге-Кутты решения уравнений Эйлера
Понятие s-чисел вполне непрерывного и ограниченного операторов
В 3.4 проведен анализ обратной задачи в терминах коэффициентов Фурье qk(x), основанный на построении сингулярных чисел оператора обратной задачи А о (к) : L2(0,L) — L2(0,L) в случае дН(х) = 1 и дискретного аналога этого оператора в случае одномерного рельефа дна Н(х). Для построения матрицы оператора AD(k) к обратной задаче применялась явная конечно-разностная консервативная схема второго порядка аппроксимации.
В 3.5 получена явная формула градиента целевого функционала, основанная на решении соответствующей сопряженной начально-краевой задачи и используемая в численных расчетах для метода простой итерации.
В 3.6 приведены результаты численных расчетов методами усеченного сингулярного разложения и простой итерации. Для зашумленных модельных данных с ошибкой є номер обрыва быстро убывающих сингулярных чисел был согласован с ошибкой є согласно 1.1.
Четвертая глава посвящена разработке, исследованию и численной реализации вариационного подхода определения начального возмущения q(x,y) из (2) по дополнительным измерениям rj(xm,ym,t) = fm(t), t Є (Tmi,Tm2), вертикального смещения свободной поверхности rj(x,y,t) в конечном числе заданных точек (хт,ут), т = 1,2,...,М (обратная задача 1), и численного алгоритма решения совмещенной обратной задачи, состоящей в определении функции q(x,y) одновременно по двум типам данных: fm(t), т = 1,2, ...,М и rj(x,y,T) = f(x,y) (определение q(x,y) по данным второго типа будем называть обратной задачей 2).
В 4.1 приведен краткий исторический обзор изучения свойств обратной задачи 1 (03 1) для уравнений мелкой воды (А.С. Запреев, ВА. Цецохо, В.Г. Романов, ТА. Воронина, ВА. Чеверда и др.).
В 4.2 приведена постановка прямой и 03 1 Ащ = F\ := (/і,..., /м)т- Компактность оператора Ах показана в работе ТА. Ворониной [27].
В 4.3 предложен вариационный подход решения 03 1, основанный на минимизации целевого функционала J\(q). Получена явная формула градиента J[q целевого функционала J\(q), опирающаяся на решение соответствующей сопряженной задачи, которая затем была использована в реализации метода сопряженных градиентов численного решения 03 1.
В 4.4 описана явная конечно-разностная консервативная схема второго порядка аппроксимации для численного решения прямой и сопряженной задач. Проведен анализ точности разностной схемы на аналитическом решении и получен уровень вычислительной ошибки.
В 4.5 исследована степень некорректности 03 1, основанная на анализе степени убывания сингулярных чисел дискретного аналога оператора 03 1 на коэффициенты Фурье qk (х), согласно подходу Т.А. Ворониной и В.А. Чеверды [29,41].
В 4.6 приведены и проанализированы результаты численных расчетов метода сопряженных градиентов, описано условие остановки метода.
4.7 посвящен исследованию, численной реализации и анализу результатов расчетов совмещенной обратной задачи для уравнений мелкой воды методом сопряженных градиентов. Получен явный вид градиента целевого функционала совмещенной обратной задачи и приведены результаты численных расчетов, демонстрирующие преимущество совмещения данных.
В пятой главе приведен алгоритм вычисления амплитуды переднего фронта волны S(x, у) в приближении мелкой воды, возбуждаемой кратковременным (импульсным) источником. Предложенный алгоритм позволяет сократить время и объем вычислений, переходя от решения задачи определения функции трех переменных rj(x,y,t) к задаче определения функции двух переменных S(x,y).
В 5.1 проанализированы методы решения возникающего в алгоритме уравнения эйконала ТІ + Ту = (дН(х,у)) 1, а именно схема С.К. Годунова (определение установившегося решения нестационарного уравнения rt + т + Ту = (дН(х,у)) 1) и метод бихарактеристик. Приведены результаты численных расчетов, демонстрирующие преимущества каждого из методов.
В 5.2 описан алгоритм определения амплитуды фронта волны в случае линейного источника q(x,y) = h(y)8(x), основная идея которого состоит в замене переменных (х,у) ь- (z,y), z = т(х,у) [46]. Получена задача Коши определения амплитуды переднего фронта волны S l\z,y), решение которой в одномерном случае находится в явном виде и совпадает с известной формулой Эри-Грина (1841 год), согласно которой амплитуда фронта волны в океане растет с уменьшением глубины.
В приложении А приведены вывод, основные свойства уравнений мелкой воды и четыре варианта граничных условий для системы уравнений мелкой воды (отражающая граница, впускающая и выпускающая границы, подвижная граница). В приложении В доказаны теоремы сходимости градиентных методов, используемых в работе: методов простой итерации и сопряженных градиентов. В приложении С изложен физический вывод формулы Эри-Грина определения амплитуды волны в одномерном случае. В приложении D приведены краткие сведения (формулы, теоремы) сингулярного разложения матриц, интегральных и компактных операторов. В приложении Е представлено свидетельство о регистрации программы для ЭВМ в Фонде алгоритмов и программ СО РАН.
Автор выражает благодарность своему научному руководителю член-корреспонденту Российской академии наук, профессору Сергею Игоревичу Кабанихину за интересную постановку задач, постоянное внимание к настоящей работе и помощь, которые способствовали успешному написанию диссертации и формированию научных взглядов. А также за руководство многих серьезных проектов, в которых автор принимал участие.
Автор выражает искреннюю и глубокую благодарность кандидату физико-математических наук М.А. Шишленину за бесценный опыт и знания, полученные в ходе выполнения работы, за всестороннюю поддержку, постоянное внимание и многочисленные обсуждения, способствовавшие успешному решению многих задач. За расширенные и подробные консультации относительно оптимизации программных комплексов автор благодарен Зятькову Николаю. Результатом полезных обсуждений стало увеличение скорости выполнения вычислений, что сделало программные инструментарии более эффективными.
За совместную работу и полученный опыт автор выражает благодарность бывшим и настоящим сотрудникам новосибирского отделения ОАО «ГеоСистема»: Грачёву Евгению, Зай-ченко Ольге, Комарову Виктору, Маринину Игорю, Романову Сергею, Семёнову Александру и Чеснокову Вадиму. Именно в указанном коллективе под техническим и функциональным руководством И.В. Маринина автор смог реализовать многие из своих идей в виде комплексов программ и отдельных компонент.
Автор выражает благодарность всем сотрудникам Института вычислительной математики и математической геофизики СО РАН и кафедры теории функций НГУ за знания, полученные на лекциях, семинарах и практических занятиях, а также за полезные научные дискуссии в ходе выполнения работы.
Автор выражает благодарность оппонентам диссертации: профессору А.Г. Яголе, профессору В.Г. Чередниченко и к.ф.-м.н. О.Б. Бочарову за ценные советы и комментарии к работе, неоценимое внимание и помощь в ходе ее выполнения. Глава 1
Сингулярное разложение в некорректных задачах
Теорема о сингулярном разложении компактного оператора дает описание некорректных задач, а именно, степень некорректности задачи Aq = / характеризуется характером убывания к нулю последовательности сингулярных чисел. Разумеется само построение сингулярного разложения какого-либо оператора может стать слишком сложной задачей для практического применения. Однако на практике линейные некорректные задачи сводятся к системам алгебраических уравнений, а для сингулярного разложения матриц может быть применено несколько стандартных приемов.
В данной главе мы приведем теорему о сингулярном разложении матрицы и ее обобщения на случай интегральных и компактных операторов (подробное доказательство всех теорем приведено в Приложении D). Изложенные результаты использованы для исследования некорректности ряда обратных задач и построения алгоритмов регуляризации.
Пусть задана система линейных алгебраических уравнений Aq = f, где А - действительная тх гг-матрица, q Є Кга, / Є Km, и требуется решить эту систему относительно неизвестного вектора q. Очевидно, что степень трудности решения такой системы зависит от структуры матрицы А. На языке линейной алгебры матрица А является представлением некоторого линейного оператора в конкретной координатной системе. Естественно попытаться так преобразовать эту координатную систему, чтобы линейный оператор в новой координатной системе был представлен матрицей более простого вида (например, матрицей блочной структуры с нулевыми блоками). Иными словами, матрицу А можно разложить на произведение матриц специального вида. В настоящее время разработан ряд алгоритмов, основанных на разложении матрицы А. Одним из достоинств этого подхода к решению системы линейных уравнений является его универсальность, то есть применимость как к хорошо обусловленным системам, так и к системам с вырожденной или плохо обусловленной матрицей.
Алгоритм сингулярного разложения матрицы А состоит в ее разложении на произведение двух ортогональных матриц и диагональной, по диагонали которой в порядке невозрастания расположены сингулярные числа. Алгоритм удобен для исследования особенностей матрицы А посредством анализа ее сингулярных чисел (Тк(А) в силу связи числа обусловленности матрицы cond(A) = атах(А)/ат[П(А) с ее сингулярными числами. Чем меньше минимальное сингулярное число 7т;п матрицы, тем больше ее число обусловленности cond(A). Краткая история развития сингулярного разложения матриц изложена в Приложении D.I. Теорема Гильберта-Шмидта для интегрального оператора (см. Приложение D.2) явилась некоторым обобщением сингулярного разложения матриц. Затем область применения теоремы о сингулярном разложении дошла и до компактных (вполне непрерывных) операторов (см. Приложение D.3). Формулировки теоремы о сингулярном разложении матриц, интегральных и компактных операторов представлены в Диаграмме 1.1.
Математическая постановка прямой и обратньгх задач
Методы томографического исследования структуры неоднородных объектов, благодаря успехам как фундаментальной и вычислительной математики, так и современному аппаратному оснащению, произвели революцию в биологических и медицинских исследованиях. Нескольким ученым, применявшим в своих исследованиях методы томографии, были присуждены Нобелевские премии по медицине (в 1979 году G.N. Hounsfield и А.Мс. Cormack получили Нобелевскую премию «за разработку компьютерной томографии»; в 2003 году Р.С. Lauterbur и S.P. Mansfield получили Нобелевскую премию «за изобретение метода магнитно-резонансной томографии»), биологии и химии. Прогресс медицинской томографии сопровождался зарождением и развитием многих других приложений этого универсального и информативного метода. Будучи эффективным средством решения обратных задач, направленных на реконструкцию внутренней структуры сложных объектов, методы вычислительной томографии стали применяться в оптике и физическом эксперименте [78], геофизике, дефектоскопии промышленных изделий и других областях науки и промышленности [79,80]. Развитие математических методов и систем сбора информации привели к новым математическим моделям томографии, таким как диффузионная томография, термотомография, томография векторных и тензорных полей [56,81,82]. Разнообразие математических постановок связано с потребностью реконструкции характеристик сред различной сложности. Более подробный обзор работ, посвященных математическим формулировкам задач восстановления векторных и тензорных полей, приведен в кандидатской диссертации А.П. Поляковой [83].
Задача термоакустики имеет важное значение в онкологии, поскольку при исследовании пораженных тканей обнаружено, что массовая доля воды в них существенно больше, чем в здоровых тканях [2-4]. Если в здоровых тканях содержание воды составляет около 65%, то в пораженных оно достигает 80%. Столь существенное отличие массовой доли воды в пораженных и здоровых клетках дает основание предполагать, что коэффициенты поглощения энергии электромагнитного излучения пораженных и здоровых клеток также существенно отличаются. Это предположение подтверждается результатами экспериментальных работ [2], в которых показано, что коэффициент поглощения электромагнитной энергии пораженных раком жировых тканей груди в 1,8 раз больше, чем у здоровых. Грудь облучалась лазером с рабочей частотой 434 МГц. Волны акустического давления фиксировались пьезодатчика-ми, размещенными на поверхности исследуемой области, имеющей форму полусферы радиуса 15 см. Пространство между грудью и датчиками заполнялось жидкостью, плотность и скорость звука которой близка к плотности и скорости звука груди. Задача определения коэффициента поглощения энергии электромагнитного излучения (источника излучения) по измерениям акустического давления на границе (или ее части) называется обратной задачей термоакустики.
Изучение свойств обратной задачи термоакустики и ее численное решение начало бурно развиваться с 90-х годов. Е.Т. Quinto, А.К. Louis и их последователи используют преобразование Радона для решения обратной задачи термоакустики в интегральной постановке [6,7,84]. В работе О. Schertzer и соавторов ( [85], 2004) численно решена обратная задача термоакустики методом обработки изображений на основе преобразования Радона для реальных данных. В работах М. Xu, L.V. Wang ( [86,87], 2002-2003) и R.A. Krager с соавторами ( [2], 1999) построены алгоритмы восстановления источника возмущений во временной области при постоянной скорости распространения акустических волн на основе функции Грина. В случае сложных областей построение функции Грина затруднительно. В работах М. Xu, L.V. Wang ( [88], 2005) и М. Ideman ( [89], 2012) получены явные формулы определения источника излучения в случае простых (две плоскости, бесконечный цилиндр, сфера) и произвольных областей, соответственно. Полученные формулы неудобны в практическом применении, так используют многократные интегралы (суммы). В работе M.V. Klibanov и соавторов ( [90], 2008) получены глобальные оценки устойчивости решения обратной задачи термоакустики как коэффициентной обратной задачи. В работе P. Stefanov и G. Uhlmann 2009 года [91] построено явное решение обратной задачи термоакустики в терминах разложения в ряд Неймана для данных, заданных на всей границе. В этой же работе получены необходимые и достаточные условия единственности и устойчивости обратной задачи при измерениях, заданных на части границы области. В случае измерений, заданных с ошибкой (зашумленные данные), обратная задача термоакустики неустойчива. В. Kaltenbacher с учениками ( [92], 2011) обосновала численные методы регуляризации некорректной задачи термоакустики: метод М.М. Лаврентьева, регуляризация путем дискретизации и метод, основанный на явной формуле решения с применением численной регуляризации. Более подробный обзор освещен в списках литературы вышеуказанных работ. В работе С. Zhang и соавторов ( [93], 2014) приближенное решение обратной задачи термоакустики получено с помощью адаптированного метода полной вариации (Total Variation). Показано преимущество использования этого метода в случае недостатка данных измерений. В данной главе проведено исследование степени некорректности обратной задачи термоакустики для данных, заданных на одной плоскости измерений, на двух и на трех. Показано, что степень некорректности обратной задачи уменьшается при увеличении числа измерений. Построены методы регуляризации обратной задачи: усеченное сингулярное разложение, итерационная регуляризация, разложение решения в ряд Фурье.
Результаты данного параграфа опубликованы в двух статьях [94,95]. В разделе 2.1.1 приведена физическая постановка задачи термоакустики. В разделе 2.1.2 сформулирована прямая и три обратные задачи термоакустики. В разделе 2.1.3 показана единственность решения задачи восстановления функции по ее сферическим средним. В разделе 2.2 описаны используемые в работе численные алгоритмы решения прямой задачи термоакустики, а именно, численное интегрирование (раздел 2.2.1), метод Фурье (раздел 2.2.2) и конечно-разностные схемы (раздел 2.2.3). В разделе 2.3 проведен анализ некорректности обратной задачи термоакустики на основе анализа убывания сингулярных чисел дискретного оператора обратной задачи. В разделе 2.3.3 сформулировано условие устойчивости обратной задачи. В разделе 2.4 приведены результаты численных расчетов обратных задач термоакустики, полученные методами простой итерации и сопряженных градиентов. В разделе 2.5 приведены результаты численных расчетов обратных задач термоакустики, полученные методом усеченного сингулярного разложения (г-решение) с использованием априорной информации о решении, согласно приему С.К. Годунова. 2.1 Физическая и математическая модели
В данном разделе приведена физическая и математическая постановки задачи термоакустики, которая заключается в определении функции и(х,у)9 у Є М2, описывающей процесс распространения волн акустического давления. В зависимости от количества измерений на части области в разделе 2.1.2 формулируются три обратные задачи термоакустики, которые заключаются в определении коэффициента поглощения электромагнитного излучения в рассматриваемой области по измерениям акустического давления на частях границы. В разделе 2.1.3 показана единственность решения задачи восстановления функции по ее сферическим средним.
Рассматриваем область Q С R3 упругой среды (Рисунок 2.1). Предположим, что начиная с момента времени t = О, область Q подвергается электромагнитному излучению интенсивностью I(t)9 которое частично поглощается средой. Поглощенная энергия переходит в тепло, что приводит к увеличению температуры среды, к ее расширению и, в итоге, к появлению волн акустического давления. Распространяясь по среде, волны акустического давления достигают границы Г области, на части которой Г і они могут быть измерены. Требуется определить коэффициент поглощения электромагнитного излучения в Q по измерениям акустического давления на части границы 1\.
Рисунок 2.6: Расчетная область ш. Черными точками изображены заданные значения функции г)1 (согласно равенствам (2.25), (2.26) и (2.27)). Значения гц\ в остальных узлах сетки определяются из схемы согласно шаблону: значение в узле {i,j + 1} (красная точка) определяется через значения функции в узлах {і + 1, j}, {і — l,j}, {i,j — 1} (желтые точки) по равенству (2.24); граничные значения на j + 1-ом слое (красные точки) определяются через значения функции в узлах {l,j}, {0, j — 1} и {Л — l,j}, {Nx,j — 1} (желтые точки) по равенствам (2.28).
В данном разделе приведем анализ некорректности обратной задачи термоакустики, основанный на анализе степени убывания сингулярных чисел оператора обратной задачи (см. Раздел D.1.2). В силу того, что поиск совокупности {uk,Vk,(Jk}, удовлетворяющей равенствам Auk = (JkVk и A Vk = JkUk, является достаточно сложной задачей, мы будем анализировать сингулярные числа дискретного аналога оператора обратной задачи, который представляет собой прямоугольную матрицу. В случае матриц существует множество алгоритмов определения ее сингулярных чисел, большинство из которых описано в монографии С.К. Годунова [61] (разложение Холесского, методы главных компонент, бисекций, обратной итерации, Годунова с итерационным уточнением и т.п.).
В разделе 2.3.1 обратные задачи 1, 2 и 3 представлены в виде систем линейных алгебраических уравнений A\k)q{k) = f y A\k)q{k) = (/(\),/ffc)) и А3{к)д{к) = (/( ,/ffc),/ffc)) , соответственно. В разделе 2.3.2 проведен анализ сингулярных чисел матриц Alky j = 1, 2, 3. В разделе 2.3.3 сформулировано условие устойчивости обратной задачи термоакустики.
Пользуясь результатами раздела 2.2.2, при любых га, п = 0,1, 2,..., определим операторы обратных задач 1 (2.17)-(2.18), 2 (2.17)-(2.19) и 3 (2.17)-(2.20):
Решения обратных задач для ограниченного числа коэффициентов Фурье уже в некотором смысле является регуляризацией обратных задач 1, 2 и 3. При фиксированном к решение задачи (2.17) «хорошее» (т.е. задача (2.17) корректна по Адамару), однако при увеличении к степень корректности уменьшается.
Используя явную конечно-разностную схему второго порядка аппроксимации, приведенную в Разделе 2.2.3, представим обратные задачи 1 (2.29), 2 (2.30) и 3 (2.31) в матричном виде.
Пусть Nt - нечетное. Представим обратную задачу 1 (2.29) в виде: A}k-.Qk = Fk. Здесь Qk = {q{k)0,q{k)1,---,q{k)Nx)T - вектор размерности Nx + 1, F = (7 , / , , f ) вектор данных задачи размерности Nt + 1, Ah-, - прямоугольная матрица размерности (Nt + векторы размерно Пусть vi= (u{k)l)u{k)2v...)u{k)2NJj KV2 = (u{k)30,u{k)3v...,u{k)l Лемма 2.1. Векторы v\ и v2 выражаются следующим образом: v1 = B1\j\ v2 = B2Vl + B3U. Здесь матрицы В\ и В3 размерности (Nx + 1) х (2NX + 2) и квадратная матрица В2 размерности Nx + 1 имеют следующий вид: В9 = /0 2г 0 0 0\ V 0 V 0 0 0 V 0 V 0 B1=(-INx+1 В2) B3=(oNx+l -lNx+l 0 0 V 0 V \о 0 0 2г О/ Здесь NX+I и INX+I - нулевая и единичная квадратные матрицы размерности Nx + 1, соответственно. Доказательство. Действительно, для получения матрицы В\ подставим j = 1 в разностную схему задачи (2.24) и граничные условия (2.28). Аналогично получим матрицы В2 и з, подставив j = 2 в эту же разностную схему и граничные условия. Идея построения матрицы А1к, заключается в следующем. Используя результаты Леммы 2.1, составим квадратную матрицу С размерности 2NX + 2: удовлетворяющую соотношению U = CU . Затем, учитывая линейность задачи и четность Обратная задача 3 (2.31) записывается в матричном виде Ah,Qk = Fk, где Fk = Анализ сингулярных чисел дискретных операторов обратных задач Для вычисления сингулярного разложения матриц Ah-,, Ah-, и Ah-, использовалась стандартная функция GESVD библиотеки Intel Math Kernel Library LAPACK [104].
На Рисунке 2.7 приведены графики сингулярных чисел матриц Ah-,, Ah-, и Ah-, дискретных обратных задач 1, 2 и 3 на коэффициенты Фурье к = 1,5,15, построенных на равномерной сетке Nx = 500, Nt = 1055 области Lx = 5, Ly = 10, Т = 3.
Сингулярные числа обратной задачи 1 на порядок меньше сингулярных чисел обратной задачи 2. В свою очередь, сингулярные числа обратной задачи 2 меньше сингулярных чисел обратной задачи 3. Здесь наблюдается известное всем утверждение: добавление дополнительных строк (столбцов) к исходной матрице не увеличит ее число обусловленности, а только уменьшит в случае, если строки (столбы) являются линейно независимыми со строками Сингулярные числа матриц Ah , (красная линия), А?к , (зеленая линия) и ААч (синяя линия) дискретных аналогов обратных задач 1, 2 и 3 при к = 1, 5,15. Ось Ох - номер сингулярного числа, ось Оу - логарифмическая шкала значений сингулярных чисел. (столбцами) исходной матрицы. На этом эффекте и строится прием С.К. Годунова, представленный в разделе D.1.3: с увеличением информации о решении обратной задачи повышается устойчивость решения обратной задачи термоакустики.
Теорема единственности
На рисунке 2.18 изображены графики точного и полученных методом сингулярного разложения решений обратных задач 1, 2 и 3 на коэффициенты Фурье при к = 4 с учетом априорной информации о решении в виде матрицы В\ с внесенной ошибкой в правых частях є = 0.01. Отметим, что точность восстановления решения увеличивается с добавлением априорной информации В\ (сравнить рисунки 2.17 и 2.18). В свою очередь, в случае обратной задачи 3 решение восстанавливается точнее, чем в случаях обратных задач 1 и 2. На рисунке 2.19 продемонстрированные двумерные графики восстановления функции q?2 (2.53) методом сингулярного разложения с использованием априорной информации в виде матрицы Вх (прием С.К. Годунова).
На рисунке 2.20 изображены графики точного и полученных методом сингулярного разложения решений обратных задач 1, 2 и 3 на коэффициенты Фурье при к = 4 с учетом априорной информации о решении в виде матрицы В2 с внесенной ошибкой в правых частях є = 0.01. Отметим, что точность восстановления решения выше, чем в случае априорной информации В\ (сравнить рисунки 2.18 и 2.20). В свою очередь, в случае обратной задачи 3 решение также восстанавливается точнее, чем в случаях обратных задач 1 и 2. А решение "rez.dath
Графики решений обратных задач: а) 1, Ь) 2, с) 3 при к = 4, полученные методом сингулярного разложения, с априорной информацией о решении в виде матрицы В2 с параметром регуляризации а = 2 6. Уровень ошибки в правых частях є = 0.01. Красная линия - коэффициенты Фурье точного решения, зеленая (синяя) линия - коэффициенты Фурье полученного решения. обратной задачи 2 ближе к точному, чем решение обратной задачи 1. То есть, чем лучше (пол нее) априорная информация о решении обратных задач (В2 более информативна, чем В{), тем точнее полученное решение к точному.
На рисунке 2.21 продемонстрированные двумерные графики восстановления функции QT2 (2.53) методом сингулярного разложения с использованием априорной информации в виде матрицы 2 (прием С.К. Годунова).
а) График точного решения QT2. Графики восстановленных решений обратных задач Ь) 1, с) 2 и d) 3 методом сингулярного разложения с помощью приема С.К. Годунова с априорной информацией в виде матрицы В2, с внесенной ошибкой в правых частях є = 0.01.
Сравнивая графики друг с другом (Рисунки 2.19, 2.21 сравнить между собой), заключаем, что при данных на трех частях границы (обратная задача 3), решение восстанавливается лучше, нежели в случае двух измерений (обратная задача 2). Также, в случае двух измерений решение восстанавливается лучше, чем при задании всего одной дополнительной информации на границе (обратная задача 1). отметим, что аналогичный результат мы получили, применяя градиентные методы (см. Раздел 2.4.2). Причем, чем точнее априорная информация (В2 точнее и обширнее, чем В\), тем лучше восстановление, даже при зашумленных данных. Глава З
Определение источника колебаний в волновом уравнении по данным о колебаниях в фиксированный момент времени
Круглый камень, брошенный вертикально вниз на гладкую водную поверхность, возбуждает концентрические волны. Если этот процесс запечатлеть на видеокамеру, а затем просмотреть запись в обратном времени [109], то можно остановить изображение в тот момент, когда камень касается воды. Примерно так, по образному сравнению Б.Я. Зельдовича, можно представить процесс обращения волнового поля в двумерном случае. Интересно, что почти аналогичный процесс происходит при возникновении волн цунами. Но в этом случае, наоборот, водная поверхность не углубляется под действием падающего сверху камня, а поднимается вверх под действием либо подъемом части дна вследствие сдвига тектонических плит (разумеется, при этом часть дна может опускаться, и водная поверхность тоже будет опускаться), либо извержения подводного вулкана. Если бы у нас была возможность (например, со спутника) записать на видеокамеру изменение водной поверхности с момента возникновения цунами вплоть до подхода волны к берегу, то мы получили бы картину, в некотором смысле аналогичную падению камня в воду. Одной из важнейших задач при исследовании цунами является оперативное определение места и формы возникновения цунами, поскольку зная источник и рельеф морского дна, можно составить оперативный прогноз, т.е. посчитать время прихода и амплитуду волны, которая придет в данную точку берега. Такие расчеты можно провести, решая систему уравнений мелкой воды, начальными данными для которой будет г](х, у, 0) = q(x, у) - форма изменения водной поверхности в начальный момент времени.
В данном разделе рассматривается задача определения скорости вертикального изменения водной поверхности r]t(x,y,0) = q(x,y) в предположении, что в момент времени t = Т измерена форма поверхности водоема / 2 (ж,г/) = rq(x,y,T) (обратная задача). Источник возмущения г](х, у, 0) = f (x, у) расположен в области П := (0, L) х (0, L). Задача определения функции г](х, у, t) (а, значит, и q(x, у)) при заданных функциях f (x, у), f 2\x, у) и однородных граничных условиях является задачей Дирихле для волнового уравнения. Как известно, задача Дирихле для волнового уравнения некорректна. Этому вопросу посвящено множество работ. Впервые неединственность задачи Дирихле для волнового уравнения была отмечена в работах Ж. Адамара (1936 г., [ПО]), А. Губера (1932 г., [111]), Д. Манжерона (1932 г., [112]).
В своей работе Д. Буржин и Р. Даффин (1939 г., [70]) рассмотрели задачу Дирихле в прямоугольнике R : {0 t Т; 0 х S} для уравнения затухающих колебаний струны
Vtt Vxx - k2r] = 0. Используя преобразование Лапласа, они доказали, что если число T/S иррационально, то имеет место единственность решения задачи в классе функций, которые непрерывно дифференцируемы и имеют в R суммируемые по Лебегу вторые производные. Они также доказали существование решения задачи Дирихле при некоторых ограничениях на параметры Т, S и к. В работе С.Г. Овсепяна (1964 г., [113]) рассмотрена задача (1 + А)0-(1-А)0 = О, V\r = ( ) О-1) для ограниченной многосвязной области D с границей Г (здесь Л - действительный параметр, Л 1). При некоторых ограничениях на границу Г установлена единственность решения задачи (3.1) в первом классе Бэра. Показано, что задача (3.1) неустойчива относительно изменения области не только для односвязных областей, но и для произвольной ограниченной многосвязной области D с кусочно-гладкой границей Г.
В работах С.Г. Овсепяна (1967 г. и 1969 г. [114,115]) изучены вопросы единственности решения и построения обобщенных собственных функций задачи Дирихле для уравнения (3.1) в классе измеримых функций.
Метод Рунге-Кутты решения уравнений Эйлера
Первые два равенства при фиксированной х0 определяют в пространстве переменных x,y,t двухпараметрическое семейство бихарактеристик. В пространстве К2 они определяют геодезические линии Г(х,хо). Проекцию бихарактеристики на пространство (х,у) называют лучом. Первые два равенства задают луч параметрически. Для однозначного задания луча, соединяющего пару точек (х,х0), необходимо выполнения условия регулярности семейства геодезических в области, ограниченной характеристическим коноидом t = т(х,у) (вопрос единственности решения уравнения эйконала (5.3) подробно изучен в работе [56]).
Лемма 5.1. [56] Пусть D - открытая область пространства Е2 с гладкой границей, с(х,у) Є Ck(D), k 2, с(х,у) Со(х,у) 0 при (х,у) Є D. Тогда для любой точки (хо, уо) Є D и любого единичного вектора и0 решение задачи (5.8) существует и единственно для всех т, при которых точка (/і(г, 0,хо), /2(1", 0,хо)) Є D, а функции f\, f2, /з, Д, fiTT, [2ТТ, hT, fiT непрерывны no совокупности аргументов вместе с производными до порядка k-l.
Решение задачи (5.8) определяет геодезическую, которая проходит через точку (ж0,Уо) в направлении вектора (о = (со sin0, с0 cos0). Введем вектор ( = т(о = ((i,(2), гДе Сі = тс0 sin 0, (2 = тс0cos 0. В работе [47] показано, что функции х = /i(,x0) и у = /2(( хо) имеют обратные функции d = #i(x, х0) и (2 = д2(х,х0), соответственно, определенные в окрестности точки х0. Тогда получим следующее представление для функции т2(х,у):
Перейдем к описанию разностного метода решения задачи Коши для системы обыкновенных дифференциальных уравнений (5.8). Пусть t Є (0,Т). Построим разбиение области
Систему (5.8) мы решаем методом Рунге-Кутты четвертого порядка аппроксимации, который состоит в определении приближенного решения системы
Здесь Nx и Ny - количество точек разбиения промежутков (О, Lx) и (0, Ly), соответственно.
Пример 1. Расчет времени прихода первых волн от точечного источника в изотропной среде, расположенного в центре квадратной области размера Lx = Ly = 2 при с(х,у) = 1. Область покрывалась сеткой Nx = Ny = 101.
Результат, полученный при помощи схемы Годунова второго порядка точности и метода бихарактеристик показан на рисунке 5.1.
Отметим, что полученные численные решения уравнения эйконала близки к явному точному решению (х— I)2 + (у— I)2 = 1. Схема С.К. Годунова имеет второй порядок точности [179], а метод Рунге-Кутты имеет четвертый порядок точности.
Пример 2. Расчет времени прихода первых волн в глубокой воде, ограниченной прямоугольной областью размера Lx = 70 км, Ly = 10 км. Область покрывалась сеткой Л = 450,
Для расчета времени прихода первых волн от линии источников методом бихарактеристик (см. Раздел 2.3.2) источник вида h(y) 8(х) в задаче (5.12) представляется в виде ряда точечных источников, для каждого из которых отдельно решается уравнение эйконала (5.3) с условием т(хо,уо) = 0. После решения системы (5.8) для гг-го точечного источника получаем массивы координат (a;ra[i][j],yra[i][j]) и значений rra[i][j] линий уровня, п = 1,...,NP, Np -количество точек разбиения линии источников. Здесь г = 0,..., 359 обозначает градус угла, под которым выпускается направляющий луч р0 (см. систему (5.8)), j = 0,..., Nt - текущий номер по времени (или, что тоже самое, номер линии уровня т{х,у)). Итоговые линии уровня r[i][j] для линии источников получаем после построения огибающей при фиксированном j для всех і и п. На рисунке 5.2 приведены расчеты времени прихода первых волн от линии источников от точки (3.5,3) до точки (3.5,7) для методов бихарактеристик и Годунова, соответственно. Время расчетов ограничивалось Т = 10 минутами. В случае метода бихарактеристик источник разбивался на Np = 60 точечных источников, Nt = 30. При реализации схемы Годунова количество точек разбиения по времени Nt = 18000 для соблюдения условия Куранта устойчивости разностной схемы, а также достижения условия установления. 10 20 30 40 50 60 70 0 10 20 30 40 50 60 70
Расчет времени прихода первых волн от линии источников: а) методом бихарактеристик; б) методом Годунова. Отметим, что вычислительное время на персональном компьютере составило 27,5 минут для метода Годунова и 28,5 секунд для метода бихарактеристик. Вычисления проводились на процессоре Intel(R) Соге(ТМ) І7-2600 CPU 3,4 GHz.
В случае рельефа дна с волнорезом, изображенным на рисунке 5.3, расчеты времени прихода первых волн от линии источников для методов бихарактеристик и Годунова приведены на рисунке 5.4.
Метод бихарактеристик больше подходит для расчета прихода первых волн от точечного источника для произвольного рельефа дна по двум причинам: (1) в 60 раз быстрее подхода С.К. Годунова и (2) представление в форме систем обыкновенных дифференциальных уравнений позволяет применять разностные схемы высоких порядков точности (схема Рунге-Кутты). Однако в случае произвольного источника (линейный, произвольной формы) подход С.К. Годунова с заданной точностью (в работе представлены расчеты со вторым порядком точности) получит картину прихода первых волн в каком-то смысле точнее метода бихарактеристик. Метод бихарактеристик в таком случае может быть применим (см. Рисунки 5.2а и 5.4а), однако громоздкость вычислений (количество точечных разбиений на источники) затрудняет определить точность расчетов, которая зависит от плотности разбиения, принципа построения огибающей и вычислительных ошибок.
В данном разделе предложен алгоритм расчета амплитуды переднего фронта волны, возбуждаемой сравнительно кратковременными воздействиями. Известно [182], что в случае источника g(x,t), х Є К3, сосредоточенного в сравнительно небольшой области (suppg С Q х (0, to), где diamU и 0 предполагаются достаточно малыми), быстроизменяющаяся «сиг-нальная часть волнового поля ?(x,t) сосредоточена в окрестности поверхности характеристического
Здесь w(x, t) - фундаментальное решение системы уравнений распространения волны в теории мелкой воды. Сказанное означает, что амплитуда переднего фронта волны с достаточно хорошей точностью может быть приближена амплитудой соответствующей сингулярной составляющей фундаментального решения. В самом деле, формула (5.11) показывает, что чем ближе g(x0, t0) к дельта-функции Дирака #(х0, t0), тем лучше г (х, t) приближает r?(x0, t0).
Подход к построению теоретических решений для уравнений мелкой воды с постоянной глубиной был предложен С.С. Войтом и Б.И. Себекиным в 1968 году [183], основанный на построении и изучении свойств передаточных функций, представляющих собой отклик океана на определенное эталонное воздействие. В работе [184] предложен асимптотический метод определения амплитуды волнового фронта, основанный на применении обобщенной конструкции канонического оператора В.П. Маслова [185]. В работе В.Г. Романова (2005, [47]) предложен подход к определению амплитуды фронта волны для произвольной глубины в случае точечного источника в многомерном варианте. В работе Ан.Г. Марчука и Г.С. Васильева (2014, [186]) представлен быстрый численный метод определения амплитуды фронта волны от точечного источника в одномерном случае, основанный на кинетическом подходе. В работе Косых B.C., Чубарова Л.Б. и соавторов [187] исследована методика расчета максимальных высот волн цунами при реальных данных, основанная на полном решении системы уравнений мелкой воды, что в большинстве случаев является сложной вычислительной задачей. В случае произвольного дна в данном разделе приведем алгоритм определения амплитуды переднего фронта волны, основанный на выявлении свойств фундаментального решения системы уравнений мелкой воды [176]. Предложенный алгоритм в 30 раз быстрее вычисляет амплитуду фронта волны в отличие от полного решения системы уравнений мелкой воды.