Параллельная реализация явного метода Эйлера с контролем точности вычислений Текст научной статьи по специальности «Математика»

УДК 519.622.852

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

Геннадий В. Ващенко*

Институт вычислительного моделирования СО РАН, Академгородок, 50/44, Красноярск, 660036,

Россия

Евгений А. Новиков^

Институт вычислительного моделирования СО РАН, Академгородок, 50/44, Красноярск, 660036,

Россия

Получена 18.05.2010, окончательный вариант 25.06.2010, принята к печати 10.09.2010 Для численного решения задачи Коши для системы обыкновенных дифференциальных уравнений рассмотрен параллельный алгоритм явного метода Эйлера с контролем точности вычислений. Приведены результаты численных экспериментов параллельной реализации алгоритма, выполненной на языке С и функций библиотеки MPI.

Ключевые слова: явный метод, параллельный алгоритм, локальная погрешность, контроль точности, кластер.

Введение

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

Во многих случаях расчеты требуется проводить с невысокой точностью — 1 % и ниже. Это связано с тем, что измерение констант, входящих в правую часть системы дифференциальных уравнений, часто проводится достаточно грубо. Иногда такая точность расчетов является удовлетворительной с точки зрения поставленной цели. Известно (см., например, [1, 2]), что порядок аппроксимации численной схемы следует сочетать с требуемой точностью расчетов. Поэтому ниже будет рассмотрен метод первого порядка точности.

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

* [email protected]

[email protected] © Siberian Federal University. All rights reserved

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

Хуже обстоит дело в алгоритмах интегрирования на основе известных безытерационных методов, к которым относятся методы типа Розенброка и их различные модификации [3]. Следует отметить, что безытерационные методы существенно более просты с точки зрения реализации на ЭВМ, чем алгоритмы на основе численных формул, в которых стадии вычисляются с применением итераций. Поэтому они более привлекательны для пользователей-вычислителей (химиков, биологов и др.). Однако в методах типа Розенброка и их модификациях матрица Якоби влияет на порядок точности численной схемы, и поэтому возникают трудности с ее замораживанием. Если вопрос об использовании одной матрицы на нескольких шагах интегрирования оставить нерешенным, то заведомо нужно ограничиться задачами небольшой размерности. Дальнейшего повышения эффективности алгоритмов интегрирования можно ожидать за счет распараллеливания вычислений.

В настоящее время многие пользователи применяют явный метод Эйлера с постоянным шагом для решения нежестких и умеренно жестких задач. Это обусловлено простотой реализации данной численной схемы, а также ее монотонностью. Отметим, что для многих задач расчеты с постоянным шагом приводят к низкой эффективности алгоритма интегрирования. В случае жестких задач весь промежуток интегрирования можно разбить на несколько отрезков. На некоторых из них, называемых переходными, решение меняется быстро (решение типа погранслойного), но такие промежутки небольшие по величине. На данных участках шаг интегрирования всегда достаточно мал, потому что производные от решения велики. Отметим, что переходный участок при определенном выборе начальных условий может отсутствовать. Однако для явных методов это не делает лучше ситуацию, потому что величина шага интегрирования все равно будет определяться максимальным собственным числом матрицы Якоби исходной системы дифференциальных уравнений.

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

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

1. Последовательный алгоритм

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

у' = f (у), y(to) = Уо, to < t < tk, (1)

где у и / — вещественные Ж-мерные вектор-функции, £ — независимая переменная. Для численного решения задачи (1) применим одностадийную схему, (п + 1)-й шаг которой задается формулой [3, 4]

Уп+1 = Уп + Нп/(у„), у(£о) = Уо, ¿о < £ < ¿к, (2)

где Нп — шаг интегрирования. Рассмотрение автономной задачи (1) не снижает общности. Введением дополнительной переменной у'м +1 = 1, уN +1 (¿о) = ¿о, всегда можно неавтономную задачу у' = /(у,£) привести к автономному виду.

