научная статья по теме СИМВОЛЬНО-ЧИСЛЕННОЕ РЕШЕНИЕ СИСТЕМ ЛИНЕЙНЫХ ОБЫКНОВЕННЫХ ДИФФЕРЕНЦИАЛЬНЫХ УРАВНЕНИЙ С ТРЕБУЕМОЙ ТОЧНОСТЬЮ Математика

Текст научной статьи на тему «СИМВОЛЬНО-ЧИСЛЕННОЕ РЕШЕНИЕ СИСТЕМ ЛИНЕЙНЫХ ОБЫКНОВЕННЫХ ДИФФЕРЕНЦИАЛЬНЫХ УРАВНЕНИЙ С ТРЕБУЕМОЙ ТОЧНОСТЬЮ»

программирование, 2013, no 3, с. 38-46

КОМПЬЮТЕРНАЯ АЛГЕБРА

У л . .'<

СИМВОЛЬНО-ЧИСЛЕННОЕ РЕШЕНИЕ СИСТЕМ ЛИНЕЙНЫХ ОБЫКНОВЕННЫХ ДИФФЕРЕНЦИАЛЬНЫХ УРАВНЕНИЙ С

ТРЕБУЕМОЙ ТОЧНОСТЬЮ

© 2013 г. H.A. Малашонок, М.А. Рыбаков

Институт математики, физики и информатики Тамбовского государственного университета им. Г. Р. Державина 392021 Тамбов, ул. Интернациональная, д. 33 E-mail: namalaschonok@gmail.com, mixail08101987@mail.ru Поступила в редакцию 02.07.2012 г.

Приводится символьно-численный алгоритм получения решения системы обыкновенных линейных дифференциальных уравнений с постоянными коэффициентами и составной правой частью. Алгоритм основан на применении преобразования Лапласа. Частью алгоритма является определение погрешности вычислений, достаточной для требуемой точности решения системы. Алгоритм эффективен для решения систем дифференциальных уравнений большого размера, позволяет выбирать методы решения алгебраической системы - образа преобразования Лапласа - в зависимости от ее типа, а тем самым оптимизировать эффективность решения исходной системы. Алгоритм входит в состав библиотеки алгоритмов системы МаЙраг [15].

1. ВВЕДЕНИЕ

Современные системы компьютерной алгебры предоставляют некоторые методы решения систем линейных дифференциальных уравнений -численных и символьно-численных. Нередкое обращение к преобразованию Лапласа в применении к решению дифференциальных уравнений в последние годы (см. например, [2, 3, 13, 14]) не случайно. Это приводит к получению аналитического решения, что особенно важно для приложений.

Мы представляем символьно-численный алгоритм, основанный на применении преобразования Лапласа.

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

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

Преимущества построенного алгоритма в следующем:

• размер исходной системы и порядок производных не играют существенной роли;

и т.д.) выбирается наиболее эффективный способ решения алгебраической системы -тем самым регулируется эффективность решения этим методом исходной системы дифференциальных уравнений;

ся выбором соответствующей погрешности численного фрагмента вычислений.

2. АЛГОРИТМ РЕШЕНИЯ СИСТЕМ ДИФФЕРЕНЦИАЛЬНЫХ УРАВНЕНИЙ МЕТОДОМ ПРЕОБРАЗОВАНИЯ ЛАПЛАСА

Рассмотрим систему

хУ = fl, l = 1,...,n, alkj e R, (1)

n N

S^S^J Jk)

х3 = Л, 1 = 1,...,п, ]=1 к=0

п дифференциальных уравнений порядка N ^ - старший по всем неизвестным функциям) п неизвестных функций с начальными условиями х^ (0) = х^, к = 0,..., N — 1. В правой части системы функции, приводимые к виду

fi(t) = fi(t), ti <t<ti+\ i = 1,...,Ih t¡ = 0,t¡l+1 = ж,

(2)

где

