Scientific journal
Fundamental research
ISSN 1812-7339
"Перечень" ВАК
ИФ РИНЦ = 1,674

HARDWARE TARGETED IMPLICIT DIFFERENTIAL EQUATIONS SOLVERS

Butusov D.N. 1 Karimov T.I. 1 Tutueva A.V. 1
1 Saint-Petersburg State Electrotechnical University
The paper describes a family of hardware-targeted solvers for ordinary differential equations based on really-implicit Runge-Kutta numerical integration methods and Aitken-Neville extrapolation formula. The mathematical notion of three-stage L-stable really-implicit method of accuracy order 3 is given and the stability region is presented. New parallel architectures of solvers for ordinary differential equations are discussed. Several test problems were considered, including Van der Pol oscillator, chaotic Sprott Case G system and autocatalytic Robertson reaction. A performance of methods was compared by plotting computational costs via obtained accuracy for each method. The comparison of studied methods with Radau IIA and linearly-implicit method of order 3 showed the superior properties of really-implicit method parallel implementation for Sprott G and Robertson problems. The implicit Euler extrapolation solver showed better performance on Van der Pol oscillator. Both of the proposed methods are suitable for implementation in parallel computing systems.
numerical integration
implicit runge kutta method
ODE Solver
parallel computing
simulation
1. Butusov D.N., Karimov A.I., Karimov T.I., Dolgushin G.K. Semejstvo apparatno-orientirovannyh metodov chislennogo integrirovanija. Sovremennye problemy nauki i obrazovanija. 2014. no. 4.
2. Alexander R. Diagonally Implicit Runge–Kutta Methods for Stiff O.D.E.’s. SIAM Journal on Numerical Analysis. 1977. no. 6 (14). рp. 1006–1021.
3. Bulirsch R., Stoer J. Numerical treatment of ordinary differential equations by extrapolation methods. Numerische Mathematik. 1966. no. 1 (8). pр. 1–13.
4. Hairer E., N?rsett S.P., Wanner G. Solving Ordinary Differential Equations I. Berlin, Heidelberg: Springer Berlin Heidelberg, 1993.
5. Richardson L.F. The Approximate Arithmetical Solution by Finite Differences of Physical Problems Involving Differential Equations, with an Application to the Stresses in a Masonry Dam. Philosophical Transactions of the Royal Society A: Mathematical, Physical and Engineering Sciences. 1911. no. 459–470 (210). pр. 307–357.

В современной вычислительной математике выделяется направление так называемых аппаратно-ориентированных численных методов интегрирования – алгоритмов численного решения обыкновенных дифференциальных уравнений, пригодных для использования как на компьютерах общего назначения, так и во встраиваемых системах. Сформулируем некоторые из требований к таким методам [1]:

1. Метод должен быть самостартующим.

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

3. Метод должен работать с равномерной выборкой точек входного сигнала.

4. Алгоритм метода должен предусматривать возможность распараллеливания решения на вычислителе с соответствующей архитектурой.

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

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

Первое требование исключает многошаговые методы, широко применяемые в настоящее время в САПР радиоэлектронных устройств. Второму требованию не соответствуют любые полностью явные численные методы. Третье и пятое требования исключают применение большинства неявных методов Рунге – Кутты. С учетом четвертого требования на роль аппаратно-ориентированных методов претендуют:

1) неявные методы семейства Рунге – Кутты особого вида;

2) экстраполяционные методы с неявными опорными методами;

3) экстраполяционные методы с полунеявными опорными методами.

В настоящей статье рассматриваются представители первых двух классов алгоритмов: экстраполяционный метод на основе неявного метода Эйлера и трехстадийный вещественно-неявный метод Рунге – Кутты третьего порядка. Проводится практическое сравнение данных аппаратно-ориентированных методов между собой и с другими численными методами на ряде тестовых задач при последовательной и параллельной реализации.

Постановка задачи. Общий вид неявных методов Рунге – Кутты

Рассмотрим решение обыкновенного дифференциального уравнения

butusov01.wmf (1)

где t – независимая переменная (время), а f(t, y) – некоторая, в общем случае нелинейная, функция, которую мы будем называть функцией правой части (ФПЧ). Переформулируем стандартную запись методов Рунге ? Кутты, введя обозначение

butusov02.wmf

Здесь zi = ki – y0. Тогда, учитывая, что butusov03.wmf [4]:

butusov04.wmf (2)

Из (2) следует общая форма записи для s-стадийного метода Рунге ? Кутты в новых переменных z:

butusov05.wmf

butusov06.wmf, (3)

где d = bA–1, Ij – единичная матрица размерности j×j, а ⊗ – умножение Кронекера. Решение уравнения (2) осуществляется с помощью итераций Ньютона, общий вид которых

butusov07.wmf

при этом отыскивается такое значение x*, при котором F(x*) = 0.

Для схемы (3) метод Ньютона может быть записан следующим образом:

butusov08.wmf (4)

где W0 – якобиан функции F0. Затем определяется новое значение вектора переменных:

butusov09.wmf

Вещественно-неявные методы Рунге – Кутты

Идея распараллеливания методов Рунге – Кутты состоит в том, чтобы вычислять s независимых компонент вектора Δzk в параллельных потоках. Действительно, если матрица метода может быть представлена в виде

butusov10.wmf

то после замены переменных

butusov11.wmf butusov12.wmf

butusov13.wmf

итерационная схема (3) преобразуется к виду

butusov14.wmf