Для построения последовательного алгоритма с контролем точности вычислений требуется уметь находить оценку локальной погрешности и определять размер шага интегрирования, гарантирующего выполнение заданной величины погрешности. Далее предполагается, что алгоритм изменения величины шага построен на основе контроля точности численной схемы, а в неравенстве для контроля точности применяется оценка локальной ошибки ¿п метода. Введем обозначения /п+1 = /(уп+1) построим неравенство для контроля точности вычислений.

Теорема 1. Пусть для численного 'решения задачи (1) используется метод (2), и пусть £ — требуемая точность расчетов. Тогда для контроля точности на шаге имеет место неравенство

0, 5Н||/п+1 - /п|| < £, (3)

где Н = Нп, /п+1 и /п — значения правой части системы (1) соответственно в точках ¿п+1 и ¿п, || • || — некоторая норма в .

Доказательство. Используя разложения точного и приближенного решений в ряды Тейлора по степеням Н до членов с Н2 включительно, получим, что для метода (2) выражение локальной погрешности ¿п имеет вид

¿п =0, 5Н2/'/ + 0(Н3),

где элементарный дифференциал /'/ вычислен на точном решении у(£п). Далее, учитывая разложение Н/(уп+1) в ряд Тейлора в окрестности приближенного решения уп до членов с Н2 включительно

получим

Н/(уп+1) = Н/п + Н2/п/п + 0(Н3),

н(/п+1 - /п) = Н2/п/п + о(Н3).

Сравнивая данное соотношение с локальной ошибкой, имеем неравенство (3). □

Отметим, что в случае успешного шага интегрирования применение значения /п+1 в неравенстве (3) к увеличению вычислительных затрат не приводит, потому что оно применяется на следующем шаге для вычисления решения по схеме Эйлера. Теперь сформулируем алгоритм интегрирования с контролем точности вычислений.

Алгоритм 1. Пусть для численного решения задачи (1) используется метод (2) и пусть известно приближенное решение уп в точке ¿п с шагом Нп. Тогда для получения приближенного решения уп+1 в точке ¿п+1 справедлив следующий алгоритм.

Шаг 1. Вычислить уп+1 = уп + Нп/(уп). Шаг 2. Вычислить /(уп+1).

Шаг 3. Вычислить норму ||£п|| = 0, 5Н||/п+1 — /п||.

Шаг 4. Вычислить значение д по формуле д = (е/||^п||)0'5.

Шаг 5. Если д < 1, то неравенство (3) не выполняется. Перевычислить Нп по формуле Нп := дЛ.„/1,1 и перейти на шаг 1.

Шаг 6. Если д ^ 1, то вычисленное решение уп+1 в точке ¿„+1 удовлетворяет требуемой точности и значение ¿„+1 принимается в качестве следующего узла.

Шаг 7. Задать следующий шаг интегрирования ^„+1 по формуле ^„+1 = дЛ.„/1,1 для дальнейшего интегрирования. Шаг 8. Перейти на шаг 1.

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

2. Параллельный алгоритм

При разработке параллельного алгоритма будем основываться на представленном выше последовательном варианте. Для определенности предположим, что вычислительная система состоит из p процессоров, N > p. Запишем параллельный аналог метода (2), имеем

yi+1) = j + hn+ifs (y(n)),

j = V. (to), (4)

1 < j < P, (j - 1) • « +1 < js < j • S,

где s = N/p, если N кратно p, или s = [N/p] + g в противном случае. Теперь сформулируем параллельный алгоритм интегрирования с контролем точности вычислений. Как и в случае последовательного алгоритма, в основу положим проверку неравенства (3), то есть контроль точности будет осуществляться с применением оценки локальной ошибки метода (4). Алгоритм 2. Пусть для численного решения задачи (1) используется метод (4), а размерность N системы (1) превосходит размерность вычислительной системы p, N > p,

