Широкую известность получили безытерационные численные схемы за счет хороших свойств точности и устойчивости, а также за счет простоты реализации. К безытерационным методам относятся схемы типа Розенброка и различные их модификации. Однако в данных методах матрица Якоби включена непосредственно в вычислительные формулы, что приводит к определенным сложностям с замораживанием матрицы Якоби, то есть с применением одной матрицы на нескольких шагах интегрирования. В [2] доказана теорема о том, что максимальный порядок точности методов типа Розенброка равен двум, если они реализуются с применением одной матрицы на нескольких шагах интегрирования. Отметим, что данные численные формулы достаточно просты с точки зрения реализации, и как следствие, привлекательны для многих пользователей. В [3-5] предложен новый класс одношаговых безытерационных численных схем. Они обладают хорошими свойствами точности и устойчивости, а также достаточно просто реализуются с замораживанием матрицы Якоби. В дальнейшем эти численные схемы стали называть (m, k)-методами, и они определяются следующим образом.
Рассмотрим задачу Коши для системы обыкновенных дифференциальных уравнений:
,
,
, (1)
где 
, 
. Пусть заданы целые числа m и k, 
. Определим множества:
,
,
,
Тогда (m,k)-методы имеют вид [3-5]:
, 
,
, 
, (2)
, 
,
где 
, 
, - векторы приращений, E - единичная матрица, 
- матрица Якоби системы (1), 
- шаг интегрирования, α, pi, αij и βij - вещественные константы. Затраты на шаг в методах (2) следующие - один раз вычисляется матрица Якоби и осуществляется декомпозиция матрицы Dn, k раз вычисляется функция f и m раз реализуется обратный ход метода Гаусса.
Здесь представлены вычислительная схема (2,2)-метода, ее последовательный и параллельный алгоритмы интегрирования, в которых изменение величины шага основывается на оценке аналога глобальной ошибки. Формирование параллельного алгоритма (2,2)-метода состоит в использовании декомпозиции на подзадачи и установлении взаимосвязи между ними [6-7].
Постановка задачи. Для решения задачи (1) будем применять (2,2)-схему вида
,
, (3)
.
Здесь 
и 
, 
, - векторы приращений, 
, матрица An представима в виде
,
а матрица Bn не зависит от размера шага интегрирования. Такое представление An позволяет применять схему (3) с замораживанием как численной, так и аналитической матрицы Якоби. В случае использования схемы (3) с замораживанием матрицы Якоби, то есть в случае применения матрицы, вычисленной k шагов назад, получим:
,
где имеет место
, 
,
- максимальное число шагов с замороженной матрицей Якоби,
.
Если матрица Якоби вычисляется численно с шагом 
, где cj, 
, - постоянные числа, то запишем:
, 
,
где элементы 
матрицы 
определяются по формулам:
, 
.
В случае замораживания численной матрицы Якоби можно записать:
,
.
Из приведенных выше формул следует, что вопрос о замораживании и численной аппроксимации матрицы Якоби можно рассматривать одновременно.
Требование второго порядка точности и учет дополнительной ошибки за счет аппроксимации матрицы Якоби приводит к коэффициентам:
, 
, 
, 
, 
.
Изменение величины шага основано на оценке аналога глобальной ошибки 
[5]. Учитывая очевидное соотношение
,
новый шаг 
будем определять по формуле 
, где значение q находится из уравнения [4]:
.
Если q<1, то есть требуемая точность не выполняется, то происходит повторное вычисление решения с шагом 
. Если имеет место неравенство 
, то вычисляется решение в следующей точке с шагом интегрирования 
.
Ниже при записи алгоритмов будем использовать обозначения LU_Decompos() и LU_Solution() для функций, реализующих LU-факторизацию и решения систем алгебраических уравнений относительно приращений 
и 
. Аналоги этих функций для параллельного алгоритма обозначим Par_LU_Decompos() и Par_LU_Solution().
Последовательный алгоритм. Последовательный алгоритм интегрирования с контролем точности вычислений формулируется следующим образом. Пусть для численного решения задачи (1) используется (2,2)-схема (3) и известно приближенное решение y(n) в точке tn с шагом hn. Тогда для получения приближенного решения y (n+1) в точке t (n+1) справедлив следующий алгоритм.
Шаг 1. Вычислить 
.
Шаг 2. Вычислить матрицу Якоби 
.
Шаг 3. Сформировать матрицу 
.
Шаг 4. Разложить матрицу Dn, то есть 
.
Шаг 5. Вычислить 
, то есть 
.
Шаг 6. Вычислить 
.
Шаг 7. Вычислить 
, то есть 
.
Шаг 8. Вычислить величину 
.
Шаг 9. Вычислить значение q по формуле 
и 
.
Если q<1, то возврат на шаг 3 с шагом интегрирования 
.
Шаг 10. Вычислить приближенное решение по формуле (3):
.
Шаг 11. Вычислить решение в следующей точке с шагом интегрирования 
.
Практическая реализация вычислений по данному алгоритму зависит от системы программирования, способа организации вычислений правой части системы (1) и способа факторизации (способа выбора ведущих элементов), а также от подходов к формированию числовых массивов и доступа к их элементам.
Параллельный алгоритм. При записи параллельного алгоритма предполагаем, что имеется вычислительная система из p процессоров 
, 
. Матрица Dn размещена строчными блоками размером s, причем s =N/p, если N кратно p, и 
в противном случае. Для контроля точности численной схемы (3) введем функцию 
, а для ее выполнения назначим процессор 
. Параллельный алгоритм вычисления приближенного решения 
с переменным шагом интегрирования формулируем следующим образом.
Пусть известно решение y(n) в точке tn с шагом hn. Тогда для получения значения y (n+1) в точке t n+1 справедлив параллельный алгоритм, в котором на каждом процессоре proc(j) формируется своя j-я часть 
вектора решения.
Шаг 1. В каждом proc(j), 1 ≤ j ≤ p ,
,
выполнить 
,
вычислить 
и матрицу Якоби Jj, 1 ≤ j ≤ p.
Шаг 2. Сформировать матрицу 
.
Шаг 3. Разложить матрицу Dn, 
.
Шаг 4. Вычислить 
, 
.
Шаг 5. В каждом proc(j), 1 ≤ j ≤ p, ,
выполнить 
и вычислить 
.
Шаг 6. Вычислить 
, то есть 
.
Шаг 7. В каждом proc(j), 1 ≤ j ≤ p ,
вычислить локальные нормы погрешности:
а затем вычислить глобальную норму:
, 1 ≤ j ≤ p .
Шаг 8. В proc(1) выполнить accur_control(): 
, 
.
Если q < 1 , то возврат на шаг 2 с шагом интегрирования 
.
Шаг 9. В каждом proc(j), 1 ≤ j ≤ p ,
вычислить:
и выполнить 
.
Шаг 10. Вычислить решение в следующей точке с шагом интегрирования 
.
Аналогичный алгоритм может быть сформулирован и для столбцовой схемы размещения матрицы . Теоретические и экспериментальные оценки показали, что подходящими топологиями для выполнения алгоритма являются полный граф и гиперкуб. Реализация алгоритма осуществлялась на 99-процессорном кластере ИВМ СО РАН [8] c использованием языка и функций библиотеки MPI.
Заключение. На примере (2,2)-схемы рассмотрены основные особенности разработки параллельных алгоритмов (m,k)-методов для многопроцессорных вычислительных систем кластерной архитектуры с использованием библиотеки MPI. Численные эксперименты на модельных задачах подтвердили, что количество пересылок и объем пересылаемых данных связан со структурой правой части системы. Подходящими топологиями для реализации являются полный граф и гиперкуб.
Работа выполнена при финансовой поддержке Российского фонда фундаментальных исследований (проект 11-01-00106).
Рецензенты:
- Богульский И.О., д.ф.-м.н., профессор кафедры сопротивления материалов и теоретической механики Института управления инженерными системами Красноярского государственного аграрного университета, г. Красногорск.
- Нужин Я.Н., д.ф.-м.н., профессор кафедры математического обеспечения дискретных устройств и систем (МОДУС) Института фундаментальной подготовки (ИФП) Сибирского федерального университета, г. Красноярск.
- Кирьянов Б.Ф., д.т.н., профессор кафедры прикладной математики и информатики Новгородского государственного университета имени Ярослава Мудрого Министерства образования и науки РФ.
Работа получена 12.11.2011