Например, Бобцов

УЛУЧШЕНИЕ СХОДИМОСТИ МЕТОДА КОНЕЧНЫХ РАЗНОСТЕЙ С ПОМОЩЬЮ ВЫЧИСЛЕНИЯ ПРОМЕЖУТОЧНОГО РЕШЕНИЯ

УЛУЧШЕНИЕ СХОДИМОСТИ МЕТОДА КОНЕЧНЫХ РАЗНОСТЕЙ …

8 СИСТЕМЫ АВТОМАТИЗИРОВАННОГО ПРОЕКТИРОВАНИЯ

УДК 519.635.4
УЛУЧШЕНИЕ СХОДИМОСТИ МЕТОДА КОНЕЧНЫХ РАЗНОСТЕЙ С ПОМОЩЬЮ ВЫЧИСЛЕНИЯ ПРОМЕЖУТОЧНОГО РЕШЕНИЯ
А.Ю. Гришенцев, А.Г. Коробейников
Рассмотрена актуальность применения метода конечных разностей для решения дифференциальных уравнений с граничными условиями. Предложен способ вычисления промежуточного решения для повышения эффективности метода конечных разностей. Приведены результаты программной реализации указанного способа. Рассмотрен принцип алгоритмической реализации способа в двумерном пространстве. Ключевые слова: метод конечных разностей, дифференциальные уравнения.
Введение
Многие процессы окружающего нас мира описываются различными дифференциальными уравнениями [1], значительную часть которых практически невозможно решить аналитически. По этой причине актуальным вопросом является разработка новых и повышение производительности имеющихся численных методов решения дифференциальных уравнений. Метод конечных разностей (МКР) широко распространен и применяется для решения дифференциальных уравнений путем их замены конечными разностями [2, 3]. Метод достаточно просто реализуется программно для решения дифференциальных уравнений в произвольном n-мерном пространстве.
Основным недостатком МКР является относительно невысокая скорость сходимости при значительном числе элементов, подлежащих итерационным операциям (расчетам). Уменьшение числа этих элементов за счет увеличения дискретности не всегда может пойти на пользу решению задачи, так как при таком подходе во многих случаях происходит снижение скорости сходимости к точному решению. Можно сделать выбор в пользу более быстрого метода конечных элементов (МКЭ), который имеет ряд преимуществ. Например, он позволяет увеличить степень дискретизации там, где необходимо повысить точность решения, а также практически исключить эффект ступенчатости границ. При этом МКЭ в большинстве случаев обеспечивает существенно более высокую скорость решения задачи по сравнению с МКР. Сложность применения МКЭ обусловлена необходимостью разбиения пространства на конечные элементы. На сегодняшний день существует ряд алгоритмов, позволяющих эффективно производить триангуляцию в n-мерных пространствах. Многие из этих алгоритмов имеются в свободном доступе в виде библиотек открытого кода. Таким образом, в некоторых случаях есть смысл отказаться от применения МКР в пользу МКЭ.
Несмотря на все преимущества МКЭ, есть ряд особенностей, ограничивающих повсеместное практическое применение МКЭ вместо МКР. Большинство алгоритмов и программ предназначено для триангуляции в двумерном пространстве. Для трехмерного пространства доступные готовые инструменты триангуляции существенно скуднее. Для пространства размерностью более трех, например, в задачах математической физики, где часто встречается четырехмерное пространство (пространственно временной континуум), таких инструментов в открытом доступе вообще не обнаружено. Разработка собственных библиотек алгоритмов для построения диаграммы Вороного в n-мерном пространстве является достаточно трудоемким процессом и не всегда оправдана с точки зрения затрат ресурсов.
Сущность предлагаемого способа
Целью работы является разработка способа сокращения времени вычислений МКР за счет повышения скорости сходимости итерационного процесса при сохранении алгоритмической простоты реализации метода для произвольного n-мерного пространства.
Известен способ сокращения времени вычислений МКР для мягких дифференциальных уравнений [3] в некоторой области W с граничными (краевыми для физических областей и начальными для времени) условиями G, сущность которого состоит в укрупнении шага h сетки U в области W+G. При этом происходит потеря точности решения вплоть до того, что оно становится непригодным. Например, в основной задаче электростатики могут быть рассмотрены электрические заряды, описанные единственной точкой. При укрупнении сетки такие граничные условия могут быть потеряны, что принципиально меняет сущность промежуточного решения.
По этой причине обычным является решение, содержащее несколько последовательно применяемых сеток. Первоначально применяется самая крупная сетка U1, позволяющая получить приближение решения. Затем решение уточняют, последовательно применяя сетки с меньшим шагом Uk, k > 1. Общее число различных сеток обычно составляет 2–3. Шаг сетки h может быть различным в разных направле-

124

Научно-технический вестник информационных технологий, механики и оптики, 2012, № 3 (79)

