А.А. Булычев, И.В. Лыгин, В.Р. Мелихов - Численные методы решения прямых задач грави- и магниторазведки, страница 19

Такая оценка нами была сделана:r p+1sup G ( s ) − G p ( s ) ≤ 2G C (δ , D ) p+ 2Rs ≥R1r1−R.В этой формуле C (δ , D ) = ∫ δ (ξ ,ζ ) dS , r – расстояние от центраDразложения σ0 до наиболее удаленной точки σ, принадлежащей области D,R – радиус круга, включающего в себя область D. Определив величинуm n (σ 0 ),числа p, можно перейти к вычислению значений членов ряда( s − σ 0 ) n+1где 0 ≤ n ≤ p .Для того чтобы вычислить эти значения необходимо определитьзначения моментов m n (σ 0 ) .

Как можно аналитически представитьмоменты нами уже рассматривалось. Так для многоугольника спостоянной плотностью δ это представление имеет вид:Niδα ν (σ ν +1 − σ 0 ) n+ 2 − (σ ν − σ 0 ) n+ 2 .m n (σ 0 ) =∑2( n + 1)( n + 2) ν =1[]Для вычисления моментов такой модели можно предложитьрекуррентную формулу:128Niδ∑ uνn ,2( n + 1)( n + 2) ν =1uν ,n+1 = uν ,n (σ ν − σ 0 ) , n = 0, 1, 2, …,m n (σ 0 ) =uν , 0 = (α ν − α ν −1 )(σ ν − σ 0 ) 2 .Как видно из приведенных формул, значение момента определяетсярасстояниями от вершин многоугольника σν до центра разложения σ0,причем эти расстояния при вычислении моментов порядка n возводятся встепень (n+2), затем суммируются или вычитаются.

Как мы уже отмечали,такие действия могут привести к потере значащих цифр и, соответственно,к потере точности в вычисленных значениях m n (σ 0 ) . Соответственно,следует таким образом организовать вычисление членов ряда, чтобыпотери точности не происходило. Например, это можно осуществитьследующим образом. Поскольку нас в данном случае больше интересуютm n (σ 0 ), то для ихне сами моменты, а значения членов ряда( s − σ 0 ) n+1вычисления можно воспользоваться следующим представлением:Nm n (σ 0 )iδ1=⋅α ν (σ ν +1 − σ 0 ) n+ 2 − (σ ν − σ 0 ) n+ 2 =∑n +1n +1(s − σ 0 )(s − σ 0 )2( n + 1)( n + 2) ν =1[]n +1N⎤⎡⎛ σ − σ ⎞ n+1⎛ σν − σ 0 ⎞iδν +10⎜⎟⎜⎟(σ ν +1 − σ 0 ) − ⎜(σ ν − σ 0 )⎥ .=∑ αν ⎢2( n + 1)( n + 2) ν =1 ⎢⎜⎝ s − σ 0 ⎟⎠s − σ 0 ⎟⎠⎥⎦⎝⎣Вычисление по приведенной формуле может оказаться болеепредпочтительным, поскольку в данном случае возводятся в степень (n+1)числа меньше единицы.

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

Если организоватьсуммирование от первых членов ряда к последующим, то вновь появляетсявозможность потери значащих цифр, что приведет к ошибочномурезультату. Для того чтобы этого избежать, члены ряда следуетсуммировать в обратном порядке: от конечного – к начальному.7.Всесделанныезамечания, связанныес возникающимивычислительными ошибками, справедливы и для трехмерных моделей.Можно вспомнить вывод аналитического выражения поля силы129притяжения горизонтальной пластины.

При выводе этой формулы мыотмечали, что могли бы остановиться на одном из этапов, и при этомполучили бы выражение, зависящее от разницы арктангенсов, т.е. отразницы некоторых углов. При этом, при достаточном удалении расчетнойточки от пластины, эти углы были бы близки по своему значению к числуπ. А так как поле будет определяться разницей таких углов, то неизбежнопроисходит потеря значащих цифр и возникает ошибка в окончательномрезультате (вспомним пример, приведенный в начале лекции:3.1415613 − 3.1415524 = 0.0000089 ).

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

