К проблеме повышения точности вычислительных схем в задаче аналитического продолжения гравитационного поля Текст научной статьи по специальности «Математика»

2024

ВЕСТНИК ПЕРМСКОГО УНИВЕРСИТЕТА

Геология

Том 23, № 1

ГЕОФИЗИКА

УДК 550.83

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

З.З. Арсанукаев, Е.Г. Рудаковская

Российский государственный химико-технологический университет им. Д.И. Менделеева 125047, Москва, ул. Миусская, 9. E-mail: zaindy@mail.ru (Статья поступила в редакцию 18 июня 2023 г.)

Рассматриваются различные методы построения вычислительных схем разностного метода с использованием метода сеток при решении задачи аналитического продолжения заданных значений гравитационного поля на поверхности Земли в нижнее полупространство. Оценивается погрешность замены непрерывного уравнения Лапласа разностным уравнением. Приводятся результаты вычислительных экспериментов аналитического продолжения заданных значений поля при различных шагах сетки. Ключевые слова: разностный метод, уравнение Лапласа, разложение в ряд Тейлора.

DOI: 10.17072/psu.geol.23.1.57

Введение

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

© Арсанукаев З.З., Рудаковская Е.Г., 2024

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

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

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

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

Обратимся теперь к основным положениям разностного метода и используемым ап-проксимационным схемам (Крылов, Бобков и др., 1975). Рассмотрим несколько методов замены дифференциальных уравнений разностными уравнениями.

В основе первого метода лежит замена производных от функции и(х,у) разностными отношениями.

Пусть в некоторой области D, ограниченной контуром Г, задано уравнение

¿(и) = а(х,у)ихх + 2Ь(х,у)иху+с(х,у)иуу + 2й(х,у)их + 2е(х,у)иу++д(х,у)и =/(х,у). (1)

Будем считать, что это уравнение в области D имеет эллиптический тип, т.е. для него выполняется условие 3(х,у) = Ь2(х,у) — а(х,у)с(х,у) < 0 всюду в D. Будем считать также, что к уравнению (1) присоединены граничные условия. Начнем с изложения понятия сетки.

Задача нахождения решения уравнения (1) в области D часто оказывается трудно разрешимой или вообще неразрешимой в конечном виде. Поэтому в практике вычислений отказываются от нахождения решения во всей области D и ставят задачу приближенно найти решение лишь на некотором множестве точек из D. Это множество точек называют сеткой. Его берут конечным и настолько густым, чтобы оно приближенно заменяло область D. При этом расположение точек сетки в принципе может быть любым. Обозначим через множество точек сетки. Пусть (XI, у'!) есть некоторая точка из Б^, а и(х1,у]) - значение функции и(х,у) в этой точке. Таких значений конечное число. Зада-

ча их нахождения принципиально возможна. Она определяется дифференциальным уравнением, граничными условиями, присоединенными к этому уравнению, областью D и контуром Г. Суть метода сеток состоит в том, что на основе этой информации строятся некоторые соотношения, отражающие свойства дифференциального уравнения и граничного условия и позволяющие приближенно вычислить значения и(х1, х{).

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

Проведем на плоскости переменных х и у совокупность прямых

Х = Хо + Ш,у = Уо +]1,

где х0, уо - координаты точки М, лежащей внутри D; h - шаг сетки по направлению х; 1 - шаг сетки по направлению у (^ 1 - положительные числа), j = 0,±1,±2,±3 ... Точки пересечения этих прямых (х = х±,у = уу ), содержащиеся внутри области D, и образуют множество Б^ точек прямоугольной сетки. Эти точки принято называть узлами сетки.

Пусть в области D задано дифференциальное уравнение (1). Рассмотрим прямоугольную сетку Бк, заменяющую область D, и в ней точку (х^, уу) будем называть внутренней, если четыре её соседние точки (Ъ-1,У]), (Ъ+1,У]), (*иУ}-д, (х„У]+1) принадлежат замкнутой области Б^ = D+Г. Поставим задачу заменить дифференциальное уравнение (1) в точке (хьУ]) сеточным уравнением. Такую замену мы выполним путем выражения каждой производной, входящей в уравнение, через значение функции в узлах сетки. Для функции и(х,у) двух независимых аргументов х и у соответствующие выражения имеют следующий вид