А.Ю. Гришенцев, А.Г. Коробейников

ниях и областях сетки U. Таким образом, ускорение процедуры сходимости к решению задачи с заданной точностью p происходит за счет более быстрого формирования некоторого промежуточного решения в области W.
В соответствии с теоремой о сходимости [4] разностного решения y(x) с шагом сетки h к точному решению u(x) для уравнения

Au(x)  f (x) ,

(1)

где A – дифференциальный оператор; x W+G (W – область решения; G – граничные условия); f(x) – за-

данная функция, из условия y(x)  u(x)  0 , при h→0 и заданном порядке точности p для

y(x)  u(x)  O(hp ) следует, что указанный метод уменьшения шага сетки h является прямым след-

ствием этой теоремы. При этом каждое уточнение решения является итерационным, имеющим вычислительную слож-
ность

 n 
O k Ni  ,  i1 

(2)

где n – размерность рассматриваемой задачи; Ni – число узлов сетки в каждом направлении; k – число итераций, обеспечивающее заданную точность на данном этапе. В многомерных задачах величина (2)

может быть очень велика даже при разреженной сетке.

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

итерационных приближениях, а на аппроксимации промежуточного решения в области W по имеющимся

данным о граничных условиях G с учетом правой части (1). При этом вычислительная сложность может

быть оценена как

 n 
O n Ni  ,  i1 
где, как правило, n  k , причем n и k имеют тот же смысл, что и в (2). Указанную задачу нахождения промежуточного решения в области W можно решить путем ап-
проксимации значений f(x) правой части уравнения (1), образующих последовательно вдоль линий однородную сетку разбиения U c граничными условиями G. Аппроксимацию f(x) можно производить различными функциями в зависимости от физической сущности задачи, например, многочленом вида

a0  a1z  a2 z2  ...  ap z p ,

(3)

где a0, a1, a2, …, ap – коэффициенты многочлена; z – переменная, вдоль которой происходит интегрирование.

В общем случае, если линии сетки не параллельны координатным осям, образующим простран-

ство W, z может не совпадать с набором переменных исходной задачи. В этом случае необходимо учиты-

вать поворот системы координат, в которых рассмотрен аргумент z , относительно исходной системы

координат. В случае n-мерной (n > 1) задачи аппроксимирование необходимо производить последова-

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

чений для каждого узла сетки. Среднее значение узла сетки вычисляется как среднее значение аппрок-

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

что для нахождения некоторого промежуточного решения задача МКР разбивается на множество одно-

мерных задач МКЭ, число которых равно числу линий сетки МКР.

Возможно совместное применение способа укрупнения сетки и предлагаемого способа, ведь даже

для укрупненной сетки процесс ускорения сходимости остается актуальным. Применение предложенно-

го способа должно производиться с учетом физической сущности процессов и условий каждой конкрет-

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

также использовать некоторые усредненные граничные условия. В этом случае фактически для получе-

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

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

сти от особенностей граничных условий), может быть «не очень» гладким. На практике возможно при-

менение некоторых дополнительных действий, повышающих гладкость. Например, при аппроксимации

возможен учет (с некоторыми весовыми коэффициентами) значений в соседних узлах, не вошедших в

качестве значения или аргумента в функцию аппроксимации, применяемую к данной линии сетки.

Апробация способа

Способ тестировался на различных задачах электротехники, описываемых уравнением Лапласа и Пуассона в пространствах с числом измерений n = 2 и n = 3:

Научно-технический вестник информационных технологий, механики и оптики, 2012, № 3 (79)

125

УЛУЧШЕНИЕ СХОДИМОСТИ МЕТОДА КОНЕЧНЫХ РАЗНОСТЕЙ …

 решение основной задачи электростатики (в том числе с зарядами, описываемыми единственным граничным элементом);
 расчет поля конденсатора сложной формы;  расчет магнитного поля в зазоре электрической машины.
Дополнительно способ тестировался на некоторых других видах задач, описываемых уравнениями в частных производных.
Во всех случаях полученное промежуточное решение позволяло существенно сократить число итераций МКР для получения решения с заданной точностью. Особый выигрыш в сокращении числа итераций наблюдается при гладких граничных условиях, вплоть до того, что при некоторой достаточной гладкости может быть получено решение только с применением указанного способа без итерационного процесса МКР.
Для разбиения области W была выбрана прямоугольная сетка, образованная линиями, параллельными координатным осям, с равным во всех направлениях шагом h. В качестве многочлена (3) было использовано уравнение прямой a0 + a1x, коэффициенты которой рассчитывались по двум элементам граничных условий из G. Полученные результаты (рис. 1) показывают существенное ускорение сходимости итерационного процесса МКР.
10000000

Логарифм суммарной величины невязки по всем узлам

1000000 100000

10000 1000

100 1

11 21 31 41 51 61 71 81 91 Номер итерации

