книги / Метод продолжения решения по параметру и наилучшая параметризация в прикладной математике и механике
..pdf3.3. Жесткие системы |
101 |
точности вычислений 10-4 и начальном шаге интегрирования Ад = 0,1. Интересно отметить, что если из программы вычислений исключить про цесс простой итерации, то время счета увеличивается, т. е. наилучший вычислительный процесс обеспечивается совместным использованием итерационных формул (2.30) и (3.24).
Также была разработана программа, учитывающая в итерационном процессе не только элементы, стоящие на главной диагонапе матрицы Якоби, но и наибольшие по абсолютной величине один или два внедиагональных элемента. Так, например, если наибольшим внедиагонапьным элементом является элемент Ьу, то i -тую формулу итерационной схемы (3.24) следует подкорректировать слагаемым
(У ),т + f c / j f a l l l . W l ) ~ |
y j , m + l ) h b i j |
(hbu - 1)(Ы>И - |
1) |
Однако, проведенные вычисления показали, что такая корректи ровка не приводит к уменьшению времени счета, хотя наименьший шаг интегрирования уравнений увеличивается. Это, по-видимому, объ ясняется тем, что полученное преимущество по шагу интегрирования компенсируется усложнением алгоритма процесса.
Кроме того, для решения непреобразованной задачи (3.39), (3.40) была использована программа DLSODE, которая является одной из по следних версий программ серии Гир, предназначенных для интегрирова ния жестких систем уравнений. Эта программа взята из пакета программ ODEPACK, собранных Хиндмаршем [100]. С ее помощью задача была решена за 13с.
Большое отличие во времени счета программ РС1М и DLSODE объясняется тем, что последняя использует формулы интегрирования до 5-го порядка точности, имеет алгоритм выбора оптимального порядка BDF-формул и оптимального шага интегрирования.
Однако, задача (3.41) была решена при помощи программ РС2М за 70с., РСЗМ за 20с., РС4М за 17с. счета на PC АТ 286/287.
Анализируя численные результаты, можно сделать выводы:
—программа РС1М, конечно, значительно уступает по времени сче та программе DLSODE, но она не требует хранения информации, полученной на предыдущих шагах. Это может оказаться решаю щим при изучении многомерных задач механики сплошной среды, которые приводятся к системам обыкновенных дифференциальных уравнений высокого порядка;
—программа РС4М по времени счета уже незначительно отличается от своего зарубежного аналога, однако, несомненным ее преимуще ством является то, что при ее помощи могут быть проинтегрированы задачи, имеющие замкнутую интегральную кривую или интеграль ную кривую, содержащую предельные точки, а также могут быть
102 |
Глава 3. Жесткие системы дифференциальных уравнений |
решены некоторые уравнения, правые части которых обращаются в некоторых точках в бесконечность.
Рассмотрим решение модельной жесткой задачи, которую предло жил И. Бабушка [80, 37, 106]. С учетом того, что один интеграл этой задачи известен, ее можно представить в виде
dt |
ym |
dy2 |
= ~аУз> |
|
yi(0) =У2(0) = 1, |
У2 |
dt |
|
|||
dyi |
(су4 ~ 9a~ ~ |
+ (c + аУз)уз ~ bj У2 + 9' |
Уз(°) = о, |
||
dt |
|||||
dm |
a |
- d(y3 + У4), |
|
У4(0) = 0,1, |
|
. dt |
|
||||
У2 |
|
|
|
|
|
|
t € [0,1] |
a = d = 100, 6 = 0,9, |
с = |
1000. |
|
В |
[37] отмечается, |
что ненулевой спектр |
матрицы Якоби этой |
системы лежит по обе стороны от нуля на достаточном расстоянии. Из-за этого шаг аппроксимации в явных и неявных методах с течением времени не будет расти, оставаясь достаточно малым, и в неявных методах возникнут дополнительные сложности при итерационном методе нахождения решения.
После применения к задаче Л-преобразования, она была решена при помощи программы РСЗМ на PC АТ 286/287 за 24с. Наиболь шая ошибка при t = 1 была равна 0,04. Заметим, что при решении непреобразованной задачи такой точности достичь не удается.
Эта же преобразованная задача решалась на PC с 486 процессором. Для достижения той же точности программе DUMKA [37, 106] потребо валось 1,4с., а программе DLSODE [100] 0,03 с. При этом, если в про грамме DUMKA правые части системы уравнений вычислялись 36074 раза, то в программе РСЗМ они вычислялись 4 223 раз, а в программе DLSODE только 320 раз.
«Точное» решение было получено при помощи программ RADAU3,
DLSODE и при t — 1 было равно |
|
yi = 0,58367615856 • 105, |
у2 = 0,171327881975 • 10“4, |
уз = 0,1903613486 •10"5, |
у4 = 0,583667158388 • 104. |
3.4. Жесткие уравнения в частных производных
Ранее Л-преобразование использовалось при решении задач малой размерности (п < 4). Теперь исследуем решение жесткой системы обык новенных дифференциальных уравнений большой размерности (п > 10), к которой сводится решение уравнений в частных производных [33].
3.4. Жесткие уравнения в частных производных |
103 |
В качестве примера исследуем жесткую задачу [97] о движении прямоугольной упругой пластины при переезде ее автомобилем. Переме щения пластины W{x> y ,t) описываются следующими соотношениями
|
..aw |
a2w |
= f(x, у, t), |
|
D V 2V 2W + G ^ - |
+ ^ |
|
|
at |
+ at2 |
|
П = о = о, |
aw |
|
= o, v 4 r |M a n . a |
feo = o, n . , ) « n |
Здесь V2V2 — бигармонический оператор; D = 100, G = 1000 — параметры жесткости пластины и трения; V2 — оператор Лапласа, Ш — граница области П = {(я, у), 0 < * < 2, 0 ^ у < 4/3}, которая разбивается на ячейки точками г,- = ih, yj = jh , h = 2/9, t = 1,8,
3 = V 5-
Нагрузка f(x ,y ,t) приближается суммой двух гессианов, которые описывают движение в направлении оси х четырех колес
/(*, У, 0 = |
200(е-5^- *-2^ + |
У = У1, |
У = У4» |
|
0, |
У Ф у1, |
У ФУ4- |
||
{ |
Бигармонический оператор представляется в стандартном конечно разностном виде, использующим центральные разности
v2vV w(o = 2j - m |
- u +20W i j - m + u +w i+ 2 J + |
+ W i j - 2 - |
- 8W ij +l + W iJ + 2 + 2 W i + \ j - i + |
+ 2Wt+ij+i + 2W i-ij-i + 2Wi_ij+i)/h4.
Согласно граничным условиям, в контурных точках рассматривае мой пластины Wij(t) = 0, а в предконтурных значение прогиба в два раза меньше, чем в предшествующих точках, лежащих по нормали к контуру.
Полученная система обыкновенных дифференциальных уравнений
интегрировалась на отрезке [0, £], t € [0,7] при помощи программы
л
РС1, в которой из-за малых значений прогиба, не превосходящих 10 , контролировалась не абсолютная, а относительная погрешность вычисле ний, полученная делением абсолютной погрешности на сумму квадратов функций Wij. Оказалось, что применение преобразования в тринадцать раз сокращает время счета и требует в 16 раз меньшего числа обращений к вычислению правой части. А шаг интегрирования увеличивается более чем в 10 раз.
Глава 4
Дифференциально-алгебраические уравнения
Среди задач, решениями которых являются гладкие непрерывные однопараметрические множества, можно выделить дифференциально алгебраические уравнения (ДАУ), которые сочетают в себе особенности нелинейных алгебраических или трансцендентных уравнений с особен ностями ОДУ. Для таких задач можно поставить задачу Коши, кор ректная постановка которой связана с необходимостью решения систем нелинейных уравнений. В этой главе мы сформулируем и исследуем алгоритм численного продолжения решения задачи Коши для различных форм дифференциально-алгебраических уравнений.
4.1. Классификация систем ДАУ
Напомним, что системой неявно заданных обыкновенных диффе ренциальных уравнений называется система вида
$ (1Л У, *) = о, |
|
|
y(t) : R1 -* К", <€[*0, 11, Ф : Е |
2п+1 - К ”, у = % |
(4Л) |
|
at |
|
Если матрица дФ/ду невырождена, то эта система может быть разрешена относительно производных и представлена в явном виде (в нормальной форме Коши)
У = <р(У, t), f = (<Pi, • • •, <Рп)Т■ |
(4.2) |
Системой дифференциально-алгебраических уравнений называется система вида
Р(У, У, *, t) = 0, |
G(y, х, t) = 0, |
|
||
F : flj2n+m+1 |
R»> |
G : ! " +m+1 - » Rm, |
(4.3) |
|
y{t) : l 1 - » 1 ”, |
a(<): l |
1 - » l m, |
t € [<0> Г], |
|
состоящая из обыкновенных дифференциальных уравнений F и не дифференциальных соотношений G, в качестве которых обычно рассма триваются нелинейные алгебраические или трансцендентные уравнения.
4.1. Классификация систем ДАУ |
105 |
При формулировании системы (4.3) предполагается, что матрица dF/dy не является вырожденной. Очевидно, что к системе уравнений такого типа может быть сведена система (4.1) в случае сингулярной матрицы дФ/ду этой системы. Поэтому наряду с установившемся в по следнее время термином дифференциально-алгебраические уравнения, для обозначения систем вида (4.3) используется термин «сингулярные системы уравнений» (см., например [81, 82, 8]).
Существует глубокая связь между дифференциально-алгебраиче скими уравнениями и сингулярно возмущенными уравнениями, которые, как мы видели ранее, относятся к классу жестких систем обыкновен ных дифференциальных уравнений. В самом деле, когда в сингулярно возмущенных уравнениях мы имеем малый параметр е, то получаем жесткую систему уравнений. Если же этот параметр е устремить к нулю, то приходим к дифференциально-алгебраическим уравнениям.
Рассмотрим уравнение Ван-дер-Поля, которое изучалось в класси ческой работе [23], записав его в виде
ez" + (z2 —1)z' + 2 = 0.
Так как
+ (*2 - 1У = I |
+ ( у - * ) ) = |
то, введя обозначение у = ez' + z3/3 - z, имеем систему
которую в общем случае можно записать в виде
У = f(y> *)> sz = g(y, z).
Полученная система сингулярно возмущенных уравнений при малых значениях параметра е будет жесткой. Если теперь параметр е устре мить к нулю, то получаем систему дифференциально-алгебраических уравнений
или в общем случае
У1= /(».*). 9(y,z) = 0-
106 Глава 4. Дифференциально-алгебраические уравнения
Задача (4.4) имеет точное решение
z2
\n\z\ - — = t + C.
Это решение реализуется, если только начальные значения ле жат на кривой L, задаваемой вторым уравнением системы (4.4). Осо бый интерес представляют точки этой кривой, имеющие координа
ты у = ±2/3, |
z = |
т 1 |
в которых производная дг = dg/dz обращается |
в нуль. Ветвь |
-1 |
< z |
< 1 кривой L является неустойчивой (дг > 0), |
поэтому в отмеченных точках происходит перескок решения на другую устойчивую ветвь (gz < 0).
Системы дифференциально-алгебраических уравнений принято ха рактеризовать дифференциальным индексом. Рассмотрим это понятие вначале на примере системы (4.1). Система уравнений (4.1) имеет диф ференциальный индекс г, если г есть минимальное число дифференци рований
Ч у,у, О = о, |
* Ч у, уЛ) |
#4v,y,t) |
|
dt |
dtr |
||
|
переводящее систему (4.1) в явную систему (4.2).
Из определения следует, что если матрица системы (4.1) дФ/ду не является сингулярной, то эта система имеет дифференциальный ин декс, равный нулю. В этом случае система приводится к каноническому виду (4.2) без дифференцирований.
Дифференциальный индекс системы (4.3) определяется следующим образом. В силу того, что матрица dF/dy не сингулярна систему можно
представить в виде |
|
У= f(y, *, 0. G(y, *»0 = 0- |
(4.5) |
Эта система имеет индекс, равный единице, если после дифферен цирования второго векторного уравнения по t получаем систему вида
GiVy + G pi + Gtt = 0, |
(4.6) |
где GjX не сингулярна, поэтому система (4.5) может быть записана явно
внормальной форме Коши
ГУ = f(x ,y ,t)
\x = - G - ,(Gy/ + G>i).
Если же матрица G^ сингулярная, то выражение (4.6) принимает
вид
G,yf + G ,t = 0 ,
4.1. Классификация систем ДАУ |
107 |
и его следует дифференцировать еще раз относительно t, и так далее до тех пор, пока матрица, стоящая перед производной а , не станет обратимой.
Таким образом, индекс системы дифференциально-алгебраических уравнений (4.3) равен наименьшему числу дифференцирований по t, позволяющих определить а как непрерывную функцию у, х, t.
Заметим, что система уравнений (4.1) с несингулярной матрицей дф/ду, имеющая индекс равный нулю, после введения новой функции у = х может быть преобразована к системе типа (4.S)
у = х, Ф(у, х, t) = 0, |
(4.7) |
которая имеет уже индекс, равный единице, т.е. преобразование (4.1)
квиду (4.7) увеличивает дифференциальный индекс на единицу [77].
Вданной работе изучается численное решение задачи Коши системы дифференциально-алгебраических уравнений, имеющей индекс, равный нулю или единице. Однако, допускается, что в некоторых точках ин тегральной кривой производная G x может вырождаться. Такие системы наиболее часто встречаются в вычислительной практике.
Решение дифференциально-алгебраических уравнений является бо лее сложной задачей по сравнению с решением обыкновенных диффе
ренциальных уравнений. В [77] отмечаются следующие трудности:
1)начальные условия должны быть согласованными с недифференци альными соотношениями;
2)система линейных алгебраических уравнений, решаемая на каждом шаге интегрирования, является плохо обусловленной для мелких шагов. В [77] показано, что обусловленность системы равна 0(Л"), где v — индекс системы, h — шаг интегрирования;
3)ошибка ограничения при выборе шага интегрирования чувствитель на к несогласованности в начальных условиях и резкому изменению решения;
4)численное решение в большей степени зависит от точности аппрок симации итерационной матрицы, чем для обыкновенных дифферен циальных уравнений.
Предложенный здесь подход построения численного решения задачи позволяет ослабить некоторые из отмеченных трудностей.
Что касается согласованности, то начальные условия для системы (4.3) можно принять в виде
у(к) = Уо, х(к) = *о.
причем векторы уо и хо должны удовлетворять системе уравнений
G(yо, ®о, *о) = 0.
108 |
Глава 4. Дифференциально-алгебраические уравнения |
По-видимому, первой работой, в которой исследуется численное решение дифференциально-алгебраических уравнений, является рабо та [94]. В ней система уравнений, описывающая процессы, протекающие в электрических сетях, интегрировалась при помощи формул диффе ренцирования назад. Проблема заключалась в решении линейных отно сительно производных дифференциальных уравнений и алгебраических соотношений типа
А(у)у + В(у) = № , |
(4.8) |
где А(у) — матрица.
Внастоящее время помимо многочисленных статей, посвященных данной тематике, имеется несколько монографий [97, 81,82,8,77,95,96].
Вработе [99] описывается программа LSODI, предназначенная
для интегрирования задачи Коши для неявных систем обыкновенных дифференциальных уравнений, а в работе [114] программа интегрирова ния задачи Коши для систем вида (4.8), возникающих в электротехни ке. Обе программы используют версию программы GEAR, основанную на формулах дифференцирования назад. Программа DASSL [77, 108] разработана для решения начальной задачи для неявной системы ви да (4.1), имеющей индекс, равный нулю или единице. В соответствии с общим подходом, предложенным в [94], производная заменяется раз ностной аппроксимацией, использующей формулы дифференцирования назад, а полученная система нелинейных уравнений решается методом Ньютона—Рафсона, в котором матрица Якоби вычисляется не на каждом шаге, а по мере необходимости.
Для решения аналогичной задачи разработана программа STIFSP [14], в основу которой положена программа STIFF и метод равных корней, предложенный в [14]. Программа предназначена для работы с разряженными матрицами Ф|У Ф^.
Программа RADAU5 [97, 96] разработана для решения начальной задачи для системы типа (4.8), имеющей индекс более двух. В отличие от предыдущих программ здесь применяются неявные методы типа Рунге—Кугга. Получающиеся при этом системы нелинейных уравнений решаются модифицированным методом Ньютона. В [97] проводится также сравнение эффективности различных программ.
В [8] описывается программа SINODE, осуществляющая решение задачи Коши для системы уравнений типа (4.8) с вырожденной матрицей некоторыми численными методами (Рунге—Кутта, Адамса).
Численное решение задачи Коши для системы дифференциально алгебраических уравнений вида
»(<) = /(»> М ), x(t)= g(y,x,t)
рассматривается в [35]. Предлагается комбинировать неявный метод Эй лера с методом простой итерации или неявный метод Адамса с методом Ньютона.
4.2. Аргумент системы дифференциально-алгебраических уравнений |
109 |
Предложенный краткий обзор численных подходов к решению диф ференциально-алгебраических уравнений отражает несомненные дости жения в данной области, однако полученные результаты не ликвидируют те трудности решения задач рассматриваемого типа, которые были отме чены выше.
Подход, предложенный в данной главе, позволяет ослабить часть из этих трудностей. Как будет видно ниже, система линейных алгебраи ческих уравнений, получающаяся на каждом шаге процесса интегрирова ния, будет наилучшим образом обусловленной, а в силу выбора аргумента задачи ошибка будет менее чувствительна к резкому изменению решения.
Вследующем разделе этой главы задача Коши для системы диффе ренциально-алгебраических уравнений рассматривается с позиций ме тода продолжения решения по параметру и изучается вопрос о выборе наилучшего параметра.
Вмонографии [77] отмечается, что при разработке программы DASSL в систему дифференциально-алгебраических уравнений вводился
параметр и для отыскания решения использовался метод продолже ния, но вопрос о выборе наилучшего параметра продолжения решения не рассматривался.
4.2. Наилучший аргумент системы дифференциально-алгебраических уравнений
Рассмотрим задачу Коши для системы дифференциально-алгебра ических уравнений
Г F (V , У, *, 0 |
= 0, |
y (to ) = Уо, |
(4 9) |
||
I G(y, х, <) = 0, |
|
я(*о) = *о> |
|
||
F = |
|
<7 = |
|
Gm)T, |
t e r n 1, |
y (t ) = (2/1 (<). • • ■. yn(t)f, |
x(t) = (x, ( t ) , . . . , xm(t))T, |
||||
. |
dy |
f dy\ |
|
dyn\ T |
|
У |
dt |
\ dt |
’ |
dt ) |
|
Уо — (j/lO) • • • >УпО)Т > |
X0 —(x,0, .. . , Xmo) . |
Причем векторы уо, xo должны быть согласованными, т.е. удовлет ворять системе уравнений G(yo, ®о> *о) = 0.
Интеграл задачи (4.9)
/(у, X, t) = 0, /(уо, Xo, to) = о, / = |
fn+mf |
(4.10) |
задает в (n + m + l)-мерном евклидовом пространстве Rn+m+l интеграль ную кривую К , процесс построения которой может быть представлен
110 |
Глава 4. Дифференциально-алгебраические уравнения |
как процесс продолжения решения у = y(t), х = x(t) по параметру t. Та кой подход позволяет поставить вопрос о выборе наилучшего параметра продолжения решения системы (4.10), а значит, и наилучшего аргумента задачи (4.9).
Будем, как и в главах 1 и 2, вводить наилучший аргумент локально, т. е. в малой окрестности каждой точки интегральной кривой К . Чтобы найти наилучший аргумент, введем в окрестности рассматриваемой точки параметр д так, что
dfi = otidyi + Pjdij + 7<ft, |
i = T7n, j = |
l,m . |
(4.11) |
Здесь a,-, pj, 7 — компоненты ранее рассмотренного единичного |
|||
вектора a = (at|........а„, p i, ... ,p m, |
i f € Дп+т+|, |
задающего на |
правление оси, по которой отсчитыается аргумент д. Предполагается суммирование в произведениях по повторяющимся индексам в огово ренных пределах.
Уравнения продолжения решения |
задачи (4.10) получим, |
если |
в предположении дифференцируемости |
функций у,(д), х; (д), |
£(д) |
равенство (4.11) разделить на dfi, а первое из соотношений (4.10) продифференцировать по д:
OiVi,li + Pj*j,/i +7<,р — 1) U Уг,11+ } , X j Xj,p + fjtjx |
— 0. |
(4.12) |
Здесь приняты обозначения у{ф = dyt/dfi, f iVi = df/dyt, . . . . |
|
|
Однако такой подход неконструктивен, так как |
интеграл |
(4.10) |
до решения задачи (4.9) неизвестен.
Уравнения продолжения могут быть получены иначе. Линеаризуем вектор-функцию F относительно производных у,- в окрестности неко торых значений & = у*, полученных, например, на предыдущем шаге итерационного процесса или процедуры итерирования. Тогда имеем
F* + F*m(l/i ~ Vi) = * = ТГп-
Здесь вектор-функции F* и F*^ вычисляются при у,- = у*.
Принимая во внимание первое уравнение системы (4.12) и равенства Vi = if* = У*фЦ*ц, а также дифференцируя вектор-функцию G относительно д, получаем следующее представление уравнений продол жения
<*т,г +Р)Ч р |
= 1, |
~ К А ) ** = 0, |
(4.13) |
G,9tVi^ + G'XjH/t + |
= 0. |
Интегральная кривая задачи (4.9) может быть построена в результа те интегрирования системы дифференциальных уравнений, полученной