В качестве примера мы уже приводили “ручной” способрасчета аномальных значений полей силы тяжести или магнитного поля спомощью палеток (палетки Гамбурцева, Юнга и т.п.). Здесь ошибка методавозникает за счет того, что область, которая может быть намипредставлена с достаточной детальностью и от которой рассчитываетсяаномальное поле, в процессе вычисления аппроксимируется набором ячеекпалетки, накрывающих ее для каждой расчетной точки.Такая же ситуация возникает и при аппроксимации заданныхобъектов при расчетах прямых задач на ЭВМ. Это связано с тем, что взависимости от вида входных данных, от поставленных задач,исследуемый объект может аппроксимироваться разным наборомэлементарных тел – многогранниками, многоугольниками, системой130призм, пластинами, материальными линиями или точечными источниками.При этом в зависимости от задачи, эти тела могут иметь постоянную илипеременную плотность или намагниченность. Попытка создать одиналгоритм на все случаи, возникающие при практической работе, неувенчались успехом, даже для двухмерных задач, для которых, как мывидели, аналитические решения получены для широкого круга моделей.Это тем более актуально при решении трехмерных задач, где количествовводимой информации на порядки больше информации, необходимой длядвухмерных задач.

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

изолированные источники. Такие задачипринято называть рудными. В то же время в практике гравиразведки имагниторазведки большую роль играют так называемые структурныезадачи, когда модель среды представлена набором субгоризонтальныхпластов и границ. В этом случае, помимо способа представлениямодельной среды возникает вопрос о радиусе учета влияния источниковполя. Это приводит к появлению ошибки, связанной с ограничениемобласти “интегрирования”, поскольку не все массы исходной моделиучаствуют в расчетах. Выбор радиуса учета влияния масс при расчететаких моделей является очень важным, поскольку от этого зависитточность расчета аномальных значений.8.Рассмотрим часть из этих вопросов на некоторых примерах. В началерассмотрим вопрос, связанный с аппроксимацией исходной моделинабором тел простой формы.

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

Теперь будем разбивать эту призму на 2, на 4на 16 равных частей и т.д. так, чтобы сечение каждой такой ячейки былозаменятькаждуюячейку(призму)представленоквадратом,соответствующим горизонтальным цилиндром и рассчитывать от нихэффект. Результаты показывают, что в данном случае уже при числецилиндров, равном 16 точное поле и поле от модели практическисовпадают.Здесь стоит отметить следующее. Выпишем интегральную формулудля вертикальной компоненты силы тяжести для прямоугольника изаменим ее приближенной интегральной суммой:132ξ2 ζ 2g z ( x , z ) = G ∫ ∫ δ (ξ , ζ )ξ1 ζ 1ζ −zdξdζ ≈(ξ − x ) 2 + (ζ − z ) 2N x Nz≈ G ∑∑ δ (ξ i , ζ j )i =1 j =1ζj−z∆ξ∆ζ ,(ξ i − x ) 2 + (ζ j − z ) 2где ∆ξ и ∆ζ – размеры ячеек (призм), на которые поделена наша исходнаяпризма, Nx и Nz – количество этих ячеек вдоль оси oX и oZ, ξi и ζj –координаты центров ячеек.Чем меньше значения ∆ξ и ∆ζ, и,соответственно, чем больше значения Nx и Nz тем лучше эта суммапредставляет исходный интеграл.

Замена исходного интеграла такойсуммой в численных методах носит название формулы прямоугольниковдля вычисления интегралов. В численных методах также показывается,какова точность вычисления интеграла по такой формуле. Для вычисленияисходного интеграла можно было бы предложить и другие способы егочисленной оценки (формула трапеций, Симпсона и т.д.), но остановимся натой формуле, которую мы записали (формула прямоугольников). Сфизической точки зрения эта формула представляет собой эффектгоризонтальных цилиндров линейная плотность которых равнаδ л = δ ∆ξ∆ζ .

Таким образом в данном случае интегрирование по методупрямоугольников совпадает с рассматриваемой аппроксимационноймоделью в виде горизонтальных цилиндров.9.Продолжим деление исходного многоугольника на большее числоячеек, которые будем представлять в виде горизонтальных цилиндров.Казалось бы, такое деление должно повысить точность вычислениязначения поля. Однако, при числе ячеек Nx=8192, Nz=4096 результатоказывается совсем не похожим на ожидаемый.133При дальнейшем увеличении детальности (Nx = 16384, Nz = 8192)расхождения в решении становятся еще более значимыми.Заметим, что расчеты проводились с 16 битными числами.

Этосоответствует типу переменной “float” в языке Cи. Появившийся дребезг врешении связан с ошибками округления, которые, накапливаясь, вносятосновной вклад в результат. Однако если те же самые расчеты произвестис 32 битными числами (“double”), то такие ошибки при данном разбиениине будут заметны.