The increase in the rate of convergence of finite difference method based on the use of middleware solutions

Korobeinikov Anatolii Grigor'evich Doctor of Technical Science professor, Pushkov institute of terrestrial magnetism, ionosphere and radio wave propagation of the Russian Academy of Sciences St.-Petersburg Filial 199034, Russia, g. Saint Petersburg, ul. Mendeleevskaya, 1 Korobeynikov_A_G@mail.ru Other publications by this author

Grishentsev Aleksei Yur'evich Doctor of Technical Science Associate Professor, St. Petersburg National Research University of Information Technologies, Mechanics and Optics 197101, Russia, St. Petersburg, Kronverkskiy prospect, d. 49 grishentcev@ya.ru Other publications by this author

Большое количество процессов окружающего нас мира описываются математическими моделями на базе дифференциальных уравнениями (ДУ) ^[1], значительную часть которых невозможно решить аналитически. По этой причине актуальным вопросом является разработка новых и повышение производительности имеющихся численных методов решения ДУ. Метод конечных разностей (МКР) – широко распространён и применяется для решения ДУ путём их замены конечными разностями ^[2,3]. Метод достаточно просто реализуется программно для решения ДУ в произвольном n-мерном пространстве. В работе рассмотрен способ повышения вычислительной эффективности метода конечных разностей.

Недостатком МКР является относительно невысокая скорость сходимости при значительном числе элементов подлежащих итерационным операциям (расчётам). Уменьшение числа этих элементов за счёт увеличения дискретов не всегда может пойти на пользу решения задачи, т.к. при таком подходе, во многих случаях, происходит снижение скорости сходимости к точному решению. Можно сделать выбор в пользу более быстрого метода конечных элементов (МКЭ), который так же позволяет решать ДУ и даже имеет ряд преимуществ. Например, можно увеличить число дискретизации там где необходимо повысить точность решения, а так же практически исключить эффект ступенчатости границ. При этом МКЭ в большинстве случаев обеспечивает существенно более высокую скорость решения задачи по сравнению с МКР. Сложность применения МКЭ обусловлена необходимостью разбиения пространства на конечные элементы. На сегодняшний день существует ряд алгоритмов позволяющих эффективно производить триангуляцию в n-мерных пространствах. Многие из этих алгоритмов имеются в свободном доступе в виде библиотек открытого кода. Таким образом, в некоторых случаях есть смысл отказаться от применения МКР в пользу МКЭ. Но несмотря на все преимущества, есть ряд практических особенностей ограничивающих замену МКР во всех случаях МКЭ. К таким особенностям можно отнести, то что, большинство алгоритмов и программ предназначены для триангуляции в двумерном пространстве. Для трёхмерного пространства доступные готовые инструменты триангуляции существенно скуднее. Для пространства размерностью более трёх, например, в задачах математической физики часто встречается четырёхмерное пространство (пространственно временной континуум), таких инструментов в отрытом доступе вообще не обнаружено. Разработка собственных библиотек алгоритмов для построения диаграммы Вороного в n-мерном пространстве является достаточно трудоёмким процессом и не всегда оправданна с точки зрения затрат ресурсов.

В работе представлен метод, позволяющий сократить время вычислений в МКР за счёт повышения скорости сходимости итерационного процесса, сохранив при этом простоту реализации метода для произвольного n-мерного пространства.

Для сокращения времени ^[3] вычислений МКР при решении ДУ в некоторой области W с граничными (краевыми для физических областей и начальными для времени) условиями G применяют укрупнение шага h сетки U в области W+G.При этом происходит потеря точности решения вплоть до того, что оно становится непригодным. Например, в основной задаче электростатики могут быть рассмотрены электрические заряды, описанные одной единственной точкой. При укрупнении сетки граничные условия могут быть потеряны, что меняет решение. Поэтому обычным является решение, основанное на использовании нескольких сеток: первоначально применяется самая крупная U₁, получая приближение решения, затем решение уточняют, последовательно применяя сетки с меньшим шагом U_k, k > 1. Общее число различных сеток обычно 2–3. Отметим, что, вообще говоря, шаг сетки h может быть различным в разных направлениях и областях сеткиU. Таким образом, ускорение процедуры сходимости, к решению задачи с заданной точностью p, происходит за счёт более быстрого формирования некоторого промежуточного решения в области W.

