Модели и алгоритмы решения прямых и обратных задач гравиразведки Рязанцев Владимир Андреевич

Содержание к диссертации

Введение

Глава 1. Постановка задачи математического моделирования прямых и обратных задач гравиразведки 13

1.1. Обзор методов решения прямых и обратных задач гравиразведки.. 13

1.2. Прямые задачи гравиразведки 18

1.3. Обратные задачи гравиразведки 26

Выводы 34

Глава 2. Аналитический и численный методы решения задачи одновременного восстановления плотности, формы и глубины залегания гравитирующего тела 35

2.1. Построение математических моделей 37

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

2.3. Построение численного метода решения обратных нелинейных задач гравиразведки 49

2.4. Газработка методики построения обобщенных обратных задач 50

2.5. Гешение модельных задач и интерпретация результатов 52

Выводы 60

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

3.1. Обзор современных результатов в области продолжения физических полей 62

3.2. Построение вычислительных схем аппроксимации потенциальных полей 63

3.3. Построение вычислительных схем аппроксимации тепловых полей 72

3.4. Гешение модельных примеров и интерпретация результатов 86

Выводы 95

Глава 4. Построение критериев устойчивости динамических систем и вычислительных алгоритмов 97

4.1. Исследование устойчивости решений дифференциальных уравнений в частных производных 97

4.2. Исследование устойчивости разностных схем 107

Выводы 117

Глава 5. Газработка программного комплекса 119

5.1. Описание алгоритма и программы восстановления физических полей в заданной области пространства 119

5.2. Описание алгоритма и программы определения формы и плотности гравитирующего тела при решении задачи потенциала 121

5.3. Выводы 124

Заключение 125

Список сокращений и условных обозначений 128

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

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

Введение к работе

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

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

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

Фундаментальное значение для современной гравиразведки играет такое свойство обратной задачи, как ее некорректность, обусловливающая значительные трудности при ее решении. Эти трудности связаны с возможной неединственностью решения обратной задачи при одинаковых входных данных, а также с чувствительностью решения задачи к малым колебаниям заданных характеристик гравитационного поля. Теория некорректно поставленных задач начала активно развиваться с исследований А. Н. Тихонова и получила дальнейшее развитие в работах В. К. Иванова и М. М. Лаврентьева. В развитие теории решения некорректно поставленных задач большой вклад внесли В. Я. Арсенин, А. Б. Бакушинский, В. Г. Васильев, В. В. Васин, В. В. Воеводин, В. Б. Гласко, А. В. Гончарский, О. А. Лис-ковец, Г. И. Марчук, В. А. Морозов, В. Н. Страхов, В. П. Танана, А. Г. Ягола и др. Важнейшие результаты в области решения обратных задач гравираз-

ведки были получены такими учеными, как 3. 3. Арсанукаев, В. М. Березкин, Ю. И. Блох, Е. Г. Булах, Г. М. Воскобойников, В. Б. Гласко, Г. Я. Голиздра, М. С. Жданов, А. А. Заморев, А. И. Кобрунов, А. К. Маловичко, П. С. Мар-тышко, Е. А. Мудрецова, П. С. Новиков, Б. В. Нумеров, С. М. Оганесян, И. Л. Пруткин, В. И. Старостенко, В. Н. Страхов, В. Г. Чередниченко, А. Ф. Шестаков, A. Bethencourt, М. Commer, О. L. Colombo, Т. Gerya, D. P. Hajela, D. J. Hofsommer, М. Н. Payne, F. Sanso и др.

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

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

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

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

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

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

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

Для достижения поставленной цели необходимо решить следующие основные задачи:

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

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

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

3. Разработать оптимальные алгоритмы аппроксимации потенциальных и тепловых полей.

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

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

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

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