Рис. 1. Значения логарифмов суммарныхБнеезвпяозиосккпаопврсеедмваурзилтаемльснеоткгои рвезшавеинсииямости от номера итерации в двумерном пространстСвеп:о—ис—комбепзрепдовиасркаитперлеьднвоагроиртешльенноигяо решения;
— с поиском предварительного решения

Применение рассмотренного способа в некоторых двумерных задачах позволило уменьшить число итераций в 10–100 раз, что говорит о существенном повышении эффективности МКР.
Устойчивость МКР для мягких уравнений, в частности, для уравнения Лапласа, показана в [2–4]. В результате нахождения некоторого промежуточного решения y(x) улучшается сходимость
МКР. Ускорение итерационного процесса зависит от значения величины y(x)  u(x) в области решения
W. Чем меньше y(x)  u(x) , тем меньше итераций требуется для уточнения решения (снижения невязки
до некоторой заданной величины). Промежуточное решение y(x) будет зависеть от способа аппрокси-
мации (выбора функции аппроксимации, ее размерности). На практике дополнительный эффект дает учет соседних значений, принадлежащих линиям сетки, соседним по отношению к той, для которой производится итерация.

u(x) В
y(x) y0(x) d d0 С
А
Рис. 2. Расчет промежуточного решения На рис. 2 показаны точки A, B, C G – элементы граничных условий; u(x) – точное решение;
y(x) – промежуточное решение и начальные значения y0(x). Области между элементами множества G принадлежат области W, в которой необходимо отыскать решение. Начальные значения y0(x) задаются

126

Научно-технический вестник информационных технологий, механики и оптики, 2012, № 3 (79)

А.Ю. Гришенцев, А.Г. Коробейников

произвольно, обычно y0 (x)  const , как в случае, представленном на рис. 2. Величина y(x)  u(x)  d
представляет разность между точным и промежуточным решениями, а y0 (x)  u(x)  d 0 – между точ-
ным решением и начальным значением. Строгое математическое доказательство возможности аппроксимации функции u(x) функцией
y(x) и, в частности, алгебраическим полиномом (3) приведено в [5].

Реализация в двумерном пространстве

Рассмотрим вариант реализации способа на примере двумерной задачи с граничными условиями,
описанной уравнением Лапласа, в случае применения равномерной сетки с шагом h , образованной ли-
ниями, параллельными осям координат пространства. Поиск промежуточного решения производится последовательно: сначала в направлении оси X , а затем в направлении оси Y , или наоборот. Совершим проход по первой строке, если встречаем элемент, являющийся элементом граничного условия, запоминаем его значение и двигаемся далее, пока не встретим следующий элемент, являющийся элементом граничного условия. Зная значение двух граничных элементов в пределах одного вектора, применяем аппроксимацию для вычисления значений элементов, расположенных между ними. Таким образом аппроксимируем все значения в пределах одного столбца. Последовательно проходим все столбцы в направлении X , а затем совершаем такой же проход по всем строкам в направлении Y . При проходе по строкам в направлении Y учитываем ранее полученные значения (при проходе по X ) и вычисляем фактическое значение как среднее между проходом в направлении X и Y .
Добавление проходов по диагоналям позволит получить промежуточное решение, более близкое к точному.

Заключение

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

Литература

1. Гришенцев А.Ю., Коробейников А.Г. Разработка модели решения обратной задачи вертикального зондирования ионосферы // Научно-технический вестник СПбГУ ИТМО. – 2011. – № 2 (72). – С. 109–
113. 2. Демидович Б.П., Марон И.А., Шувалова Э.З. Численные методы анализа. Приближение функций,
дифференциальные и интегральные уравнения. – СПб: Лань, 2010. – 400 с. 3. Бахвалов Н.С., Воеводин В.В. Современные проблемы вычислительной математики и математическо-
го моделирования: в 2 томах. – Т. 1. Вычислительная математика. – М.: Наука, 2005. – 343 с. 4. Формалёв В.Ф., Ревезников Д.Л. Численные методы. – М.: Физматлит, 2004. – 400 с. 5. Самарский А.А., Гулин А.В. Численные методы: Учебное пособие для вузов. – М.: Наука. Гл. ред.
физ.–мат. лит., 1989. – 432 с.

Гришенцев Алексей Юрьевич Коробейников Анатолий Григорьевич

– Санкт-Петербургский национальный исследовательский университет информационных технологий, механики и оптики, кандидат технических наук, доцент, tigerpost@ya.ru
– Санкт-Петербургский филиал Института земного магнетизма, ионосферы и распространения радиоволн им. Н.В. Пушкова Российской Академии наук (СПбФ ИЗМИРАН), доктор технических наук, профессор, зам. директора, Korobeynikov_A_G@mail.ru

Научно-технический вестник информационных технологий, механики и оптики, 2012, № 3 (79)

127