К2 ,д3и

и(х1+1'У] ~)-и(х1-1'У] ) 2К

^-^(Я+О КУ) -1 < 0 <! (2)

(ду)(х1-У])

_ и(хьу]+1 )-и(х1,у]-1) 21

ду

I2 ,д3и Л

7(-ф-)(?<1'У+ О -1 < т <1 (3 }

д2и _

и(*1+1>У] )-2и(х1,у] )+и(х-1,у] ) К2

К2 ,д4и . ,

12 ( дх4 Хъ+т'Ку})' 1<т <1

,д2и _

^ Я,,2 )

(4)

ду2 (Ъ,У¡) '

— — (д4и \

—12 ( ду4 )(хьу ¡+(т1)

и(хьу]+1)-2и(х1,у])+и(х1,у]-1)

1<а' <1

(5)

д2и

( дхду ду( дх)

д ди

ХЬУ])

+ 0(12) -

'дхду 1ду дх

(д^)(хУУ ]+1)-(д^)(хуУ]-1) 21

1 и(Х1 + 1,у]+1 ))+и(Х1-1,у] + 1 ) 2

21[ 2Н + и(1г

и(х1+1,у]-1 У)-и(х1-1,у]-1)

2Н

+0(12)

-о (к2)]

(6)

в предположении, что y=const Отметим, что

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

д2и

дхду

по формуле (6) используются точки в 4 углах большого квадрата и отмеченные на рис. 4 черными кружочками (формулы (2)-(5) получены из разложения функции и(х^у^ в окрестности точки (х^у^) по формуле Тейлора с остаточным членом в форме Лагран-жа.

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

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

отсутствует смешанная производная, т.е. Ь(х, у)= 0. Обозначим

а(х1,у]) - ац .,Г(х0у]) - [ц.

Используя формулы (2)-(5) и уравнение (1), получим

[Ь(и)] х.,у.) = аи(-

д 2 и

х,,у,)~ дх2

д2и ди

(хьУ ¡)

+

+С1)(~д^)(хьУ1) + 2ЛЧ(д$(хьУ]) + 2ечф(ХЬУ]) + д1]и(Х1,У}) -

и(х1+1,У] )-2и(хьУ] )+и(х1-1,У] ) . = О-И--7 ---+

1} ь2

и(хьу]+1 )-2и(хьу} )+и(хьу}-1 ) Сц--^---+

I] I2

. и(х1+1,У] )-и(х1-1,У]) .

-7ь-+

2еи-1—--1—++д1]и(х1,у]) +

0(к2+12) -

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

[ Ь(и)] (ХиУ.) - 1к (и(х1,у])) +0(к2+12) - Гц,

(7)

+

где положено

(и(х1,у])У) - Аци{М]) + Вииц-1П+ ац о _ аи _ г) _ сч | еч

ь2 + ь ,оц - ь2 ь, и1}- г2 + 1 ,

Р -са^ел р — 2ач 2сч | \ п - 12 ^ пу - ь2 12 + ь +

и(1,]) и т.д. - точные значения функции и(х,у) в точке (х„ уу). Если решение уравнения (1) имеет в области D непрерывные производные до четвертого порядка включительно, то в формуле (7) величиной 0(к2+12) при достаточно малых h и 1 можно пренебречь. Тогда окончательное сеточное уравнение, заменяющее дифференциальное уравнение (7), может быть записано в виде:

где под Щ] понимается приближенное значение функции и(х, у) в точке (х±,у]). Переход к этим приближенным значениям вызван тем, что в формуле (7) отбрасывается член О (к +12). Этот член имеет смысл погрешности, с которой сеточный оператор Ьк (щ¡) заменяет дифференциальный оператор L(u) во внутренней точке сеточной области. Для замены дифференциального оператора сеточным в этом случае привлекается схема точек (рис. 1), получившая название «крест». Сеточное уравнение (8) имеет место для всех внутренних точек сеточной области.

Рассмотрим теперь случай, когда в уравнение (1) входит смешанная производная

дд™ду и Ь(х, у), следовательно, не равняется