Соответствие паспорту специальности. Диссертация выполнена в соответствии с требованиями специальности 05.13.18 - Математическое моделирование, численные методы и комплексы программ. Области исследования: 1) разработка новых математических методов моделирования объектов и явлений; 2) развитие качественных и приближенных аналитических методов исследования математических моделей; 3) разработка, обоснование и тестирование эффективных вычислительных методов с применением современных компьютерных технологий и 5) комплексные исследования научных и технических проблем с применением современной технологии математического моделирования и вычислительного эксперимента.

Научная новизна работы заключается в следующем:

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

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

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

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

5. Разработаны устойчивые разностные схемы восстановления потенциальных и тепловых полей, положенные в основу построения оптималь-

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

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

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

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

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

Основные результаты, выносимые на защиту:

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

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

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

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

10. Критерии устойчивости разностных схем с неравномерными сетками узлов.

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

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

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

Диссертационные исследования были проведены на кафедре высшей и прикладной математики факультета вычислительной техники ФГБОУ ВПО «Пензенский государственный университет». Выбранная тема исследований является частью научной работы, которая проводится на кафедре в рамках научно-исследовательских работ, выполняемых по государственному заданию Минобрнауки РФ:

- «Оптимальные методы вычисления гиперсингулярных интегралов, решения гиперсингулярных интегральных уравнений и их применение к задачам аэродинамики, электродинамики и геофизики» (2005-2009), регистрационный номер 0120.0502705;

«Аналитические и численные методы исследования динамических процессов в биосистемах и физической кинетике» (2010-2011), регистрационный номер 0120105992;

«Численные методы анализа прямых и обратных задач переноса излучения на наноструктурах» (2012-2013), регистрационный номер 1.656.2011.

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

Разработан программный комплекс (свидетельства № 2014663162 и № 2015610120 о государственной регистрации программы для ЭВМ) решения задачи восстановления потенциального поля на заданную глубину по заданным значениям потенциала на поверхности Земли, а также на некоторой высоте от поверхности. Указанный программный комплекс, использованный в исследовательской, производственной и проектно-конструкторской деятельности при исследовании и разработке перспективных для поисков полезных ископаемых районов и при оценке залежей углеводородного сырья, позволяет улучшить методы расчета и анализа источников возмущения полей силы тяжести.

Апробация диссертации. Основные положения и результаты работы докладывались и обсуждались на таких научных конференциях, как: X научная конференция «Дифференциальные уравнения и их приложения в математическом моделировании» с участием зарубежных ученых (г. Саранск, 2012); VII Международная научно-техническая конференция «Аналитические и численные методы моделирования естественно-научных и социальных проблем» (г. Пенза, 2012); VII Международная научно-техническая конференция молодых специалистов, аспирантов и студентов «Математическое и компьютерное моделирование естественно-научных и социальных проблем» (г. Пенза, 2013); VI Международная математическая школа-семинар им. Е. В. Воскресенского «Математическое моделирование, численные методы и комплексы программ» (г. Саранск, 2013); VIII Международная научно-техническая конференция «Аналитические и численные методы моделирования естественно-научных и социальных проблем» (г. Пенза, 2013); 41-я сессия Международного семинара им. Д. Г. Успенского «Вопросы тео-

рий и практики геологической интерпретации гравитационных, магнитных и электрических полей» (г. Екатеринбург, 2014); VIII Международная научно-техническая конференция молодых специалистов, аспирантов и студентов «Математическое и компьютерное моделирование естественно-научных и социальных проблем» (г. Пенза, 2014); Международная научная конференция «Математические методы в технике и технологиях» (г. Тамбов, 2014); XI научная конференция «Дифференциальные уравнения и их приложения в математическом моделировании» с участием зарубежных ученых (г. Саранск, 2014); IX Международная научно-техническая конференция «Аналитические и численные методы моделирования естественно-научных и социальных проблем» (г. Пенза, 2014); IX Международная научно-техническая конференция молодых специалистов, аспирантов и студентов «Математическое и компьютерное моделирование естественно-научных и социальных проблем» (г. Пенза, 2015); ежегодные научные конференции профессорско-преподавательского состава Пензенского государственного университета (2012-2015).

