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

Широкую известность получили безытерационные численные схемы за счет хороших свойств точности и устойчивости, а также за счет простоты реализации. К безытерационным методам относятся схемы типа Розенброка и различные их модификации. Однако в данных методах матрица Якоби включена непосредственно в вычислительные формулы, что приводит к определенным сложностям с замораживанием матрицы Якоби, то есть с применением одной матрицы на нескольких шагах интегрирования. В [2] доказана теорема о том, что максимальный порядок точности методов типа Розенброка равен двум, если они реализуются с применением одной матрицы на нескольких шагах интегрирования. Отметим, что данные численные формулы достаточно просты с точки зрения реализации, и как следствие, привлекательны для многих пользователей. В [3-5] предложен новый класс одношаговых безытерационных численных схем. Они обладают хорошими свойствами точности и устойчивости, а также достаточно просто реализуются с замораживанием матрицы Якоби. В дальнейшем эти численные схемы стали называть (m, k)-методами, и они определяются следующим образом.

Рассмотрим задачу Коши для системы обыкновенных дифференциальных уравнений:

, , ,    (1)

где , . Пусть заданы целые числа m и k, . Определим множества:

,

,

  ,

Тогда (m,k)-методы имеют вид [3-5]:

, ,

, ,     (2)

, ,

где , , - векторы приращений, E - единичная матрица,  - матрица Якоби системы (1),  - шаг интегрирования, α, p_i, α_ij и β_ij - вещественные константы. Затраты на шаг в методах (2) следующие - один раз вычисляется матрица Якоби и осуществляется декомпозиция матрицы D_n,  k  раз вычисляется функция f и m раз реализуется обратный ход метода Гаусса.

Здесь представлены вычислительная схема (2,2)-метода, ее последовательный и параллельный алгоритмы интегрирования, в которых изменение величины шага основывается на оценке аналога глобальной ошибки. Формирование параллельного алгоритма (2,2)-метода состоит в использовании декомпозиции на подзадачи и установлении взаимосвязи между ними [6-7].

Постановка задачи. Для решения задачи (1) будем применять (2,2)-схему вида

,

,   (3)

.  

Здесь  и , , - векторы приращений, , матрица A_n представима в виде

,

а матрица B_n не зависит от размера шага интегрирования. Такое представление A_n позволяет применять схему (3) с замораживанием как численной, так и аналитической матрицы Якоби. В случае использования схемы (3) с замораживанием матрицы Якоби, то есть в случае применения матрицы, вычисленной k шагов назад, получим:

,

где имеет место

, ,

 - максимальное число шагов с замороженной матрицей Якоби,

.

Если матрица Якоби вычисляется численно с шагом , где c_j, , - постоянные числа, то запишем:

, ,

где элементы  матрицы  определяются по формулам:

, .

В случае замораживания численной матрицы Якоби можно записать:

,

.

Из приведенных выше формул следует, что вопрос о замораживании и численной аппроксимации матрицы Якоби можно рассматривать одновременно.

Требование второго порядка точности и учет дополнительной ошибки за счет аппроксимации матрицы Якоби приводит к коэффициентам:

,   ,   ,   ,   .

Изменение величины шага основано на оценке аналога глобальной ошибки  [5]. Учитывая очевидное соотношение

,

новый шаг  будем определять по формуле , где значение q находится из уравнения [4]:

.

Если q<1, то есть требуемая точность не выполняется, то происходит повторное вычисление решения с шагом . Если имеет место неравенство , то вычисляется решение в следующей точке с шагом интегрирования .

Ниже при записи алгоритмов будем использовать обозначения  LU_Decompos() и LU_Solution() для функций, реализующих LU-факторизацию и решения систем алгебраических уравнений относительно приращений  и . Аналоги этих функций для параллельного алгоритма обозначим  Par_LU_Decompos() и Par_LU_Solution().

Последовательный алгоритм. Последовательный алгоритм интегрирования с контролем точности вычислений формулируется следующим образом. Пусть для численного решения задачи (1) используется (2,2)-схема (3) и известно приближенное решение y⁽ⁿ⁾ в точке t_n с шагом h_n. Тогда для получения приближенного решения y ⁽ⁿ⁺¹⁾ в точке t _(n+1) справедлив следующий алгоритм.

Шаг 1. Вычислить  .

Шаг 2. Вычислить  матрицу Якоби .

Шаг 3. Сформировать матрицу .

Шаг 4. Разложить матрицу D_n, то есть .

Шаг 5. Вычислить , то есть .

Шаг 6. Вычислить .

Шаг 7. Вычислить ,  то есть .

Шаг 8. Вычислить величину .

Шаг 9. Вычислить значение q по формуле  и .

Если q<1, то возврат на шаг 3  с шагом интегрирования .

Шаг 10. Вычислить приближенное решение по формуле (3):

.

Шаг 11. Вычислить решение в следующей точке с шагом интегрирования .

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

Параллельный алгоритм. При записи параллельного алгоритма предполагаем, что имеется вычислительная система из p процессоров , . Матрица D_n размещена строчными блоками размером s, причем s =N/p, если N кратно p, и  в противном случае. Для контроля точности численной схемы (3) введем функцию , а для ее выполнения назначим процессор . Параллельный алгоритм вычисления приближенного решения  с переменным шагом интегрирования формулируем следующим образом.

Пусть известно решение y⁽ⁿ⁾ в точке t_n с шагом h_n. Тогда для получения значения y ⁽ⁿ⁺¹⁾ в точке t _n+1 справедлив параллельный алгоритм, в котором на каждом процессоре proc(j) формируется своя j-я часть  вектора решения.

Шаг 1. В каждом proc(j), 1 ≤ j ≤ p , ,

выполнить ,

вычислить  и матрицу Якоби J_j, 1 ≤ j ≤ p.

Шаг 2. Сформировать матрицу .

Шаг 3. Разложить матрицу D_n, .

Шаг 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