тождественно нулю. Замену этой производной осуществим по формуле (6) таким же путем, как и ранее. Сеточное уравнение, заменяющее дифференциальное уравнение (1) в точке к (Х{,уу), в этом случае будет таким

М (Щ¡) = ^.

(9)

часть узлов к внутренним, а часть узлов к граничным. Возьмем какую-нибудь внутреннюю точку с координатами х0, у0 и обозначим ее через М0 . Наряду с точкой М0 рассмотрим некоторую совокупность других внутренних точек М1, М2 ..., . Всего точек будет N+1. Обозначим через и^) точное значение функции и(х, у) в точке М^. Будем решать, как и прежде, задачу о замене уравнения (1) некоторым сеточным уравнением. С этой целью рассмотрим разностный оператор

1к(и(0)) = !£=о

Соки(к)

(10)

При записи сеточного уравнения (9) использовалась схема точек (рис. 4), получившая название «ящик». Уравнение (9) имеет место во всякой внутренней точке сеточной области, при этом точку (х^, у^) будем называть в случае схемы «ящик» внутренней тогда, когда восемь окружающих ее точек принадлежат области . Изложенный выше метод замены дифференциального уравнения сеточным прост и удобен, но он существенно использует то обстоятельство, что сетка прямоугольная и равномерная (в ней узлы по каждому из направлений взяты равноотстоящими). Если сетка такими свойствами не обладает, то указанный способ замены не может быть осуществлен.

Рассмотрим теперь второй метод, основанный на аппроксимации дифференциального оператора в целом. Будем по-прежнему рассматривать уравнение (1). Предположим, что уже построена некоторая сеточная область Бк. Узлы сетки будем считать расположенными произвольно (это может быть и неравномерная и не прямоугольная сетка). Пусть по некоторому правилу мы отнесли

и в нем неизвестные параметры Сок определим таким образом, чтобы значение дифференциального оператора L(u) в точке М0 и значение Ьк(и(0)) совпадали бы возможно более точно. Предположим, что все точки М1, М2 ..., Мн лежат в достаточно малой окрестности точки Мо. Тогда разности хк -Хо и Ук -Уо (1— к—N будут малыми числами и значение функции и(к) можно разложить по формуле Тейлора в окрестности точки Мо(хо.,уо) (до производной г-го порядка включительно и с остаточным членом в форме Лагранжа). Этим разложением мы воспользуемся для того, чтобы добиться упомянутого выше совпадения. Подставив это разложение в формулу (10), получаем:

^(и(о)) = Цр+Я=о Ар^Ц^дуч )(о) +

%к=о сокК(Ю

где через (—)(о) обозначено значение производной ^ в точке Мо, а через R(k) - остаток

разложения.

Явные значения для коэффициентов Арс{ слишком громоздкие, поэтому ограничимся только указанием на их смысл:

1) эти коэффициенты линейно зависят от величин сок (к=0,1,...№);

2) они содержат множители вида

(хк — Хо)р(ук—уо)4.