где Ji = ∂f ? ∂y(t0 + cih, y0 + zi) – якобиан функции правой части в точке(t0 + cih, y0 + zi). В конце каждой итерации потоки собираются, преобразуются в старую систему координат. Такие методы будем называть вещественно-неявными методами Рунге ? Кутты, ВНРК (англ. Really Implicit Runge-Kutta, RIRK). Заметим, что λi ≠ λj, i ≠ j, в противном случае матрица T сингулярная и T–1 не существует. С помощью подходов многокритериальной оптимизации авторами был синтезирован L-устойчивый метод 3-го порядка, имеющий таблицу Бутчера

butusov15.wmf (5)

собственные векторы матрицы A которого порождают матрицу

butusov16.wmf

Область устойчивости метода (5) приведена на рис. 1.

pic_1.wmf

Рис. 1. Область устойчивости вещественно-неявного метода (5)

pic_2.tif

Рис. 2. Схема параллельного решателя на основе трехстадийного метода ВНРК

На рис. 2 представлена архитектура параллельного решателя ОДУ на основе ВНРК. В случае если применяемый вычислитель поддерживает работу с комплексными числами, подобная организация вычислений теоретически возможна и для методов с комплексными элементами матрицы Λ. Однако на практике оказывается, что точность решения такими методами существенно ниже по сравнению с использованием непреобразованной матрицы A.

Экстраполяционные методы

Идея экстраполировать решения, взятые с разным шагом интегрирования, для повышения порядка точности метода лишь немного моложе классических разностных схем Рунге – Кутты. Предложенный Ричардсоном [5] способ экстраполяции получил свое развитие в работах Грэгга, Булирша, Штёра [3], а также многих других авторов. По сравнению с методами Рунге – Кутты, экстраполяционные методы, с одной стороны, требуют большего числа вычислений функции правой части на шаге. С другой стороны, такие методы обладают способностью организации вычислительных потоков с единственным обменом данных на шаге, а не с двумя, как в случае ВНРК, что сокращает накладные расходы.

Суть экстраполяционных методов заключается в следующем. Выберем ряд натуральных чисел

n1 < n2 < n3... (6)

и определим соответствующие длины шагов, относящиеся к исходному шагу H как

butusov17.wmf

После этого c помощью так называемого опорного метода вычислим значения функции y в точке t0 + H, используя различные длины шагов:

butusov18.wmf (7)

Затем применим алгоритм рациональной интерполяции Эйткена – Невилла, результатом которого является так называемая экстраполяционная таблица:

butusov19.wmf

Первый столбец таблицы формируется по формуле (7). Последующие колонки заполняются с помощью рекурсивной формулы [4]:

butusov20.wmf

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

butusov21.wmf

где N1 – число итераций Ньютона на 1-й стадии экстраполятора.

Экспериментальная оценка свойств численных методов

Оценим свойства предложенных аппаратно-ориентированных методов, рассмотрев решение популярных тестовых задач.

Задача VDPL – уравнение осциллятора Ван дер Поля

butusov22.wmf

butusov23.wmf

При моделировании используется значение μ = 4.

BRUS – модель циклической химической реакции, предложенная Лефевером и Николисом в 1971 г. и известная в литературе [4] как «брюсселятор»:

butusov24.wmf

butusov25.wmf(8)

Параметры системы (8) выбраны следующим образом: A = 2, B = 8,553.

SPRG – система уравнений Case G, предложенная Спроттом в 1994 г. и описывающая хаотическую нелинейную динамическую систему:

butusov26.wmf

butusov27.wmf

butusov28.wmf

ROBER – разноскоростная система уравнений, описывающая автокаталитическую реакцию Робертсона:

butusov29.wmf

butusov30.wmf

butusov31.wmf

pic_3.tif

Рис. 3. Организация параллельных вычислений экстраполяционным методом третьего порядка в двух или трех потоках

pic_4_1.tif

а б

pic_4_2.tif

в г

pic_4_3.tif

д е

Рис. 4. Сравнение эффективности численных методов интегрирования на задачах: BRUS при последовательной (а) и параллельной (б) реализациях, VDPL при последовательной (в) и параллельной (г) реализациях, SPRG (д) и ROBER (е) при параллельной реализации

Результаты решения тестовых задач аппаратно-ориентированными методами сравним с решениями, полученными неявным методом Эйлера, двухстадийным методом Радо IIA 3-го порядка и двухстадийным методом ОДНРК 3-го порядка, описанным Р. Александером [2]

Оба исследуемых аппаратно-ориентированных алгоритма разделены на три вычислительных потока. Метод ВНРК исследовался при 10 итерациях Ньютона, а ЭПНЭ3 – на последовательности итераций [10, 5, 3] для различных значений шага. Для оценки погрешности использовалось эталонное решение, рассчитанное методом высокого порядка с шагом в 100 раз меньшим, чем шаг тестируемых методов. Представлены графики отношения числа вычислений в единицу времени к достигаемой точности решения.

Заключение

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

Анализ экспериментальных данных показывает, что при параллельной реализации метод ВНРК почти всегда оказывается самым экономичным по числу вычислений функции правой части f(t, y). При этом на задачах Спротта и Робертсона он показывает лучшую точность. Несмотря на то, что неявный экстраполяционный метод Эйлера имеет в большинстве случаев более высокую вычислительную сложность (кроме задачи VDPL), он в то же время наименее подвержен негативному влиянию машинного шума коэффициентов.

Дальнейшим направлением исследований является исследование решения аппаратно-ориентированными методами задач высокого порядка, а также разработка методов ВНРК более высокого порядка точности с аналитической формулировкой коэффициентов.

Работа выполнена при финансовой поддержке Российского фонда фундаментальных исследований в рамках договора № НК 14-01-31277\14 от 06.02.2014 г.