Личный вклад автора. Все основные результаты, представленные в диссертационной работе, сформулированы и получены автором самостоятельно. Работы [1-8, 12, 13] опубликованы в соавторстве с научным руководителем, которому принадлежат формулировка решаемой проблемы и концепция ее решения. В работе [9] автором предложен подход к анализу устойчивости решения систем уравнений с дробными производными. В работе [10] разработанный в цикле работ [3, 4, 13] адаптивный разностный метод продолжения физических полей распространен автором на решение уравнения Гельмгольца. В работах [11, 14] автором проведено обобщение методики исследования устойчивости систем дифференциальных уравнений в частных производных, изложенной в работах [1, 2, 5, 6, 8, 9]. В программном комплексе [15, 16] автором разработаны основные алгоритмы и составлены программные коды.

Публикации. По материалам диссертационного исследования опубликовано 16 работ, включая 4 работы в изданиях, рекомендованных ВАК РФ, и 2 свидетельства о регистрации программ для ЭВМ; 4 работы опубликовано без соавторов.

Структура и объем работы. Диссертация состоит из введения, пяти глав с выводами, заключения, списка литературы и приложения. Общий объем работы составляет 159 страниц, из них 127 страниц основного текста, включая 19 рисунков. Список литературы содержит 224 наименования.

Прямые задачи гравиразведки

Гравиразведка является методом разведочной геофизики, основанным на изучении вариаций поля силы тяжести Земли, вызываемых плотностными неод-нородностями земной коры [15]. Анализ и интерпретация гравитационного поля позволяет изучать распределение неоднородных по плотности масс в земной коре.

Следуя работе [16], дадим краткое описание исторического развития математических моделей геофизики. Современные математические модели геофизики уходят корнями в XVII-XIX века, когда были выведены основные физические законы, составляющие фундамент современной физики, разработаны основы теории гравитационных, магнитных и электрических полей, а математическая физика выделилась в самостоятельную ветвь математики. В период с середины XIX века по 20-е годы XX века возникает и развивается общая геофизика, и вместе с ней появляются такие дисциплины, как гравиметрия, сейсмология, учения о магнитном и электрическом полях Земли и т. д. [17].

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

Одним из наиболее разработанных и широко распространных геофизических методов в настоящее время является гравиметрический метод (гравираз-ведка). Теоретическая база гравиметрического метода начала формироваться в конце XVI века с опытов Г. Галилея, впервые вычислившего значение ускорения свободного падения на поверхности Земли. Изучение сил тяготения затем продолжилось в работах X. Гюйгенса, И. Кеплера и И. Ньютона. Последнему принадлежит честь открытия закона всемирного тяготения, составляющего фундамент современной гравиметрии. Позднее М. В. Ломоносов в 1753 году выдвинул гипотезу о том, что вариации силы тяжести на поверхности Земли обусловлены неоднородностями ее внутренней структуры и сконструировал измерительный прибор для исследования упомянутых вариаций. Дальнейшее интенсивное развитие гравиметрии происходило во многом благодаря всесторонней разработке такими учеными, как, например, К. Гаусс, П.-С. Лаплас, К. Маклорен, М. В. Остроградский, С. Пуассон и Дж. Стоке, методов математической физики и математической теории поля. Это, а также разработка и распространение к концу XIX века эффективных приборов для определения величины силы тяжести, позволили в дальнейшем перейти к систематическому изучению строения земной коры на основе данных гравиметрических измерений [17]. В Госсии активное развитие гравиметрических исследований следует связать в первую очередь с именами Ф. А. Слудского и П. К. Штернберга [18], которые провели большую работу по измерению сил тяжести для изучения гравитационного поля Земли на территории страны. Важным этапом в истории отечественной гравиметрии является исследование под руководством П. П. Лазарева Курской магнитной аномалии, в ходе которого полученные разведочные данные были впервые с успехом применены для решения одной из важнейших задач геофизики - определению основных характеристик рудного тела, являющегося источником аномального гравитационного поля [19]. В советскую эпоху гравиметрическая наука нашла дальнейшее развитие в трудах таких ученых как А. Д. Архангельский [20], Б. К. Балавадзе [21], Ю. И. Блох [22], Г. А. Гам-бурцев [23], В. Н. Дахнов [24], А. И. Заборовский [25], А. А. Михайлов [26], Б. В. Нумеров [27], [28], Н. Н. Парийский [29], Л. В. Сорокин [30], В. И. Старостенко [31], В. Н. Страхов [32], В. В. Федынский [17], Э. Э. Фотиади [33] и др.