Из теоремы о сходимости ^[4] разностного решения y(x) с шагом сетки h к точному решению u(x) уравнения:

`Au(x) = f(x)` (1)

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

и условия ||y(x) - u(x)||`->` 0, следует, что если при `h->0` и заданном порядке точности p выполнено ||y(x) - u(x)||` ` `0(h^p)` , то указанный метод уменьшения шага сетки h, является прямым следствием этой теоремы. При этом каждое уточнение решения является итерационным, имеющее вычислительную сложность:

`0(k prod_(i=1)^n N_(i))` (2)

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

аппроксимации промежуточного решения в области W, по имеющимся данным о граничных условиях G, с учётом правой части (1). При этом вычислительная сложность может быть оценена:

`0(n prod_(i=1)^n N_(i))` (3)

где: обычно n<<k, причём n и k имеют тот же смысл, что и в (2).

Указанную задачу нахождения промежуточного решения, в области W, можно решить путём аппроксимации значений интеграла `f(x)` – правой части уравнения (1), последовательно вдоль линий образующих однородную сетку разбиения U, c граничными условиям G. Аппроксимацию в зависимости от вида (1), можно производить различными функциями в зависимости от физической сущности задачи, например, многочленом вида:

`a_(0) + a_(1)z + a_(2)z^2 + ... + a_(p)z^p` (4)

В общем случае, если линии сетки не параллельны координатным осям образующим пространство W, z может не совпадать с набором переменных исходной задачи. В этом случае необходимо учитывать поворот системы координат, в которых рассмотрен аргумент z, относительно исходной системы координат. В случае n – мерной (n > 1) задачи, интегрирование необходимо производить последовательно для всех линий образующих сетку в каждом направлении с последующей оценкой средних значений для каждого узла сетки. Среднее значение узла сетки вычисляется как среднее значений аппроксимации в точках линий сетки принадлежащих данному узлу. В определённом смысле можно говорить, что для нахождения некоторого промежуточного решения задача МКР разбивается на множество одномерных задач МКЭ, равное числу линий сетки МКР.

Предлагаемый метод исследовался на примере решения уравнений Лапласа ∆u = 0 в декартовой системе координат с числом измерений n = 2 и n = 3 (рис. 1) .

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

В качестве разбиения области W была выбрана прямоугольная сетка (образованная линиями параллельными координатным осям) с равным во всех направлениях шагом h. В качестве многочлена (4) было использовано уравнение прямой `a_(0) + a_(1)x` , коэффициенты которой рассчитывались по двум элементам граничных условий из G.

Представленные результаты (рис. 1) показывают существенное ускорение сходимости итерационного процесса МКР.

Применение представленного метода^[5] в некоторых задачах для достижения заданной точности решения позволило уменьшить число итераций в 10 – 100 раз, что говорит о существенном повышении эффективности МКР.

Устойчивость МКР, для ДУ, в частности для уравнения Лапласа, показана в ^[2-4].

Представленный метод способствует сходимости МКР, путём нахождения некоторого промежуточного решения y'(x). Ускорение итерационного процесса будет зависеть от значения в области решения W величины ‖y'(x) – u(x)‖. Чем меньше эта величина, тем ближе мы к точному решению, и тем меньше итераций требуется для уточнения решения (снижения невязки до некоторой заданной величины). Промежуточное решение y'(x) будет зависеть от способа аппроксимации. На Рис. 2 показаны точки A,B,C `in` G – элементы граничных условий, u(x) – точное решение, y'(x) – промежуточное решение и начальные значения y⁰(x). Области между элементами множества G, принадлежат области W – в которой необходимо отыскать решение. Начальные значения y⁰(x) задаются произвольно, обычно y⁰ = const, как в случае на изображении (рис. 2). Величина ‖y'(x) – u(x)‖=d' между точным решением и промежуточным, ‖y⁰(x) – u(x)‖=d⁰ между точным решением и начальным значением.

Строгое математическое доказательство возможности аппроксимации функции u(x) функцией y'(x), и в частности, алгебраическим полиномом (4) приведено в ^[6].

Рис. 2. Расчёт промежуточного решения

Строгое математическое доказательство возможности аппроксимации функции u(x) функцией y'(x), и в частности, алгебраическим полиномом (4) приведено в ^[6].

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