fi(t) = E pis(ty

к t

i = 1,... ,Il, l = 1,... ,n,

s=l

тж Г>г __Лг +т

и Г1в(ь) = 2^,т=0 свть ■

Здесь обозначены через х^= 1,... ,п неизвестные функции аргумента > 0, через производные этих х^ к = 0,..., N порядка к.

Предлагаемый нами алгоритм представим в трех основных этапах.

Этап 1. Преобразование Лапласа данной системы дифференциальных уравнений.

Обозначим образы преобразования Лапласа функций х^ (Ь) и Л (Ь) соответственно через X^ (р) и ^ (р).

Этап 1.1. Преобразование Лапласа левой части системы (1) с учетом начальных условий выполнятся формально согласно свойствам преобразования Лапласа. Получим

n N n N-l

E E akj pk Xj (P)—E E djk (p)x0j, l = 1,...,n j=1 k=0 j=1 k=0

где

djk (P) =

N -1

ai

/ j ™i+1,j P i=k

i-k

(3)

I - номер уравнения.

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

Хевисайда принимающей значение 1 для

0

остальных.

Представим /¡(Ь) в виде суммы функций

т = фКЬ—ьк) = е —ьк вь1(—).

s=1

Здесь ris(t - tk) = Pís(t) and ^\s(t — tf) =

Sj=o Yu¡j (t — tk)j, а коэффициенты y™* ВЬ1ЧИС~

j=0 llsj

ляются по формуле

Mks-j

YÜ3 = E

r=0

lk I i + r ,k\r cs,i+r\ m l(tl ) .

Окончательно, функция fl(t) представляется в виде

Ii-1

fl (t) = E m—n)v(t—n)—m—t+1)v(t—t+1)]+

i=1

ф1 (t — til )n(t — til).

Этап 1.3. Обозначим через Ф},г(р) образ функции ф1(t—t¡)n(t—Ц). Так как образ преобразования Лапласа функции (t—t*)nea(t-t )n(t — t*)

n!

Равен (p-n)„+1

fl (t) приводит к следующему:

e t*p, то преобразование Лапласа

fi(p) = е [*Г(р) — Ф

i=1

(р) +<PÍl,li (р). (4)

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

Этап 2. Решение линейной алгебраической системы с полиномиальными коэффициентами.

п

нейных уравнений с полиномиальными коэффициентами относительно неизвестных функций хз,3 = 1,...,п

n N n N-1

E E alkj pk Xj (p) = EE djk (p)xkj+fi (p), j=1 k=0 j=1 k=0

l = 1, . . . , n

Этап 2.1. Обозначим

п N-1 £ £ =

j= 1 0=0

Систему (5) в этих обозначениях приводим к виДУ

са для функций (7) записываются формально:

п N

£ £ а^ркX(р) = ЗД + ВД, 1 = 1, j=l 0=0

, п.

А,

тк

0—атр

т к

(Р - Pjk)втк

(7)

Атк

?? (втк - 1)Г

X № =

(6)

Для каждого 1 = 1,..., п выражения в правых частях системы (6) приводим к общему знаменателю.

Этап 2.2. Напомним, что в системе (6) полиномы основной матрицы имеют коэффициентами действительные числа (что видно из (1)). Метод решения выбирается в зависимости от вида системы (например, р-адический [9], модулярные методы, методы, основанные на детерминантных тождествах [7, 8] в случае, если коэффициенты рациональные, и другие).

Этап 3. Обратное преобразование Лапласа.

Этап 3.1. Обозначим определитель основной матрицы системы (6) через ^(р). Решение системы (6), т.е. функции Xj(р), 3 = 1,...,п, пред-

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

Этап 3.2. Подготовка к обратному преобразованию Лапласа сводится к разложению рациональных коэффициентов перед экспонентами на простые дроби А/(р — р*)^, р* е С. Для этого требуется найти корни многочлена ^(р). Это именно тот шаг алгоритма, который требует численных методов. На этом этапе делается оценка точности этих вычислений (см. Главу 2). В итоге каждая функция Xj (р) представляется в виде

3 = 1,...,п.

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

3. ОЦЕНКА ТОЧНОСТИ

На Этапе 3.2 алгоритма численно находятся корни многочлена ^(р). Задача - определить, какая погрешность вычисления корней достаточна, чтобы получить требуемую точность решения системы дифференциальных уравнений.

Будем рассматривать все функции и делать вычисления и оценки на отрезке [0,Т], где Т > для всех 1 = 1,...п, и достаточно велико для поставленной задачи. Обозначим через Хг(£), 1 = 1,...п, решение системы (1), полученное с использованием приближенных корней ^(р). Будем говорить, что решение исходной системы на [0, Т] имеет точность е, если

тах*е[о,т]|х(¿) - X(¿)| < е, 1 = 1,... ,п. (8)

Требуется найти погрешность А вычисления корней ^(р), достаточную для получения е

Обозначим через Т(р) минор (¿,1) матрицы системы (6). Решение XI(р) системы (6) может быть записано в виде

XI (р) =

А(р) ДР):