Сегодня гравиметрический метод наряду с сейсмическим методом является одним из самых развитых и распространенных методов прикладной геофизики. В круг задач, решаемых в рамках гравиразведки, входит: 1. исследование внутреннего строения земной коры и разведка участков, предположительно богатых полезными ископаемыми с тем, чтобы изучить эти участки более тщательно в том числе и с использованием иных геофизических методов; 2. решение ряда задач геологии и геодезии, в частности, дальнейшее уточнение фигуры Земли; 3. изучение планетарного строения Земли, а также других планет [18].

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

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

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

Замечание 2.1.1. При реализации предлагаемых ниже методов к исследованию аномалий поля силы тяжести необходимо провести предварительную обработку гравитационных данных [192]. Это связано с необходимостью исключить влияние боковых источников и источников, залегающих выше или ниже рассматриваемого тела.

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

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

Предположим вначале, что известны результаты измерения силы тяжести на поверхностях z = 0, z = —hi, z = —hi- В этом случае, связь между источнками гравитационного поля и полями силы тяжести описывается при z 0 уравнениями (2.5). Полагая в уравнении (2.5) z = 0, z = —hi, z = —h2, приходим к системе уравнений

Замечание 2.2.1. В случае, если значения функций f(x), fi(x), f2(x) определены только на сегменте [—/, /] и их граничные значения /(±/), /i(±/), f2(±l) не все равны нулю, то эти функции могут быть продолжены на числовую ось следующим образом. Пусть /(/) ф 0. Функцию /(ж) можно продолжить на интервал (/, оо) функцией д+(х) = f(l)e x l\ Если известны значения /(/) ф О и / (/) ф 0, то в качестве продолжения можно взять функцию

Замечание 2.2.2. Как будет показано ниже при рассмотрении модельного примера, в выражении Н сомножители, зависящие от переменной ш, сокращаются, и Н является константой. В случае, если значения F( x ), Р\{ш), F2{uS) в аналитическом виде неизвестны и получены в результате приближенного вычисления преобразования Фурье функций f(x), fi(x), f2(x), то величина \Н(ш) — Н\ определяется погрешностью вычисления преобразования Фурье. В этом случае представляется естественным вычислить Н(ш) на некоторой последовательности точек Wk = —В + Ц -, к = О, N, где [—В, В] - сегмент, в котором вычисляются преобразования F( x ), Р\{ш), F2{uS), и положить

Рассмотрим теперь случай, когда на поверхности z = 0 известны значения поля силы тяжести и его первой и второй производных по переменной Z.

Замечание 2.2.4. Значения производных могут быть вычислены по разностным схемам с использованием формулы Пуассона для получения разгонных значений.

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

Применяя к (2.24) обратное преобразование Фурье, находим функции wi(x), w2(x) и параметр Я. Замечание 2.2.5. В случае, если функции f(x), fi(x), /2(00) заданы не в аналитическом виде, и их преобразования Фурье не вычисляются точно, то применение обратного преобразования Фурье к функциям WI(UJ), И ( ), выраженным формулами (2.18), (2.24) может привести к большим погрешностям и, в ряде случаев, может быть неосуществимо. В этом случае необходимо применение приближенных методов. Соответствующие приближенные методы будут изложены ниже применительно к трехмерным задачам.