Если считать, что хк — Хо и ук — у о являются малыми величинами одного порядка малости, то из последнего свойства следует, что коэффициенты Арс{ являются величинами, порядок малости которых относительно

(Хк - Хо) или (ук - Уо) есть p+q. Потребуем, чтобы равенство

Г

p+q=0

А (■

дР+Чи

pqK дхРдуЧ

■)(0) = [L(u)]

(0)

Аг0 =

А

г-11

А

1г-1

(р + q = г).

Система уравнений возникает потому, что, вообще говоря, для коэффициентов Арч верно представление

А = YN

"■pq L¡

к=0

apq^ )(-0к,

(11)

где а

p q

(к) _

этой системы являются числа соо, со1*,.. ст*. Таким образом, можно утверждать, что равенство

[L(u)]

(0)

N

= Lk=o

Сок*и(к)

имело место для любой функции u(x,y), имеющей в области D непрерывные производные до порядка r включительно. Для этого необходимо, чтобы совпали коэффициентами .

ты правой и левой частях при ( dxPg q )(о)

для всех p и q, удовлетворяющих неравенству 0 < р + q < г. Это дает такую систему уравнений для определения величин с00, Со!^Con':

Аоо = д(хо ,Уо) (p + q = о)

А10 = 2d(xo, у о), А01 = 2е (Хо ,Уо) (р + q = 1)

А20 = а(хо,уо), Ап = 2Ь(хо,уо) А02 =

с (Хо ,Уо) (p + q = 2)

А30 =, А21 = Ап = Аоз = о (р + q = 3)

выполняется с погрешностью R *=ZN=0 С0к*и(к).

= Аог = 0

известные величины. Для определения величин Сок получим, таким образом, 1 (г + 1)(г + 2) уравнений. Эти уравнения линейны. Встает вопрос о разрешимости системы линейных алгебраических уравнений. Если число уравнений системы будет больше числа неизвестных, то система может быть несовместной. Будем поэтому считать, что выполняется неравенство

1 (г + 1)(г + 2) <Ы + 1.

Рассмотрим случай, когда 1 (г + 1)(г +

2) = N + 1. Предположим, что система уравнений для определения величин Сок имеет ранг, равный N+1, и, следовательно, она однозначно разрешима. Пусть решением

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

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

Определение параметров шаблонов 1-й крест», «2-й крест» для дискретной

аппроксимации уравнения Лапласа

Применим приведенные выше сведения из теории разностного метода при замене непрерывного уравнения Лапласа (уравнение

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

Пусть в некоторой области D, ограниченной контуром Г, задано уравнение Лапласа

Ли =

2 о и

дх2

+

2 о и

ду2

- 0

(12)

и к нему присоединено граничное условие (— С()1 + С()3) и=о,

(ди)(0) и т.д. с аналогичными коэффициентами в выражении (Ли)(0), получим систему линейных алгебраических уравнений для определения величин С0{ :

с00+ с01+с02+с03+с04=0,

ь

и\г -ф(М).

(13) (-С02 + С04) ]Т0,

Поставим задачу: заменить оператор Ли сеточным оператором. Будем рассматривать случай квадратной сетки, когда h = 1. Координаты узла сетки будем определять по формулам

хк - х0 + кь, ук-у0 + кь, к - 0, ±1,...

Заменим область Б сеточной областью Бь и в ней выделим йь* внутренних узлов и множество Г граничных узлов. Узел будем называть внутренним, если четыре его соседних узла принадлежат области Б.=0+Г. Пусть узел с координатами (х^,у^~) является внутренним. Для замены в этом узле дифференциального оператора сеточным изберем схему точек «крест» и в ней пронумеруем точки в таком порядке, как показано на рис. 1.

3

Т |} * 1

2

N

4

(С01 + Соз) 1,

( с 02 + с 04) "¡р1,

' ь3

(- С01 + С03) ^=0,

33

(-С02 + С04)-^0.

Решение этой системы легко выписывает-

ся:

С01 - с02 - с03 - с04~ _ 4

с°°-- р.

(15)

Таким образом, сеточный оператор, заменяющий дифференциальный оператор Ли в точке (х^ ,уу) на схеме точек «крест», будет иметь вид

1ь(и(0)) = -(и(1)+и(2)+и(3) +и(4)-4и(0)).

(16)

Рис. 1. Схема точек «крест»

Составим выражение для сеточного оператора Ьь(и(0)) = ^{=0 ^ки(к) (14) и в нем величины С0к определим так, чтобы равенство ((Ли)(0) - Ьь(и(0)) выполнялось с возможно большей точностью. Используя разложение в ряд Тейлора выражения (14) и сравнивая, как и выше, коэффициенты, стоящие в этом выражении при и(0), (—)(0),

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

(Ли)(0) - Ьь(и(0)) + 0(ь2).

Таким образом, замена дифференциального оператора (12) разностным оператором (14) приводит непрерывное уравнение Лапласа Ли - 0 , рассматриваемое в узле ((х1 ,У]') к дискретному уравнению Ьь(и(0)) =

(и( 1)+и(2)+и(3)+и(4)-4и(0))=0;

или после

сокращения на ^ Ф 0 к и(1)+и(2)+и(3)+и(4)-4и(0)-0.

2

1

2

ь

Рис. 2. Схема точек «1-й крест»

i-1j + 1

1-1 i-1

1- 2 0 -< 3 Л

1 i- 4

i + 1¡ + i

í+ij-1

Рис. 3. Схема точек «2-й крест»

U-1

Рис. 4. Схема точек «ящик»

Рассматривая дискретное уравнение Лапласа (17) во всех внутренних узлах заданной сеточной области при схеме точек «крест» с присоединением граничных условий, придем к системе линейных дифференциальных уравнений:

Ах = f, (18)

где А - матрица системы (18) - состоит из коэффициентов уравнения (17)

С01 = с02 = с03 = с04=1, с00=-4;

х - вектор неизвестных значений поля; f -вектор заданных значений поля (граничные условия). Следует отметить, что матрица А

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

Решение задачи аналитического продолжения гравитационного поля с использованием дискретного уравнения Лапаласа на модельном примере

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

дУЛ

циала —) вычисляется в конечном виде, т.е.

дг'

решение прямой задачи является точным. Решения прямой задачи используются далее в качестве «входных» значений в задаче аналитического продолжения в модельных условиях. Помимо этого, точные решения прямой задачи позволяют оценивать точность расчетных значений в нижнем полупространстве, полученных в результате аналитического продолжения, путем сравнения в соответствующих узлах с точными решениями прямой задачи. Одним из таких тел является прямоугольный параллелепипед бесконечного простирания в направлении одного из своих измерений (в данном примере в направлении оси Оу, т.е. перпендикулярно плоскости рисунка) и имеющего сечение 2.0х 2.4 км). Таким образом, решается двумерная задача, избыточная плотность тела принимается равной 1 г/см 2. Заданные значения поля (решения прямой задачи) располагаются на двух уровнях (профилях): на уровне Земли z = 0 и на уровне шага сетки у над уровнем Земли z = (ось Oz направлена вниз). Длина каждого профиля 16 км, шаг сетки равен 0.4 км, сеточная область представляет собой прямоугольник с основанием, равным длине профиля (16 км) и высотой,

равной расстоянию от z = 0 до верхней кромки параллелепипеда, находящейся в нижнем полупространстве на глубине 4 км.

Отметим сразу, что задача является некорректной (в классической постановке это задача о решения уравнения Лапласа в прямоугольнике (сеточной области) - так называемая задача Дирихле, где заданы граничные условия по всему контуру прямоугольника), поскольку почти всегда измеренные значения поля находятся на верхней стороне прямоугольника на уровне z = 0 или выше, иногда на боковых сторонах прямоугольника (при наличии данных шахтной гравиразвед-ки) и никогда на нижнем основании прямоугольника.

Матрица системы (18) в условиях данного модельного примера имеет размерность 390x410, т.е. число уравнений равно числу внутренних узлов: 39 узлов на одном профиле умножено на число уровней, равное 10, всего 390; а число неизвестных значений поля равно числу неизвестных значений поля на одном уровне, умноженному на число уровней, т.е. 41 х 10 = 410. Решение системы (18) проводилось методом итерации (Страхов В.Н., Страхов А.В., 1999), и оно оказалось неустойчивым, что объяснимо, когда система является недоопределенной, т.е. число уравнений меньше числа неизвестных. В нижнем полупространстве только на отметке, равной одному шагу сетки от поверхности Земли, ^ = 0) относительная погрешность (в среднеквадратичной норме) расчетных значений была порядка 10~ (т.е. разница в расчетных значениях поля и точных решений прямой задачи составила около 10 %), а ниже этой отметки относительная погрешность была уже более 100 %, т.е. начался распад поля. Было решено обеспечить устойчивость решения системы (18), рассматривая совместно две дискретные аппроксимации уравнения Лапласа: к схеме точек «крест» (которая дальше будет называться «1-й крест» (рис. 2)) добавляется схема точек «2-й крест» (рис. 3). Для решения предыдущей задачи о замене дифференциального оператора Ли в точке (х; ,уу) изберем иную схему точек (рис. 3). В этом случае оператор имеет вид

Lh'(u(0))= 12 (и(1)+и(2)+ +и(3) +и{4) - 4и(0)).

Соответствующие этому вычисления мы опускаем, они аналогичны предыдущим. Имеет место также формула

(Ли)(о) = Lh* (и(0)) + O(h2).

Операторы Lh(v.(0)) и Lh* ((и(0)) одинаковы по виду, однако качество их различно, так как различными являются схемы точек, на которых эти операторы заменяют дифференциальный оператор Ли.

Как и в случае схемы точек «1-й крест», записывая дискретное уравнение Лапласа для схемы точек «2-ый крест»

L¿(u(0)) = —2 (и(1) + и(2)+ и(3) + и(4) - 4и(0)) = 0, после сокращения на ^ Ф 0 получим вторую дискретную аппроксимацию уравнения Лапласа

и(1) + и(2) + и(3) + и(4) - 4и(0)=0 (19),

где С01 = с02 = с03 = с04=1, с00=-4

являются элементами матрицы A вместе с коэффициентами первой дискретной аппроксимации уравнения Лапласа.

Полученная система (18) теперь будет переопределенной (т.е. число уравнений будет больше числа неизвестных) с матрицей, имеющей размерность 780х 410. Решение ее является уже устойчивым - результаты расчетов, представленные с помощью специальной программы в графической форме (Арсанукаев З.З., Арсанукаев И.З., 2015), приведены на рис. 5 а (на оси ординат расположены значения — поля в мГал, на оси абс-

dz

цисс - координаты узлов сетки на профиле). Чтобы увеличить точность вычисляемых значений поля, расчеты проводились и на более густых вариантах сетки, т.е. с меньшими шагами сетки (рис. 5 б, в). Отметим, что на рис. 5 а, б, в самая верхняя кривая построена для расчетных значений поля в узлах уровня, расположенного на расстоянии шага сетки от поверхности Земли (z = h) в нижнем

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

двум дискретным аппроксимациям «1-й крест» и «2-й крест» девятиточечной схемы «ящик» (рис. 4). Имеет место следующая формула:

(Ли)Ш) = -[4Ьк(и((х1,у])) +

Рис. 5 а. Аномальные кривые при шаге сетки 400 м (-0 кривых)

Рис. 5 б. Аномальные кривые при шаге сетки 200 м (20 кривых)

Рис. 5 в. Аномальные кривые при шаге сетки 100 м (40 кривых)

Анализ результатов расчета показывает,

что на глубине, равной шагу сетки от уровня

z = 0, относительная погрешность на шагах

сетки 400, 200,100 м имеет значения, поря-

-6

док которых соответственно -0 3 ,-0 5,

-0'

2К

2К

4

к (и(( Х(,у]*)~)~\ —(А2и)(ц) — [Л3и +

2д4

Ли\ (Х + °ь,у} +т К) , где -1< а, т<1 (Канторович, Крылов, 1962).

Аналогично предыдущему, с точностью до бесконечно малых 2-го порядка, получим для схемы «ящик» на равномерной сетке дискретное уравнения Лапласа:

4Ьн(и((х1,у])) + V (и((х1,у]))=0, или

4(и(-)(1)+и(-)(2)+и(-)(3)+и(-)(4)-4и(0))+ (и (2)( -)+и (2)(2) +и (2)(3)+и (2)(4) -4и(0)=0,

или

0,8(и(-)(1)+и(-)(2)+и(-)(3)+и(-)(4))+0.2( и(2 (-)+и(2(2)+и(2) (3)+и(2 (4)) -4и(0) = 0,

где и(--) (-), и(--) (2), и(-) (3), и(-) (4) - значения поля, стоящие в узлах схемы «1-й крест» внутри схемы «ящик» (кроме узла 0, рис. 4); и(2)(-),и(2)(2),и(2)(3),и(2)(4) - значения поля, стоящие в узлах схемы «2-й крест» внутри схемы «ящик» (кроме узла 0, рис. 4); и(0) - значения поля, стоящие в центре схемы «ящик» (рис. 4). Коэффициенты схемы «ящик»

(т.е. восстановленные значения отличаются от точных, например для шага 200 м в среднем на 1/1000 %), а на глубине, расположенной выше отметки верхней кромки тела на 1 шаг сетки, погрешность имеет порядок с0(==—4

^ = с02^1- = С03(-) = С01(-) = 08, с0/ ) = с02( ) = с03( ) = с01< 2 = 02,

_- _2 _3

10 ,-0 , -0 (т.е. для шага 200 м восстановленные значения отличаются от точных в среднем на несколько процентов).

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

служат далее элементами матрицы, а вместе с коэффициентами двух дискретных аппроксимаций «1-й крест» и «2-й крест» уравнения Лапласа. Матрица А системы линейных алгебраических уравнений, полученной в результате использования совместно трех дискретных аппроксимаций уравнения Лапласа:

0706100807111107090313060504075511080620373611050609040609083619030609033721050407318107060437071103

«1-й крест», «2-й крест» и «ящик», будет переопределенной, и ее размерность для данного модельного примера равна 1170 X 410. Как можно видеть из табл., точность расчетных значений поля для уровня, расположенного на 1 шаге сетки в нижнем полупространстве от уровня Земли ^ = 0), увеличивается на порядок по сравнению с вычислительной схемой «1-й крест» + «2-й крест».

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

Таблица

Относительная погрешность

Глу- точ х - X зыч Е

бина, / точ

км Е

«1-й к» + «2- «1 -й к» + «2-й

й к» к» + I- «ящ.»

0,4 1,38440Е-03 4,92070Е-04

0,8 3,90422Е-03 3,492533Е-03

1,2 9,60364Е-03 4,769514Е-03

1,6 1,85948Е-02 1,069080Е-02

2,0 2,96881Е-02 2,459834Е-02

2,4 4,46757Е-02 4,197743Е-02

2,8 7,05478Е-02 6,123032Е-02

3,2 1,15451Е-01 8,713836Е-02

3,6 1,84324Е-01 1,312850 Е-01

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

плотности масс и таким образом решить обратную задачу.

Заключение

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

Библиографический список

Арсанукаев З.З., Арсанукаев И.З. Программа аналитического продолжения заданных на поверхности Земли значений гравитационного поля в нижнее полупространство с использованием дискретных схем и решения систем линейных алгебраических уравнений больших порядков на модельных примерах. Свидетельство о государственной регистрации программ для ЭВМ (SOFT). #2015661026, 15 октября 2015 г.

Канторович Л.В., Крылов В.И. Приближенные методы высшего анализа. М., 1962. 531 с.

Крылов В.И., Бобков В.В., Монастыр-ный П.И. Вычислительные методы высшей математики. Т. 2. Минск: «Вышэйшая школа», 1975. 672 с.

Страхов В.Н., Страхов А.В. Основные методы нахождения устойчивых приближенных решений систем линейных алгебраических уравнений, возникающих при решении задач гравиметрии и магнитометрии. II М.: ОИФЗ РАН, 1999. 51 с.

On the Problem of Improving the Accuracy of Computational Schemes in Analytical Continuation of Gravitational Field

Arsanukaev Z.Z., Rudakovskaya E.G.

Russian State University of Chemical Technology named after D. I. Mendeleev, 9 Miusskaya Str., Moscow 125047, Russia. E-mail: zaindy@mail.ru

The article discusses the variants of constructing the computational schemes of the difference grid method in solving the problem of analytical continuation of the gravitational field observed on the Earth's surface into the lower half-space. The error of replacing the continuous Laplace equation by a difference equation is estimated. The results of computational experiments of analytical continuation of given field values at different grid sizes are presented.

Key words: difference method; Laplace equation; Taylor series expansion.

References

Arsanukaev Z.Z., Arsanukaev I.Z. 2015. The program for the analytical continuation of the values of the gravitational field given on the surface of the Earth into the lower half-space using discrete schemes and solving systems of linear algebraic equations of large orders on model examples. Certificate of state registration of computer programs (SOFT).#2015661026, October 15, 2015.

Kantorovich L.V., Krylov V.I. 1962. Pribli-zhonnye metody vysshego analiza [Approximate methods of higher analysis]. Moskva, p. 531. (in Russian)

Krylov V.I., Bobkov V.V., Monastyrnyy P.I. 1975. Vychislitelnye metody vysshey ma-tematiki [Computational methods of higher mathematics]. T.2. Vysshaya shkola, Minsk, p. 672. (in Russian)

Strakhov V.N., Strakhov A.V. 1999. Osnov-nye metody nakhozhdeniya ustoychivykh prib-lizhonnykh resheniy system lineynykh algedraicheskikh uravneniy, vosnikayushchikh pri reshenii zadach gravimetrii i magnitometrii [Basic methods for finding stable approximate solutions to systems of linear algebraic equations arising in solving problems of gravimetry and magnetometry. II Moskva, OIFZ RAN, p. 51. (in Russian)