где

А(р) = £ И(р)+ $ (р)] Т(р).

г=1

Обозначим через рг точные корни и через Этап 3.4. Оригиналы преобразования Лапла- р*,г = 1,... ,п, приближенные корни полинома

В(р). Если\рГ — р*\ < А, то погрешность вычислений меньше А.

Теорема Для любого е существует А такое, что при условии \рГ — р*\ < А, имеет место (8) Доказательство

Рассмотрим многочлен В(р+ Аега), а € [0, 2п]. Обозначим

X (р) =

вы

(9)

Б(р + Аега)' А

ся (8).

Оценим оригинал разности

Рг(р) Рг(р) Б(р) Б(р + Аега).

Линейность преобразования Лапласа позволяет оценить независимо оригиналы

£— £Ш- ™

г=1

г=1

Б(р + Аега)

(10)

Теперь оценим (¿). Функция 0\(Ь) может быть записана в форме

Я I "т

ош) = Т. |Е в

Г=1 \"=1

еРт *

(16)

где рг - корн и В (р) кратн ос ти ц,Г, Я - количество различных корней.

Коэффициенты Вгг^ можно вычислить и оценить как Тейлоровские коэффициенты функции (р — Рг )"т В(рр) в 5-окрестности точки рг. Пусть 5 < 1/2тт 1рг — р3\,т, в = 1,... Я, 5 < 1. Мы должны получить А < 5, в противном случае А=5

Из неравенств Коши для тейлоровских коэффициентов следует:

1К»1 < (т*хр-рт\=* (р — рг Г Щ)) /5".

Рассмотрим многочлен В(р) = ^т=о сири ■ Так как

£ ШЦ—£ т- Т'м

Б(р + Аега)'

(11)

г=1 у г=1

Запишем выражение (10) в форме:

£ Рг(р)

г=1

Т(р)

ТМ

Б(р) Б(р + Аега)

(12)

Обе функции - (р) и разность, стоящая в скобках - имеют оригиналы. Оригинал Ф} (Ь) выражения (10) можно найти как сумму сверток Л (Ь) и оригинал а Щ(Ь) функции

Т (р)

Т (р)

Б(р) Б(р + Аега)' Представим выражение (13) в виде

(13)

Т (р) Т} (р + Аега) + Т} (р + Аега) — Т'} (р)

Б(р) Б(р + Аега)

Б(р + Аега)

~ (14)

Обозначим через 0\(Ь) и 0\(Ь) оригиналы функ-

- ТКр) Т(р) т!(р+Ае™) тт

цИИ ч ~Щр) — в(р+Ае*а) соответственно. Из

свойств преобразования Лапласа следует, что:

?(Ь) = (1 — е ) 0}(Ь).

(р — рг

Т}(р)

В(р)

\т (р)1

\Cmn\Il 1^=г 1р — рг

то нам нужно оценить \Т\(р)\ в 5-окрестности рГ. Обозначим через р радиус круга, содержащего все нули В(р). В качестве та кого р можно взять число

р = тах < 1;

\

Для дальнейшего прочтения статьи необходимо приобрести полный текст. Статьи высылаются в формате PDF на указанную при оплате почту. Время доставки составляет менее 10 минут. Стоимость одной статьи — 150 рублей.

Показать целиком