Содержание к диссертации
Введение
Глава 1. Выход в (n+1)-мерное пространство 8
1.1. Постановка задачи 8
1.2. Теоретическое обоснование
1.3 Алгоритм 11
1.4 Моделирующая программа
1.5. Тестирование метода 17
1.6. Выводы 27
Глава 2. Некоторые сведения из топологии 28
Глава 3. Использование свойств неподвижной точки для решения систем линейных уравнений 34
3.1. Постановка задачи 34
3.2. Теоретическое обоснование 34
3.3. Алгоритм 38
3.4. Моделирующая программа 42
3.5. Тестирование метода
3.5.1. Тестирование на СЛАУ, полученных с помощью датчика случайных чисел 46
3.5.2. Тестирование на СЛАУ, коэффициенты которых представляют собой матрицу Гильберта 65
3.6. Выводы 74
Глава 4. Использование свойств неподвижной точки для решения систем нелинейных уравнений 76
4.1. Постановка задачи 76
4.2. Теоретическое обоснование 76
4.3. Алгоритм для случая n=2 82
4.4. Моделирующие программы 85
4.5. Тестирование метода 4.
5.1. Тестирование на системах квадратных уравнений 88
4.5.2. Тестирование метода с построением таблицы углов 97
4.5.3 Тестирование метода на системе, возникающей при расчете пространственной траектории горизонтальной скважины 99
4.6. Выводы 106
Заключение 107
Список литературы
- Моделирующая программа
- Теоретическое обоснование
- Тестирование на СЛАУ, коэффициенты которых представляют собой матрицу Гильберта
- Тестирование метода с построением таблицы углов
Моделирующая программа
Метод, рассматриваемый в этой главе, представляет собой алгоритм отыскания области, содержащей решение системы уравнений. В данной работе под уравнением от п независимых переменных х1,...,хп понимается связь между двумя функциями от этих переменных, определенная знаком равенства. Если имеется п уравнений, то можно говорить о решении системы этих уравнений, называя решением набор координат (х10, х20,…, хп0) некоторой точки л-мерного пространства, которые при подстановке в каждое из этих уравнений обращает эти уравнения в тождества. Таких решений может не быть совсем, а может быть несколько и даже бесконечное число решений.
Задача ставится следующим образом: разбить область поиска решений на подобласти и выделить те подобласти, внутри которых лежат решения. Таким образом, взяв любую точку подобласти, можно говорить о приближенном решении, которое отличается от точного не больше, чем на диаметр данной подобласти. Уточнять решение можно уменьшая размер подобласти.
Каждое из уравнений системы (1.1) описывает л-мерное многообразие, которое представляет собой нулевую линию уровня для (л+1)-мерной поверхности, заданной соответствующим уравнением из системы (1.2).
Каждое i-ое многообразие обладает тем свойством, что делит л-мерное пространство на две области: область, где Fi(x1,x2,...,xn) 0 и область, гдеFj(x1,x2,...,xn) 0. В окрестности многообразия находятся как точки одной области, так и точки другой. Имея л многообразий, получим 2л соответствующих областей.
Графики двух поверхностей в трехмерном пространстве У этих поверхностей нулевые линии уровня описываются уравнениями исходной системы. Точка плоскости (х,у), являющаяся решением системы, должна лежать как на одной, так и на другой линии уровня, то есть она лежит на их пересечении (рисунок 1.2).
Графики нулевых линий уровня График функции z=F(x,y) делит пространство на области, где z = F(x,y) 0 и z= =F(x,y) 0, поэтому в окрестности точки-решения найдутся и точки, где F1(x,y) 0 и точки, где F1(x,y) 0, а также точки, где F2(x,y) 0 и точки, где F2(x,y) 0 (рисунок 1.3). Таким образом, область, содержащая набор всевозможных сочетаний плюсов и минусов (в двухмерном случае этот набор состоит из четырех типов «+,+», «+,-», «-,+», «-,-»), и будет областью, содержащей решение.
Представим рассматриваемую область гиперплоскости в виде набора ячеек, размеры которых определяются шагом dx{dx\, dx2,..., dxn). 1. В каждом из 2" узлов ячейки вычисляем значение функции Fi(xi,x2,...,xn) и обращаем внимание на знаки F\ в узле. 2. Если во всех узлах ячейки знаки одинаковые, переходим к следующей ячейке, так как через данную ячейку многообразие Fi(xi,x2,...,xn)=0 не проходит. 3. Если в некоторых узлах ячейки встречается как знак «плюс», так и знак «минус», то это означает, что через данную ячейку проходит многообразие Fx{xx,x2,...,xn)=0. Аналогичным образом проверяем данную ячейку на содержание многообразий F2{x\,x2,...,хп)=0, Рз(х\,х2,...,хп)=0 ... Fn(x\,x2,...,xn)=0. 4. Если окажется, что ячейка содержит в себе все многообразия F1{x1,A 2,...,x/;)=0, i=0..п, то это ячейка становится «подозрительной» на решение, т.е. кандидатом на содержание внутри себя решения системы. При этом стоит заметить, что расстояние между любыми двумя функциями системы в любой точке «подозрительной» ячейки не превосходит диаметра этой ячейки. В некоторых задачах такого приближенного решения достаточно, но для поиска области, содержащей точное решение, перейдем к пункту 5. 5. Внутри «подозрительной» ячейки рассмотрим более детальную сетку, в узлах которой снова будем вычислять значения Fi(x1,x2,...,xn)=0, i=0..n. И если среди этих узлов найдутся такие, в которых будут содержаться все комплекты знаков (всевозможные сочетания «плюсов» и «минусов» всех функций), то это будет означать, что в данной ячейке находится решение системы.
Программная реализация включает в себя две подпрограммы. Первая отыскивает область, через которую проходят все многообразия (так называемая «подозрительная» область), вторая же проверяет, содержит эта область решение или нет.
В основе первой подпрограммы лежат пункты 1-4 вышеизложенного алгоритма. А именно, если в узлах ячейки нашлись два разных знака («+» и «-») функции F1, то многообразие, представляющее из себя нулевую линию уровня поверхности xn+1=F1(x1,...xn), содержится в данной ячейке. Таким образом, проверяются все ячейки на содержание внутри них нулевых линий уровня для всех поверхностей xn+1=Fi(x1,...xn) (i=1..n) и выбираются те ячейки, в которых найдутся все нулевые линии уровня. Блок-схема первой подпрограммы изображена на рисунке 1.4.
Блок-схема первой подпрограммы Рассмотрим подробнее алгоритм второй подпрограммы, реализующей пункт 5 из вышеизложенного алгоритма. На её вход поступают сами функции Fi (i=1..n), границы области поиска x0 (x10, x20,…, xn0) и x (x1, x2,…, xn), шаг детализации dx (dx1, dx2,… , dxn). Для определения все ли сочетания знаков имеются, заводится массив булевых переменных (переменные, принимающие два значения – false или true) размерности 2n, переменная bin_sign, принимающая значения в двоичной системе счисления, а также счетчик, который считает, сколько сочетаний знаков уже имеется на данном этапе подпрограммы.
На первом этапе вычисляется значение F1(x), где x – это рассматриваемый узел. В начале программы x=x0, затем далее по узлам (порядок обхода всей области будет приведен ниже). Если F1(x)0, тогда bin_sign=1. Затем вычисляется F2(x) и если F2(x)0, тогда bin_sign=10+bin_sign. Так вычисляем все Fi(x) (i=1..n) и если Fi(x)0, тогда bin_sign=10i-1+bin_sign. Тем самым переменная bin_sign характеризует сочетание знаков Fi(x) (i=1..n). Например,
Чтобы запомнить такое сочетание знаков, переведем bin_sign в десятичную систему счисления m=convert(bin_sign, dec, bin) и в специально заведенном нами булевом массиве m-ому элементу присвоим значение true, mas[m]=true.
При этом введем счетчик, который считает, сколько различных вариантов на данном этапе программы. Как только вариантов будет 2n, это будет означать, что найдены всевозможные сочетания знаков, следовательно, точное решение находится внутри рассматриваемой области. Блок-схема второй подпрограммы изображена на рисунке 1.5.
Теоретическое обоснование
Пусть С есть аффинное отображение пространства К1 на себя, при котором некоторая (л-І)-мерная плоскость отображается на себя. Аффинное отображение плоскости F?-\ порожденное отображением С обозначим через С.
Обозначим через є число, равное +1, если при отображении С каждое из полупространств отображается на себя, и равное -1, если при отображении С полупространства меняются местами. Тогда для знаков отображений Си С имеется соотношение sgnCsgn С =є.
Ориентировать л-мерное евклидово пространство - значит определить на всех его упорядоченных остовах функцию Деь,еі,...,еп)=±1, принимающую на всех остовах одного класса число, равное +1, а на всех остовах другого класса число, равное -1.
Ориентация остова симплекса.
Рассмотрим различные упорядочивания остова симплекса. Их (л+1)!. Два упорядоченных остова принадлежат к одному классу, если их нумерации переводятся друг в друга четной перестановкой. Ориентировать данный остов -это значит определить ориентацию каждого класса. То есть, у данного остова может быть только две ориентации. Для определения ориентации данного остова или симплекса достаточно задать ее на одном каком-нибудь упорядоченном остове этого симплекса (e0,e1,…,en).
Ориентация, принимающая на упорядоченном остове (e0,e1,…,en) значение +1, обозначается (e0,e1,…,en), а ориентация, принимающая значение -1, обозначается -(e0,e1,…,en).
Симплекс (остов) с определенной ориентацией называется ориентированным и обозначается (e0,e1,…,en) или tn, -(e0,e1,…,en) или n. При этом tn=Tn, n=Tn. Пусть tn - ориентированный симплекс данного Rn. Выпуклое замыкание остова симплекса называется телом ориентированного симплекса tn. Неподвижная точка Пусть дано отображение f:XX множества X в себя. Неподвижной точкой f называется любой такой элемент xX, для которого f(x)=x. Иначе говоря, неподвижная точка остается на месте при отображении f. Степень отображения
Пусть заданы две сферы S1 и S2 в n-мерном пространстве. Пусть поверхность S1 отображается непрерывным преобразованием f в поверхность S2. Определим степень отображения f. На сфере S1 выберем конечное число точек, расположенных достаточно плотно, представляющих собой -сеть. Рассмотрим триангуляцию сферы, т.е. разобьем сферическую поверхность, которая представляет собой (n-1)-мерное многообразие, на криволинейные (n-1)-мерные симплексы (вершины которых - точки -сети).
У каждого симплекса возможны две ориентации: положительная и отрицательная. Пусть у рассматриваемой триангуляции все ориентации положительны. Отобразим непрерывным образом триангуляцию сферы S1 в сферу S2.
Получим на сфере S2 набор симплексов, каждый из которых имеет либо положительную, либо отрицательную ориентацию. Составим сумму которая представляет собой сумму (n-1)-мерных объемов i симплексов сферы S2, причем каждый объем входит в сумму со знаком ориентации того симплекса, которому он принадлежит. В этом случае получится, что некоторые слагаемые будут положительными, а некоторые отрицательными.
Предельное значение суммы (2.1), когда измельчаются симплексы с коэффициентом обратно пропорциональным площади образа сферы, носит название степени отображения сферы S1 в сферу S2.
Пусть в евклидовом пространстве задано непрерывное преобразование F, имеющее конечное число неподвижных точек. Рассмотрим преобразование Ф, отображающее сферу S1 в единичную сферу S2 c центром в начале координат, определенное следующим образом. К каждой точке x0S1 применим преобразование F. Соответствующие этому преобразованию векторы перенесем в начало координат и пронормируем до длины равной единице. На сфере S2 получим точки, которые и будут образом Ф(x).
Пусть теперь сфера S1 cодержит одну неподвижную точку преобразования F. Тогда степень отображения Ф будет равна либо +1, либо -1.
Если же сфера S1 не содержит неподвижных точек преобразования F, то степень отображения преобразования Ф равна нулю. Если сфера S1 содержит несколько неподвижных точек, то степень отображения преобразования Ф этой сферы в единичную сферу складывается из степеней отображений соответствующих каждой из этих неподвижных точек и окружающих их сфер в случае, если эти сферы содержат только одну неподвижную точку.
В случае линейного преобразования F, когда возможна всего одна неподвижная точка преобразования F, степень отображения преобразования может быть либо нулем, либо плюс единицей, либо минус единицей.
Тестирование на СЛАУ, коэффициенты которых представляют собой матрицу Гильберта
Таким образом, получаем, что значение минимакса после размножения с нормировкой для симплекса (3.9) равно 0.4634, а для симплекса (3.10) 0.4824. Следовательно, значение разности в случае симплекса (3.9) равно 0.2284, для симплекса (3.10) равно 0.2277. Это происходит из-за того, что до нормировки было несколько векторов, длина которых была меньше 1. После нормировки эти векторы «растянулись» и оказали некоторое влияние на исследуемую разность.
Отсюда следует сделать вывод, что случай, когда векторов с длиной меньше единице получается немногим меньше п (где п - размерность) следует рассматривать более внимательно, при этом обращая внимание на то, в каких координатных частях гиперплоскости лежат эти векторы. Так как если они лежат в одной координатной части гиперплоскости, то можно говорить о том, что все рассматриваемые векторы лежат в одном л-мерном полупространстве, а значит решение снаружи.
Стоит отметить, что при многочисленных генерациях матриц, случай, когда разность положительна, а решение снаружи, встретилось только один раз, а именно в предложенном примере.
Для исследования поведения метода на плохо обусловленных системах были взяты системы уравнений, коэффициенты которых представляют собой матрицу Гильберта, известную своей плохой обусловленностью [51,56]. То есть в СЛАУ ЛА=, А=Н.= 1 . Коэффициенты вектора Ъ подбирались таким образом, чтобы решением системы была точка ль(1.1,1.1,..., 1.1). Симплексы строились аналогично тому, как это делалось в п.3.5.1. Исследование проводилось на размерностях от п=2 до л=100. Для сравнения данные системы также решались с помощью встроенных в математический пакет Maple методов LU-разложения и PLU-разложения [44, 46, 48], а также был реализован один из методов крыловского типа BiCGStab [34], для которого в качестве начального приближения была выбрана точка, удаленная от решения на расстояние 49-Jn. Все методы были ограничены 5000 итерациями. В таблице 3.1 в первой колонке указана размерность рассматриваемой СЛАУ, во второй - определитель матрицы Гильберта этой размерности, в третьей - число обусловленности, в четвертой, пятой и шестой - нормы разности между точным решением и решением полученным с помощью метода LU-разложения, PLU-разложения и BiCGStab соответственно, в седьмой и восьмой - разность минимаксов в соответствии с предлагаемом алгоритмом, в случае когда решение находится внутри рассматриваемого симплекса без нормировки после «размножения» и с нормировкой, в девятой и десятой - разность минимаксов в соответствии с алгоритмом, в случае когда решение не содержится внутри рассматриваемого симплекса без нормировки после «размножения» и с нормировкой. Исследование метода на матрице Гильберта, сравнение с другими методами размерность определитель число обусловленности Норма невязки для LU-разложения Норма невязки для PLU-разложения Норма невязки для BiCGStab r1, без r1, безнормировки нормировки,решение внутри решение снаружи r2, с нормировкой r2, с нормировкой решение внутри решение снаружи
Конечно, говорить о том, что в случае, когда решение не содержится внутри симплекса, разность всегда равна нулю, неверно. Однако если варьировать размер симплекса и вращать его, то во многих случаях разность получается таковой. Как видно из таблицы, с увеличением размерности системы, методы, заложенные в математический пакет Maple, дают очень большую погрешность, BiCGStab дает каждый раз разную погрешность от относительно маленькой до большой, это связано с зависимостью данного метода от точки, с которой начинается итерационный процесс.
В таблице 3.2 приведены результаты исследования разности минимаксов без нормировки после «размножения» для различных размеров симплекса, в том случаях, когда решение находится внутри симплекса. По строкам меняется размерность СЛАУ (в таблице обозначена как N), по столбцам расстояние от точки, являющейся решением системы, до вершин симплекса (в таблице обозначено как S). Заметим, что разность до четвертого знака не меняется. Это происходит из-за того, что вершины симплекса равноудалены от точки-решения.
В таблице 3.3 и таблице 3.4 приведены результаты исследования разности минимаксов для случаев, когда решение находится внутри симплекса, но при этом вершины симплекса не равноудалены от точки решения. По строкам меняется размерность системы (в таблице обозначена через N) , по столбцам сдвиг от центра симплекса до точки-решения (в таблице обозначен H). Сдвиг при этом происходит по прямой соединяющей точку-решение и n+1 вершину симплекса (напомним, что данная вершина строится по диагонали). Как видно из таблиц, чем ближе точка-решение к вершине, тем меньше разность. Это связано с тем, что в случае совпадения вершины и точки решения, вектор, полученный после преобразования F, будет просто точкой начала координат.
Тестирование метода с построением таблицы углов
С помощью алгоритма из пункта 4.3. и соответствующей программы также была решена система, возникающая при расчете пространственной траектории горизонтальной скважины. Задача взята из диссертационной работы Иткина В. Ю. «Математические модели пространственных траекторий при проектировании кустовых скважин» [33].
Для начала распишем, как получается необходимая система. Рассматривается модель траектории скважины вида "дуга окружности -отрезок - дуга окружности". Рассматривается скважина, наклонно-направленная часть которой состоит из 4-х участков (рисунок 4.19): - вертикальный участок от устья до точки зарезки – отрезок длины L0; - участок набора зенитного угла – дуга в вертикальной плоскости радиуса R1; - участок стабилизации – наклонный отрезок длины L1; - участок изменения зенитного угла и азимута – дуга в наклонной (вообще говоря) плоскости радиуса R2.
Направление вертикального участка – строго вниз, определяется вектором V0=[0,0,1]. Направление участка стабилизации определяется единичным вектором V1. Направление в точке входа в пласт определяется единичным вектором V2.
Исследование показало, что метод хорошо работает в случаях, когда детальность сетки позволяет разделить несколько решений между собой. Метод обладает свойством естественного параллелизма, что позволяет эффективно использовать многопроцессорные системы. А именно, можно поделить узлы между процессорами и поручить каждому процессору работу со своей группой узлов. Тем самым время счета уменьшается пропорционально количеству процессоров[55].
Разработан новый метод отыскания области, содержащей решение системы уравнений, основанный на свойствах нулевых линий уровня функций, входящих в систему. Метод включает в себя: алгоритм нумерации узлов сетки; алгоритм отыскания «подозрительной» на решение подобласти; алгоритм проверки «подозрительной» подобласти на содержание внутри себя решения.
Разработан новый метод отыскания области, содержащей решение системы линейных алгебраических уравнений, основанный на свойствах неподвижной точки и степени отображения непрерывного преобразования. Метод включает в себя: алгоритм построения непрерывного преобразования, неподвижная точка которого является решением исходной системы; алгоритм проверки расположения векторов, а именно получение ответа на вопрос лежат ли n+1 вектор n-мерного пространства в одном n мерном полупространстве.
Предложен новый метод отыскания области, содержащей решение системы нелинейных уравнений, основанный на свойствах неподвижной точки и степени отображения непрерывного преобразования. Метод включает в себя: алгоритм построения непрерывного преобразования, неподвижная точка которого является решением исходной системы; алгоритм определения ориентации расположения точек; алгоритм подсчета степени отображения.
Написан комплекс моделирующих программ, реализующих разработанные методы. Проведено исследование методов на линейных системах различных размерностей, с различным числом обусловленности, на системах, коэффициенты которых представляют собой матрицу Гильберта; на нелинейных системах различной размерности; проиллюстрировано изменение степени отображения в случае системы, состоящей из двух квадратных уравнений. Также проведено тестирование одного из методов для решения задачи расчета пространственной траектории горизонтальной скважины.
Результаты исследования показали, что методы хорошо работают при достаточно мелком шаге сетки, на которую делится область. В нелинейном случае, это означает, что детальность сетки позволяет разделить несколько решений между собой. Также следует отметить, что методы обладают свойством естественного параллелизма, что позволяет эффективно использовать многопроцессорные системы.