Теперь обратимся к построению аналитического метода решения задачи одновременного восстановления плотности, формы и глубины залегания грави-тирующего тела в трехмерной контактной задаче. Вначале рассмотрим случай, когда значения поля силы тяжести известны на поверхностях z = 0, z = — h\, z = —hi и обозначим функции /(ж, 2/,0), f(x,y, —/її), f(x,y, —/12) соответственно через f(x,у), /і(ж,2/), /2(2,3/).

Построение вычислительных схем аппроксимации тепловых полей

Для простоты обозначений положим А = 0, = 1, Т = 1, а = 2. Обозначим через Г2о множество точек (x,t): удовлетворяющих неравенствам {0 x l,0 t 2_дг}. Через Qk обозначим множество точек, удовлетворяющих неравенствам 0 х 1, \г- fw \ Каждую область Г2& покроем прямоугольниками Д , у которых длина ребер, параллельных оси ОХ, равна hk, к = О, Х-1. Положим h0 = r,hk= Nl/\ (f). A; = 1,X-1.

Пусть /(/:) Є 0[ 2, &]. Обозначим через Ls(f, [a, 6]) интерполяционный полином, построенный по s узлам, расположенным в сегменте [a, Ъ]. В качестве узлов интерполяции можно взять равноотстоящие точки, а также узлы ортогональных многочленов, линейно отображенных с сегмента [—1,1] на [a, Ъ]. В случае, когда необходимо построить непрерывное приближение, достаточно воспользоваться приемом, описанным в [203], [204].

Пусть f(t,x) Є C(D), D = [a, b; c, d]. Через Ls s(f,D) обозначим интерполяционный полином, действующий по формуле LSyS(f,D) = = L (Lf (/, [с, d\), [a, 6])), т.е. вначале оператор Щ действует на функцию f(t, х) в сегменте [с, d] по переменной ж, а затем оператор L\ действует на функцию L (f, [с, d]) по переменной t в сегменте [а, Ъ].

В каждой области А , к = 0, X — 1, функцию u(t,x) будем аппроксимировать интерполяционным полиномом Ls s(u,Af), где s = [aN] + 1. Сплайн, составленный из полиномов LS;S(w, Af), обозначим через S(t,x).

Оценим точность аппроксимации функций из множества РГ;7;а(Г2,1) сплайном S(t,x). Оценку начнем с областей Д9. Рассмотрим для определенности область до = [0, 2 N; 0, ho]. Так как функция u(t, х) по переменной t удовлетворяет условию Ге ль дера с показателем а, то

Так как функция u(t, х) по переменной х имеет ограниченные единицей производные до г-го порядка включительно, то

Отметим, что построенный выше сплайн S(t,x) может иметь линии разрывов на границах областей А , к = О, N — 1. По аналогии с рассуждениями, приведенными в работах [191], [218], можно построить непрерывный в области Q сплайн, имеющий погрешность (3.26).

Отсюда вытекает следующая теорема. Теорема 3.3.2. Пусть Q = [0, 1; 0, Т]. Справедлива оценка Замечание 3.3.2. Здесь положено а = 2, так как эта константа возникает при решении задачи Коши (3.17)-(3.18) Замечание 3.3.3. Используя метод оценки снизу поперечников Бабенко компактного множества РГ;7(Г2, М) (см. [204]) можно показать, что

Теперь на основе принадлежности тепловых полей введенному классу РГ;7;а(Г2, М, а) перейдем к построению адаптивных разностных методов решения задачи (3.17)-(3.18). Введем декартову систему координат Oxt. Пусть в области

Замечание 3.3.4. Заметим, что на практике значения функции u(t,x) при t = 0 известны не во всем пространстве R, ав некоторой конечной области G, которую с целью упрощения дальнейших выкладок будем считать сегментом:

Опишем неравномерную сетку узлов, используемую для построения предлагаемого разностного метода. Пусть начальное условие щ(х) задано в виде множества значений мо,/г}, Щ,к = щ(%к) в узлах XQ = (О, —A + kti), к = О, N, h = построенной в квадрате G равномерной сетки узлов. Будем вычислять значения теплового поля в области

