CC BY
33
УДК 519.688
ПАРАЛЛЕЛЬНОЕ ВЫЧИСЛЕНИЕ ОБЩЕГО РЕШЕНИЯ НЕОДНОРОДНОЙ СИСТЕМЫ ОБЫКНОВЕННЫХ ЛИНЕЙНЫХ ДИФФЕРЕНЦИАЛЬНЫХ УРАВНЕНИЙ С ПОСТОЯННЫМИ КОЭФФИЦИЕНТАМИ
© М.А. Рыбаков
Ключевые слова: неоднородная система обыкновенных линейных дифференциальных уравнений с постоянными коэффициентами; преобразование Лапласа; параллельный алгоритм.
Обсуждается параллельный алгоритм получения символьного аналитического решения неоднородной системы обыкновенных линейных дифференциальных уравнений с постоянными коэффициентами. Решение системы находится в аналитическом виде и может быть найдено с требуемой точностью. Полученный алгоритм эффективен для решения систем дифференциальных уравнений большого размера. Алгоритм входит в состав библиотеки алгоритмов системы МаШраг.
1 Постановка задачи
Пусть задана неоднородная система обыкновенных линейных дифференциальных уравнений с постоянными коэффициентами:
п т ¿к
^Б^х^) = ¡¿г),В.ц(г) = ^ а^ —к ,г = 1,.. .,т,акц е К,п,т Е Н, (1)
]=\ к=О
где ак^{ — действительные числа, ^(¡), х() — ограниченные на К+ функции, имеющие конечное число точек разрыва I рода и удовлетворяющие условиям: ^(¡) = 0 при г < 0, Цг{г)1 < Ме/0* при г > 0, где М > 0, з0 ^ 0 — некоторые действительные постоянные. Обозначим А(г) = (В^(¡)) систему (1), и запишем в матричном виде:
А(г)Х(г) = В (г),Х(г) = [х1(г),...,хп(г)]т, В(г) = [fl(t),..., fn(t)]T.
Пусть х(к) (¡) обозначает к-ую производную функции х(), а числа х^ определяют начальные условия:
х?’(0) = хЬ, (2)
где к = 1, 2,... ,т — 1, г = 1, 2,... ,т, х^ Е К.
Будем полагать, что в общем случае каждая из функций ^(¡) в правой части может иметь вид конечной суммы:
^(¡) = ^2Рз(¡)е6гз* (7ч¡) (вцt)UnгtStep(t — а^),
з
где аз ,вз ,^гз ,&з Е К,^ц ,Щ е Н,рз (^ — полином переменной ^ ^ ^такцпя UnгtStep(t) принимает значение 1 для неотрицательных аргументов и значение 0 для остальных. Требуется найти решение X(^ системы (1), удовлетворяющее условиям (2), в аналитическом виде.
2 Алгоритм Лапласа
Алгоритм состоит из трёх этапов [1-3].
Этап I. Прямое преобразование Лапласа системы дифференциальных уравнений.
Преобразования Лапласа для функции f (t):
F(p)= f (t)e-ptdt, p e C. (3)
J 0
В результате преобразования функций, стоящих в левой и правой частях системы дифференциальных уравнений (1), и начальных условий (2) по формуле (3) получаем систему алгебраических уравнений:
A(p)X(p) -B(p) = F(p). (4)
Здесь A(p) — образ левой части системы (1), F(p) — образ правой части системы (1), B(p) — вектор, появляющийся при введении начальных условий (2).
Этап II. Решение полученной алгебраической системы.
Решение алгебраической системы (4) ищем в виде
X (p) = A(p)-1(F (p) + B(p)). (5)
p
Этап III. Обратное преобразование Лапласа.
F( p)
1 Г ^
f (t) = L-1{F} = — eptF(p)dp.
™ J-**,
Искомое решение системы (1) задается вектором X(t):
1 f ^
X (t) = ^ <?X (p)dp. (6)
2m J _iX}
3 Параллельный алгоритм
Пусть имеется кластер с п процессорами с ном ерами 0,1, 2,..., (п — 1),
Будем считать, что процессор с номером 0 — корневой.
Разработанный параллельный алгоритм состоит из следующих шагов.
Шаг 1.
Прямое преобразование Лапласа исходной системы, дифференциальных уравнений.
Пусть Т(р) — образ функции f (^ при прямом преобразовании Лапласа, В результате прямого преобразования Лапласа производной п-ого порядка функции f (^ получим полином степени п.
В результате прямого преобразования Лапласа левая часть системы дифференциальных уравнений преобразуется в матрицу полиномов А(р) = (Рз(р)) одной действительной переменной р.
Результатом прямого преобразования Лапласа (3) правой части будет вектор функций,
Т (р) = [11(р),...,1и(р)]Т.
Каждая из функций (р) является суммой произведений экспоненциальных и дробнорациональных функций.
Результатом преобразования начальных условий будет В(р) — сумма произведений р
В результате преобразования Лапласа системы обыкновенных линейных дифференциальных уравнений получаем алгебраическую систему линейных уравнений
А(р) •Х(р) = Т(р) + В(р).
Все вычисления на этом шаге могут быть выполнены на нулевом процессоре.
Шаг 2.
Решение алгебраической системы линейных уравнений в кольце полиномов — алгоритмически сложная задача, для решения которой будут использованы все процессоры кластера.
Для матрицы полиномов А(р) Е С[р]пхт вычисляем определитель д,еЛ(А(р>)) и обратную матрицу А-1 (р) = А*(р)/д,еЛ(А(р)), где А*(р) — присоединённая матрица. Воспользуемся параллельными алгоритмами [12],
Матрицу А-1(р) получим па пулевом процессоре.
Шаг 3.
Полипом д,еЛ(А(р>)) раскладываем па свободные от квадратов множители.
Задача факторизации полинома распараллеливается и используются все узлы кластера.
Область поиска комплексных корней разбивается на сектора, количество секторов равно количеству процессоров. Нулевой процессор рассылает ¿^(А(р)) всем процессорам. Каждый процессор берет случайную точку в своем секторе и по методу Ньютона начинает спуск к корню уравнения. Таким образом, когда нулевой процессор получает к различных комплексных корней, все процессоры получают сообщение о прекращении вычислений [13],
В результате разложения полинома д,еЛ(А(р)) на линейные сомножители в области С получим:
Г
<!&(А(р)) = (р — рк)£к, где г ^ пт, рк Е С, £к Е N. (7)
к=0
Шаг 4.
Раскладываем дробь 1/det(A(p)) в сумму простых дробей в области C
¿вг{А{р)) ^=0 - РкУк’
На нулевом процессоре собираем матрицу полученную в результате решения данной задачи с помощью метода неопределенных коэффициентов. Матрицу Ы^ формируем следующим образом: количество столбцов равно к + 1, а количество строк равно старшей степени д,еЛ(А(р)), в первых к столбцах будут записаны коэффициенты полинома ГС=0(Р - Рк Ук для к = ] , в последнем столбце записаны коэффициенты д,еЛ(А(р)).
Решаем параллельно задачу и полученные на каждом процессоре неопределенные коэффициенты (к собираем па пулевом процессоре.
Шаг 5.
Обозначим через Н(р) = Т(р) + В(р) и Н(р) = [к\,..., кп]т . Тогда, получим:
n Q ( ) p n m-1 m-1
hi = Z ljцe— + Z Z dikip) = Z ai+np
j=l i=l k=0 i=k
Найдем решение
X(р) = А*(р) • Н • (1/авг(А(р)). (9)
Пусть х (р) = [х\,..., Хп]Т .Тогда хг(р) = Ек=о п-р^к , гДе Пк е С.
Для каждого элемента вектора X(р) = [х\(р),... ,Хп(р)]Т находим обратное преобразование Лапласа (6) X(I) = [х\(Ь),... ,хп(1)]т :
хг(Ь) = Ь~1{хг(р)}, где г = 1, 2,...,п. (10)
Умножение матрицы на вектор (9) и обратное преобразование Лапласа (10) выполним
к
матрицы А* (р), вектор Н (р) и сумму дробей (8), где к = , п — число строк матрицы
А*(р), а N — количество процессоров. Для полученных па каждом из процессоров к элементов вектора Хг (р) = [х^ (р),... ,Х<к (р)]Т находим обратное преобразование Лапласа (6), Затем нулевой процессор объединяет полученные результаты в искомое решение системы дифференциальных уравнений.
ЛИТЕРАТУРА
1. Микусинский Я. Операторное исчисление. М.: ИЛ, 1956.
2. Диткин В.А., Прудников А.П. Интегральные преобразования и операционное исчисление. М.: Физ-матгиз, 1974.
3. Дёч Г. Руководство к практическому применению преобразованию Лапласа. М.: Наука, главная редакция физико-математической литературы, 1965.
4. Malaschonok N.A. An Algorithm for Symbolic Solving of Differential Equations and Estimation of Accuracy // Computer Algebra in Scientific Computing. LNCS 5743. Berlin: Springer, 2009. P. 213-225.
5. Малашонок Г.И., Бетин А.А., Рыбаков M.A., Смирнов P.А. Параллельная компьютерная алгебра. Часть 3. Учебное пособие. Тамбов: Издательский дом ТГУ им. Г.Р. Державина, 2012.
6. Malaschonok G.I. Project of Parallel Computer Algebra // Tambov University Reports. Series: Natural and Technical Sciences. Tambov, 2010. V. 15. Issue 6. P. 1724-1729.
7. Малашопок Г. И. Компьютерная математика для вычислительной сети // Вестник Тамбовского университета. Серия Естественные и технические науки. Тамбов, 2010. Т. 15. Вып. 1. С. 322-327.
8. Малашонок Г. И. О проекте параллельной компьютерной алгебры // Вестник Тамбовского университета. Серия Естественные и технические науки. Тамбов, 2009. Т. 14. Вып. 4. С. 744-748.
9. Малашонок Г. И. О вычислении ядра оператора, действующего в модуле // Вестник Тамбовского университета. Серия Естественные и технические науки. Тамбов, 2008. Т. 13. Вып. 1. С. 129-131.
10. Рыбаков М.А. Решение систем линейных дифференциальных уравнений с кусочно-непрерывными правыми частями с помощью преобразования Лапласа // Вестник Тамбовского университета. Серия Естественные и технические науки. Тамбов, 2010. Т. 15. Вып. 4. С. 339-341.
11. Рыбаков М. А. Решение систем линейных дифференциальных уравнений с постоянными коэффициентами с помощью преобразования Лапласа // Вестник Тамбовского университета. Серия Естественные и технические науки. Тамбов, 2009. Т. 14. Вып. 4. С. 791-792.
12. Бетин А.А. Параллельный рекурсивный алгоритм вычисления присоединенной матрицы // Вестник Тамбовского университета. Серия Естественные и технические науки. Тамбов, 2009. Т. 14. Вып. 1. С. 265-269.
13. Бетин А.А. Параллельный алгоритм вычисления комплексных корней уравнения // Вестник Тамбовского университета. Серия Естественные и технические науки. Тамбов, 2013. Т. 18. Вып. 1.
14. Рыбаков М.А. О нахождении общего и частного решений неоднородной системы обыкновенных линейных дифференциальных уравнений с постоянными коэффициентами / / Вестник Тамбовского университета. Серия Естественные и технические науки. Тамбов, 2012. Т. 17. Вып. 2. С. 552-565.
БЛАГОДАРНОСТИ: Работа выполнена при поддержке гранта РФФИ № 12-07-00755-а,
Поступила в редакцию 20 декабря 2012 г.
Rybakov М.A. PARALLEL COMPUTATION OF THE GENERAL SOLUTION OF THE INHOMOGENEOUS SYSTEM OF ORDINARY DIFFERENTIAL EQUATIONS WITH CONSTANT COEFFICIENTS.
We discuss a parallel algorithm to obtain numerical-analytical solutions of the inhomogeneous system of ordinary differential equations with constant coefficients. The solution of the system can be found with the required accuracy. This algorithm is effective for solving large systems of differential equations. The algorithm is part of the library of algorithms Mathpar.
Key words: inhomogeneous system of ordinary differential equations with constant coefficients, Laplace transform, parallel algorithm.
