CC BY
14
УДК 681.326
ГИБРИДНЫЙ МЕТОД РЕШЕНИЯ СИСТЕМ ОБЫКНОВЕННЫХ ДИФФЕРЕНЦИАЛЬНЫХ УРАВНЕНИЙ С АВТОМАТИЧЕСКИМ ВЫБОРОМ АЛГОРИТМА ИНТЕГРИРОВАНИЯ
ЧАДОВ С.Н., асп.
Рассмотрен подход к решению систем обыкновенных дифференциальных уравнений, основанный на комбинировании классических алгоритмов. Приводится сравнение скорости работы и результатов расчетов на модельной задаче для предложенной схемы и классических алгоритмов.
Ключевые слова: адаптивные алгоритмы, модельная система, обыкновенные дифференциальные уравнения.
HYBRID METHOD OF SYSTEM SOLUTION OF DIFFERENTIAL EQUATIONS WITH AUTOMATIC SELECTION OF INTEGRATION ALGORITHM
S.N. CHADOV, Post Graduate Student
The article considers the approach of solving the systems of common differential equations. This approach is based on classical algorithms combination. The author compares the operating speed and model results for the given scheme and classical algorithms.
Key words: adaptive algorithm, model system, common differential equations.
Введение. Известно достаточно большое количество методов решения систем обыкновенных дифференциальных уравнений. Большинство из них основаны на одном из «классических» методов, приведенных в большинстве учебников по численному анализу. Однако эти «классические» методы часто сами по себе не могут обеспечить решение достаточного качества для сложных систем, поэтому в «промышленных» пакетах для решения систем дифференциальных уравнений используются часто достаточно сложные адаптивные алгоритмы с автоматическим изменением порядка и шага. Так, широко используемый пакет MATLAB в качестве основного использует алгоритм Дорманда-Принса [1].
Ниже предлагается другой подход к созданию качественных алгоритмов решения систем обыкновенных дифференциальных уравнений, который может использоваться как сам по себе, так и дополнить существующие методы. Предлагаемый подход заключается в том, чтобы использовать не автоматическую адаптацию порядка и шага метода, а автоматический выбор алгоритма интегрирования. Получившийся метод будем называть гибридным.
Модельная система уравнений. Теоретическое исследование гибридных алгоритмов затруднено, поскольку их поведение сильно зависит от характера интегрируемых функций. Поэтому будем изучать их поведение на примере модельной системы. В качестве модельной возьмем систему из [2], поскольку она достаточно проста, но при этом содержит участки, сложные для решения «классическими» методами. Также интерес представляет исследование получившегося метода решения способом, описанным в [2].
Итак, система уравнений имеет следующий вид:
% = a"V + Xn=o2nfi sin(2nf-f),y(0) = 0. (1)
Для данной системы легко получить аналитическое решение:
У = e
t / а'
i=0
Qi 2a2 qi 2a2 +1
I (-
i=0
Qi 2a2 qi 2a2 +1
cos qt
Qi
qi 2a2 +1
sin qt), (2)
График решения данной системы (на рис. 1) сочетает в себе как участки с большим значением производной, так и маленькие по амплитуде колебания.
Рис. 1. График функции (2)
Алгоритм. Рассмотрим решение модельной системы несколькими классическими методами:
1) методом Рунге-Кутты 4-го порядка [3];
2) явным методом Эйлера [3];
3) явным методом Эйлера-Коши [3];
4) неявным методом Эйлера [3];
5) с применением правила трапеций [3];
6) методом Bader и Deuflhard [4] (полунеявным).
Для каждого решения приведем время вычислений t, RMS(root mean square, см., например, [5]) разности между численным и аналитическим решением, максимальную ошибку вычислений max err. В табл. 1 приведены средние данные, округленные до целого числа, полученные на основании 8 испытаний.
Таблица 1. Результаты решения системы (1) классическими методами
2nfj.
Метод (h= 256 , te= 512) Время t RMS max err
РК-4 4302 174 569
ЯМЭ 1531 31 102
ЯМЭК 3600 30 97
НЯМЭ 12334 29 100
Правило трапеций 18346 29 95
БД 4376 97 403
Анализ полученных данных свидетельствует, что неявные методы показывают более точный результат, однако требуют значительно больше времени на вычисления (поскольку требуют решения системы уравнений на каждом шаге). Поэтому очевидна идея объединить лучшие черты использованных алгоритмов. В результате должен получиться метод, который на каждом участке использует более подходящий алгоритм. Трудность заключается в выборе алгоритма для каждого участка. Этот выбор может быть сделан на основе множества факторов, таких как жесткость системы, спектральный портрет интегрируемой функции и т.д. Рассмотрим простейший вариант. Из общих соображений понятно, что наибольшую проблему представляют собой жесткие участки (т. е. для модельной системы - просто участки с большим значением производной по времени). Поэтому на этих участках требуется более точный и стабильный метод интегрирования. На нежестких участках можно без сильной потери точности воспользоваться более быстрым методом. Таким образом, в самом общем виде алгоритм вычислений следующий:
1. for t е (ti | ti = to+ih, ti < te)
a. d = приближенное значение производной на нескольких предыдущих шагах
b. if (d > D) использовать Алгоритм1 иначе использовать Алгоритм2
здесь (to, te) - интервал, на котором производится интегрирование; h - шаг; D - некоторое число.
Исходя из данных табл. 1, а также некоторых теоретических соображений, наиболее выгодной представляется комбинация быстрого явного метода (метод Эйлера или метод Рунге-Кутты) и устойчивого неявного (правило трапеций). Рассмотрим оба варианта. В более развернутом виде алгоритм представлен ниже.
//n - номер итерации, D - некоторая константа for (real_t t = to;t< te;t+=h)
{
//если не накоплено достаточно истории для
аппроксимации производной, используем «точный»
алгоритм
if(n<k)){
nmethod = STABLE;
}
else{
k
ave d = -1 dn-i ; //k — длина промежутка ап-i=0
проксимации производной
if( ||ave d|| > D) //производная велика, использовать «точный» алгоритм nmethod = STABLE; else //производная невелика,
использовать «быстрый» алгоритм nmethod = FAST;
}
//вычислить следующее значение по выбранному алгоритму
if(nmethod==STABLE){
Yn=StableMethod();
}
else if(nmethod==FAST){ yn=FastMethod();
}
dn = ||yn - yn-i||;
}
Очевидно, что приведенный алгоритм может быть усовершенствован. Например, количество используемых методов может быть увеличено, усовершенствована процедура выбора конкретного алгоритма и т.д. Однако и в таком варианте результаты вычислений по двум алгоритмам, построенным по приведенной выше схеме, значительно лучше результатов вычислений по «классическим» алгоритмам, применяемым по отдельности.
Таблица 2. Результаты решения системы (1) двумя гибридными методами
Метод t RMS max err
(t = 1/(32*8)), tmax = 512
Эйлера-Трапеций 2389 17 59
РК-Трапеций 5686 7 33
Параллельная реализация. Параллельная реализация предложенных алгоритмов во многом опирается на качественную параллельную реализацию использованных «классических» методов. Это делает возможным использование множества существующих оптимизированных параллельных реализаций. Используемая схема выбора метода интегрирования не требует ресурсоемких операций и специального распараллеливания.
Спектральные характеристики. Характер спектральных искажений [2], вносимых предлагаемыми алгоритмами, трудно определить однозначно. Он зависит в первую очередь от используемых «классических» алгоритмов, характера интегрируемых функций, выбора констант О,к (непосредственно влияющих на выбор алгоритма на конкретном участке). Для модельной системы, решаемой с параметрами из табл. 1, график 1одН [2] имеет вид, показанный на рис. 2:
Рис. 2. График 1одН для метода РК-Трапеций
Анализ графика показывает, что в решении присутствует высокочастотный шум, который, однако, убывает с ростом частоты.
Заключение
Рассмотренный альтернативный подход к созданию алгоритмов численного решения систем обыкновенных дифференциальных уравнений достаточно перспективен и, несомненно, требует дальнейшей разработки. Проведенные численные эксперименты показывают, что даже в простейшей реализации можно получить более качественные результаты (по соотношению скорость/точность), чем при использовании «классических» алгоритмов.
Список литературы
1. Hairer Ernst; Norsett Syvert Paul & Wanner Gerhard. Solving ordinary differential equations I: Nonstiff problems. - Berlin, New York: Springer-Verlag, 2008.
2. Чадов С.Н. О спектральных характеристиках методов численного решения обыкновенных дифференциальных уравнений / Применение многопроцессорных суперкомпьютеров в исследованиях, наукоемких технологиях и учебной
работе: Сб. мат-лов регион. науч.-техн. конф. - Иваново: ИГТА, 2008.
3. Амосов А.А., Дубянский Ю.А., Копченов Н.В.
Вычислительные методы для инженеров: Учеб. пособие. -М.: Высш. шк., 1994.
4. William H. Press ... [et al]. Numerical recipes in C: the art of scientific computing . 2nd ed.
5. Lyons Richard G. Understanding digital signal processing. - Prentice Hall, 2001.
Чадов Сергей Николаевич,
ГОУВПО «Ивановский государственный энергетический университет имени В.И. Ленина», аспирант кафедры высокопроизводительных вычислительных систем, e-mail: sergei.chadov@gmail.com
