Спектрально-согласованные сетки для моделирования волновых процессов тема автореферата и диссертации по математике, 01.01.07 ВАК РФ
Лисица, Вадим Викторович
АВТОР
|
||||
кандидата физико-математических наук
УЧЕНАЯ СТЕПЕНЬ
|
||||
Новосибирск
МЕСТО ЗАЩИТЫ
|
||||
2007
ГОД ЗАЩИТЫ
|
|
01.01.07
КОД ВАК РФ
|
||
|
На правах рукописи
ЛИСИЦА Вадим Викторович
СПЕКТРАЛЬНО-СОГЛАСОВАННЫЕ СЕТКИ ДЛЯ МОДЕЛИРОВАНИЯ ВОЛНОВЫХ ПРОЦЕССОВ
01.01.07 - вычислительная математика
АВТОРЕФЕРАТ диссертации на соискание ученой степени кандидата физико-математических наук
НОВОСИБИРСК 2007
003053150
Работа выполнена в Институте нефтегазовой геологии и геофизики им. A.A. Трофимука Сибирского отделения Российской академии наук и Новосибирском государственном университете
Научный руководитель:
кандидат физико-математических наук, доцент В.И.Костин
Официальные оппоненты:
доктор физико-математических наук В.Д. Лисейкин
кандидат физико-математических наук Г.В. Решетова
Ведущая организация:
Институт прикладной математики им. М.В Келдыша РАН (г. Москва)
Защита состоится 21 февраля 2007г. в 15 час. на заседании диссертационного совета Д 003.061.01 при Институте вычислительной математики и математической геофизики СО РАН, в конференц-зале.
Адрес: 630090, Новосибирск-90, пр-т Ак. Лаврентьева, 6а Факс: (383) 330 87 83
С диссертацией можно ознакомиться в библиотеке ИВМиМГ СО РАН
Автореферат разослан 20 января 2007г.
ОЫЦЛЯ ХЛРЛК'1 СГИС I IПСА РАБОТ Ы
Объектом исследовании в данной работе являются сетки для конечно-разностных схем на предмет их оптимизации для моделирования волновых процессов в упругих средах, в том числе анизотропных.
Актуальность. При моделировании волновых процессов достаточно распространенным требованием является знание волнового поля не внутри расчетной области, а лишь на некоторой ее границе. Наиболее показательные в этом смысле геофизические эксперименты Применение стандартных конечно-разностных методов для численного решения таких задач приводит к гигантским требованиям на вычислительные ресурсы и временным затратам на реализацию.
Причина этого явления в том, что отображение совокупности параметров среды и особенности волновых процессов в среде в "данные наблюдений" (сейсмограммы) происходит с катастрофической потерей информации В частности, в трехмерных динамических задачах число параметров, входящих в описание задачи, оценивается величиной 0(И4) , где N пропорционально линейному размеру области, в которой развивается волновой процесс и обратно пропорционально шагу сетки Размерность пространства образов - число отсчетов на сейсмограммах - имеет порядок
не вышеО(А^2). По-видимому, даже эта оценка является завышенной, если принять во внимание гладкость зондирующего сигнала Однако уже этого сопоставления достаточно, чтобы представить насколько сильна потеря информативности при действии оператора, отображающего параметры среды в отсчеты на сейсмограммах.
Стандартные способы аппроксимации этого оператора посредством построения разностной схемы для начально-краевой задачи, соответствующей выбранной модели волновых процессов, никак не учитывают этой потери информативности В определенном смысле оператор оказывается слишком плохо аппроксимирован - каждая значащая цифра в ответе появляется как результат выполнения огромного числа арифметических операций над входными данными, лежащими в пространстве очень большой размерности Эту размерность нельзя уменьшить, уменьшая так как это приведет к увеличению локальной ошибки в промежуточных и окончательных результатах Повышение порядка аппроксимации разностной схемы хотя п может приводить к некоторым улучшениям, но они не принципиальны
Аналогичная ситуация возникает и при реализации "неотражающих" граничных условий, использование которых также необходимо для качественного моделирования геофизических экспериментов. Наиболее употребимым является Идеально Согласованный Слой (англ. - PML: Perfectly Matched Layer) Основной идеей метода является окружение расчетной области специально сконструированным слоем, таким, что волна проходит через границу без отражений и затухает внутри этого слоя. Несложно видеть, что PML является исключительно лишь вспомогательной конструкцией и необходимо лишь знание решения на границе PML - расчетная область, однако стандартные методы реализации эгого подхода, не учитывающие этого требования, приводят к тому, что затраты на вычисление решения внутри PML доходят до 75% от общих затрат на моделирование эксперимента.
Иных результатов удается достичь, если, вводя определенную свободу в строение схемы во внутренних точках, пытаться оптимизировать погрешность в отдельных точках. В задачах, упоминавшихся выше, возможна экспоненциальная и даже сверхэкспоненциальная сходимость.
Практическая реализация идей, лежащих в основе упомянутого подхода, наталкивается на ряд трудностей, связанных с необходимостью обеспечения высокой точности в решении целого ряда вспомогательных задач. Понятно, что для достижения оптимального результата, погрешности во всех промежуточных задачах должны быть минимизированы, в частности важным аспектом является построение количественной оценки скорости сходимости метода для волновых задач, что представлено в данной работе. Использование спектрально-согласованных сеток для моделирования волновых процессов в упругих средах является актуальной задачей и, как следствие, требует умения их построения и теоретического обоснования возможности их применения для таких задач.
Цель исследований - повышение эффективности расчета волновых полей в упругих средах за счет использования спектрально-согласованных сеток для конечно-разностных схем, ориентированных на аппроксимацию решения лишь в заданных, определяемых условиями эксперимента, точках расчетной области.
Научные задачи: построение спектрально-согласованных сеток для волнового уравнения с переменными коэффициентами и получение количественной оценки скорости сходимости численного решения для волновых задач, полученного с помощью данного подхода; построение идеально согласованного слоя, основанного на применении спектрально-согласованных сеток, для уравнений теории упругости в изотропном слу-
чае; моделирование волновых процессов в анизотропных упругих средах с использованием спектрально-согласованных сеток.
Основные этапы исследовании:
1. получить количественную оценку точности расчета волновых полей с использованием конечно-разностных схем на спектрально-согласованных сетках, основанную на теории сходимости аппроксимаций Паде-Чебышева;
2. построить оптимальный идеально согласованный слой, основанный на использовании спектрально-согласованных сеток для схемы Вирье для системы уравнений динамической теории упругости в случае изотропной среды;
3. разработать алгоритм построения спектрально-согласованных сеток для схемы Лебедева для системы уравнений динамической теории упругости в случае анизотропной среды,
4. провести серию численных экспериментов для верификации полученных теоретических результатов.
Научные методы исследований.
Теоретической основой решения поставленной научной задачи являются
• современная теория аппроксимаций - рациональные аппроксимации, аппроксимаци Паде, Паде-Чебышева;
• численные методы линейной алгебры - решение прямой и обратной спектральной задач для симметричных трехдиагональных матриц с использованием последовательностей Штурма, алгоритма Ланцоша, метод спектрального исчерпывания, основанный на цепочках двумерных вращений Якоби,
• результаты современной теории конечно-разностных методов для решения уравнений математической физики - схемы Вирьё, схемы Лебедева, схемы на повернутых сетках;
• спектральная теория линейных дифференциальных операторов; современная теория упругости
Рафабо1анные алгортмы и программы использовались для проведения численного моделирования волновых полей в средах различной степени сложности Результат!,1 подвергались сравнительному анализу с результатами, полученными с использованием других методов - аналитическими и конечно-разностными на равномерных сетках
Защищаемые научные результаты:
1. Доказана экспоненциальная скорость сходимости конечно-разностного решения волнового уравнения с переменными коэффициентами вычисленного с помощью спектрально-согласованных сеток.
2. Разработан алгоритм построения оптимального идеально согласованного слоя, основанный на применении спектрально-согласованных сеток, для системы уравнений динамической теории упругости в случае изотропной среды.
3. Разработан и реализован алгоритм построения спектрально-согласованных сеток для решения системы динамической теории упругости для случая анизотропной упругой среды.
Новизна работы. Личный вклад.
• Получена оценка скорости сходимости конечно-разностного решения при использовании спектрально-согласованных сеток для волнового уравнения с переменными коэффициентами, основанная на теории сходимости аппроксимаций Паде-Чебышева для мероморфных функций.
• Построен и реализован оптимальный идеально согласованный слой для изотропной упругой среды с использованием модификации схемы Вирьё на спектрально-согласованных сетках.
• Разработан, теоретически обоснован и реализован алгоритм построения спектрально-согласованных сеток для схемы Лебедева для системы уравнений динамической теории упругости в случае линейных анизотропных сред, и проведены серии тестов для верификации полученных результатов.
Теоретическая и практическая значимость результатов.
Получена количественная оценка сходимости метода оптимальных сеток для волновых задач, в том числе с переменными, включая разрывные, коэффициентами. Построен оптимальный идеально согласованный слой для уравнений теории упругости в случае изотропной среды, позволяющий существенно уменьшить вычислительные затраты при моделировании волновых процессов в неограниченных областях. Разработан алгоритм построения спектрально-согласованных сеток для анизотропных упругих задач, позволяющий эффективно моделировать распространение
волн в таких средах и расчитывать волновые ноля в заданных точках среды (приемниках)
Апробация работы и публикации.
Основные положения и результаты докладывались на седьмой международной конференции, посвященной математическим и численным аспектам распространения волн "Waves 2005" (Провидэнс США, 2005), ежегодном форуме Общества Индустриальной и Прикладной Математики "SIAM Annual Meeting 2006" (Бостон, США, 2006), 12 Международном симпозиуме по сейсмической анизотропии 12IWSA (Пекин, Китай, 2006); II Всероссийской конференции, посвященной памяти академика А Ф Сидорова «Актуальные проблемы прикладной матемашкп и механики» (Абрау-Дюрсо, 2004), V Международной научно-пракшческой геоло-го-геофи шческой конференции-конкурсе молодых ученых и спецмалиаов "Геофизика 2005"(Санкт-Петербург, 2005); Второй Сибирской международной конференции молодых ученых по наукам о Земле (Новосибирск, 2004)
Результаты исследований по теме диссертации изложены в 6 опубликованных работах Из них одна статья в Сибирском Журнале Вычислительной Математики, 2 работы - это материалы Международных конференции, 3 работы - материалы российских международных конференций
Структура и объем работы.
Диссертация состоит из введения, 4 глав, заключения к списка литературы из 61 наименования Работа содержит 90 страницу основного текст
благодарности.
Диссертация выполнена в лаборатории вычислительных методов геофизики Института нефтегазовой геологии и геофизики СО РАН Автор выражает искреннюю признательность научному руководителю, к ф -м н В И Костину за всестороннюю поддержку и постоянное внимание, а 1ак-же В А Чеверде. В Л Друскину п С В Гольдину за регулярные обсуждения
ОСНОВНОЕ СОДЕРЖАНИЕ РАБОТЫ.
Во введении обосновывается актуальность, научная новизна, важное и» полученных в диссертации результатов. Формулируются основные положения, выносимые на защиту
В первой I лаве приведен обзор существующих на ссюдня конеч-но-рашосшых меюдов расчета волновых нолей в упрушх средах и мею-
дов ограничения расчетной области при таком моделировании. Дано краткое описание этих методов и краткий анализ их достоинств и недостатков. Кроме того, указано место развиваемого в диссертации подхода среди многообразия других методов
Во второй главе приводится алгоритм построения спектрально-согласованных сеток для волнового уравнения с переменными коэффициентами Рассматривается одномерное волновое уравнение с переменными - кусочно-непрерывными ограниченными коэффициентами в ограниченной по пространству области.
В разделе 2 1 ставится задача построения спектрально-согласованных сеток и приводится краткое описание алгоритма решения. Для построения сетки удобно перейти в спектральную область, применив преобразование Фурье по времени. В этом случае задача сводится к построению сеток для соответствующей задачи Штурма-Л иувилля
а{х)1ь[х{х)%) = °
Х(х)~г(х = */) = -и и(х = хг) = О ах
(1)
Нас интересует решениеи{х = X/) . Для данной задачи оно представляется в виде:
«(*,) =/(л) = £-
<р,
~ л * ' (2)
где - собственные числа, а <р} - нормированные значения (при х = х)) квадратов собственных функций задачи (1).
С другой стороны, если ввести некоторую сетку, где и определена
в точках х , при ] = 1,..., к + 1, а ее производная в точках Х/, при У = 0,..., к , можно построить конечно-разностную схему
1
м/+1
и, -«,_,
У,
У,-1
Лг/; = О,
м, - и.
\ у I
1
Лм, = —, ик+] = 0. У\
Веса схемы естественно определить по правилу:
тогда конечно-разностную схему можно записать в виде: Ни - Л//7 = —
Ж"
(4)
где С, - первый базисный вектор, I - единичная матрица, Н - симметричная трехдиагональная матрица с положительными внедиагональными коэффициентами, однозначно определяемая параметрами (3) Решение задачи (4) можно представить в виде.
где у1 - нормированные квадраты первых компонент собственных векторов, - собственные числа матрицы Н.
Задача построения спектрально-согласованной сетки для решения уравнения (I) с использованием схемы (4) сводится к решению следующей задачи минимизации
В рассматриваемых задачах спектральный интервал [Л^О], зависящий от временной частоты, определяется носителем Фурье-образа функции источника и предполагается конечным
Суммируя представленные выше рассуждения, алгоритм построения сетки может быть разбит на следующие подзадачи-
• построение точного решения /(Л), представленного формулой (2) для задачи (I);
• решение задачи (6) с использованием теории рациональных аппроксимаций с получением функции /„'„(Л) вида (5);
• восстановление матрицы системы (4) посредством решения обратной спектральной задачи для симметричных трехдиагональных
(5)
(6)
матриц;
• отыскание весов у и у из матрицы Н, восстановление узлов и шагов сетки с использованием (3)
В разделе 2.2 приводится решение задачи (6) с использованием аппроксимаций Паде-Чебышева. В разделе 2.2.1 излагается алгоритм построения такой аппроксимации. В разделе 2.2.2 формулируется и доказывается ключевая теорема об оценке скорости сходимости аппроксимации.
Теорема 1. При использовании аппроксимации Паде-Чебышева для решения задачи (б) верна оценка
тах /(Л) -/* (Л) < С
-с2(4(*-И)+1)
Ле|Д,.0]1....."Р' 'I " 4(*-и) + 1
где к - оби/ее число сетки (полюсов функции (5)), П - число полюсов функции (2), лежащих внутри спектрального интервала, константы С | и С , -
положительные не зависящие от к
Следует заметить, что построение точного решения не всегда возможно в явном виде (2), по этой причине в качестве точного решения для волновых задач разумным является рассмотрение в качестве точного решения - конечно-разностное решение на некоторой достаточно мелкой
Н
равномерной сетке, т.е. У(Л) = / (Л) . В этом случае во временной области получим:
е-( ,(4Ц-и)+|) £
4{к -п) + \
где N - число точек исходной равномерной сетки.
В разделе 2.3 представлено оригинальное решение обратной спектральной задачи для симметричных трехдиагональных матриц с положительными внедиагонапьными элементами или матриц Якоби. Обратной спектральной задачей мы называем задачу восстановления матрицы Якоби по набору всех ее собственных чисел и первых компонент собственных векторов. Предлагаемый в работе алгоритм решения основывается на цепочках двумерных вращений Якоби и является, фактически, обратным к алгоритму спектрального исчерпывания для матриц Якоби.
Раздел 2.4 содержит описание конечно-разностной схемы во временной области. Данная схема нспользуося с равномерным шаюн но времени и спектрально-согласованной сеткой по пространству Численные эксперименты представлены в разделе 2.5. Целью экспериментов было
вычисление решения волнового уравнения на заданной границе с использованием спектрально-согласованных сеток и сравнение результатов с расчетами на равномерных сетках Эксперименты иллюстрируют, что использование спектрально-согласованных сеток с 3 точками на длину волны обеспечивают ту же точность решения, что и равномерные с 40 точками на длину волны Выбор значения 40 точек на длину волны для схемы второго порядка аппроксимации на равномерных сетках не случаен, а определяется ее дисперсионными свойствами и является наиболее употребимым для таких схем. Примеры одномерных сеток для волнового уравнения с линейным возрастанием скорости. с(д:) = X представлены на рис. I, по оси абсцисс представлена координата узла по пространству.
О К
* О *
О х о > О * О
*ЮхО * О к О * О
*00«0 «0*0*0* о» о
1 1.2 1.4 1.6 1 >5 2
Рис I
В третье» главе представлено построение оптимального идеально согласованного слоя (РМЬ) для системы уравнений динамической теории упругости в случае изотропной среды Идеально согласованные слои используются для ограничения расчетной области без отражений от внешних границ. Стандартным подходом к построению таких слоев является применение комплексной замены переменных в спектральной области
~ ( Р\
внутри этого слоя Л' = (X Н--Л", здесь предполагается, что РМЬ сгро-
V ко)
ится по переменной х при х > 0 Несложно видеть, что такая замена пе-
ременных переводит распространяющиеся волновые решения в экспоненциально затухающие. В основе оптимального РМЬ лежит использование
~ х
мнимой замены переменных: X =— в области Х> I). Раздел 3.1 по-
1С0
священ постановки задачи построения оптимального РМЬ для системы уравнений динамической теории упругости в случае изотропной среды, что включает в себя как применение замены переменных к системе уравнений, так и постановку граничных условий.
В разделе 3.2 приводится алгоритм построения спектрально-согласованной сетки для системы уравнений, описывающей оптимальный РМЬ. В разделе 3.2.1 представлено построение точного решения дифференциальной задачи, соответствующей оптимальному идеально согласованному слою. Для построения спектрально-согласованных сеток нам необходимо знание решения исключительно лишь при X = 0. В следующем разделе приводится построение решения соответствующей конечно-разностной задачи, как функции шагов сетки В данном случае нас вновь интересуют лишь первые компоненты искомых векторных переменных. Используя полученные представления решений, в разделе 3.2.3 формулируется теорема.
Теорема 2. Решения дифференциальной и конечно-разностной задач для оптимального РМЬ для системы уравнений динамической теории упругости представляются в виде
Ч(оЛ , /МО)
(«л),
= {со, к_, И/, И))
("г) о
где и = (мд,Мг)/ - вектор скорости смещений, СГ(1, СТ.,, 0\_ - компоненты тензора напряжений, СО -временная частота, а к. - соответствующая компонента волнового вектора Функции
и
((О, к.,И / ,И ) имеют представление
Е{со,кг) = Я> {со,к:)/(А,) + Я"(О),кг )/(Л„),
,И,) = Л'\со,к2)/*(Л,,И/,И,) + Я"{соЛ__ )/*(Л,,И,,И,),
где матричные функции Я' {&>, к: ) и Я'' {(0, к_ ) - совпадают для обеих задач и равномерно ограниченны, как функции своих аргументов Скачяр-
пая функция / (А) = .— . а / * (Л) задается формулой (5) При том л/Л '
Л
/ 1- е
V СО' V?
и верпа оценка ошибки 0)
\ Г , чА \
(»V):
I /
».-(0)
по числу
при на-
При использовании подходящей рациональной аппроксимации квадратного корня, возможно достижение экспоненциальной сходимости
узлов сетки, тах |/(Л) - (к (Л)| < Се'( Далее,
личин построенной рациональной функции, удовлетворяющей приведенной оценке, используя результаты главы 2 можно построить сетку для поставленной задачи Примеры одномерных сеток для оптимального РМЬ представлено на рис. 2 По оси ординат указано число точек сетки, по оси абсцисс - координата узла сетки по пространству.
О х С >
Рис 2
В разделе 3.4 представлен численный эксперимент расчета волнового поля в некоторой однородной среде в ограниченной области при наличии и отсутствии окружающего его оптимального РМЬ. На рис. 3 и 4 представлены волновые поля на границе расчетной области (сейсмограммы) при отсутствии РМЬ и его наличии По оси абсцисс отложен номер
приемника, tío ори ординат номер слоя по времени. Как несложно видеть, отсутствие ¡'Ml. приводит к появлению отражений от граним расчетной области, il то время как использование оптимального PML на трех точках позволяет эффективно бороться с этими отражениями и на рис. 4 й'идны лишь прямые волны.
Рис.3
Рис.4
Более детальный анализ, проведенный в разделе 3,4, показывает, что оптимальный ]JML с 3 точками позволяет уменьшить отражения до 1% от тех, которые появляются при отсутствии РМI Мри использовании оптимального PML С четырьмя т очками отражения составляют порядка О, I %.
Для достижения тех же результатов с использованием стандартной реализации РМЬ необходимо использовать от 30 до 50 точек.
В четвёртой главе рассматривается задача построения спек-трапьно-согласованной сетки для системы динамической теории упругости в случае произвольной линейной анизотропной среды в двумерной постановке
В разделе 4 1 приводится постановка задачи. Для анизотропных задач, прежде всего, необходимо подобрать подходящую конечно-разностную схему После чего поставить задачу построения спектрально-согласованной сетки для выбранной схемы
Выбору и исследованию схемы посвящены разделы 4 2 и 4.3. В разделе 4.2 приводится и обосновывается схема Лебедева, в разделе 4 3 -хорошо известная схема на повернутых сетках Далее проводится их детальный сравнительный анализ: исследование устойчивости - раздел 4 3 2, дисперсионный анализ - 4 3 3 и сравнение затрат на реализацию - раздел 4.3 4. В результате этих исследований показано, что использование схемы Лебедева эффективнее с вычислительной точки зрения, чем схемы на повернутых сетках.
Раздел 4 4 посвящен построению спектрально-согласованных сеток для схемы Лебедева для системы уравнений динамической теории упругости в случае анизотропной среды. Здесь формулируется и доказывается следующая ключевая теорема
Теорема 3. При рассмотрении системы уравнений упругости в спектральной области (х,СО,к_) и испотьзовании произвольной сетки по направлению X , решение соответству юи/ей конечно-разностной задачи представляется в виде
Ы )> = К К )/< (Л„,) + /V {со, к. )/* (Л¥Л )0,
где (и'7 )0 - входные параметры, зависящие от граничных усчовий дифференциальной задачи, (и'*), - вектор искомых переменных, ГцГ {со, к ^ ) и Т7^ {со, к_ ) - .матрицы, равномерно ограниченные, как функции своих аргументов, в области изменения аргументов, рациональная функция j (Л) задается формулой (5) При том
Л„.. ,л = , -cos2 (а), а = arceos
и Ач1, s е [Л, max(v~2(а)),0], где А, - определяется функцией источника Согласно результатам главы 2 верпа оценка
. <,(4Ц_„)+|)
и >г - к )Г'|| * с, |/>|+к. Ib )„II4(^_/7)+i
Данная теорема обосновывает возможность построения спектрально-согласованных сеток схемы Лебедева для рассматриваемой системы уравнений и доказывает более чем экспоненциальную скорость сходимости при их применении. В силу рассуждений, приведенных в главе 2, в качестве точного решения может быть использовано конечно-разностное решение на некоторой достаточно мелкой равномерной сетке.
Раздел 4.6 содержит численные эксперименты расчета волновых полей в анизотропных средах. На рис. 5 представлено волновое поле в момент времени Í = 0.2 с указанием размеров расчетной области Физический размер представлен на осях. Размер области в точках указан на самом рисунке. В области слева всюду использовались равномерные сегки, справа - в половине области вводились спектрально-согласованные сетки по направлению X. Видно, что в данном эксперименте число точек спектрально-согласованной сетки в шесть раз меньше, чем равномерной.
На рис. 6 представлена сейсмограмма - поле, зафиксированное
при х = х^ внутри левых подобластей, для обоих экспериментов. По оси
абсцисс - расчетное время, по оси ординат - номер точки по вертикали. Относительная норма ошибки на данной сейсмограмме составляет порядка 1%.
Из представленных экспериментов можно заключить, что использование спектрально-согласованных сеток для схемы Лебедева для системы уравнений динамической теории упругости в случае анизотропной среды является эффективным инструментом, позволяющим существенно, более чем в 5 раз, сократить вычислительные затраты на проведение моделирования
о
О 2 О.-4 Об
0 8 "1
1.2
1 -1
1.6 1 8
1=0,2
N =425 V N =^25 *
S p
II
Г* ^
,вС
t=0 2 N ~425 N =70
0.5
О 5
Рис.5
s . > . .J --i.'-..,,
ЩМтш
t
о
0.5
1.5
ЗАКЛЮЧЕНИЕ.
Основным результатом работы является развитие, реализация и тестирование подхода спектрально-согласованных сеток для конечно-разностных схем, который позволяет существенно сократить вычислительные затраты при моделировании волновых процессов. Его отличительной особенностью является локальная минимизация ошибки численного решения в заданных точках расчетной области связанных с нуждами конкретного эксперимента (модели). В данной работе показано, что использование таких сеток для численного решения волновых задач позволяет добиться более чем экспоненциальной сходимости, что в свою очередь обеспечивает существенное сокращение, более чем в пять раз, вычислительных затрат на расчет решения, в сравнении с равномерной сеткой Следует также отметить, что затраты на построение сетки пренебрежимо малы в сравнении с затратами на моделирование самого процесса (эксперимента), что позволяет вычислять их каждый раз непосредственно перед моделированием.
В работе представлено обоснование возможности использования подхода спектрально-согласованных сеток для расчета волновых полей в упругих, в том числе анизотропных, средах. Задача моделирования распространения упругих волн в анизотропных средах с использованием конечно-разностных схем на сдвинутых сетках сама по себе является весьма актуальной. Как показано в работе, использование схем Лебедева является эффективным методом решения этой задачи. Более того, в работе представлена оптимизация сеток для схемы Лебедева, обеспечивающая сверх-экспоненциапьную сходимость. Приведенные эксперименты демонстрируют работоспособность и высокую эффективность применения схемы Лебедева на спектрально-согласованных сетках для решения этой задачи.
Построение оптимального идеально согласованного слоя для системы уравнений динамической теории упругости является, пожалуй, наиболее значимым результатом. Представленный в работе оптимальный РМЬ, основанный на методе спектрально-согласованных сеток позволяет сократить число точек, необходимое для достижения нужной точности более чем в десять раз, 3-4 точки вместо 30-50 - в стандартной реализации РМЬ. Если учесть тот факт, что идеально согласованный слой окружает всю расчетную область, то при стандартной реализации он требует более 50% от общего числа вычислительных ресурсов, использование оптимального РМЬ сокращает эти затраты до 5-10%
Основные результаты диссертации опубликованы в следую-
щих работах:
1. Лисица В.В. Оптимальные сетки для численного решения волнового уравнения с переменными коэффициентами.// Сибирский журнал вычислительной математики, 2005, т.8(3), с.219-229
2. Лисица В.В. Оптимальные сетки для численного решения волнового уравнения.// Тезисы докладов II Всероссийской конференции, посвященной памяти академика А.Ф. Сидорова "Актуальные проблемы прикладной математики и механики", 2004, с. 71-72.
3. Лисица В.В. Оптимальный идеально согласованный слой для моделирования волновых геофизических полей в неограниченных областях.// Тезисы докладов. V международная научно-практическая геолого-геофизическая конференция молодых ученых и специалистов "Геофизика 2005", 2005, с.163-166.
4. Lisitsa V. Optimal Unsplit Perfectly Matched Layer (PML) for 2D Elasticity.// Proceedings of the 7th International Conference on Mathematical and Numerical Aspects of Wave Propagation "Waves 2005", 2005, pp.97-99.
5. Lisitsa V. Lebedev Schmes for Simulation of Waves' propagation in Anisotropic Elastic Media.// Expanded Abstracts of The 12th International Workshop on Seismic Anisotropy, 2006, pp.122-124.
Технический редактор О М Вараксина Подписано к печати 22 12 2006 Бумага 60x84/16 Бумага офсет № 1 Гарнитура «Тайме» Печать офсетная _Печ л 0,9 Тираж 100 Заказ № 11_
НП АИ «Гео» 630090, Новосибирск, пр-т Ак Коптюга, 3
ВВЕДЕНИЕ
Глава 1. ИЗУЧЕННОСТЬ ВОПРОСА.
Глава 2. ПОСТРОЕНИЕ СПЕКТРАЛЬНО-СОГЛАСОВАННЫХ СЕТОК ДЛЯ ВОЛНОВОГО УРАВНЕНИЯ С ПЕРЕМЕННЫМИ КОЭФФИЦИЕНТАМИ
2.1 Постановка задачи.
2 1.1 Аналитическое решение.
2 12 Конечно-рашостная задача.
21.3 Построение решения конечно-разностной задачи
2 14 Построение рациональной аппроксимации.
2 15 Построение сетки.
2 2 Аппроксимации Паде-Чебышева
2 2 1 Построение аппроксимаций
2 2 2 Порядок аппроксимации
2 3 Обратная спектральная задача.
2 3 1 Алгоритм, основанный на вращениях Якоби.
2 4 Решение во временной области
2 4 1 Задача со смешанными краевыми условиями
2 5 Эксперименты.
2 5 1 Линейное возрастание скорости с глубиной
2 5 2 Слоистые среды.
Глава 3. ОПТИМАЛЬНЫЙ ИДЕАЛЬНО СОГЛАСОВАННЫЙ СЛОЙ ДЛЯ
СИСТЕМЫ УРАВНЕНИЙ ТЕОРИИ УПРУГОСТИ.
31 Постановка задачи
3 2 Существование оптимальной сетки.
3 2 1 Построение импедансной функции
3 2 2 Построение конечно-разностной импедансной функции.
3 2 3 Порядок сходимости.
3 2 4 Восстановление шагов сетки.
3 3 Конечно-разностные схемы
3 3 1 Схема Вирье
3 3 2 Схема для PML.
3 4 Численный эксперимент.
Глава 4. МОДЕЛИРОВАНИЕ ВОЛНОВЫХ ПРОЦЕССОВ В АНИЗОТРОПНЫХ УПРУГИХ СРЕДАХ
41 Постановка задачи
4 2 Схема Лебедева
4 2 1 Свойства системы уравнений теории упругости для анизотропной среды.
4 2 2 Модификация системы уравнений
4 2 3 Конечно-разностная схема.
4 3 Схемы на повернутых сетках.
4 3 1 Построение схемы.
4 3 2 Исследование устойчивости.
4 3 3 Дисперсионный анализ.
4 3 4 Затраты на реали *ацию
4.4 Построение оптимальной сетки
4 4 1 Построение импедансной функции
4 4 2 Рациональная аппроксимация
4 5 Построение схемы с применением оптимальных сеток.
4 6 Численные эксперименты
4 6 1 Однородная среда
4 6 2 Слоистая среда.
Объектом исследований в данной работе являются сетки для конечно-разностных схем на предмет их оптимизации для моделирования волновых процессов в упругих средах, в том числе анизотропных
Актуальность. При моделировании волновых процессов достаточно распространенным требованием является знание волнового поля не внутри расчетной области, а лишь на некоторой ее границе. Наиболее показательные в этом смысле 1еофизические эксперименты Применение стандартных конечно-разностных методов для численною решения таких задач приводит к 1шантским требованиям на вычислительные ресурсы и временным затратам на их реализацию
Причина этою явления в том, что отображение совокупности параметров среды и особенности волновых процессов в среде в "данные наблюдений"(сейсмо1раммы) происходит с катастрофической потерей информации В частности, в трехмерных динамических задачах число параметров, входящих в описание задачи, оценивается величиной 0(N4), 1де N пропорционально линейному рашеру области, в которой развивается волновой процесс и обратно пропорционально шагу сетки Размерность пространства образов - число отсчетов на сейсмограммах - имеет порядок не выше 0{N) По-видимому, даже эта оценка является завышенной если принять во внимание гладкость зондирующего сигнала Однако уже этого сопоставления достаточно, чтобы представить насколько сильна потеря информативности при действии оператора, отображающего параметры среды в отсчеты на сейсмограммах
Стандартные способы аппроксимации этого оператора посредством построения разностной схемы для начально-краевой задачи, соответствующей выбранной модели волновых процессов, никак не учитывают этой потери информативности. В определенном смысле оператор оказывасчся слишком плохо аппроксимирован - каждая значащая цифра в ответе появляется как результат выполнения огромного числа арифметических операций над входными данными, лежащими в пространстве очень большой размерности Эту размерность нельзя уменьшить, уменьшая N, так как это приведет к увеличению локальной ошибки в промежуточных и окончательных результатах. Повышение порядка аппроксимации разностной схемы хотя и может приводить к некоюрым улучшениям, но они не принципиальны.
Аналогичная ситуация возникает и при реализации "неотражающих"!раничных условий, использование которых также необходимо для качественного моделирования геофизических экспериментов Наиболее употребимым является Идеально Согласованный Слой (PML) от ашлийскою Perfectly Matched Layer. Основной идеей метода является окружение расчетной области специально сконструированным слоем, таким что волна проходит через границу без отражений и затухает внутри этого слоя. Несложно видеть, что PML является исключительно лишь вспомогательной конструкцией и необходимо лишь знание решения на границе PML - расчетная область, однако стандартные методы реализации эгою подхода, не учитывающие этого требования, приводят к тому, что затраты на вычисление решения внутри PML доходят до 75% от общих затрат на моделирование эксперимента
Иных результатов удается достичь если, вводя определенную свободу в строение схемы во внутренних точках, пытаться оптимизировать погрешность в отдельных точках. В задачах, упоминавшихся выше, возможна экспоненциальная [32] и даже сверхэкпоненци-альная сходимость.
Практическая реализация идей, лежащих в основе упомянутого подхода, наталкивается на ряд трудностей, связанных с необходимостью обеспечения высокой точности в решении целого ряда вспомоттельных задач Понятно, что для достижения оптимальною результата, погрешности во всех промежуточных задачах должны быть минимизированы, в частности важным аспектом является построение количественной оценки скорости сходимости метода для волновых задач, чю представлено в данной работе Использование оптимальных сеток для моделирования волновых процессов в упругих средах является актуальной задачей и, как следствие, требует умения их построения и теоретического обоснования возможности их применения для таких задач.
Цель исследований - повышение эффективности расчета волновых нолей в упругих средах за счет использования спектрально-согласованных сеток для конечно-разностных схем, ориентированных на аппроксимацию решения лишь в заданных, определяемых условиями эксперимента, точках расчетной области
Научные задачи - построение спектрально-согласованных сеток для волновою уравнения с переменными коэффициентами и получение количественной оценки скорости сходимости численного решения для волновых задач, полученно! о с помощью данного подхода, построение идеально согласованного слоя, основанного на применении спектрально-со1ласованных сеюк, для уравнений теории упругости в изотропном случае, моделирование волновых процессов в анизотропных упругих средах с использованием спектрально-согласованных сеток
Основные этапы исследований:
1 получить количественную оценку точности расчета волновых полей с использованием конечно-разностных схем на спектрально-согласованных сетках, основанную на теории сходимости аппроксимаций Паде-Чебышева,
2 построить оптимальный идеально согласованный слой, основанный на использовании спектрально-согласованных сеток для схемы Вирье для системы уравнений динамической теории упругости в случае и зотропной среды;
3 разработать алгоритм построения спектрально-согласованных сеток для схемы Лебедева для системы уравнений динамической теории упругости в случае анизотропной среды;
4 провести серию численных экспериментов для верификации полученных теоретических результатов.
Научные методы исследований.
Теоретической основой решения поставленной научной задачи являются- современная теория аппроксимаций - рациональные аппроксимации, аппроксимации Паде, Паде-Чебышева, численные методы линейной ал1ебры - решение прямой и обратной спектральной задач для симметричных трехдиагональных матриц с использованием последовательностей Штурма, алюритма Ланцоша, метод спектральною исчерпывания, основанный на цепочках двумерных вращений Якоби, результаты современной теории конечно-разностных меюдов для решения уравнений математической физики - схемы Вирье, схемы Лебедева, схемы на повернутых сетках, спектральная теория линейных дифференциальных операторов, современная теория унруюсти
Разработанные алгоритмы и программы использовались для проведения численного моделирования волновых полей в средах различной степени сложности. Результаты подвергались сравнительному анализу с результатами, полученными с использованием других методов - аналитическими и конечно-разностными на равномерных сетках
Защищаемые научные результаты:
• Доказана экспоненциальная скорость сходимости конечно-разностного решения волнового уравнения с переменными коэффициентами вычисленного с помощью спектрально-со1ласованных сеток
• Разработан алгоритм построения оптимального идеально согласованного слоя (PML), основанный на применении спектрально-согласованных сеток, для системы уравнений динамической теории унруюсти в случае изотропной среды.
• Разработан и реализован алгоритм построения спектрально-согласованных сеток для решения системы динамической теории упругости для случая анизотропной упругой среды.
Новизна работы. Личный вклад.
• Получена оценка скорости сходимости конечно-разностною решения при использовании спектралыю-со1ласованных сеток для волнового уравнения с неременными коэффициентами, основанная на теории сходимости аппроксимаций Паде-Чебышева для мероморфных функций
• Построен и реализован оптимальный идеально согласованный слой для изотропной упругой среды с использованием модификации схемы Вирье на спектрально-согласованных сетках
• Разработан, теоретически обоснован и реализован алюритм построения спектрально-согласованных сеюк для схемы Лебедева для системы уравнений динамической теории упругости в случае линейных анизотропных сред, и проведены серии тестов для верификации полученных результатов
Теоретическая и практическая значимость результатов.
Получена количественная оценка сходимости метода оптимальных сеток для волновых задач, в том числе с переменными, включая разрывные, коэффициентами. Построен оптимальный идеально согласованный слой для уравнений теории упругости в случае изотропной среды, позволяющий существенно уменьшить вычислительные затраты при моделировании волновых процессов в neoi раниченных областях Разработан алгоритм построения спектрально-согласованных сеток для анизотропных ynpyiих задач, позволяющий эффективно моделировать распространение волн в таких средах и рассчитывать волновые ноля в заданных точках среды (приемниках).
Апробация работы и публикации.
• Основные положения и результаты докладывались на седьмой международной конференции, посвященной математическим и численным аспектам распространения волн "Waves 2005"(Провидэнс США, 2005) [56], ежегодном форуме Общества Индустриальной и Прикладной Математики "SIAM Annual Meeting 2006"(Бостон, США, 2006), 12 Международном симпозиуме по сейсмической анизотропии 12IWSA (Пекин, Китай, 2006) [55], II Всероссийской конференции, посвященной памяти академика А Ф Сидорова «Актуальные проблемы прикладной математики и механики» (Абрау-Дюрсо, 2004) [20], V Международной научно-практической геолою-1еофизической конференции-конкурсе молодых ученых и специалистов "Геофизика 2005"(Санкт-Петербург, 2005) [23], Второй Сибирской международной конференции молодых ученых но наукам о Земле (Новосибирск, 2004) [22], Научной конференции "Трофиму-ковские чтения"(Новосибирск, 2006)
• Результаты исследований по теме диссертации изложены в 6 опубликованных работах Из них одна статья в Сибирском Журнале Вычислительной Математики [21], 2 работы - это материалы Международных конференции [56], [55], 3 работы [20], [23], [22]- эго материалы российских международных конференций
Диссертация выполнена в Лаборатории вычислительных методов геофизики Института нефтегазовой 1еологии и геофизики СО РАН и Новосибирском Государственном Университете
Автор выражает искреннюю признательность научному руководителю к.ф-м н. В И Костину за всестороннюю поддержку и постоянное внимание.
Заключение
Основным результатом работы является развитие, реали зация и тестирование подхода построения сеюк для конечно-разностных схем, коюрый позволяет существенно сократить вычислительные затраты при моделировании волновых процессов, как в изотропных, так и в анизотропных упругих средах Его отличительной особенностью является локальная минимизация ошибки численного решения в заданных точках расчетной области, связанных с нуждами конкретного эксперимента(модели). В данной работе показано, что использование таких сеток для численного решения волновых задач позволяет добиться более чем экспоненциальной сходимости, что, в свою очередь, обеспечивает существенное сокращение, более чем в пять раз, вычислительных затрат на расчет решения. Следует также отметить, что затраты на построение сетки пренебрежимо малы в сравнении с затратами на моделирование самого процесса(эксперимента), что позволяет вычислять их каждый раз непосредственно перед моделированием.
В работе представлено обоснование возможности использования подхода оптимальных сеток для расчета волновых полей в ynpyiHX, в том числе анизотропных, средах Задача моделирования распространения упругих волн в анизотропных средах с использованием конечно-разностных схем на сдвинутых сетках сама по себе является весьма актуальной Как показано в работе, использование схем Лебедева является весьма эффективным методом решения этой задачи. Более того, в работе представлена оптимизация сеток для схемы Лебедева, обеспечивающая сверхэкспоненциальную сходимость Приведенные эксперименты демонстрируют работоспособность и высокую эффективность применения схемы Лебедева на оптимальных сетках для решения этой задачи
Построение оптимального идеально согласованного слоя для системы уравнений динамической теории упругости является, пожалуй, наиболее значимым результатом Представленный в работе оптимальный PML, основанный на методе оптимальных сеток, позволяет сократить число точек, необходимое для достижения нужной точности более чем в десять раз, 3-4 точки вместо 30-50 в стандартном. Если учесть тот факт, что идеально согласованный слой окружает всю расчетную область, то при стандартной реализации он требует более 50% от общею числа вычислительных ресурсов, использование оптимального PML сокращает эти затраты до 5-10%
1. Ахиезер П.И Классическая проблема моментов и некоторые вопросы анализа, связанные с нею // М : Физматгиз, 1961, 310 с
2. Бабенко К.И. Основы численною анализа.// М : Наука, 1986, 744 с
3. Бахвалов Н С., Жидков II П , Кобельков Г.М. Численные методы.// М • Наука, 1987, 542 с.
4. Бейкер Дж., Грейвис-Моррис П Аппроксимации Паде // М. Мир, 1986, 502 с.
5. Березин И.С , Жидков Н П., Методы вычислений.// М/ Государственное издательство физико-математической литературы, 1962, 464 с
6. Владимиров В С Обобщенные функции в математической физике //М/Наука, 1976, 280 с.
7. Владимиров В С. Уравнения математической физики.// М • Наука, 1972, 395 с8j Гангмахер Ф.Р Теория матриц // М • Наука, 1988, 552 с
8. Годунов С К Современные аспекты линейной алюбры.// Новосибирск- Научная книга, 1997, 390 с
9. Годунов С К., Антонов А Г., Кирилюк О П , Костин В И Гарантированная точность решения систем линейных уравнений в евклидовых пространствах // Новосибирск, Наука, 1990, 352 с
10. Годунов С К , Рябенький В.С Разностные схемы, введение в теорию // М.: Наука, 1973, 400 с
11. Голуб Дж ,Ван Лоун Ч. Матричные вычисления // М.: Мир, 1999, 548 с
12. Друскин В Л., Книжнерман Л А. Оценки ошибок в простом процессе Ланцоша при вычислении функций от симметричных матриц и собственных значений // Журнал вычислительной математики и математической физики, 1991, т 31(7), с 970-983
13. Зорич В А Математический анализ ч.1 // М,- Наука, 1981, 544 с
14. Зорич В.А Математический анализ ч 2 // М.: Паука, 1984, 640 с
15. Колмоюров А.Ф , Фомин С В. Элементы теории функции и функциональною анализа // М : Наука, 1989, 624 с
16. Коновалов А Н. Введение в вычислительные методы линейной алгебры // Новосибирск Наука, 1993, 159 с
17. Книжнерман JT А Качество аппроксимации к хорошо отделенному собственному значению и расположение "чисел Ритца"в простом процессе Ланцоша // Журнал вычислительной математики и математической фишки, 1995, т35(10), с.1459-1475
18. Книжнерман Л А. Простой процесс Ланцоша оценка погрешности гауссовой квадратурной формулы и их приложения // Журнал вычислительной математики и математической физики, 1996, т36(11), с 5-19
19. Лисица В В Оптимальные сетки для численного решения волнового уравнения // Тезисы докладов II Всероссийской конференции, посвященной памяти академика А Ф. Сидорова "Актуальные проблемы прикладной математики и механики", 2004, с 7172.
20. Лисица В В Оптимальные сетки для численною решения волновою уравнения с неременными коэффициентами // Сибирский журнал вычислительной математики, 2005, т8(3), с 219-229
21. Лисица В В. Оптимальные сетки для численною решения волнового уравнения с переменными коэффициентами // Тезисы докладов Второй Сибирской международной конференции молодых ученых по наукам о Земле, 2004, с 110
22. Наймарк М А. Линейные дифференциальные операторы // М. Паука, 1969, 528 с
23. Никишин ЕМ, Сорокин ВН Рациональные аппроксимации и ортоюнальность// М • Наука, 1988, 256 с
24. Парлетт Б. Симметричная проблема собственных значений.// М Мир, 1983, 384 с
25. Пашквоский С. Вычислительные применения многочленов и рядов Чебышева // М Наука, 1983, 327 с.
26. Самарский А А Теория разностных схем.// М Наука, 1983, 616 с.
27. Сеге Г. Оргоюнальные мноючлены // М.: Государственное издательство физико-математической литературы, 1962, 500 с.
28. С>етин П.К. Классические ортогональные многочлены // М : Наука, 1976, 328 с.
29. Тихонов А Н , Самарский А А Уравнения математической фишки // М.: Паука, 1972, 736 с
30. Asvadurov S , Druskin V, Guddati M.N., Knizhncrman L On optimal finite-difference approximation of PML // SAIM Journal of Numerical Analysis, 2003, n 41, pp 287-305
31. Asvadurov S., Druskin V., Knizhnerman L. Application of the difference Gaussian rules to solution of hyperbolic problems // Journal of Computational Physics, 2000, n. 158, pp 116-135.
32. Asvadurov S , Druskin V , Knizhnerman L Application of the difference Gaussian rules to solution of hyperbolic problems II Global expansion // Journal of Computational Physics, 2002, n 175, pp 24-29.
33. Asvadurov S , Druskin V. and Moskow S Optimal Grids for Anisotropic Problems // Electronic Transactions on Numerical Analysis (ETNA), accepted
34. Becache E., Fauqueux S., Joly P. Stability of Perfectly Matched Layers, Group Velocities and Anisotropic Waves // Rapport de recherche n° 4304, Novembre 2001, 35p.
35. Berenger J -P. Perfectly matched layer for the absorption of electromagnetic waves// J. Comput Phys , 1994, v 114, pp 185-200.
36. Borcea L , Druskin V. Optimal finite-difference grids for direct and inverse Sturrn-Liouville problems // Inverse Problems, 2002, n 18, pp 979-1001.
37. Collino F , Conditions d'orde eleve pour des modeles de propagation d'ondes dans de domaines rectangulaires.// INRIA report de recherche n° 1790, 1993
38. Collino F., Tsogka C. Application of PML absorbing layer model to the linear elastodynamic problem in anisotropic heterogeneous media // Geophysics, 2001, v.66, pp 294- 307.
39. Davydycheva S , Druskin V., Habashy T. An efficient finite-difference scheme for electromagnetic logging in 3D anisotropic inhornogeneous media // Geophysics, 2003, v 68(5), pp 1525-1536
40. Engquist В., Majda A Radiation boundary conditions for acoustic and elastic wave calculations// Comm Pure Appl Math., 1979, v.32, pp. 313-357
41. Golub G.IL, Welsch J.H Calculation of Gaussian Quadrature Rules// Mathematics of Computation, 1969, v.23, no 106, pp 221-230.
42. Gragg W.B , Ilarrod W J The Numerically Stable Reconstruction of Jacobi Matrices from Spectral Data // Numenshe Mathematik, 1984, v.44, pp 317-335
43. Guddati M., Tassoulas J. Continued-fraction absorbing boundary conditions for the wave equation // J. Comput Acoust, 1998, v 8, pp 139-156
44. Hagstrom T , Goodrich J. Accurate Radiation Boundary Conditions for the Linearized Euler Equations in Cartesian Domains // SIAM J Sci Comput, 2003, v 24, pp.770-795.
45. Hald H.O. Inverse Eigenvalue Problems for Jacobi Matrices.// Linear Algebra and Its Applications, 1976, v 14, pp 63-85
46. Ingerman D., Druskin V, Knizhnerman L. Optimal finite-difference grids and rational approximation of square root I Elliptic problems.// Comm of Pure and Applied Mathematics, 2000, n 53, pp. 1039-1066.
47. Levander A R Tourth-ordr finite-differrence p-sv seismorgems.// Geopthysics, 1988, v.53, pp. 1425-1436.
48. Lisitsa V. Lebedev Schmes for Simulation of Waves' propagation in Anisotropic Elastic Media // Expanded Abstracts of The 12th International Workshop on Seismic Anisotropy, 2006, pp.122-124.
49. Lisitsa V. Optimal Unsplit Perfectly Matched Layer (PML) for 2D Elasticity.// Proceedings of the 7th International Conference on Mathematical and Numerical Aspects of Wave Propagation "Waves 2005", 2005, pp 97-99
50. Metodes numeriques d'ordre eleve pour les ondes en regime transitoire. Editeur Cohen G // INRIA, Collection Didactique, 1994, p 523
51. Saenger E H , Gold N., Shapiro S A. Modeling the propagation of the elastic waves using a modified finite-difference grid // Wave Motion, 2000, v 31(1), pp 77-92.
52. Simon H D. The Lanczos Algorithm with Partial Reorthogonalization // Mathematics of Computation, 1984, v 42(165), pp 115-142
53. Vacus О Mathematical analysis of absorbing boundary conditions for the wave equation: the corner problem // Mathematical Computations, 2005, v.74, pp 177-200
54. Virieux J. P-SV wave propagation in heterogeneous media1 Velocity-stress finite-difference method // Geophysics, 1986, v 51(4), pp 889-901