j € Comp(j), 1 < j < p, (j - 1) • s + 1 < js < j • s.

Пусть известно решение y(n) в точке tn с шагом hn. Тогда для получения значения y(n+1) в точке tn+i справедлив следующий параллельный алгоритм.

Шаг 1. Переслать y(0) из Comp(1) всем Comp(j), 2 ^ j ^ p.

n) 1 и , •

Шаг 2. В каждом Сотр(з), 1 ^ з ^ р, вычислить у(„+1') = у^ + (у(п)). Шаг 3. Переслать у(„+1) всем от всех. Шаг 4. В каждом Сотр(з), 1 ^ з ^ р, вычислить /^ (у(п Шаг 5. В каждом Сотр(з), 1 ^ з ^ р, вычислить 0, 5Л.„||/^ (у(п+1))-/¿а (у(и))||^ .Вычисление нормы по всем компонентам - глобальная операция. Шаг 7. В Сотр(1):

a) вычислить значение д по формуле д = (е/||^п||)05,

b) если д < 1, то есть неравенство (3) не выполняется, перевычислить шаг по формуле := д^„/1,1,

c) переслать из Сотр(1) всем компьютерам Сотр(з), 2 ^ з ^ р, и перейти на шаг 2, ^ если д ^ 1, то вычисленное решение у(п+1) в точке ¿п+1 удовлетворяет требуемой

точности. Тогда значение ¿„+1 принимается в качестве следующего узла, следующий шаг ^„+1 задается по формуле ^„+1 = д^„/1,1 и пересылается всем Сотр(з), 2 ^ з ^ р. Шаг 8. В каждом Сотр(з), 1 ^ з ^ р, положить = у(п+1) и /^ (у(п)) = /¿а (у(п+1)).

Шаг 9. Переслать всем от всех.

Шаг 10. Выполнить следующий шаг интегрирования.

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

3. Результаты расчетов

При проведении вычислительных экспериментов ограничимся иллюстрацией принципиальной возможности параллельной реализации явной схемы первого порядка с контролем точности вычислений и построением экспериментальной зависимости времени выполнения для достижения заданной точности е от числа процессоров p и правой части исходной системы дифференциальных уравнений, то есть t = t(p, f). Это дает возможность оценить значения ускорения Sp параллельного алгоритма [5]. Кроме того, дадим справку по суммарному количеству шагов интегрирования isa, числу повторных вычислений решения (возвратов) при нарушении требуемой точности расчетов iwo и числу вычислений правой части исходной системы ifu.

Приведем описание используемой тестовой системы и исходных данных. Пусть процесс многостадийного синтеза вещества без ветвления [6] описывается следующей системой:

где

£ = F ">■

F (x)

(5)

( g(xN) - xi \ c(x 1 - X2)

\ CXN _1 — ©xN

FX =

-c

c

V 0

0

-c 0

0 q'(Xn) \ 00

c —©/

где FX — матрица Якоби системы (5), N — размерность системы, с = (N — 1)/т, т и в — параметры. Численное решение задачи (5) рассмотрим на отрезке t G [0, 9; 1] с точностью е = 0,1. Размерность системы N составляет 1 000 000 уравнений, значения параметров имеют вид т =1, в = с. Тогда из вида матрицы Якоби следует, что задача (5) является жесткой. Неравенство для контроля точности имеет вид (3), где норма ||£n|| оценки локальной ошибки Sn вычисляется по формуле ||£n|| = max |^П|/(|уП| + г). Если по i-й компоненте

l^i^N

решения выполняется неравенство | yln | < г, то контролируется абсолютная погрешность ег, иначе относительная ошибка е. Значение параметра г принималось равным 1, при котором выполняется требуемая точность расчетов.

Вычисления выполнялись на 99-процессорном кластере ИВМ СО РАН [7] с использованием языка C и функций библиотеки MPI. На рис. 1 представлены графики зависимости времени выполнения от количества используемых процессоров: по оси ординат — время в сек., по оси абсцисс — число процессоров. Здесь F 1(x) соответствует случаю, когда в

правой части исходной системы функция ) имеет вид ) = 2/(1 + 3xN), F2(x) для случая g(xN) = 10/(1 + 300xN) и F3(x) при g(xN) = 100/(1 + 30 000xN). Начальные условия: yi(0, 9) = 100, y,(0, 9) = 0,1, i = 3, 5,...,N - 1 и y¿(0, 9) = 0, 2, i = 2,4,..., N.

