книги / Метод продолжения решения по параметру и наилучшая параметризация в прикладной математике и механике
..pdf4.2. Аргумент системы дифференциально-алгебраических уравнений |
111 |
после разрешения уравнений продолжения (4.13) относительно произ водных с учетом начальных условий
У<(0) = У«о, * # ) = *,(,, <(0) = *0- |
(4.14) |
Здесь предполагается, что аргумент ц отсчитывается от начальной точки задачи (4.9). Обусловленность системы (4.13) зависит от выбора аргумента ц, который определяется вектором а. Структура этой систе мы полностью совпадает со структурой системы (1.37), рассмотренной при доказательстве теоремы 1 о необходимых и достаточных условиях выбора наилучшего параметра продолжения решения системы нели нейных алгебраических или трансцендентных уравнений, содержащих параметр. Таким параметром, доставляющим наилучешую обусловлен ность системе продолжения, является длина дуги А, отсчитываемая вдоль кривой множества решений системы (4.10), которая в данном случае является интегральной кривой К задачи (4.9), а параметр про должения системы (4.10) — аргументом задачи (4.9). Аргумент р = А, обеспечивающий системе уравнений продолжения (4.13) наилучшую обусловленность, будем называть наилучшим. Напомним, что в качестве меры обусловленности принимается величина определителя системы, деленная на произведение квадратичных норм^ его строк. При доказа тельстве теоремы 1 было показано, что в случае выбора наилучшего аргумента ошибки, возникающие при численном решении задачи, будут оказывать на решение наименьшее влияние.
Согласно правилу Крамера, решение системы (4.13) в данном случае можно представить в виде
^У» _ |
Д | |
dxj _ A n +j |
dt |
A n +m + 1 |
|
||
dX ~ |
Д ’ |
dX ~ |
A ’ |
dX ~ |
Д |
’ |
(4.15) |
|
|
i = TTn, |
j = |
1, m, |
|
|
|
где A — определитель |
системы, |
A*. |
= ( - 1)*+1б* |
(к = |
1, п+тп+1), |
6k — определитель, получающийся при вычеркивании в матрице пос
ледних п + тп уравнений системы А:-го столбца. |
Эти определители |
удовлетворяют равенству типа (1.46), которое можно записать в виде |
|
Д2 = Д*Дк, к = 1, п+тп+1. |
(4-16) |
Это равенство показывает, что квадратичная норма правой части си стемы обыкновенных дифференциальных уравнений (4.1S) всегда равна единице. Если аргумент А отсчитывать от начальной точки задачи (4.9), то начальные условия примут внд (4.14).
Таким образом, опираясь на результаты главы 1, доказана следующая теорема.
112 |
Глава 4. Дифференциально-алгебраические уравнения |
Теорема 4. Для того, чтобы задачу Коши (4.9) для системы диф ференциально-алгебраическихуравнений сформулировать относительно наилучшего аргумента, необходимо и достаточно выбрать в качестве такого длину дуги А, отсчитываемую вдоль интегральной кривой за дачи. При этом задача (4.9) формулируется в виде (4.15), (4.14),
аправые части системы (4.15) удовлетворяют равенству (4.16)j.
Вследующих разделах, основываясь на этой теореме, сформулируем алгоритм и опишем программы решения некоторых систем дифферен циально-алгебраических уравнений. Эффективность работы программ изучим при решении некоторых тестовых задач для дифференциально алгебраических уравнений. Отметим, что в [66, 30, 67] рассматривались некоторые подходы к решению таких задач, но здесь описывается более совершенный алгоритм.
4 .3 . Явно заданные дифференциально-алгебраические уравнения
Под явно заданными уравнениями будем понимать систему, определяющую задачу Коши вида
| ! = /< » ,* ,« > , |
y(to) = У0, |
(4.17) |
|
1 G(y, х, t) = 0, |
x(to) = Xo, |
||
|
|||
у : R1 — * R", х : R1 — ►R1m, / : Rn+m+l — ►R", |
|
||
G : Rn+m+l —+ Rm, |
G(yo, XQ, to) —0. |
|
Сформулируем задачу (4.17), которая, очевидно, является частным случаем задачи (4.9) относительно наилучшего аргумента А. Пусть функ ции у = у(А), х = х(А), t = <(А) являются дифференцируемыми. Введем обозначения
<*У= у |
’ d \ |
= х |
— = Т |
(4.18) |
d \ |
|
’ dX |
||
Y = (Yu ... ,Y a)T, |
X |
= (X u . . . , X m)T. |
|
|
После дифференцирования по А |
вектор-функции G, |
принимая |
во внимание соотношения (4.18) и смысл наилучшего аргумента, запи шем систему (4.17) в виде
Yi ~ fiT —°. |
|
< G,yY i + G ^X j + GrtT = 0, |
(4.19) |
YiYi + X jX j + T T ^ X , |
j ~ 1,m- |
4 3 . Явно заданные дифференииально-алгебраические уравнения |
113 |
Из-за последнего уравнения эта система является нелинейной отно сительно функций Y, Х , Т . Однако, ее можно представить в линейном виде, если использовать решение, найденное на предыдущем (к - 1)-ом шаге. Для этого перепишем систему (4.19) в виде
' у}к) - |
= О, |
|
|
|
< а д (4) + G ^ x f + |
= 0, |
к = 1 ,2 ,.... |
(4.20) |
|
yfr-OyW + |
|
+ т(* _1)т(*) = |
1, |
|
Обозначим через Z ^ = (У ^ , Х^к\ Т ^ ) т вектор размерности n+m+1. В силу структуры системы (4.20) этот вектор является касатель ным к интегральной кривой К задачи (4.17) в точке, соответствующей к-му шагу, поэтому последнее уравнение системы (4.20) представля
ет собой скалярное произведение векторов и Z^k~l\ касательных к интегральной кривой на к-ом и (А:-1)-ом шагах. Это уравнение утвер
ждает, что проекция вектора Z ^ на направление вектора Z^~1^ рав на единице (см. рис. 4.1).
Заменив в (4.20) неизвест ный вектор Z ^ на извест
ный z(k~l\ мы обеспечива ем такой локальный выбор аргумента, который являет ся близким к наилучшему.
Очевидно, что вектор
Z ^ , удовлетворяющий си стеме линейных алгебраи ческих уравнений (4.20), вообще говоря, не будет
единичным, как этого требует система (4.19), поэтому после решения
системы (4.20) найденный вектор Z ^ |
следует нормировать по формулам |
|
Z)(*) |
hJ 1, n + m + 1. |
(4.21) |
Z-(i) = |
У?
При этом мы получаем решение системы (4.19). Далее звездочку в формуле (4.21) опускаем.
Так как начальная точка обычно не бывает предельной, то в качестве
начального приближения вектора Z можно взять вектор |
|
Я(0) = (0, . . . , 0, 1)г . |
(4.22) |
114 |
Глава 4. Дифференциально-алгебраические уравнения |
Полагая, что аргумент А отсчитывается от начальной точки задачи (4.17) , можно предложить следующий алгоритм ее решения.
Находим решение системы дифференциальных уравнений (4.18), удовлетворяющее начальным условиям
У(0) = уо, х(0) = х0, Щ = t0. |
(4.23) |
Правые части системы (4.18) определяем из решения методом Гаус са системы линейных уравнений (4.20) с последующей нормировкой по формулам (4.21).
При таком подходе можно обходить не только ситуации, возникаю щие из-за обращения в нуль якобиана G ^ , но и решать системы (4.17), у которых правые части / дифференциальных уравнений обращаются в некоторых точках в бесконечность. Чтобы преодолеть возникающие при этом трудности, достаточно переписать, если это возможно, пер вые п уравнений системы (4.19) в виде
QaYa - РаТ = 0.
Здесь а = 1, я и по индексу нет суммирования, а функции Q„, Ра не обращаются в бесконечность.
Если же правые части / дифференциальных уравнений системы (4.17) не обращаются в бесконечность, то размерность системы (4.20) можно уменьшить на п единиц, записав ее в виде
|
G,Xix f + (Git + <7,ft/,)T(*> = 0, |
|
Зная решение системы (4.24) значения |
' определяем по формулам |
|
у/** = |
после чего следует произвести нормировку вектора |
согласно соотношениям (4.21) и использовать полученные значения как правые части системы (4.18).
В соответствии с этими двумя алгоритмами были разработаны на ал горитмическом языке ФОРТРАН программы DA1EXG, в которой решает ся система (4.20), и DA1EXP, решающую систему (4.24). Интегрирование системы дифференциальных уравнений (4.18) осуществляется при помо щи программы РС1, а система линейных уравнений решается методом
Гаусса. |
___ |
|
При разработке программ неизвестные функции j/i, Xj, t, |
i = |
l,n , |
j = l,m , задавались единым массивом уи, к = l,n + m + l; |
yt = |
2h, |
yn+j = Xj, yn+m+i = t, а коэффициенты системы линейных уравнений (4.20) как элементы а« матрицы порядка n+m+1.
4.3. Явно заданные дифференциально-алгебраические уравнения |
115 |
|
|
Примеры |
|
В качестве первого примера рассмотрим задачу [82] |
|
|
dy _ |
2 |
|
< й ~ У + Х > у(0) = 4, х(0) = 2, |
(4.25) |
|
{i/y - х |
— 0, |
|
имеющую точное решение
у —4е2*, х = 2е‘.
Интегрирование задачи (4.25) на отрезке t € [О, I] при помощи программы DA1EXG требовало на 25% большего времени счета по срав нению с программой DA1EXP. Точность вычислений была практически одинаковой.
В качестве второго примера исследуем задачу [77, 96]
Г У|+У*У2 + (1+»?)У2 = 0, |
. _ * у |
|
1у1+7<У2 = У(0> |
У~ |
(4.26) |
d t’ |
||
где g(t) — дифференцируемая функция, у — число. |
|
|
Задача имеет точное решение |
|
|
Г Vi = 9 { t ) + n t g ( t ) , |
|
(4.27) |
I У2 = |
|
|
|
|
|
Погрешность вычислений оценивалась по формулам |
||
Дщ —Ут ~ Ут> тп = I, 2t |
(4.28) |
где yin — численное решение, ут — точное решение (4.27).
В [77, 96] отмечается, что эта задача является трудной для любого численного метода. Она имеет индекс 2, а численные методы на высоких индексах часто неустойчивы. Так, решение задачи неявным методом Эйлера неустойчиво, если у < -0,5.
Применим для решения задачи (4.26) вышеизложенный алгоритм, понизив предварительно ее индекс и используя метод j -продолжений (см. [8]). Для этого второе уравнение (4.26) дважды продифференцируем по t. С учетом первого уравнения получаем
Г У1 +17*У2 + |
0 + 1 7 )У 2 = 0, |
(4.29) |
I «2 = -§(!)■ |
|
|
|
|
|
Решим эту задачу при начальных условиях |
|
|
У!(0) = 9(0), |
уг(О) = -д (0). |
(4.30) |
116 |
Глава 4. Дифференциально-алгебраические уравнения |
|
Задача (4.29), (4.30) решалась при g(t) = е* и rj = 1,77 = -2 на От |
резке t Е [0,1]. При практически одинаковой точности вычислений программа DA1EXG требовала на 30% большего времени счета. При по мощи программы DA1EXP в момент t — 1 были получены следующие результаты (tC4 — время счета)
>7 = |
1, |
£сч —57с, |
Д[ |
=0,3- 10-2 |
Д2 = |
-0,73 |
• 10~3, |
V - |
“ 2, |
—- 68с, |
Д[ |
= -0 ,2 8 -10"2, |
Д2 = |
-0,57 |
• 10,-3 |
Таким образом, если нет необходимости, то лучше проводить вычи сления при помощи программы DA1EXP.
Рассмотрим решение задачи, в дифференциальное уравнение кото
рой производная входит нелинейно, |
|
|
y - t - y + lny = 0, у(1) = е. |
(4.31) |
|
Ее аналитическое решение у = е*. Точность вычислений оцениваем |
||
по формуле |
|
|
Д = у - е ‘, |
|
(4.32) |
где рассматриваются численные значения функций. |
|
|
Уравнение (4.31) имеет также особое решение у = t + 1, |
именно |
|
поэтому решение отыскивается при 1 ^ t ^ |
2, где интегральные кривые |
|
не соприкасаются. |
|
|
Преобразуем задачу (4.31) к виду (4.17) |
|
|
dy |
|
|
dt = х, |
У(1) = е, |
|
{ y - t - х + \пх = 0, |
ж(1) = е. |
|
Эта задача решалась при помощи программы DA1EXP при началь ном шаге интегрирования 0,001 с точностью 10"5 на PC АТ 286/287 31с.
Вмомент t = 2 ошибка, подсчитанная по формуле (4.32) была равна
Д= 0,28 • 10-2. Начальный шаг интегрирования трижды удваивался.
4.4. Неявно заданные обыкновенные дифференциальные уравнения
Рассмотрим задачу
ЖУь- -,Уп, Уи---,Уп) = 0,
Vi(to) = Ко> * = !.» . f |
= (fu • • • >fn)T, Ш= |
^ ^ |
Если матрица Якоби df/ dy |
не является сингулярной, то рассма |
триваемая система неявных уравнений является системой дифферен циально-алгебраических уравнений индекса нуль.
4.4. Неявно заданные обыкновенные дифференциальные уравнения |
117 |
Введением новых переменных х,- = у{ эта задача преобразуется в задачу с расширенным пространством решений (4.17). Кажется, что это снимает проблему решения задачи (4.33), так как алгоритм решения задачи (4.17) разработан и подробно описан в предыдущем разделе.
Однако, как показывает пример, приведенный ниже, такое преобра зование не всегда оправдано и может даже привести к непреодолимым вычислительным трудностям.
Рассмотрим задачу Коши для вырожденного уравнения Ван-дер-
Поля (4.4), записанного в неявном |
виде |
|
|
( l - V ) ^ f ~ У |
= °’ |
у(0) = 2‘ |
(4,34) |
Частный интеграл этой задачи |
можно представить в виде |
|
|
I n f - ^ |
+ 2 - f |
= 0. |
(4.35) |
График этой непрерывной функции в осях (у, t) показан на рис. 4.2 сплошной линией.
Если же задачу (4.34) преобразовать к виду (4.17), в каком она приводится в [96], то получим
dy _ т |
У(0) = 2, |
|
< dt |
||
(4.36) |
||
( 1 - у 2) х - у = 0, |
x(0) = - f . |
118 |
Глава 4. Дифференциально-алгебраические уравнения |
Интегральная кривая этой задачи в пространстве переменных t, у, х, |
|
пунктирная кривая |
на рис. 4.2, имеет разрыв вдоль прямой у -р 1, |
t = 3/2 - In 2 и, естественно, при интегрировании задачи (4.36) любым численным методом в окрестности этой прямой возникнут непреодоли мые вычислительные трудности.
Приведенный пример показывает, что следует быть осторожным при переходе от задачи (4.33) к задаче с расширенным пространством решений (4.17), хотя в отдельных случаях такой переход допустим и оправдан, несмотря на то, что он увеличивает индекс системы на еди ницу.
Рассмотрим алгоритм численного решения задачи (4.33) без пре образования ее к виду (4.17). Ясно, что задача (4.33) является частным случаем задачи (4.9). Сформулируем задачу (4.33) относительно наилуч шего аргумента.
Пусть величины у,-, t являются функциями наилучшего аргумента Л, отсчитываемого от начальной точки задачи (4.33). Введем обозначения
dyj |
Yh |
i = l,n, |
(4.37) |
|
d\ |
||||
|
|
|
правые части которых, учитывая смысл наилучшего аргумента, удовле творяют равенству
YiYi+T2 = L |
(4.38) |
Линеаризуем систему (4.33) относительно функций j/y, а равенство (4.38) — относительно функций и Г. Тогда с учетом соотношений у,- =
Yi/Т получаем систему линейных уравнений относительно функций Y ^
и Т^к\ вычисленных на fc-ом шаге итерационного процесса
( Y}k~l)Y^k) + т(*_1)г(*) = 1,
(4.39)
\ т (к_1)/*-.у;(*) + (/*г(‘_1) - f*yy}k~l))T^k) = о.
Здесь звездочкой помечены векторные функции, вычисленные при
* = ! ? - ^/Т^к *). Систему (4.39) рекомендуется записывать таким обра зом, чтобы, по возможности, отсутствовали соотношения, обращающие ся в бесконечность, и деление на Т.
Если начальная точка не является предельной, то начальное значе ние вектора Z — (Y\ , . . . , Yn, Т)т можно взять в виде
# (0) = (0,..., 0,1). |
(4.40) |
Таким образом, проблема заключается в интегрировании систе мы дифференциальных уравнений (4.37), удовлетворяющей начальным условиям
V i( 0) = y i о, i(0) = fo. |
(4.41) |
4.4. Неявно заданные обыкновенные дифференциальные уравнения |
119 |
Правые части системы (4.37) определяются из решения системы (4.39) с последующей нормировкой по формулам типа (4.21).
Понятно, что последние п уравнений системы (4.39) определя ют процедуру Ньютона—Рафсона, поэтому систему уравнений (4.39) следует решать до тех пор, пока итерационный процесс не сойдется с принятой точностью е : | | —Z^k~x\ < е. Заметим, что при ис пользовании при интегрировании системы ОДУ (4.37) программы РС1, это условие реализуется благодаря методу прогноза и коррекции. В ка честве начального приближения итерационного процесса принимается решение, вычисленное на предыдущем шаге процедуры интегрирования.
Очевидно, система (4.39) будет иметь наиболее простой вид в том случае, когда система дифференциальных уравнений (4.33) линейна относительно функций у,-, т. е. имеет вид
При этом, как ранее отмечалось, несмотря на приближенность последнего уравнения системы (4.39), после нормировки получается решение нелинейной системы
не зависящее от Yfk
В соответствии с изложенным алгоритмом была разработана вычи слительная программа DE1ILN, в которой интегрирование задачи (4.37), (4.41) осуществлялось с помощью программы РС1, а система (4.39) решалась методом Гаусса.
Для решения нелинейной задачи (4.33) можно предложить дру гой алгоритм, не требующий линеаризации уравнений, но требующий дополнительного дифференцирования. Рассмотрим его.
Принимая во внимание равенства гц = Y i / T , перепишем уравнения (4.33) таким образом
чтобы, по возможности, отсутствовали соотношения, обращающиеся в бесконечность, и деление на Y i, Т .
После дифференцирования выражений (4.38), (4.42) по Л получаем систему линейных уравнений
(4.43)
120 |
Глава 4. Дифференциально-алгебраические уравнения |
относительно функций
(4.44)
и проблема заключается в интегрировании системы дифференциальных уравнений (4.37), (4.44), удовлетворяющей начальным условиям (4.41) и условиям
Yi(0) = Yi0, Г(0) = Г0) t = 1, п. |
(4.45) |
Правые части уравнений (4.44) определяются в результате решения линейной системы (4.43), а начальные значения (4.4S) функций Y$, 3"о находятся из решения следующей системы
F(t0, 2/Ю) ■■■>2/п0) -ГЮ) ■• • >-РпО) 2о) —0,
(4.46)
ВД о + 2о = 0.
Всоответствии с изложенным алгоритмом разработана вычисли тельная программа DE1IN, в которой интегрирование систем диффе
ренциальных уравнений производилось посредством программы РС1, а система линейных алгебраических уравнений решалась методом Гаусса.
Примеры
В качестве примера рассмотрим решение задачи (4.34). Погрешность вычислений подсчитывается по формуле, стоящей слева в равенстве (4.3S), в которой используется численное решение. Задача интегрируется от точки А до точки В (рис. 4.2), соответствующей значению t —-2 .
Решение задачи при помощи программ DE1ILN было получено за 28с. Ошибка А при t = - 2 была равной А —0,33 ■10~2. Начальное значение вектора Z принималось в виде (4.40).
Эта же задача при тех же условиях решалась при помощи програм мы DE1IN 43 с. Погрешность А составила значение А = 0,27- 10—1. Начальные значения YQ = —2/у/\Ъ, То = З/уДЗ определялись из реше ния нелинейной системы уравнений
(1 - yl)Yo - уоТо = 0,
{ Y i + T i = l.
Отметим, что уравнение (4.34) может быть записано в виде, разре шенном относительно производной. Если применить А-преобразование,