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

Горбань А.И.
Кафедра ПМиИ ДонНТУ
nastya_gorban@mail.ru

Abstract
Gorban’ A.I. The analysis of parallel numerical methods of solving the Cauchy problem for ordinary differential equations and stability of algorithms. In represented article the parallel methods of the decision of solving the Cauchy problem for ordinary differential equations are described. In their basis the so-called block method lays which will be coordinated to architecture of parallel computing systems. The Coefficients of blocks for difference amount of points are determine by the medium of packet Mathematica®. The approaches to research of stability of methods of solving the Cauchy problem for systems of the ordinary differential equations are described.

Введение

Концепция параллелизма привлекает внимание специалистов потенциальными возможностями повышения производительности и надежности вычислительных систем, поэтому в современной науке количество работ по теории параллельных вычислений постоянно растет. 
В течение последних десятилетий пиковая производительность суперЭВМ возросла на несколько порядков, изменилась технологическая база и архитектура. Однако отсутствие параллельных алгоритмов является основным препятствием к внедрению многих параллельных архитектур [1]. 
Становится очевидным, что в ближайшем будущем успехи в развитии многопроцессорных вычислительных систем будут зависеть в большей степени от достижений вычислительной математики, а не от технического прогресса. Поэтому наука уделяет большое внимание модернизации существующих алгоритмов с целью применения широкого параллелизма [2], а также разработке совершенно новых методов и алгоритмов. 
Что касается теории численных методов, и решения дифференциальных уравнений в частности, то многие алгоритмы, успешно используемые в последовательных вычислениях, оказываются не применимыми при проведении параллельных вычислений. Отсюда вытекает необходимость построения специальных методов, адаптированных к вычислительным системам с параллельной архитектурой.
При этом следует учитывать, что метод должен обладать численной устойчивостью и иметь оптимальную для сохранения эффективности структуру. Данным требованиям удовлетворяют блочные методы [3].
В данном методе для блока из k точек новые k значений функции вычисляются одновременно, что вполне согласуется с архитектурой параллельных вычислительных систем. Такой принцип позволяет вычислять коэффициенты разностных формул не в процессе интегрирования, а на этапе разработки метода, что значительно увеличивает эффективность счета. В то же время для блочных методов отсутствуют оценки погрешностей и условия сходимости.
При моделировании систем, описываемых системой дифференциальных уравнений, важнейшим вопросом является проблема устойчивости. Существует целый ряд методов, реализованных на практике, которые используют последовательные алгоритмы решения жестких систем. Однако разработка алгоритмов и методов для их использования в системах с параллельной архитектурой является открытым вопросом.
Данная работа посвящена исследованию сходимости и получению оценок погрешностей параллельных блочных методов решения задачи Коши для обыкновенных дифференциальных уравнений, а также вопросам устойчивости алгоритмов.
Рассмотрим решение задачи Коши
f(t,x),  x(t0) = x0                                                                            (1)
k точечным блочным разностным методом. Множество M точек равномерной сетки tm, m=1,2,…,M и tM=T с шагом разобьем на N блоков, содержащих по k точек, при этом kN ≥ M. В каждом блоке введем номер точки i =0,1,…,k  и обозначим через tn,i точку n –го блока с номером i. Точку tn,0  назовем началом блока n, а tn,k – концом блока. Очевидно, что имеет место tn,k = tn+1,0 . Условимся начальную точку в блок не включать. Если для расчета значений в новом блоке используется только последняя точка предшествующего – блочный метод считается одношаговым, а если все точки предшествующего блока – многошаговым.
В общем случае уравнения многошаговых разностных методов для блока из k точек с учетом введенных обозначений можно записать в виде
N                                       (2)
где .
Для одношаговых разностных методов уравнения могут быть представлены как 
 N                                                   (3)

Погрешность аппроксимации многошаговых блочных методов
Выражение для погрешности аппроксимации рассматриваемого разностного метода (2) на решении уравнения (1) запишем следующим образом
                                (4)
где , ,

Воспользуемся подходом, предложенным в работе [5]. Разлагая , и  в ряды Тейлора в окрестности точки , получим

,

Подставляя эти разложения в выражение (4) для погрешности аппроксимации и учитывая, что конечная точка предыдущего блока совпадает с начальной точкой следующего блока, будем иметь

                                    (5)
Сгруппируем члены с одинаковыми производными и изменим порядок суммирования в последнем выражении, тогда получим