Вычислительные затраты, полученные в процессе численного решения исходной системы, имеют вид:

для F1 (x) - isa = 141 450, iwo = 26, ifu = 141 576; для F2(x) - isa = 141 460, iwo = 32, ifu = 141 492; для F3(x) - isa = 141 468, iwo = 34, ifu = 141 502.

Особенность правой части исходной системы позволила организовать равномерное распределение данных по процессорам, а также свести к минимуму количество пересылаемых данных и обменов при вычислении вектора правой части F (x). Пересылка N-й компоненты функции g(x_N) осуществлялась в каждом tn узле с помощью функций MPI_Send(... ) и MPI_Recv(... ).

5000 4000 3000 2000 1000 0

Рис. 1. Изменение времени выполнения при увеличении числа процессоров

Заключение

На примере явной схемы метода Рунге-Кутты первого порядка рассмотрены основные особенности разработки параллельных алгоритмов с контролем точности для решения жестких задач на вычислительных кластерах с использованием библиотеки MPI.

Теоретические и экспериментальные оценки показали, что алгоритм с контролем точности для явной схемы первого порядка обладает достаточным потенциалом для распараллеливания и может служить основой для построения алгоритмов с переменным шагом и более высоким порядком.

Численные эксперименты подтвердили, что количество пересылок и объем пересылаемых данных связан со структурой правой части системы и требуемой точностью вычислений.

Работа выполнена при финансовой поддержке Российского фонда фундаментальных исследований (грант 08-01-00621) и Президента (грант НШ-3431.2008.9).

F1(x) —♦— F2(x) F3(x)

Список литературы

[1] Е.А.Новиков, Алгоритм интегрирования переменной структуры для решения жестких задач на основе явного и L-устойчивого методов, Вестник СибГАУ, 1(2008), №18, 75-78.

[2] Е.А.Новиков, Явные методы для жестких систем, Новосибирск, Наука, 1997.

[3] Э.Хайрер, Г.Ваннер, Решение обыкновенных дифференциальных уравнений. Жесткие и дифференциально-алгебраические задачи, М., Мир, 1999.

[4] Г.В.Ващенко, Е.А.Новиков, Параллельная реализация явных методов типа Рунге-Кутты, Вестник КрасГАУ, (2010), №2, 14-18.

[5] В.В.Воеводин, Вл.В.Воеводин, Параллельные вычисления, СПб., БХВ, 2002.

[6] Г.В.Демиденко, В.А.Лихошвай, Т.В.Котова, Ю.Е.Хропова, Об одном классе систем дифференциальных уравнений и об уравнениях с запаздывающим аргументом, Сиб. журн. индустр. математики, 17(2006), №1, 58-68.

[7] С.В.Исаев, А.В.Малышев, В.В.Шайдуров, Развитие Красноярского центра параллельных вычислений, Вычислительные технологии, 11(2006), 28-33.

Parallel Algorithm Explicit Euler Method with Accuracy Control

Gennadiy V. Vashchenko, Evgeniy A. Novikov

Parallel algorithm explicit Euler method with accuracy control for numerical solving of Cauchy problem

for differential equation systems is considered. Parallel implementation of the algorithm is written in C

and MPI- functions. The numerical results obtained are presented.

Keywords: explicit method, parallel algorithm, local error, accuracy control, clusters.