Зафиксируем шаг г = - по переменной t) где М 0 целое число. Для аппроксимации уравнения (3.17) в точке (tj,Xk), tj = jr, Xk = —A + kh используется следующий четырехточечный шаблон.

Из рисунка видно, что область, в которой рассчитываются значения теплового поля, имеет трапециедальную форму. При этом если на J -OM слое (j ф 0) по переменной t имеется N значений Ujk, то на {j + 1)-ом слое будем иметь [уJ — 1 значений. При этом если на начальном слое по переменной t значения числового поля заданы в интервале [—-В, В]: то на слое j ф 0 значения будут вычислены в интервале [—В + Sjh(0), В — Sjh(0) \, где Sj - сумма первых j членов геометрической прогрессии с первым членом Ъ\ = 1 и знаменателем 9 = 2:

Вычислительный эксперимент показывает, что при неудачном выборе числа N разбиений промежутка [—-В, В]: на котором задаются начальные значения Uok теплового поля, вышеописанные рассуждения в общем случае не приводят к построению нужной разностной сетки; в частности, сетка может оказаться несимметричной относительно прямой х = 0, либо на некоторых слоях по переменной х могут возникнуть системы неравноотстоящих узлов.

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

Зафиксируем число М слоев по переменной t и требуемое число узлов N = = N(M) на слое j = М по переменной t (т. е. при t = t ). Тогда начальные шаги т(0) и /г(0) по переменным t и h соответственно вычисляются по формулам

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

Замечание 3.3.5. Поскольку реализация описанной вычислительной схемы требует задания начальных значений на интервале [—-В, В] D [— 4, А]: то требуется доопределить функцию щ{х) на множестве [—-В, В] \ [—А, А] значений аргумента х. Сделать это можно, в частности, одним из следующих способов:

Описание алгоритма и программы определения формы и плотности гравитирующего тела при решении задачи потенциала

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

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

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

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

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

В рамках программного комплекса пользователь определяет решаемую задачу (адаптивный алгоритм восстановления потенциальных или тепловых полей) и размерность решаемой задачи (двух- или трехмерный случай). Поскольку для работы программы требуется массив значений потенциального поля в узлах равномерной сетки, то для работы программы эти значения необходимо загрузить, нажав кнопку "Данные..." и выбрав файл с данными. Файл с данными должен иметь расширение .txt и последовательно содержать границу А содержащего известные значения потенциала интервала [—Л, Л], число разбиений интервала [—А, А] узлами равномерно сетки, а также (по одному на каждой строке) значения потенциала в узлах сетки на поверхности Земли. Значение полуинтервала определяет границу В интервала [— , ], на котором в программе фиксируется сетка с начальными значениями поля; при этом требуется, чтобы выполнялось неравенство В А] если В А: то известные значения потенциала доопределяются на множество [—В, В] \ [—А, А] нулями. Кроме того, пользователем фиксируется необходимое число расчетных слоев по переменной z (переменной t): а также значение (глубина), до которого поле должно быть восстановлено. Структура используемой адаптивной сетки определяется следующим образом: шаги адаптивной сетки увеличивается через каждые "Период увеличения" слоев в "Кратность увеличения" раз. Результатом работы программы является текстовый файл, содержащий последовательность значений восстановленного потенциального поля в узлах адаптивной сетки, а также (по выбору пользователя) двух- или трехмерный график восстановленного потенциального поля на последнем расчетном слое. Пример такого графика показан на рис. 5.2.

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

Из соображений простоты алгоритм излагается для двухмерного случая при заданных значениях потенциального поля на повер-хости Земли, а также на высоте 5 0 от поверхности. Однако этот алгоритм может быть с успехом распространен как на трехмерный случай, так и на случай заданных производных потенциального поля по переменной z на поверхности Земли.