.                                  (6)
Отсюда следует, что погрешность аппроксимации имеет порядок p, если выполнены условия 
            (7)
Система уравнений (7) для каждого фиксированного i содержит p уравнений и 2k неизвестных  . Потребуем, чтобы p=2k, тогда из системы (7) можно будет определить все неизвестные коэффициенты. Отсюда следует, что наивысший порядок аппроксимации многошагового k- точечного блочного метода равен 2k. Его погрешность в соответствии с (6) определяется формулой
              (8)
Для одношаговых блочных методов аппроксимацию порядка p обеспечивают следующие условия [4]
                                         (9)
Система уравнений (9) для каждого фиксированного i содержит p уравнений и k+1 неизвестных  . Если потребовать, чтобы p=k+1, тогда из системы (9) можно будет определить все неизвестные коэффициенты. Отсюда следует, что наивысший порядок аппроксимации одношагового k- точечного блочного метода равен k+1. Его погрешность определяется формулой 
                             (10)

Примеры нахождения погрешности разностных методов
Элементы bij и ai,j ,i,j=1,2,…,k матриц B и A соответственно можно найти в зависимости от выбранного метода, решая системы (7) или (9) для конкретной размерности блока k. Погрешности аппроксимации одношагового блочного метода могут быть получены по формулам (8) или (10).
Решение систем выполним с помощью пакета Mathematica . Получившиеся при этом коэффициенты блочных методов и погрешности в соответствующих точках блока будут иметь следующий вид:
-для двухточечного блока двухшагового метода

,	

- для трехточечного блока одношагового метода

- для четырехточечного блока одношагового метода

Сходимость и оценка погрешности многошаговых блочных методов
Обозначим через un - вектор значений приближенного решения в точках n- го блока с компонентами un,i, i=1,2,…k и через  вектор, компоненты которого и  равны значениям правых частей уравнения (1) в точках, соответствующих значениям приближенного решения. Запишем систему (2) в векторной форме
                                                      (11)
где D =(dii) – диагональная матрица с  элементами dii=i, i=1,2,…,k,
e- единичный вектор столбец, B и A –матрицы размерности  с компонентами и  cоответственно,  i=1,2,…,k, j=1,2,…,k.
Обозначим через xn - вектор значений точного решения задачи (1) в точках блока n. Получим уравнение, которому удовлетворяет вектор погрешностей в блоке zn=un-xn. Подставим в левую часть уравнения (11) выражение un= zn +xn, добавим к правой части и вычтем из нее , тогда уравнение для погрешности примет вид 
+.
Входящее в правую часть выражение
                                                      (12)
представляет собой вектор невязок разностных уравнений (2) на точном решении уравнения (1). Поскольку разностные уравнения (11) аппроксимируют исходное уравнение в точках блока с порядком , то имеет место оценка
                                                                                      (13)
Оставшиеся члены правой части уравнения для погрешности обозначим через 
                                                             (14)
Тогда уравнение для погрешности запишется короче
                                                           (15)
Вектор функция  от погрешности  зависит нелинейно. Вид этой зависимости определяется функцией f(t,x). В дальнейшем будем предполагать, что f(t,x) удовлетворяет условию Липшица по второму аргументу, т. е. 
                                                       (16)
для всех t, x1, x2 из рассматриваемой области. Согласно формуле конечных приращений Лагранжа для компоненты i имеем
 
где  
Подставляя последние в (14), получим 
,
где Dl – диагональная матрица с элементами d li,i =ln,i , i=1,2,…,k и ее норма имеет в силу (13) следующую оценку 
L                                                                                                  (17)
Заменим  в уравнении (15) полученным для него выражением, запишем его в виде

Введем нормы
и 
Далее, учитывая (13), (17) и то, что , получим неравенство
,
которое преобразуем к виду 
(1-kL)
Если на   наложить ограничение 
,                                                                                  (18)
то оценка, связывающая нормы погрешностей соседних блоков примет вид
                                                          (19)
т. к. в силу условия (18) имеет место   1- kL >0. Подставляя последовательно в (19) значения погрешностей для блоков n-1, n-2,…,1, получим

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

Последнее неравенство справедливо для любого n, , из этого следует: 
                                        (20)
Теорема. Пусть правая часть уравнения (1) f(t,y) удовлетворяет условию Липшица по второму аргументу с константой L и rn –невязка многошагового к- точечного блочного метода (2) определенная согласно (3) с оценкой (13). Тогда при
  и   для погрешности метода имеет место оценка (20).
