close

Вход

Забыли?

вход по аккаунту

?

МироненкоВКР

код для вставкиСкачать
 МИНИСТЕРСТВО ОБРАЗОВАНИЯ И НАУКИ РОССИЙСКОЙ ФЕДЕРАЦИИ
Федеральное государственное бюджетное образовательное учреждение
высшего профессионального образования
"БАШКИРСКИЙ ГОСУДАРСТВЕННЫЙ УНИВЕРСИТЕТ"
(БашГУ)
Факультет математики и информационных технологий
Кафедра информационных технологий
Мироненко Кирилл Олегович
Бакалаврская работа
Численное решение интервальных СЛАУ методом вращений в задачах компьютерного моделирования электрических полей катодной защиты протяженных сооружений
Направление 010500 - Прикладная математика и информатика К защите допущено:
Заведующий кафедрой,
д.ф.-м.н., профессор ___________ А. М. Болотнов "___"_____________ 2013г. Научный руководитель:
д.ф.-м.н., профессор
__________А. М. Болотнов
"___"____________2013 г.
Дата защиты: "___" _________ 2013 г.
Оценка: _________________________
Уфа - 2013
Оглавление
Глава 1. Интервальные вычисления.- 2 -
Введение- 2 -
История интервального анализа- 3 -
Интервальная арифметика.- 6 -
Основные обозначения- 8 -
Существующее ПО- 9 -
Библиотека интервальных вычислений- 11 -
Немного об IDE.- 11 -
Структура и принцип работы библиотеки- 11 -
Некоторые проблемы- 14 -
Тестирование библиотеки- 17 -
Пример кода для работы с библиотекой- 18 -
Глава 2. Моделирование электрических полей.- 19 -
Постановка задачи- 19 -
Сведение краевой задачи (8)-(11) к ГИУ- 21 -
Переход к интервальной СЛАУ- 23 -
Численные методы решения СЛАУ- 24 -
Метод Гаусса- 24 -
Метод вращений- 26 -
Проблема численного решения интервальной СЛАУ.- 29 -
Неоднородность ширины интервала решения- 30 -
Заключение- 32 -
Использованная литература- 33 -
Приложения- 35 -
Приложение А. Исходный код библиотеки интервальных вычислений- 35 -
Приложение Б. Исходный код программы.- 42 -
Приложение В. Результаты.- 49 -
Глава 1. Интервальные вычисления. Введение
Первоначально интервальные методы возникли как средство автоматического контроля ошибок округления на ЭВМ и впоследствии превратились в один из разделов современной прикладной математики. В основе лежала идея двусторонней аппроксимации, которая при учете погрешностей приводит к необходимости обобщения понятия вещественного числа до интервального числа. Непосредственное применение интервальных методов в вычислительных процессах позволяет заключить в интервалы решения задач, о входных данных которых известно лишь то, что они лежат в определенных пределах. При этом в получаемые интервалы включаются и встречающиеся в процессе вычислений ошибки округлений. При точно определенных входных данных задачи получаемые интервалы содержат точное решение исходной задачи, и интервальный метод служит для учета ошибок аппроксимации и округлений.
Основная идея интервального анализа состоит в замене арифметических операций и вещественных функций над вещественными числами интервальными операциями и функциями, преобразующими интервалы, содержащие эти числа. Ценность интервальных решений заключается в том, что они содержат точные решения исходных задач. [2] (Добронец Б.С.)
Двусторонние методы численного анализа известны раньше интервальных, и их аппарат до последнего времени не использовал понятия интервального анализа. Для получения двусторонних оценок применяются различные приемы и методы, как то: операторные неравенства, априорные и апостериорные оценки погрешности, оценки остаточных членов, и т.п. Следует отметить, что двусторонние методы несколько проще интервальных в реализации, однако обладают некоторыми недостатками. В частности, в них, как правило, используются точные входные данные, не учитывается влияние погрешностей, связанных с машинным округлением, некоторые приемы гарантируют включение точного решения только при достаточно малых шагах сетки. Однако, двусторонние методы довольно просто обобщаются с обыкновенных дифференциальных уравнений на уравнения эллиптического и параболического типа. Применение интервального анализа к уравнениям в частных производных встречает большие трудности, а результаты его часто неудовлетворительны. [4] (Шарый С.П.)
Проблемы интервального анализа можно условно разделить на три группы:
Исследование самого множества интервальных чисел как некоторой математической структуры;
Применение интервальных методов к различным задачам прикладной математики;
Программирование интервальных методов.
История интервального анализа
Истоки интервального анализа, как и многих других глубоких и плодотворных идей, могут быть прослежены задолго до фактического оформления соответствующего научного направления. Наиболее курьёзной является на этот счёт точка зрения, высказанная изобретателем термина "интервальный анализ" американцем Р.Э. Муром [9] и состоящая в том, что первым "интервальщиком" следует считать Архимеда, широко использовавшего в своих расчётах двусторонние приближения, в частности, для определения границ числа π - отношения длины окружности к её диаметру. Так или иначе, развитие "интервальной идеи" состоялось лишь в XX веке, причём оно оказалось тесно связанным с развитием и распространением практических вычислений. А оформление интервального анализа в самостоятельную научную дисциплину вообще стало возможным лишь с появлением ЭВМ. В 1931 году англичанка Розалинда Янг [10] разработала арифметику для вычислений с множествами чисел. В 1951 году П. Двайер [11] в США рассматривал специальный случай замкнутых интервалов (числовые диапазоны) в связи с необходимостью учёта погрешностей в численном анализе. В 1956-58-м годах появились работы Мечислава Вармуса в Польше [12] и Теруо Сунаги в Японии [13], предлагавшие классическую интервальную арифметику и намечавшие её приложения. При этом впервые были использованы и современные термины "интервал", "интервальный". Кроме того, Т. Сунага заложил основы интервального алгебраического формализма и дал весьма нетривиальные примеры применений новой техники, к примеру, в численном решении алгебраических уравнений и задачи Коши для обыкновенных дифференциальных уравнений. 1959-м годом датируется начало многосторонней деятельности Раймона Э. Мура [14], написавшего в 1966-м году первую систематическую монографию по интервальному анализу [9]. Ясный и свежий язык этой книги, новые интересные постановки задач, поучительные применения интервальной техники - всё это оказало громадное влияние на становление и развитие нового научного направления во всём мире. Перу Р.Э. Мура принадлежат также современные книги по интервальному анализу [15][16], причём издание последней было задумано как переработка на современный лад классической книги [9]. В России и Советском Союзе "интервальную" историю можно отсчитывать с 20-х годов прошлого века, и связана она с именем замечательного русского советского математика и педагога В.М.Брадиса. Владимир Модестович Брадис известен своими вычислительными математическими таблицами. Но мало, кому известна другая, гораздо более серьёзная сторона его научно-педагогической деятельности. С середины 20-х годов прошлого века он работал над так называемым методом границ - способом организации вычислений, приводящим к достоверным двусторонним границам точного значения вычисляемого результата, фактически аналогичный интервальной арифметике. Работая в Тверском Педагогическом институте, он выступал с научно-методическими произведениями на эту тему[17][18][19].
В докомпьютерную эпоху "метод границ" имел мало шансов на широкое воплощение в вычислительной практике, так как требовал увеличения числа выполняемых операций и скрупулёзного их выполнения. Тем не менее, идеи В.М.Брадиса были замечены и его статья "Устный и письменный счёт. Вспомогательные средства вычислений" была внесена в "Энциклопедию элементарной математики". Параграф 8 этой статьи содержит чёткое изложение учёта погрешностей вычислений по "методу границ", который В.М.Брадис рекомендовал даже для изучения и использования в средней школе. "Энциклопедия элементарной математики" была переведена на немецкий язык в Германской Демократической Республике [20], и, кроме того, переведена и издана в Японии, став, таким образом, известной за пределами нашей страны.
В 1962-м году в одном из первых выпусков "Сибирского математического журнала" появилась статья Леонида Витальевича Канторовича [6], обозначившего эту тематику как приоритетную для нашей вычислительной науки. Она написана чрезвычайно ясно и чётко, будучи кратким конспектом лекций, прочитанных автором в Ленинграде, Москве и Киеве. Примечательно, что в этой статье, где новое научное направление ещё не называется явно, но рельефно очерчивается, акцент в приложениях нового подхода делается как на повышении точности и надёжности численных алгоритмов, так и на перспективах развития аппарата для оперирования с ограниченными неопределённостями. Создание настоящей школы исследователей по интервальным вычислениям в СССР стало заслугой других людей - выдающегося советского математика и механика, многолетнего директора Института теоретической и прикладной механики Сибирского Отделения Академии Наук академика Николая Николаевича Яненко и его ученика Юрия Ивановича Шокина. Сам Н.Н.Яненко имел немного работ по этой тематике, но лично создал интервальную группу, поддерживая новое направление у себя в Институте теоретической и прикладной механики в Новосибирске. Н.Н.Яненко был редактором первой на русском языке книги по интервальному анализу, написанной Ю.И.Шокиным. Интервальная арифметика.
Интервалом называют замкнутый отрезок вещественной оси, а интервальная неопределенность - это состояние неполного (частичного) знания об интересующей нас величине, когда известна лишь её принадлежность некоторому интервалу, т.е. когда представляется возможным лишь указать границы возможных значений величины (пределы её измерения) [Шарый. Конечномерный интервальный анализ]. Соответственно, интервальный анализ - это отрасль математического знания, исследующая задачи с интервальными неопределённостями и методы их решения. Более развернуто, можно выделить предмет интервального анализа и его специфический метод. Интервальный анализ - это раздел математики, предметом которого является решение задач с интервальными (или, более общо, ограниченными) неопределённостями и неоднозначностями в данных, возникающими либо в условиях задачи, либо в процессе её решения, чьей характеристической особенностью является рассмотрение множеств неопределённости как самостоятельных целостных объектов, посредством установления арифметический, аналитических и т.п. операций и отношений медлу ними.
Интервальный анализ и его специфичные методы имеют, таким образом, наивысшую ценность в задачах, где неопределённости и неоднозначности возникают с самого начала и являются неотъемлемой частью постановки задачи. Хотя это никоим образом не исключает других плодотворных применений интервального анализа, в частности, в задачах, формулируемых вообще без привлечения понятия интервала. Например, в последние десятилетия интервальный анализ получил широчайшее распространение в качестве основы для так называемых доказательных (достоверных, надёжных) вычислений на ЭВМ, вычислений с гарантированной точностью и т. п., несмотря на то, что в этих приложениях интервальные методы являются всего лишь вспомогательным средством для решения задач, не интервальных по своей природе. Можно выделить три основные сферы применения методов интервального анализа: • Решение практических задач, имеющих интервальную неопределённость в данных (ввод данных как интервальных величин);
• Строгий учёт ошибок округления при вычислениях с числами с плавающей точкой на цифровых ЭВМ (направленное округление границ интервала при проведении операций); • Новые подходы к решению традиционных математических задач (таких, к примеру, как задача глобальной оптимизации, глобальное доказательное решение систем нелинейных уравнений и т.п.).
Основные обозначения
Введём некоторые основные понятия и обозначения интервального анализа. [Добронец. Интервальная математика]
а - интервальное число [ а, a̅ ] с границами а ≤ a̅ , а, a̅ ∈ R, a ∈ I(R)
|а| = max (|a|, |a̅|), mid a = (a + a̅)/2, wid a = a̅ - a.
b = (bi) = (b1, ..., bn) - вектор с элементами bi, i = 1, ..., n.
b = (bi) = (b1, ..., bn) - интервальный вектор с элементами bi, i = 1, ..., n, |b| = (|bi|), mid b = (mid bi), wid b = (wid bi). [далее: Алефельд, Херцбергер: введение в интервальные вычисления]
Пусть * ∈ {+,-,*, :} - бинарная операция на множестве вещественных чисел. Если A, B ∈ I(R), то A*B = {z = a*b | a ∈ A, b ∈ B} определяет бинарную операцию на I(R). В определении полагается, что в случае деления 0 ∉ B. Результат операции над интервалами A = [a, a̅] и B = [b, b̅] может быть получен явно с помощью формул:
A + B = [a + b, a̅ - b̅];
A - B = [a - b̅, a̅ - b];
A · B = [min{a·b, a·b̅, a·b̅, a·b̅}, max{ a·b, a·b̅, a·b̅, a·b̅}];
A : B = [a·a̅]·[1/b̅, 1/b];
Их обоснованием служит тот факт, что z = f(x, y) = z * y, где * ∈ {+, -, ·, : } - непрерывная функция на компактном множестве. Следовательно, f(x, y) принимает как наименьшее и наибольшее значения, так и все прочие значения между ними. Таким образом, A * B - также замкнутый вещественный интервал.
Существующее ПО
На данный момент существует уже достаточно много решений использования интервальной математики в вычислениях на различных языках. Среди них: Интервальная библиотека boost (http://www.boost.org/libs/numeric/interval/doc/interval.htm)
Интервальная библиотека filib++ (http://www2.math.uni-wuppertal.de/wrswt)
Интервальная библиотека C-XSC (http://www2.math.uni-wuppertal.de/wrswt/xsc-sprachen.html)
Интервальная библиотека Gaol (http://sourceforge.net/projects/gaol)
проект JInterval (http://kenai.com/projects?q=JInterval) Int4Sci - интервальное расширение Scilab'а (http://www-sop.inria.fr/coprin/logiciels/Int4Sci/)
Scilab - система компьютерной математики, созданная как европейский аналог Matlab'а. Её особенностью является то, что это свободно распространяемый некоммерческий продукт с открытыми исходными кодами. Существуют его версии для Gnu/Linux, Windows 2000/XP/Vista и других платформ. INTLAB - интервальное расширение MATLAB'а (http://www.ti3.tu-harburg.de/rump/intlab/index.html)
Помимо реализации низкоуровневых операций и отношений с интервальными типами данных INTLAB содержит ряд процедур для решения интервальными методами стандартных задач численного анализа. Например, процедура VERIFYLSS предназначена для решения систем линейных уравнений (как обычных, так и интервальных), и выдаёт гарантированные двусторонние оценки решения. Крупным недостатком пакета INTLAB с некоторых пор является то, что он не бесплатен. intpakX - интервальное расширение для Maple (http://www2.math.uni-wuppertal.de/~xsc/software/intpakX)
Система компьютерной математики Sage (http://www.sagemath.org/)
- это свободно распространяемая система компьютерной математики, во многом напоминающая Maple, Mathematica и Matlab, в которой реализована интервальная арифметика переменной точности. На сайте Sage есть русскоязычная ветка с обучающими видеороликами. Программы И.Рона (http://www.cs.cas.cz/~rohn/matlab/index.html)
- эти программы, написанные в среде MATLAB / INTLAB, предназначены для решения традиционных и интервальных задач линейной алгебры и оптимизации, в частности, для внешнего оценивания множеств решений интервальных линейных систем и доказательного решения точечных задач линейного и квадратичного программирования.
Программа MCN для построения линейной регрессионной модели по данным с интервальной неопределённостью (http://www.nsc.ru/interval/Programing/MCN/mcn.tar)
- написана Н.М.Оскорбиным и А.В.Максимовым в Алтайском госуниверситете ещё в конце 80-х годов XX века для 32-разрядной MS DOS. [материалы взяты с сайта "Interval Computations": http://www.cs.utep.edu/interval-comp/]
Библиотека интервальных вычислений
Приложение написано на языке C++ в IDE wxDev-C++ (http://wxdsgn.sourceforge.net/), компилятор - g++.
Немного об IDE.
Dev-C++ - свободная интегрированная среда разработки приложений для языков программирования C/C++. В дистрибутив входит компилятор MinGW. Сам Dev-C++ написан на Delphi. Распространяется согласно GPL.
Проект поддерживается SourceForge. Основатель проекта Колин Лаплас, компания Bloodshed Software. В настоящее время актуализирована только Windows-версия.
На настоящий момент не разрабатывается, вместо него активно разрабатывается порт интерфейса Dev-C++ на wxWidgets - wxDev-C++.
wxDev-C++ является развитием проекта Dev-C++, но дополнительно содержит дизайнер форм для библиотеки разработки wxWidgets. WxDev-C++ включает все свойства Dev-C++, а также новейшую версию wxWidgets необходимую дизайнеру форм для среды быстрой разработки приложений.
Структура и принцип работы библиотеки
Для разработки выбрана парадигма объектно-ориентированного программирования, что обеспечивает создание достаточно простой и понятной архитектуры, и позволяет позднее расширять функционал без изменения уже существующего кода. Работа выполнена в виде подключаемого файла библиотеки. Сама по себе библиотека включает один основной класс - класс interval, обеспечивающий работу с интервалами теми же способами, что и с вещественными числами. Класс интервала содержит два скрытых поля вещественного типа - правую (верхнюю) и левую (нижнюю) границы интервала, и функции для работы с ними. Функции библиотеки.
Элементарные операции реализованы через перегрузку бинарных и унарных операторов (+, -, *, /, =, ==, >, <), что упрощает работу с объектами класса interval. В частности, арифметические операции реализованы таким образом, что в выражениях могут любым образом комбинироваться элементы вещественного и интервального типов. Результатом будет интервал. Кроме того, реализованы некоторые функции интервалов:
interval isin() - синус
interval icos() - косинус
interval ilog() - натуральный логарифм
В данной реализации они требуют вызова через экземпляр класса; в будущем планируется вынести их за рамки класса для приведения к привычному виду, применяющемуся с вещественными числами.
Несколько информационных функций, которые могут оказаться полезными при использовании библиотеки: long double left() - левая граница
long double right() - правая граница
long double mid() - середина интервала
long double rad() - радиус интервала
long double wid() - ширина интервала
long double mag() - магнитуда (наибольшее отклонение значения от 0)
long double mig() - мигнитуда(наименьшее отклонение значения от 0)
Полный код класса представлен в приложении А.
Библиотека работает с правильными интервалами (значение правой границы всегда больше, чем левой); это гарантируют функции записи значений интервала: interval() - конструктор по умолчанию; значения интервала инициализируются нулем;
interval(long double s) - конструктор интервала нулевой ширины; левая и правая границы интервала принимают значение s;
interval(long double sl, long double sr) - классический конструктор интервала; при этом не требуется задавать значение sr больше sl; функция сравнит их и расположит в правильном положении;
interval(long double mid, mong double wid, int flag) - конструктор, также позволяющий задать границы интервала, но другим способом: здесь указываются середина интервала, и его ширина. Третье передаваемое значение служит лишь для возможности перегрузки данной функции и не используется;
void set(long double sl, long double sr) - функция ручного задания значений интервала; так же, как и конструктор, контролирует положение границ;
void set(long double mid, mong double wid, int flag) - аналогично конструктору, с заданием ширины и середины интервала.
Такая схема работы библиотеки введена для последующего упрощения некоторых операций, в частности, умножения: операция сравнения вещественных чисел менее ресурсозатратна, по сравнению с вычислением произведения тех же чисел, часто и не одного.
Некоторые проблемы
Многие операции с интервалами могут быть выполнены несколькими способами, как, например, вычисление частного, которое, как и разность, может быть вычислено через умножение. На данном этапе работы были использованы простейшие варианты операций с интервалами, которые были упомянуты выше. В дальнейшем планируется заменить их более эффективными и точными. Во многих задачах требуется сравнение чисел, и в зависимости от результата, принимается определённое решение. С введением интервалов встаёт вопрос: как сравнивать накладывающиеся друг на друга интервалы? Однозначного ответа на этот вопрос дать невозможно, и подобные ситуации придется обходить средствами методов. Поэтому с подобными вопросами стоит разбираться уже в рамках конкретной задачи, а не интервальной математики. В данной библиотеке сравнение может даль однозначный ответ только для интервалов, значения которых не пересекаются
Стоит упомянуть также о проблеме округления, которая возникает при попытке реализации библиотеки интервальных вычислений на ЭВМ. Вещественных чисел, как известно, континуум, а память ЭВМ ограничена, и при попытке представить вещественные числа в памяти, а, тем более, выполнять операции с ними, происходит округление последнего числа мантиссы. Для обеспечения более высокой точности часто происходит округление к ближайшему числу. Интервальная арифметика призвана, прежде всего, гарантировать попадание ответа в интервал решения, поэтому такой метод округления при использовании методов интервального анализа не подходит. С другой стороны, слишком дорого в плане ресурсов для каждой арифметической операции переключать режимы округления процессора: скорость вычислений падает в десятки, а то и в сотни раз, по сравнению с округлением "к ближайшему". Для расширения интервала в обе стороны при работе с данной библиотекой должен быть предварительно установлен режим округления вниз (в стандарте C99 это позволяет сделать команда fesetround(FE_DOWNWARD)). Покажем, как это работает на примере сложения. Обозначим A = [a, a̅], B = [b, b̅] - слагаемые интервалы, Z = [z, z̅] - сумма интервалов. Используем указанную выше формулу для сложения: Z = A + B = [a + b, a̅+b̅]
Нижняя граница вычисляется, как и положено. Это позволяет сделать округление вниз. Верхняя же граница вычисляется в два этапа, согласно следующей формуле:
z̅ = -(-a̅ - b̅); Под первым этапом подразумевается значение внутри скобок. Это значение будет округлено "наружу" интервала, поскольку сам интервал по сути, отражен относительно нуля. Это хорошо видно на рис. 1.
Рис. 1. Интервалы
Здесь интервал A отмечен зеленым цветом, интервал B - синим; результирующий интервал Z - красным, а "отраженный", полученный в результате выполнения операции в скобках (обозначим Z̅') - оранжевым. Для простоты интервалам заданы конкретные значения. Максимальное для интервала Z значение 8 является для обращенного Z' минимальным. Таким образом, мы дважды находим минимум значений, а затем один из них - отражаем. Таким образом вычисляется верхнее значение суммы интервалов. Натуральный логарифм - монотонная функция, поэтому с его реализацией проблем не возникло. Левая граница интервала I переходит в левую границу результирующего интервала Ln(I). В отличие от него, синус и косинус - периодические функции. Таким образом, нельзя гарантировать того, что границами результирующего интервала будут функции границ исходного интервала. Это становится очевидно при взгляде на рис. 2. Рис. 2. Интервальные синус и косинус.
Для интервала A=[x1,x2] значение синуса должно вычисляться не по формуле [sin(x1),sin(x2)] (что, однако, верно для интервала [x3,x4]), поскольку его нижней границей будет 0. То есть sin(A) = [max(x1,x2),0]. Каждая функция, таким образом, разбивается на пять возможных вариантов: Включение максимума и минимума функции на отрезке;
Монотонное возрастание на отрезке;
Монотонное убывание на отрезке;
Включение максимума функции;
Включение минимума функции;
Первый случай отслеживается без вычисления значений синуса и косинуса при слишком большой ширине интервала: больше 2π. Отслеживать остальные случаи можно через производные: синус возрастает на отрезке, где косинус положительный, и убывает там, где косинус отрицательный. Тестирование библиотеки
Для проверки корректности работы библиотеки были проведены тесты роста ширины интервалов при выполнении элементарных операций с вырожденными интервалами. Сложение (вычитание из) интервала с числом 0.00000001 и умножение (деление) интервала на число 1.00000001.
Рис 3. Рост относительной ширины интервала при выполнении элементарных операций. По оси абсцисс - число выполненых операций, по оси ординат - относительная ширина интервала.
Такой рост ширины интервалов вызван использованием направленного округления. Однако, количество вычислений, необходимых для того, чтобы рост был заметен, весьма велико (более миллиона), а расширение интервала возникает в 12 знаке, так что такой рост погрешности не будет сильно сказываться на решении. Пример работы с библиотекой
#include "mk_interval.h" //подключение библиотеки
<...>
interval A(1.0,2.0), B=2.0, C(4.0, 1.0, 1); //В интервале А задаются обе границы; интервал В - вырожден. Интервал С задается через середину и радиус.
A = 2*A+B;
cout<<"A = "<<A.toStr()<<", C = "<<C.toStr()<<endl;
На экран будет выведено сообщение "A = [3.5;5.5], C = [3.0;5.0]".
Пример показывает простоту использования библиотеки. Для достаточно простых вычислений от пользователя не требуется никаких дополнительных знаний, что позволяет гораздо шире и эффективнее применять аппарат интервального анализа. Синтаксис максимально приближен к стандартному синтаксису языка С для других типов данных.
Глава 2. Моделирование электрических полей.
Постановка задачи
Рис. 4. Область поиска решения
Задача для потенциала u(p) в области ("полукруг" с двумя внутренними окружностями):
, , (1)
, (2)
. (3)
. (4) Решение задачи производится в три этапа:
Сведение к граничному интегральному уравнению;
Переход к интервальной системе линейных алгебраических уравнений методом конечных сумм;
Решение полученной ИСЛАУ одним из стандартных методов.
Обезразмеривание
R - характерный размер; - характерный потенциал. Безразмерные величины задаются следующими формулами:
, , ; (5)
, , ; (6)
, . (7)
Задача в безразмерных величинах (знак "тильда" опущен):
, , (8)
, (9)
, (10)
. (11)
Оразмеривание
Переход к размерным величинам осуществляется по следующим формулам:
, , ; (15)
, ; (16)
, . (17)
Задача (1) - (3) описывает распределение потенциала в двумерном сечении труб при применении катодной защиты протяженных сооружений от коррозии.
Сведение краевой задачи (8)-(11) к граничному интегральному уравнению
Для двух функций u, v∈C^2 (Ω), непрерывных и имеющих непрерывные первые производные вплоть до границы Г= ∂Ω конечной области Ω, вторая формула Грина применительно к оператору Лапласа имеет вид:
∫_Ω▒〖(v∆u- u∆v)dΩ= ∫_Г▒〖(v ∂u/∂n-u ∂v/∂n )dГ〗〗(18) Где n - внешняя нормаль к границе Г.
Возьмём в качестве v функцию E - фундаментальное решение уравнения Лапласа:
E(x, ε)=ln 1/r(19) Где r=|x- ε| расстояние между точками x и ε.
Тогда получаем:
ωu(x)= ∫_Г▒〖u(ε)(∂E(x, ε))/∂v〗 dГ- ∫_Г▒〖E(x, ε) ∂u(ε)/∂v〗 dГ+ ∫_Ω▒E(x, ε)∆u(ε)dΩ
(20)
Где величина ω определяется формулой Гаусса
ω= ∫_Г▒〖∂E(x, ε)dГ= {█(-|S_1 |, x∈Ω^+@-|S_1 |/2, x∈Г@ 0, x∈Ω^- )┤ 〗
(21)
Интегральная формула Грина для двумерной области с учетом (8), нормаль - внешняя:
. (22)
Из (9) - (11) выразим производные по нормали:
, . (23) Для граничной задачи (8) - (11) на основе формулы (22) и с учетом соотношений (23) построим граничное интегральное уравнение (далее p и q опущены):
.
Или:
.
Последнее уравнение можно записать в следующем виде:
, (24)
Где (25)
Переход к интервальной системе линейных алгебраических уравнений
От полученного интегрального уравнения (24) можно перейти к ИСЛАУ с помощью метода конечных сумм. Метод основывается на приближенном вычислении интеграла с помощью некоторой квадратурной формулы. ∫_a^b▒〖F(x)dx= ∑_(i=1)^n▒〖A_i F(x_i )+R[F]〗〗(26)Где x_i - абсциссы точек отрезка [a, b], A_i - числовые коэффициенты, не зависящие от выбора функции F(x).
Выберем точки p_(i ), q_i∈Г и введём обозначения:
u(p_i )=u_i, K(p_i,q_j )=K_ij, F(p_i )= F_i(27)Тогда, с учётом квадратурной формулы из (22) получаем систему алгебраических уравнений:
πu_i+ ∑_(i=1)^N▒〖u_j A_i K_ij 〗= F_i(26)Полученная система далее решается одним из стандартных математических методов.
Численные методы решения СЛАУ
Метод Гаусса
Метод Гаусса - классический метод решения системы линейных алгебраических уравнений. Это метод последовательного исключения переменных, когда с помощью элементарных преобразований система уравнений приводится к равносильной системе треугольного вида, из которой последовательно, начиная с последних (по номеру), находятся все переменные системы. Рассмотрим метод подробнее. Он условно разделается на два этапа: прямой ход, когда формируется верхнетреугольная матрица коэффициентов, и обратный, в течение которого находятся, собственно, решения системы. Пусть исходная система выглядит следующим образом
Матрица называется основной матрицей системы, - столбцом свободных членов.
Тогда, согласно свойству элементарных преобразований над строками, основную матрицу этой системы можно привести к ступенчатому виду (эти же преобразования нужно применять к столбцу свободных членов):
При этом будем считать, что базисный минор (ненулевой минор максимального порядка) основной матрицы находится в верхнем левом углу, то есть в него входят только коэффициенты при переменных .
Тогда переменные называются главными переменными. Все остальные называются свободными.
Если хотя бы одно число , где , то рассматриваемая система несовместна, т.е. у неё нет ни одного решения.
Пусть для любых .
Перенесём свободные переменные за знаки равенств и поделим каждое из уравнений системы на свой коэффициент при самом левом (, где - номер строки):
, где Если свободным переменным системы (2) придавать все возможные значения и решать новую систему относительно главных неизвестных снизу вверх (то есть от нижнего уравнения к верхнему), то мы получим все решения этой СЛАУ. Так как эта система получена путём элементарных преобразований над исходной системой (1), то по теореме об эквивалентности при элементарных преобразованиях системы (1) и (2) эквивалентны, то есть множества их решений совпадают.
Метод вращений
Метод вращений решения систем линейных уравнений, в отличие от метода Гаусса, не допускает роста элементов матрицы. В методе Гаусса при решении систем уравнений, последовательно сокращается количество переменных (при приведении матрицы к треугольному виду), тем самым увеличивается разрядность элементов. При большой размерности системы возникают заметные погрешности округления. Метод вращений позволяет сократить погрешности вычислений путём нормирования элементов матрицы. Рассмотрим систему уравнений
(1) Преобразуем первые два уравнения системы (1) следующим образом:
Первое: (с_1*A_11+s_1*A_21 )*x_1+(c_1*A_12+s_1*A_22 )*x_2+...+(c_1*A_1n+s_1*A_2n )*x_n=c_1*b_1+s_1+b_2
Второе: (-s_1*A_11 + c_1*A_21 )*x_1 + (-s_1*A_12+c_1*A_22 )*x_2+...+(-s_1*A_1n+c_1*A_2n )*x_n=-s_1*b_1+c_1*b_2
Где s1- синус угла φ, c1-косинус угла φ. Их можно посчитать по формулам:
(2) На значения с1, s1 накладываются следующие ограничения: 1. -s1*A11+c1*A21=0 - исключения переменной x1 из второго уравнения 2. c12+s12=1 - нормировка элементов. Заменяем выражение (c1*A11+s1*A21) на A11(1), а (c1*A12+s1*A22) на A12(1) и т.д. Система уравнений (1) примет следующий вид:
(3) После того как мы сократили переменную x1 из второго уравнения, сократим переменную x1 из третьего уравнения. Изменим значения с, s согласно формулам (4).
(4) Преобразуем первое и третье уравнения системы (рис. 3) таким образом: первое: (c_2*A_11^1+s_2*A_31 )*x_1+(c_2*A_12^1+s_2*A_32 )*x_2+...+(c_2*A_1n^1+s_2*A_3n)*x_n= c_2*b_1^1+s_2*b_3
Третье:
(-s_2*A_11^1 + c_2*A_31 )*x_1 + (-s_2*A_12^1+c_2*A_32 )*x_2+...+(-s_2*A_1n^1+c_2*A_3n)*x_n=-s_2*b_1^1+c_2*b_3.
Заменим выражение (c2*A111+s2*A31) на A112 ,а (c2*A121+s2*A32) на A122 и т.д. Система уравнений (3) примет следующий вид:
(5) Далее подобным образом сократим переменную x1 из всех остальных уравнений. Затем переходим ко второму уравнению и сокращаем переменную x2 с 3-его по n-ое уравнение. Продолжаем так до тех пор, пока не получим треугольный вид системы уравнений:
(6) Далее, используя обратный ход метода Гаусса, можно получить значение неизвестных xn, xn-1,..x1.
Проблема численного решения интервальной СЛАУ.
Характерной чертой численного решения интервальных систем линейных алгебраических уравнений является экспоненциальный рост ширины интервала в зависимости от размерности системы. Проведены несколько тестов, в ходе которых выявлена данная зависимость.
Рис. 5. Зависимость ширины интервала от размерности матрицы.
Данная проблема не позволяет решать системы больших размерностей. Интервалы, разумеется, содержат решение, но их ширины делают такое "решение" бессмысленным. В задачах, где допускается малое число разбиений, можно применять классический метод Гаусса, который показал хорошие результаты для задачи катодной защиты в области с числом разбиений границы менее 120.
Неоднородность ширины интервала решения
Рассмотрим график напряжения на аноде для 52 точек разбиения. Уже здесь явно заметно сильное расширение интервала вначале, и его сужение к концу (остальные графики и данные можно найти в приложении В). Рис. 6. Полученный график распределения потенциала на аноде
Это можно объяснить издержками метода (в данном случае использован метод Гаусса) вкупе с особенностями интервалов. После выполнения прямого и обратного хода с элементами матрицы, стоящими ближе к началу, производится гораздо большее число операций, чем с теми, что расположены ниже их. Все элементы матрицы являются интервальными величинами, а интервалы, как было отмечено выше, сильно расширяются, если участвуют в большом количестве операций. Решена данная проблема может быть заменой метода - на другой, либо его оптимизация в пользу уменьшения числа операций с элементами матрицы.
Заключение
Целью данной работы была поставлена разработка программного комплекса, использующего средства интервального анализа, для последующего использования в различных задачах, требующих аппарата интервального анализа, и решения с его помощью задачи моделирования электрических полей катодной защиты протяженных сооружений. В рамках данной работы разработана и протестирована библиотека интервальных вычислений, пригодная для решения достаточно широкого класса задач; были выявлены и решены некоторые проблемы использования интервалов в вычислениях, как ошибки округления и рост ширины интервалов при большом числе операций.
Возможности библиотеки продемонстрированы на примере решения задачи для потенциала в двумерной области. Использованная литература:
Алефельд Г., Херцбергер Ю. Введение в интервальные вычисления. - Москва: Мир, 1987. Добронец Б.С. Интервальная математика. - Красноярск: Издательство КГУ, 2004. Калмыков С.А., Шокин Ю.И., Юлдашев З.Х. Методы интервального анализа. - Новосибирск: Наука, 1986. Шарый С.П. Конечномерный интервальный анализ. - XYZ, 2010. Шокин Ю.И. Интервальный анализ. - Новосибирск: Наука, 1981. Амосов А.А., Дубинский Ю.А., Копченова Н.В. Вычислительные методы для инженеров - Москва: ВШ, 1994.
Болотнов А.М. Методы граничных элементов в расчетах электрических полей электрохимических систем. Уфа: РИО БашГУ, 2002.
Демидович Б.П., Марон И.А., Шувалова Э.З. Численные методы анализа. Москва: Наука, 1967. Moore R.E. Interval analysis. - Englewood Cliffs: Prentice Hall, 1966. Young R.C. Algebra of many-valued quantities // Mathematische Annalen. - 1931. - Bd. 104. - S. 260-290
Dwyer P.S. Linear Computations. - New York: John Wiley & Sons, 1951. Warmus M. Calculus of approximations // Bull. Acad. Polon. Sci. - 1956. - Cl. III, vol. IV, No. 5. - P. 253-259.
Sunaga T. Theory of an interval algebra and its application to numerical analysis // RAAG Memoirs. - 1958. - Vol. 2, Misc. II. - P. 547-564. Moore R.E. Automatic error analysis in digital computation // Technical report LMSD84821 of Lockheed Missiles and Space Division. - Sunnyvale: Lockheed Corp., 1959. Moore R.E. Methods and applications of interval analysis. - Philadelphia: SIAM, 1979.
Moore R.E., Kearfott R.B., Cloud M.J. Introduction to interval analysis. - Philadelphia: SIAM, 2009.
Брадис В.М. Опыт обоснования некоторых практических правил действий над приближёнными числами // Известия Тверского педагогического института. - 1927. - Вып. 3. Брадис В.М. Теория и практика вычислений. Пособие для высших педагогических учебных заведений. - Москва: Учпедгиз, 1937.
Брадис В.М. Средства и способы элементарных вычислений. - Москва: Издательство Академии педагогических наук РСФСР, 1948. Enzyklopaedie der Elementarmathematik. Band I "Arithmetik". Dritte Auflage / Grell H., Maruhn K., Rinow W., eds. - Berlin: VEB Deutscher Verlag der Wissenschaften, 1966. Канторович Л.В. О некоторых новых подходах к вычислительным методам и обработке наблюдений // Сибирский Математический Журнал. - 1962. - Т. 3, No. 5. - С. 701-709. Приложения
Приложение А. Исходный код библиотеки интервальных вычислений
#include <math.h>
#include <iostream>
#include <sstream>
#include <string>
#pragma once
#ifndef MK_Interval
#define MK_Interval
//all intervals is right
#define min(a,b) (a>b)?b:a
#define max(a,b) (a>b)?a:b
using namespace std;
class interval
{
protected:
long double l;
long double r;
public:
interval(){
l=0; r=0;
};
interval(long double s){
l=s; r=s;
};
interval(long double sl, long double sr){
l = min(sl,sr);
r = max(sl,sr);
};
interval(long double mid, long double wid, int flag){
l = mid-wid;
r = mid+wid;
};
void set(long double sl, long double sr){
l = min(sl,sr);
r = max(sl,sr);
}
void set(long double mid, long double wid, int flag){
l = mid-wid;
r = mid+wid;
};
//====================
long double left(){
return l;
}
long double right(){
return r;
} string toStr(){
ostringstream s;
s<<"["<<l<<";"<<r<<"]";
return s.str();
}
//====================Operators
friend interval operator +(const long double a, const interval i)
{interval tmpi;
long double tmpd = -a - i.r;
tmpi.l = a + i.l;
tmpi.r = -tmpd; return tmpi;};
interval operator +(const interval i)
{interval tmpi;
long double tmpd = -r - i.r;
tmpi.l = l + i.l;
tmpi.r = -tmpd; return tmpi;};
interval operator +(const long double n)
{interval tmpi;
long double tmpd = -r - n;
tmpi.l = l + n;
tmpi.r = -tmpd;
return tmpi;};
friend interval operator -(const long double a, const interval i)
{interval tmpi;
long double tmpd = -a + i.l;
tmpi.l = a - i.r;
tmpi.r = -tmpd;
return tmpi;};
interval operator -(const interval i)
{interval tmpi;
long double tmpd = -r + i.l;
tmpi.l = l - i.r;
tmpi.r = -tmpd;
return tmpi;};
interval operator -(const long double n)
{interval tmpi;
long double tmpd = -r + n;
tmpi.l = l - n;
tmpi.r = -tmpd;
return tmpi;};
interval operator -()
{interval tmpi;
tmpi.l = -r;
tmpi.r = -l;
return tmpi;};
interval operator *(const interval i)
{interval tmpi;
long double m[16];//multiplicate results
m[0] = l * i.l;
m[1] = l * i.r;
m[2] = r * i.l;
m[3] = r * i.r;
m[4] = l * (-i.l);
m[4] = -m[4];
m[5] = l * (-i.r);
m[5] = -m[5];
m[6] = r * (-i.l);
m[6] = -m[6];
m[7] = r * (-i.r);
m[7] = -m[7];
long double L1, L2, R1, R2;
L1 = min(m[0], m[1]);
L2 = min(m[2], m[3]);
R1 = max(m[4], m[5]);
R2 = max(m[6], m[7]); tmpi.l = min(L1, L2);
tmpi.r = max(R1, R2);
return tmpi;
};
interval operator *(const long double a)
{interval tmpi;
long double m[4];
m[0] = l*a;
m[1] = r*a;
m[2] = l * (-a);
m[2] = -m[2];
m[3] = r * (-a);
m[3] = -m[3];
tmpi.l = min(m[0],m[1]);
tmpi.r = max(m[2],m[3]);
return tmpi;};
friend interval operator *(const long double a, const interval i)
{interval tmpi;
long double m[4];
m[0] = i.l*a;
m[1] = i.r*a;
m[2] = i.l * (-a);
m[2] = -m[2];
m[3] = i.r * (-a);
m[3] = -m[3];
tmpi.l = min(m[0],m[1]);
tmpi.r = max(m[2],m[3]);
return tmpi;};
interval operator /(const interval i)
{interval tmp, i1, i2;
i1.l = l;
i1.r = r;
i2.l = 1/i.r;
long double tmpd = 1/(-i.l);
i2.r = -tmpd;
tmp = i1 * i2;
return tmp;
};
interval operator /(const long double a)
{interval tmp, i1, i2;
i1.l = l;
i1.r = r;
i2.l = 1/a;
long double tmpd = 1/(-a);
i2.r = -tmpd;
tmp = i1 * i2;
return tmp;
};
friend interval operator /(const long double a, const interval i)
{interval tmp, i1, i2;
i1 = a;
i2.l = 1/i.r;
long double tmpd = 1/(-i.l);
i2.r = -tmpd;
tmp = i1 * i2;
return tmp;
};
interval operator =(const interval i)
{l=i.l;
r=i.r;
};
interval operator =(const long double d)
{l=d; r=d;
};
//====================Сompare
bool operator >(const interval i)
{
if (l>i.r) return true;
else return false;
}
bool operator <(const interval i)
{
if (r<i.l) return true;
else return false;
}
bool operator >=(const interval i)
{
if (l>=i.r) return true;
else return false;
}
bool operator <=(const interval i)
{
if (r<=i.l) return true;
else return false;
}
bool operator ==(const interval i)
{
if ((r==i.r)&&(l==i.l)) return true;
else return false;
}
interval isin(){
interval result;
long double diff=r-l, R, L;
if (diff>=2*M_PI){
R=1; L=-1;}
long double cr = cos(r);
long double cl = cos(l);
long double sr = sin(r);
long double sl = sin(l);
if (diff>=M_PI){
if (cl>=0 && cr<=0){//top
R=1;
L=min(sl,sr);
}
if (cl<=0 && cr>=0){//bottom
R=max(sl,sl);
L=-1;
}
if ((cl>=0 && cr>=0)||(cl<=0 && cr<=0)){//both
R=1; L=-1;
}
}
if (diff<=M_PI){
if (cr>=0 && cl>=0){
L=sl;
R=sr;
}
if (cr<=0 && cl<=0){
L=sr;
R=sl;
}
if (cl>=0 && cr<=0){
R=1;
L=min(sl,sr);
}
if (cl<=0 && cr>=0){
L=-1;
R=max(sl,sl);
}
}
result.set(L,R);
return result;
};
interval icos(){
interval result;
long double diff=r-l, R=0, L=0;
if (diff>=2*M_PI){
R=1; L=-1;}
long double cr = cos(r);
long double cl = cos(l);
long double sr = sin(r);
long double sl = sin(l);
if (diff>=M_PI){
if (sl>=0 && sr<=0){//top
R=max(cl,cr);
L=-1;
}
if (sl<=0 && sr>=0){//bottom
R=1;
L=min(cl,cr);
}
if ((sl>=0 && sr>=0)||(sl<=0 && sr<=0)){//both
R=1; L=-1;
}
}
if (diff<=M_PI){
if (sr>=0 && sl>=0){
L=cr;
R=cl;
}
if (sr<=0 && sl<=0){
L=cl;
R=cr;
}
if (sl<=0 && sr>=0){//top
L=min(cl,cr);
R=1;
}
if (sl>=0 && sr<=0){//bottom
R=max(cl,cr);
L=-1;
}
}
result.set(L,R);
return result;
};
interval ilog(){
interval result(log(l),log(r));
return result;
};
long double mid(){
return (r + l)/2.0;
};
long double rad(){
return ((r - l) /2.0);
};
long double wid()
{
return (r - l);
};
long double mag(){
long double L=fabs(l), R=fabs(r);
return max(L, R);
};
long double mig(){
long double L=fabs(l), R=fabs(r);
if (((l>=0)&&(r<=0))||((l<=0)&&(r>=0)))
return 0;
return min(L, R);
};
};
#endif
Приложение Б. Исходный код программы.
#include <cstdlib>
#include <iostream>
#include <iomanip>
#include <math.h>
#include <fstream>
#include <string.h>
#include "mk_interval.h"
using namespace std;
typedef long double Rt;
typedef interval It;
// CONSTANTS
const int Ser = 1; const char war[] = "01";
It Xar( 10.0, 1e-12, 1), Yar(-10.0, 1e-12, 1), Rar(1.0, 1e-12, 1), Ca(1.0, 1e-12, 1), Uar=1.0; const int Na=14;
It Xkr(-10.0, 1e-12, 1), Ykr(-10.0, 1e-12, 1), Rkr(1.0, 1e-12, 1), Ck(1.0, 1e-12, 1), Ukr=0.0; const int Nk=14;
It Sigma(0.8, 1e-9, 1), Rir=200.0; const int N0=10, Ni=14;
const int
N1=Na, N2=N1+Nk, N3=N2+N0, N=N3+Ni;
It V0=Uar-Ukr, Ua=1.0, Uk=0.0, Por=Rir*100,
Xa=Xar/Por, Ya=Yar/Por, Ra=Rar/Por, K1=Ca*Sigma/Por,
Xk=Xkr/Por, Yk=Ykr/Por, Rk=Rkr/Por, K2=Ck*Sigma/Por,
Ri=Rir/Por,
_Pi=3.14159265358979323846,One=1.0,
_Pi2=_Pi*2, _1_2Pi=One/_Pi2;
// Коэффициенты и координаты Гаусса
It Ag[] = {0.41795918367346938776,
0.12948496616886969327, 0.27970539148927666790,
0.38183005050511894495, 0.38183005050511894495,
0.27970539148927666790, 0.12948496616886969327};
It Kg[] = {0.0,
-0.94910791234275852453, -0.74153118559939443986,
-0.40584515137739716691, 0.40584515137739716691,
0.74153118559939443986, 0.94910791234275852453};
// VARIABLES
It U[N],J[N];
It A[N][N],B[N];
It ds[N],Gradus[N];
It X[N][7];
It Y[N][7],Cnx[N][7],Cny[N][7];
char Tgran[N];//Char
unsigned char Ngran[N];//byte
ofstream F1,F2,F3,F4,F5;
unsigned char code;
It Ia,Ik,Ua0,Uk0;
void Geometr();
It Green2(It Xp,It Yp,It Xq,It Yq);
It dGreen2(It Xp,It Yp,It Xq,It Yq,It CosNX,It CosNY);
void OutGreen(int p,int q, It *rG,It *rDg);
void Matriza();
void SumCurrent();
void Output();
void Gauss(It A[N][N], It B[N], It X[N]);
using namespace std;
int main(int argc, char *argv[])
{
char str[20];
strcpy(str, "Da_"); strcat(str, war); strcat(str, ".dat");
F1.open(str, fstream::out);
strcpy(str, "UJa"); strcat(str, war); strcat(str, ".dat");
F2.open(str,fstream::out);
strcpy(str, "UJk"); strcat(str, war); strcat(str, ".dat");
F3.open(str,fstream::out);
strcpy(str, "UJ0"); strcat(str, war); strcat(str, ".dat");
F4.open(str,fstream::out);
strcpy(str, "UJI"); strcat(str, war); strcat(str, ".dat");
F5.open(str,fstream::out);
cout << "Files opened."<<endl;
Geometr();
cout << "Geometr() done"<<endl;
Matriza();
cout << "Matriza() done"<<endl;
SumCurrent();
cout << "SumCurrent() done"<<endl;
Output();
cout << "Out1() done"<<endl;
cout << "Ready"<<endl;
F1.close();
F2.close();
F3.close();
F4.close();
F5.close();
cout << "Press the enter key to continue ...";
cin.get();
return EXIT_SUCCESS;
}
void Geometr(){
int m,k;
It Fi,Fik,Fm,dx,dm,dim;
Fi =_Pi2/Na; Fm =Fi*0.5; Fik =Fi; // Анод в грунте
for (m=0; m<N1; m++){ Gradus[m]=(Fik*180)/_Pi;
ds[m]=Ra*Fi;
for (k=0; k<7; k++){
It tmp = Fik+Kg[k]*Fm;
X[m][k]=Xa+Ra*tmp.icos();
Cnx[m][k]=(Xa-X[m][k])/Ra;
Y[m][k]=Ya+Ra*tmp.isin();
Cny[m][k]=(Ya-Y[m][k])/Ra;
}
Tgran[m]= 'A'; Ngran[m]=1; Fik=Fik+Fi;
}
Fi=_Pi2/Nk; Fm=Fi*0.5; Fik=Fi; // Катод в грунте
for (m=N1; m<N2; m++){
Gradus[m]=(Fik*180)/_Pi;
ds[m]=Rk*Fi;
for (k=0; k<7; k++){
It tmp = Fik+Kg[k]*Fm;
X[m][k]=Xk+Rk*tmp.icos();
Cnx[m][k]=(Xk-X[m][k])/Rk;
Y[m][k]=Yk+Rk*tmp.isin();
Cny[m][k]=(Yk-Y[m][k])/Rk;
}
Tgran[m]= 'K'; Ngran[m]=2; Fik=Fik+Fi;
}
dx=(Ri*2)/N0; dm=dx*0.5; dim=Ri-dm; // Верхняя граница
for (m=N2; m<N3; m++){
Gradus[m]=dim;
ds[m]=dx;
for (k=0; k<7; k++){
X[m][k]=dim+Kg[k]*dm;
Cnx[m][k]=0.0;
Y[m][k]=0.0;
Cny[m][k]=1.0;
}
Tgran[m]='0'; Ngran[m]=0; dim=dim-dx;
}
Fi=_Pi/Ni; Fm=Fi*0.5; Fik=_Pi+Fm; // Фиктивная граница
for (m=N3; m<N; m++){
Gradus[m]=(Fik*180)/_Pi;
ds[m]=Ri*Fi;
for (k=0; k<7; k++){
It tmp = Fik+Kg[k]*Fm;
X[m][k]=Ri*tmp.icos();
Cnx[m][k]=X[m][k]/Ri;
Y[m][k]=Ri*tmp.isin();
Cny[m][k]=Y[m][k]/Ri;
}
Tgran[m]='I'; Ngran[m]=3; Fik=Fik+Fi;
}
}
It Green2(It Xp,It Yp,It Xq,It Yq){
It Rx,Ry,R2;
Rx=Xp-Xq; Ry=Yp-Yq;
R2=Rx*Rx+Ry*Ry;
return (R2.ilog())*(-0.5);
}
It dGreen2(It Xp,It Yp,It Xq,It Yq,It CosNX,It CosNY){
It Rx,Ry,R2;
Rx=Xp-Xq; Ry=Yp-Yq;
R2=Rx*Rx+Ry*Ry;
return (CosNX*Rx+CosNY*Ry)/R2;
}
void OutGreen(int p,int q, It *rG,It *rDg){
It Sum; int k;
if (p==q){
It tmp = ds[q]*0.5;
*rG=ds[q]*(One-tmp.ilog());
}
else {
Sum=0.0;
for (k=0; k<7; k++)
Sum=Sum+Ag[k]*Green2(X[p][0],Y[p][0],X[q][k],Y[q][k]);
*rG=ds[q]*Sum*0.5;
}
if (Ngran[p]==Ngran[q]){
if (Tgran[q]=='A') *rDg=(ds[q]*0.5)/Ra;
else if (Tgran[q]=='K') *rDg=(ds[q]*0.5)/Rk;
else if (Tgran[q]=='I') *rDg=(ds[q]*(-0.5))/Ri;
else *rDg=0.0;
}
else {
Sum=0.0;
for (k=0; k<7; k++) Sum=Sum+Ag[k]*dGreen2(X[p][0],Y[p][0],X[q][k],Y[q][k],Cnx[q][k],Cny[q][k]);
*rDg=ds[q]*Sum*0.5;
}
}
void Gauss(It A[N][N], It B[N], It X[N]){
It C,r,g;
for(int i = 1; i < N; i++){
for(int k = i; k < N; k++){
C = A[k][i-1]/A[i-1][i-1];
for(int j = 0; j < N; j++){
A[k][j] = A[k][j] - A[i-1][j]*C;
}
B[k] = B[k] - B[i-1]*C;
}
}
It tmp(0, A[N/2][N/2].wid()/2);
for (int k = N - 1; k >= 0; k--){
r = 0;
tmp = 0;
for(int j = k+1; j < N; j++){
g=A[k][j]*X[j];
r=r+g;
}
X[k]=(B[k]-r)/A[k][k];
}
}
void Matriza(){
It G,dG,SumA;
int p,q;
for (p=0; p<N; p++){
SumA=0.0;
for (q=0; q<N; q++){
OutGreen(p,q,&G,&dG);
if (Tgran[q]=='A'){
A[p][q]=G/K1+dG;
SumA=SumA+G;
}
else if (Tgran[q]=='K') A[p][q]=G/K2+dG;
else A[p][q]=dG;
}
A[p][p]=A[p][p]+_Pi;
B[p]=SumA*Ua/K1;
}
Gauss(A,B,U);
}
void SumCurrent(){
int k;
It Sr,U0;
for (k=0; k<N; k++)
U[k]=U[k]*V0+Ukr; // потенциал - в размерный вид
Ia=0.0; // Ток на аноде
for (k=0; k<N1; k++){
J[k]=(Uar-U[k])/Ca;
Ia=Ia+J[k]*ds[k]*Por;
}
Ik=0.0; // Ток на катоде
for (k=N1; k<N2; k++){
J[k]=(Ukr-U[k])/Ck;
Ik=Ik+J[k]*ds[k]*Por;
}
for (k=N2; k<N; k++)
J[k]=0.0; // Ток на внешней границе
}
void Output(){
int i;
It Ugol,Xr0;
F1<<"D02: Серия \t"<<Ser<<"\tВариант \t"<<war<<endl; F1<<"Xa =\t"<<Xar.toStr()<<"\tXk =\t"<<Xkr.toStr()<<endl;
F1<<"Ya =\t"<<Yar.toStr()<<"\tYk =\t"<<Ykr.toStr()<<endl;
F1<<"Ra =\t"<<Rar.toStr()<<"\tRk =\t"<<Rkr.toStr()<<"\tRi =\t"<<Rir.toStr()<<endl;
F1<<"Ca =\t"<<Ca.toStr()<<endl;
F1<<"Ck =\t"<<Ck.toStr()<<endl;
F1<<"Sg =\t"<<Sigma.toStr()<<endl;
F1<<"Ua =\t"<<Uar.toStr()<<"\tUk =\t"<<Ukr.toStr()<<"\tPor =\t"<<Por.toStr()<<endl;
F1<<"Na =\t"<<Na<<"\tNk =\t"<<Nk<<"\tNi =\t"<<N0<<" /"<<Ni<<endl;
F1<<"Анодный ток = "<<Ia.toStr()<<endl;
F1<<"Катодный ток = "<<Ik.toStr()<<endl;
F1<<endl;
F1<<"N = "<<N<<endl;
F1<<"Wid A = "<<Ia.wid()<<endl;
F1<<"Wid K = "<<Ik.wid()<<endl;
F1.close();
F2<<"#Ugol \t U \t\t\t J*1000"<<endl;
F2<<0.0<<" \t "<<U[N1-1].left()<<" "<<U[N1-1].right()<<" \t "<<J[N1-1].left()*1000<<" "<<J[N1-1].right()*1000<<""<<endl;
for (i=0; i<N1-1; i++)
F2<<Gradus[i].mid()<<" \t "<<U[i].left()<<" "<<U[i].right()<<" \t "<<J[i].left()*1000<<" "<<J[i].right()*1000<<" "<<endl;
F3<<"#Ugol \t U \t\t\t -J*1000"<<endl;
F3<<0.0<<" \t "<<U[N2-1].left()<<" "<<U[N2-1].right()<<" \t "<<-J[N2-1].right()*1000<<" "<<-J[N2-1].left()*1000<<" "<<endl;
for (i=N1; i<N2-1; i++)
F3<<Gradus[i].mid()<<" \t "<<U[i].left()<<" "<<U[i].right()<<" \t "<<-J[i].right()*1000<<" "<<-J[i].left()*1000<<" "<<endl;
F4<<"#X \t\t\t U"<<endl;
for (i=N2; i<N3-1; i++){
Xr0=Por*Gradus[i];
if ((Xr0>-100)&&(Xr0<100))
F4<<Xr0.mid()<<" \t "<<U[i].left()<<" "<<U[i].right()<<" "<<endl;
}
F5<<"#Ugol \t U"<<endl;
for (i=N3; i<N; i++)
F5<<Gradus[i].mid()<<" \t "<<U[i].left()<<" "<<U[i].right()<<" "<<endl;
}
Приложение В. Результаты.
Проведено три контрольных запуска программы при различном числе разбиений границы. Остальные начальные данные оставались неизменными:
Xa =[ 10; 10]Xk =[-10;-10]
Ya =[-10;-10]Yk =[-10;-10]
Ra =[1;1]Rk =[1;1]Ri =[200;200]
Ca =[1;1] Ck =[1;1] Sg =[0.8 - 1e-9;0.8 + 1e-9]
Ua =[1;1]Uk =[0;0]Por =[20000;20000]
Серия 1
Na =10Nk =10Ni =14 /26
Анодный ток = [0.596148;0.598293]
Катодный ток = [-0.597409;-0.597032]
Серия 2
Na =12Nk =12Ni =16 /30
Анодный ток = [0.595973;0.601384]
Катодный ток = [-0.599078;-0.598279]
Серия 3
Na =14Nk =14Ni =20 /34
Анодный ток = [0.597109;0.609158]
Катодный ток = [-0.603899;-0.602369]
Серия 1.
Серия 2. Серия 3.
- 55 -
Документ
Категория
Разное
Просмотров
199
Размер файла
462 Кб
Теги
мироненковкр
1/--страниц
Пожаловаться на содержимое документа