Следствие. Если разностное уравнение (2) аппроксимирует исходное уравнение (1), то решение разностной задачи (2) сходится при  к решению исходной задачи (1), причем порядок точности совпадает с порядком аппроксимации [4].
После проведения соответствующих преобразований можно получить оценку погрешности одношаговых блочных методов. При тех же ограничениях на длину шага она будет иметь следующий вид [3]
                                                               (21)

Устойчивость разностных методов

Пусть для численного решения задачи Коши
                                                                                       (11)
используется линейный m-шаговый разностный метод
.                                                       (12)
Pассмотрим также однородное разностное уравнение
                                                                                                              (13)
и будем искать решения (24), имеющие вид , где  число, подлежащее определению. Подставив  в формулу (13) и разделив на n-m+1  получим  для нахождения   уравнение  всех решений однородного разностного уравнения (13).
                                                                             (14)
которое называется характеристическим уравнением разностного метода (12). Устойчивость метода (12) определяется расположением корней характеристического уравнения (14). Требуется, чтобы все корни (14) удовлетворяли условию
                                                                                                                         (15) 
причем корни, для которых  не должны быть кратными. Заметим, что условия (15) обеспечивают  всех решений однородного разностного уравнения (13).
Говорят, что разностный метод (12) устойчив в смысле Дальквиста, если все корни характеристического уравнения (14) лежат внутри или на границе единичного круга комплексной плоскости, причем на границе единичного круга нет кратных корней.
Существуют определенное ограничение на порядок аппроксимации устойчивого метода [5]. 
Теорема Пусть метод (12) устойчив  в смысле Далквиста и имеет порядок аппроксимации p. Тогда  при m четном и  при m  нечетном. Для явных m - шаговых устойчивых методов порядок аппроксимации не превосходит m.      
Теорема  Пусть метод (12) устойчив  в смысле Далквиста и, 
при . Тогда при ,  и всех достаточно малых  выполнена оценка
                                                          (16)
где -погрешность аппроксимации, - погрешность в задании начальных условий и С - константа, зависящая от L,T  и не зависящая от n.
Из оценки (16) следует, что если начальные погрешности  и погрешности аппроксимации  являются величинами  то и  при , т. е. метод сходится и имеет p-й порядок точности. Таким образом, исследование сходимости метода (12) сводится к анализу погрешности аппроксимации и проверке условия устойчивости  в смысле Далквиста.
Исследование устойчивости используемых разностных схем при решении задачи Коши для систем обыкновенных дифференциальных уравнений чаще всего проводится на модельном уравнении (17). 
								(17)
Приведем для примера определение условий устойчивости для явного метода Рунге-Кутты третьего порядка точности.
Найдем с помощью пакета Mathematica для модельного уравнения с правой частью  выражения для коэффициентов  и приближенного решения .
In[1]:=

Out[1]=

Запишем полученное выражение для  в виде  где . Построим график  функции .
In[2]:=

Out[2]=

Рисунок 1 – График функции  .

Из рисунка 1 видно, что условие устойчивости метода () выполняется при . Если же , метод Рунге-Кутта третьего порядка точности условно устойчив при 
Аналогичным образом определяются области устойчивости явных методов и более высоких порядков точности.

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

Литература

1.	Воеводин В.В.. Информационная структура алгоритмов. - М.:Изд-во МГУ,1997.
2.	Забродин А.В., Елизаров Г.С., Каратанов В.В. и др. Возможности и ближайшие перспективы создания высокопроизводительных вычислительных систем.//Тезисы докладов конференции "Высокопроизводительные вычисления и их приложения" 30 октября - 2 ноября 2000 года, Черноголовка. – М.:МГУ, 2000. С.8-9.
3.	Системы параллельной обработки: Пер. с англ./ Под. ред. Ивенса Д.- М.: Мир,1985-416с.
4.	Фельдман Л.П. Сходимость и оценка погрешности параллельных одношаговых блочных методов моделирования динамических систем с сосредоточенными параметрами // Научные труды ДонНТУ. Серия: Информатика, моделирование и вычислительные методы, (ИКВТ-2000), выпуск 15: Донецк: ДонНТУ, 2000, с. 34-39.
5.	Самарский А.А., Гулин А.В. Численные методы.- М.: Гл. ред. физ.-мат. лит., 1989.-432 с.