CC BY
32
Вычислительные технологии
Том 10, № 2, 2005
НЕОДНОРОДНЫЙ МЕТОД ВТОРОГО ПОРЯДКА ДЛЯ ЖЕСТКИХ СИСТЕМ*
E. A. Новиков Институт вычислительного моделирования СО РАН
Красноярск, Россия e-mail: [email protected]
A method of the second order of accuracy is constructed for stiff additive systems of ordinary differential equations. The error estimation and control inequalities for computation accuracy and stability of the numerical scheme are obtained. The results of calculations are presented.
Введение
С целью численного решения задачи Коши для жестких систем обыкновенных дифференциальных уравнений
У = f (t,y), y(to)= Уо, to < t < tk, (1)
обычно применяются A-устойчивые или L-устойчивые методы. Здесь y и f — гладкие вещественные N-мерные вектор-функции; t — независимая переменная. В случае большой размерности задачи (1) для методов с неограниченной областью устойчивости общие вычислительные затраты фактически полностью определяются временем вычисления и обращения матрицы Якоби системы (1). Во многих алгоритмах используется "замораживание" матрицы Якоби, т. е. применение одной матрицы на нескольких шагах интегрирования. Это позволяет значительно уменьшить вычислительные затраты. Наиболее естественно это осуществляется в итерационных методах решения обыкновенных дифференциальных уравнений, где данная матрица не влияет на порядок точности численной схемы, а только определяет скорость сходимости итерационного процесса. Такой подход широко применяется при реализации полуявных и неявных методов типа Рунге — Кутты, многошаговых методов типа Адамса и Гира (см., например, [1]). Однако для безытерационных методов, к которым относятся методы типа Розенброка [2] и их различные модификации [3, 4], вопрос о замораживании или какой-либо другой аппроксимации матрицы Якоби значительно более сложный. В таких методах матрица Якоби влияет на порядок точности численной схемы и поэтому какие-либо ее возмущения могут приводить к потере порядка точности численной формулы. Следует отметить, что безытерационные методы просты с точки зрения реализации на ЭВМ и, как следствие, привлекательны для многих вычислителей.
* Работа выполнена при финансовой поддержке Российского фонда фундаментальных исследований (грант № 02-01-00523).
© Институт вычислительных технологий Сибирского отделения Российской академии наук, 2005.
С другой стороны, задачу (1) можно записать в виде [5]
у' = [/(¿,у) - Ву] + Ву, у(£о) = Уо, ¿0 < £ < ^, (2)
где В есть некоторая аппроксимация матрицы Якоби. Предполагая, что вся жесткость сосредоточена в слагаемом Ву, выражение в квадратных скобках можно интерпретировать как нежесткую часть. Если при построении безытерационных методов учитывать этот факт, то такая постановка задачи позволяет, в частности, использовать в алгоритмах интегрирования замораживание матрицы Якоби, которая может вычисляться как аналитически, так и численно. Для некоторых задач в качестве матрицы В можно использовать симметричную часть матрицы Якоби или ее диагональную аппроксимацию.
В настоящей работе построен четырехстадийный метод второго порядка точности, допускающий различные виды аппроксимации матрицы Якоби. Получены оценка ошибки и неравенство для контроля точности вычислений. Приведены результаты расчетов, подтверждающие работоспособность и эффективность алгоритма интегрирования.
1. Численная схема для автономных задач
Рассмотрим задачу Коши для автономной системы обыкновенных дифференциальных уравнений следующего вида:
У' = ^(У)+ ^У), У(£о) = Уо, ¿о < £ < ^, (3)
где У, и д — гладкие вещественные Ж-мерные вектор-функции; £ — независимая переменная. Далее будем полагать, что вся жесткость сосредоточена в функции д(У), а <^(у) есть нежесткая часть. Для численного решения (3) рассмотрим четырехстадийный метод вида
4
Уп+1 = Уп + ^рЛ, ^ = Е - аЛ,дП, = к^(уп),
г=1
= ЛМУп) + д(Уп)], ^кэ = к4 = + в41&1 + $42^2 + $43^), (4)
где Е — единичная матрица; дП = дд(уп)/ду; к — стадии метода; а, р^, ^4?, 1 < г < 4, 1 < ] < 3, — числовые коэффициенты, определяющие свойства точности и устойчивости (4). Для исследования схемы (4) разложим стадии А^, 1 < г < 4, в ряды Тейлора по степеням к до членов с к2 включительно. Получим
А1 = к^п, к2 = + кд„ + ак2дП + ак2дПд„ + 0(к3),
кэ = + кд„ + 2ак2дП^п + 2ак2дПдга + 0(к3), (5)
к4 = + (в41 + в42 + в43)к2^П^п + ($42 + ^^^п + 0(к3),
где <р„ = ^(Уп), дп = д(Уп), < = д^(у„)/ду, дП = дд(у„)/ду. Подставляя (5) в первую формулу (4), получим
Уп+1 = Уп + (Р1 + Р2 + Р3 + Р4 )к^п + (Р2 + Р3 )кдп +
+ ($41 + $42 + в43)Р4к2^п^п + ($42 + $43)^4^^^ +
+ а(р2 + 2р3)к2дп^п + а(Р2 + 2р3)к2дпдп + 0(к3), (6)
где элементарные дифференциалы вычислены на приближенном решении уп. Представление точного решения в окрестности точки ¿п в виде ряда Тейлора по степеням Н имеет вид
у(*п+1) = у(0 + + £) + 0.5Н2(^ + + ^ + £'£) + 0(Н3), (7)
где элементарные дифференциалы вычислены на точном решении у(£п).
Сравнивая (6) и (7) до членов с Н2 включительно при условии уп = у(£п), получим условия второго порядка точности схемы (4), т.е. р1 + р2 + р3 + р4 = 1, р2 + р3 = 1, а(р2 + 2рз) = 0.5, (в41 + в42 + в4з)р4 = 0.5, (в42 + в4з)р4 = 0.5. Отсюда сразу следует
4а — 1 1 — 2а
в41 = 0, Р2 = —-, Рз = —-, Р1 + Р4 = 0, (в42 + в43)р = 0.5. (8)
2а 2а
Исследуем устойчивость схемы (4). Применение тестового уравнения у' = Ау с комплексным Л (И,е(А) < 0) в данном случае неправомерно, поскольку теряется смысл в разделении правой части системы дифференциальных уравнений на жесткую и нежесткую части. Поэтому в (3) положим <^(у) = А1у, $(у) = А2у, где А1 и Л2 есть произвольные комплексные числа, причем И,е(Л2) < 0. Смысл А1 и Л2 — некоторые собственные числа матриц Якоби функций <^(у) и ^(у) соответственно. Применяя (4) для решения задачи
у' = А1у + А2у, у(0) = уо, I > 0, (9)
и обозначая х = НА1, г = НА2, имеем
, , х + г х + г
к1 = ху„, к2 = --уп, &3 = Т"-г^
1 — аг (1 — аг )2
, = х + (в42 + в43 — 2а)хг + (^42 + в43)х2 — а^42Х2г + а(а — ^42)хг2 к4 = (1 — аг)2 уп.
Учитывая (8), подставим полученные соотношения в первую формулу (4). Получим уп+1 = Ф(х,г)уп, где
ф(х, г) = {1 + (1 — 2а)г + х + [—2ар1 — ар2 + (в42 + в43 — 2а)р4]хг + 0.5х2 — — ав42р4х2г + [а2р1 + а2р4 — ав42р4]хг2 + (а2 — ар2)г2 }/(1 — аг )2.
Необходимым условием ¿-устойчивости численной формулы (4) относительно функции ^(у) = А2у является выполнение соотношения ф(х,г) ^ 0 при г ^ —то. Из вида ф(х,г) следует, что это требование будет выполнено, если р2 = а и в42 = 0. В результате, учитывая (8), получим набор коэффициентов схемы (4) второго порядка точности, т.е.
в41 = в42 = 0, Р2 = а, Р3 = 1 — а, Р4 = —Р1 = тг1", (10)
2Р43
где в43 — свободный параметр, а коэффициент а есть корень квадратного уравнения
а2 — 2а + 0.5 = 0. (11)
Тогда функция устойчивости ф(х, г) схемы (4) имеет вид
1 + х + 0.5х2 + (1 — 2а)г + (1 — 2а)хг
^(х,г) =-у,-та-. (12)
(1 — аг)2
Заметим, что если <^(у) = 0, то схема (4) с коэффициентами (10) и (11) совпадает с ¿-устойчивым (2,1)-методом (см., например, [3], с. 22)
Уп+1 = Уп + ак2 + (1 - а)к3, (13)
функция устойчивости ф(0,£) которого имеет вид ф(0,г) = [1 + (1 — 2а)г]/(1 — аг)2, а локальная ошибка ¿п — следующий вид:
¿п = ^ - 3) к3д'2д + 6к3д''д2 + 0(к4).
Уравнение (11) имеет два вещественных корня: а = 1 - 0.5л/2 иа= 1 + 0.5^. Выберем а =1 - 0.5\/2, потому что в этом случае меньше коэффициент в главном члене локальной ошибки численной схемы (13).
Если д(у) = 0, то формула (4) вырождается в явный двухстадийный метод типа Рун-ге — Кутты вида
Уп+1 = Уп + ( 1 — 7П1— ) + 771" к4, (14)
2в43У 2^43
функция устойчивости 0) которого имеет вид 0) = 1 + ж + 0.5ж2. Учитывая (10), нетрудно видеть, что приближение к решению уп+1 в этом случае представимо в виде
Уп+1 = Уп + к^п + 0.5к2^п + 0.25в43к3^п^п + 0(к4).
Из сравнения полученного выражения с рядом Тейлора для точного решения вытекает, что локальную ошибку ¿п схемы (14) можно записать в виде
¿п = 1 кУ V +Ц - + 0(к4).
6 " V6 4 у
Отсюда следует, что локальная ошибка численной формулы (14) будет минимальной, если в43 = 2/3. Теперь окончательно имеем набор коэффициентов схемы (4) второго порядка точности, т. е.
/2 2 3 3
а =1--—, в41 = в42 = 0, в43 = 3, Р1 = -4, Р2 = а, Р3 = 1 - а, Р4 = 4. (15)
2. Контроль точности и устойчивости
Контроль вычислений численной схемы (4) будем осуществлять с помощью метода первого порядка точности. С использованием стадий (4) можно построить семейство численных формул первого порядка вида
5
Уп+1,1 = Уп + ^ (16)
г=1
где к5 = кд(уп), 6^, 1 < г < 5, — числовые коэффициенты. Применяя разложения стадий в ряды Тейлора (5), видим, что (16) будет иметь первый порядок точности, если
61 + 62 + 63 + 64 = 1, 62 + 63 + 65 = 1.
Тогда оценку ошибки еп схемы (4) можно вычислить по формуле
£п = уп+1 — уп+1,1. (17)
При выборе коэффициентов 1 < г < 5, можно руководствоваться различными соображениями. Если, например, в функции ^(у) полностью сосредоточена жесткость задачи (3), что для многих задач (2) имеет место при В = д/(у)/ду, то имеет смысл выбрать набор коэффициентов 61 = 63 = 64 = 65 = 0, 62 = 1 или 61 = 62 = 64 = 65 = 0, 63 = 1. Такой способ оценки ошибки успешно использовался при реализации метода (13) с аналитическим вычислением матрицы Якоби. Однако если, например, в задаче (2) используется диагональная аппроксимация матрицы Якоби, то для многих задач (3) функцию ^(у) нельзя считать нежесткой частью. В такой ситуации данная оценка может приводить к потере точности вычислений из-за возникающей неустойчивости явной части численной формулы (4). Исходя из этих соображений в (16) выбраны коэффициенты 61 = 65 = 1, 62 = 63 = 64 = 0. В этом случае (16) преобразуется к виду
уп+1,1 = уп + %(уп) + #(уп)].
Как показывают расчеты, применение данной схемы в оценке (17) приводит к более надежному контролю точности вычислений. Подчеркнем важную особенность построенной оценки ошибки. В силу ¿-устойчивости из (12) следует, что для функции устойчивости ) выполняется ) ^ 0 при г ^ —то. Так как для точного решения у(^п+1) = ехр(х+г)у(£п) задачи (9) выполняется аналогичное свойство, естественным будет требование стремления к нулю оценки ошибки (17) при г ^ —то. Однако для построенной оценки имеем еп = О (г). Поэтому с целью исправления асимптотического поведения вместо еп рассмотрим оценки еп(^п) вида
£пС?п) = £п, 1 < ^п < 3. (18)
Нетрудно видеть, что оценки еп и еп(^п) как главные члены при разложении ошибок в ряды Тейлора по степеням Н совпадают при любом значении причем еп(3) ^ 0 при г ^ —то. Теперь для контроля точности вычислений можно применять следующее неравенство:
1кп(7п)|| < £п, 1 < ^п < 3. (19)
Отметим, что применение еп(^п) вместо еп не приводит к существенному увеличению вычислительных затрат. При г ^ 0 оценка еп(1) = еп правильно отражает поведение ошибки и нет смысла проверять (19) при других значениях При резком увеличении шага поведение еп может оказаться неудовлетворительным, что проявляется в неоправданном уменьшении шага и повторных вычислениях решения. Поэтому при реализации алгоритма интегрирования неравенство (19) используется следующим образом. При каждом фиксированном п выбирается наименьшее значение при котором выполняется неравенство (19). Если оно не выполняется ни при каком то шаг уменьшается и решение вычисляется повторно.
Теперь перейдем к построению неравенства для контроля устойчивости явной части схемы (4). Заметим, что переход от задачи (1) к (2) в общем случае не гарантирует, что функция ^(у) = / (у) — Ву есть нежесткая часть системы обыкновенных дифференциальных уравнений. В некоторых случаях это не оказывает сильного влияния на эффективность алгоритма интегрирования, поскольку формула (4) обладает хорошими свойствами
устойчивости. Однако для многих задач неравенство для контроля устойчивости явной части численной схемы (4) может существенно повысить эффективность расчетов. Отметим также, что в отличие от явных численных схем, где на жестких задачах контроль устойчивости всегда дает положительный эффект [3], в случае схемы (4) столь обременительное ограничение на выбор шага может приводить и к существенному понижению быстродействия алгоритма интегрирования, вплоть до катастрофически низкого. Поэтому право на включение или выключение данного неравенства следует передать пользователю через посредство формального параметра подпрограммы, реализующей алгоритм интегрирования.
Неравенство для контроля устойчивости построим по аналогии с [3]. Однако здесь, в отличие от [3], нельзя воспользоваться ранее вычисленными стадиями k^, 1 < i < 4, метода (4) в силу специфики задачи (3). Поэтому поступим следующим образом. Введем в рассмотрение стадии di и d2 вида
di = h^(y„ + a2i ki), d2 = h^(y„ + a3iki + o^di).
Пусть ^(y) = Ay + b, где A и b — матрица и вектор с постоянными коэффициентами соответственно. Тогда имеют место соотношения
ki = h(Ay„ + b), di = ki + a2ihAki, d2 = ki + (a3i + a32)hAki + a2i a32h2A2ki. Полагая a2i = a3i + a32, получим
d2 — di = a2ia32h2A2ki, di — ki = a2ihAki.
Тогда с использованием степенного метода максимальное собственное число V^ = h|An max | матрицы hA можно оценить по следующей приближенной формуле:
лг 1 |d2» — dii| Vn = -г max
|«32| Ин - Ы'
В общем случае полученная оценка является грубой, потому что в степенном методе используется мало итераций и плюс дополнительные искажения вносит нелинейность функции <^(у). Поэтому для контроля устойчивости будем применять неравенство как некоторый ограничитель на рост величины шага интегрирования.
Приведем формулу для выбора шага по точности и устойчивости. Пусть приближение к решению уп в точке ¿п вычислено с шагом кп. Для оценки ошибки еп(^п) имеет место соотношение еп(^п) = О(кп). Прогнозируемый шаг касс по точности вычислим по формуле касс = 91 кп, где есть корень уравнения
З^Ы^Ц = е.
Учитывая, что ^Л = 0(кп), шаг к.^ по устойчивости определим по формуле к.^ = д2кп, где д2 есть корень уравнения = 2, а числом 2 ограничен интервал устойчивости явного двухстадийного метода типа Рунге — Кутты второго порядка точности. Тогда прогнозируемый шаг кп+1 по точности с ограничением по устойчивости можно вычислить по формуле
кп+1 = тах[кп, тт(к„сс, к^)].
Оценка максимального собственного числа матрицы Якоби приводит к двум дополнительным вычислениям функции <^(у). Если есть уверенность, что вся жесткость сосредоточена в функции д(у) или <^(у) = 0, то с целью экономии вычислительных затрат данная оценка может быть отключена.
3. Численная схема для неавтономных задач
Рассмотрим задачу Коши для неавтономной системы обыкновенных дифференциальных уравнений следующего вида:
у' = р(*,у)+ д(г,у), у(*о) = уо, ¿0 < г < 4, (20)
где у, ^ и д — гладкие вещественные N-мерные вектор-функции; г — независимая переменная. Далее снова будем предполагать, что вся жесткость сосредоточена в функции д(г, у), а у) есть нежесткая часть. Для численного решения (20) рассмотрим четырехстадийный метод вида
4
уп+1 = уп + ^ РЛ, Б = Е — аНдп, ^ = Н^п, уп),
г=1
= Н[^(^п,уп)+ д(^п + сН,уп)], = &2, (21)
^4 = Н^(£п + (в41 + в42 + в43)Н, уп + ^41^1 + ^42^2 + ^43^),
где Е — единичная матрица; дп = дд(гп, уп)/ду; к — стадии метода; а, Д^?, 1 < г < 4, 1 < ^ < 3, — числовые коэффициенты. Для исследования схемы (21) разложим стадии кг, 1 < г < 4, в ряды Тейлора по степеням Н до членов с Н2 включительно. Получим
к1 = Н^п,
к2 = Н^п + Ндп + сН2д^п + аН2дП^п + аН2дпдп + О(Н3),
к3 = Н^п + Ндп + сН2д£п + 2аН2дП^п + 2аН2 дпдп + О(Н3),
^4 = Н^п + (в41 + в42 + в43)Н2^[п + (в41 + Д42 + ^43)Н2 ^П^п +
+ (в42 + в43)Н2^пдп + О(Н3),
(22)
где ^п = ^(гп,yn), дп = дС^Уп^ ^ = д^(гп,уп)/дг, д!п = ^С^Уп^^ ^П = ^С^Уп^^
дп = дд(гп,уп)/ду. Подставляя (22) в первую формулу (21), запишем Уп+1 = Уп + (Р1 + Р2 + Р3 + Р4 )Н^п + (^2 + Р3 )Ндп +
+ (в41 + ^42 + в43)Р4Н2^п + с(Р2 + Р^Кп + + (^41 + в42 + ^43)Р4Н2^П^П + (в42 + ^^Н^дп +
+ а(р2 + 2р>3)Н2дП^п + а(р + 2р3)Н2дПдп + О(Н3),
(23)
где элементарные дифференциалы вычислены на приближенном решении уП. Представление точного решения у(гп+1) в окрестности точки ¿П в виде ряда Тейлора по степеням Н имеет вид
у(гп+1) = у(гп) + Н(^ + д) + 0.5Н2(^ + д! + ^ + р'д + д> + д'д) + О(Н3), (24)
где элементарные дифференциалы вычислены на точном решении у(гп).
Сравнивая (23) и (24) до членов с Н2 включительно при условии уП = у(гп), получим условия второго порядка точности схемы (21), т.е. р1 + р2 + р3 + р4 = 1, р2 + р3 = 1, (в41 + в42 + в43)р4 = 0.5, с(р2 + Р3) = 0.5, (^42 + ^43)^4 = 0.5, а(р + 2р?) = 0.5. Отсюда сразу следует, что с = 0.5. Теперь, рассуждая по аналогии с исследованием схемы (4), получим
коэффициенты численной формулы (21), которые имеют вид
/2 2
а = 1 - —, в41 = в42 = 0, С = 0.5, в43 = 3,
33 Р1 = -4, Р2 = а, Р3 = 1 - а, Р4 = 4. (25)
Для контроля точности вычислений будем применять неравенство (18), где в оценке еп используется приближение к решению, полученное методом второго порядка точности (21), и приближенное решение, вычисленное методом первого порядка следующего вида:
Уп+1,1 = Уп + к[^(£п,Уп) + д(^п + 0.5к,уп)].
Контроль устойчивости и выбор шага по точности и устойчивости осуществляются точно так же, как и в случае автономной системы.
4. Анализ результатов расчетов
Далее построенный алгоритм будем называть АБОВЕ. Замораживание матрицы Якоби, т. е. применение матрицы О = Е - акдп на нескольких шагах интегрирования, проводилось по следующему правилу. Если матрица Якоби не пересчитывалась, то для сохранения устойчивости численной схемы величина шага интегрирования тоже не менялась. Попытка применения старой матрицы осуществлялась после каждого успешного шага интегрирования. Три причины приводили к размораживанию: 1) нарушение неравенства для контроля точности вычислений; 2) превышение количества шагов с замороженной матрицей числа /; 3) превышение прогнозируемого шага интегрирования последнего успешного шага в ^ раз. Параметры / и ^ можно использовать для настройки метода на решение конкретных задач. Если 5/ ^ то и ^ ^ то, то число шагов с одной матрицей Якоби возрастает. Если 5/ = 0 и ^ = 0, то замораживания матрицы не происходит. Поэтому в случае большой размерности системы обыкновенных дифференциальных уравнений имеет смысл выбирать 5/ и ^ достаточно большими. В приведенных ниже расчетах / = 20 и ^ = 2.
Все рассматриваемые ниже примеры приводились к виду (2). Расчеты выполнялись с
требуемой точностью е = 10-2, на 1ВМ РС АШоп(1;т)ХР 2000+--с двойной точностью.
Схема (4) имеет второй порядок точности, и поэтому проводить с ее помощью расчеты с более высокой точностью нецелесообразно. В конкретных расчетах левая часть неравенства (19) вычислялась по формуле
||епЫН = ЙТТ1, (26)
|уп1 + г
где г — номер компоненты; г — положительный параметр. Если по г-й компоненте решения выполняется неравенство |уп| < г, то контролируется абсолютная ошибка ег, в противном случае — относительная ошибка е. В расчетах параметр Г выбирался таким образом, чтобы по всем компонентам решения фактическая точность была не хуже задаваемой.
Ниже через idec, 1эо1 обозначены соответственно суммарное число шагов инте-
грирования, правых частей системы (1), декомпозиций матрицы Якоби и число обратных ходов в методе Гаусса.
Пример 1 [4].
у' = — 0.013у1 — 1000у1у3,
у2 = — 2500У2Уз,
у3 = — 0.013у1 — 1000у1уз — 2500у2уз,
г е [0,50], У1(0) = 1, У2(0) = 1, Уз(0) = 0, Но = 2.9 ■ 10-4.
(27)
Решение задачи (27) осуществлялось методом (4) с контролем устойчивости и диагональной аппроксимацией матрицы Якоби, т. е. в численной формуле (4) использовалась диагональная матрица Б с элементами йц на диагонали вида йц = —0.013 — 1000у3, й22 = —2500у3, й33 = — 1000у1 — 2500у2. Так как в этом случае вычислительные затраты метода (4) практически такие же, как и в явных методах, сравнение эффективности проводилось с известным явным методом Мерсона [5] четвертого порядка точности. Для вычисления приближенного решения с точностью е = 10-2 алгоритмом АБОВЕ потребовалось 687 шагов (остальные затраты вычисляются из вида схемы (5)). Для решения данной задачи методу Мерсона без контроля устойчивости [5] потребовалось 400 627 вычислений правой части, а с контролем устойчивости (см. [3, с. 81]) — = 311 715.
В случае использования полной матрицы Якоби системы (27) алгоритму АБОВЕ без замораживания матрицы Якоби для решения задачи (27) потребовалось 38 шагов, 38 декомпозиций матрицы Якоби и 108 обратных ходов в методе Гаусса (остальные затраты вычисляются из вида схемы (5)). При расчетах с замораживанием матрицы Якоби вычислительные затраты следующие: = 98, idec = 15, 1бо1 = 288.
Пример 2 [4].
Решение задачи (28) осуществлялось методом (4) с контролем устойчивости и диагональной аппроксимацией матрицы Якоби, причем й11 = —55 — у3, й22 = —0.07 85у3, й33 = 0. Приближенное решение с точностью е = 10-2 алгоритмом АБОВЕ вычислено за 4953 шага. Для решения данной задачи методу Мерсона без контроля устойчивости [5] потребовалось 80 713 вычислений правой части, а с контролем устойчивости [3] — if = 62 299.
В случае использования полной матрицы Якоби системы (28) алгоритму АБОВЕ без замораживания матрицы Якоби потребовались 81 шаг, 81 декомпозиция матрицы Якоби и 388 обратных ходов в методе Гаусса. При расчетах с замораживанием матрицы вычислительные затраты следующие: is = 338, idec = 24, isol = 1124.
у1 = — 55У1+ 65У2— у^
У2 = 0.0785(у1 — у2),
У3 = 0^Уъ
г е [0, 500], У1 (0) = 1, У2(0) = 1, Уз(0) = 0, Но = 2 ■ 10-2.
(28)
Пример 3 ([1, с. 167]).
у1 = у2,
у2 = ((1 — у2)у2 — уОЯ * = ю-6,
г е [0,11], У1(0) = 2, У2(0) = 0, Но = 10-6.
(29)
Решение задачи (29) осуществлялось методом (4) с контролем устойчивости и полной матрицей Якоби. При расчетах без замораживания матрицы Якоби is = 7533, idec = 8879, isol = 28182, в случае замораживания матрицы Якоби затраты следующие: is = 36 300,
idec = 5660, iso1 = 115 324. Аппроксимация матрицы Якоби диагональю приводит примерно к двум миллионам шагов, т. е. для решения такого типа задач применение диагонали практически неприемлемо. Однако отметим, что затраты метода Мерсона примерно в 10 раз больше, и поэтому в случае необходимости применения явных численных схем построенный алгоритм эффективнее.
Пример 4 ([1, с. 167]).
у1 = -0.04у1 + 104У2У3,
У2 = 0.04у1 - 104У2У3 - 3 ■ 107у2, (30)
У3 = 3 ■ 107у2,
г е [0,1011 ], У1 (0) = 1, У2(0) = 0, У3(0) = 0, ко = 10-6.
Решение задачи (30) осуществлялось методом (4) без контроля устойчивости и с полной матрицей Якоби. При расчетах без замораживания матрицы Якоби is = 7182, idec = 7183, iso1 = 17487, в случае замораживания матрицы Якоби затраты следующие: is = 45 303, idec = 2336, iso1 = 134684. Заметим, что на данной задаче контроль устойчивости приводит к слишком обременительному ограничению на размер шага. Если задачу рассматривать на небольшом интервале, то контроль устойчивости приводит к положительному результату. Однако когда шаг порядка 103, рост шага ограничивается неравенством для контроля устойчивости. Отметим также, что свойства устойчивости метода (4) позволяют проводить расчеты со значительно большим шагом. В частности, при расчетах без контроля устойчивости с полной матрицей Якоби величина шага в конце интервала интегрирования порядка 108.
Аппроксимация матрицы Якоби диагональю тоже приводит к значительным вычислительным затратам. Дело в том, что в задаче (30) как только компонента численного решения у2 случайно становится отрицательной, она стремится к -то ([1, а 167]). Для проверки качества работы численных методов именно из этого соображения выбран такой большой интервал интегрирования. В случае диагональной аппроксимации матрицы Яко-би, которая грубо показывает направление изменения решения, параметр г в формуле (26) приходится выбирать порядка 10-14, т.е. порядка у2 в конце интервала интегрирования. Это естественно приводит к большим вычислительным затратам.
Для полноты картины приведем еще результаты расчетов задачи (30) методом (4) при <^(у) = 0, т. е. когда в схеме (4) отключена явная часть. В этом случае для алгоритма АБОБЕ без замораживания матрицы Якоби is = 2686, idec = 2689, iso1 = 5950, а для алгоритма с замораживанием — is = 12 061, idec = 641, iso1 = 35 531.
Пример 5 ([1, с. 168]).
у1 = 77.27(у1 (1 - 8.375 ■ 10-6у1 - У2) + У2),
У2 = (У3 - (1+ У1)у2)/77.27, (31)
у3 = °.161(У1 - У3)
г е [0,360], у1(0) = 1, У2(0) = 2, у3(0) = 3, ко = 10-6.
Решение задачи (31) осуществлялось методом (4) с диагональной аппроксимацией матрицы Якоби, причем ^п = 77.27(1 -1.675- 10-7у1 -у2), ^22 = -(1+у1)/77.27, ^33 = -0.161. Приближенное решение с точностью е = 10-2 алгоритмом АБОБЕ с контролем устойчивости вычислено за 27 385 шагов, а без контроля устойчивости — за 19 964 шага. Для решения
данной задачи методу Мерсона без контроля устойчивости [5] потребовалось 23 700 664 вычислений правой части, а с контролем устойчивости [3] — if = 19 037830.
В случае использования полной матрицы Якоби системы (31) алгоритму ASODE без замораживания матрицы Якоби потребовалось 2449 шагов, 2652 декомпозиции матрицы Якоби и 6964 обратных хода в методе Гаусса. При расчетах с замораживанием матрицы вычислительные затраты следующие: is = 19 807, idec = 3431, isol = 50 924.
При решении задачи (31) методом (4) при <^(у) = 0, т.е. когда в схеме (4) отключена явная часть, для алгоритма ASODE без замораживания матрицы Якоби is = 953, idec = 1112, isol = 3054, а для алгоритма с замораживанием — is = 2895, idec = 347, isol = 7944.
Пример 6 ([1, с. 168]).
yi = -1.71yi + 0.43y2 + 8.32уз + 0.0007,
у2 = +1.71yi - 8.75y2,
y3 = -10.03уз + 0.43y4 + 0.035y5,
у4 = +8.32y2 + 1.71y3 - 1.12y4,
y5 = -1.745y5 + 0.43y6 + 0.43у7, (32)
y6 = -280y6ys + 0.69y4 + 1.71y5 - 0.43y6 + 0.69у7,
у7 = +280y6y8 - 1.81y7,
у8 = ^
t G [0, 421.8122], yi (0) = 1, у2(0) = ... = У7(0) = 0, уз(0) = 0.0057, ho = 10-6.
Решение задачи (32) осуществлялось методом (4) с диагональной аппроксимацией матрицы Якоби, причем dii = -1.71, d22 = -8.75, d33 = -10.03, d44 = -1.12, d55 = -1.745, d66 = -280y8 - 0.43, d77 = -1.81, d88 = -280y6. Приближенное решение с точностью е = 10-2 алгоритмом ASODE с контролем устойчивости вычислено за 5167 шагов, а без контроля устойчивости — за 7079 шагов. Для решения данной задачи методу Мерсона без контроля устойчивости [5] потребовалось 73 308 вычислений правой части, а с контролем устойчивости [3] — if = 57912.
В случае использования полной матрицы Якоби системы (32) алгоритму ASODE без замораживания матрицы Якоби потребовалось is = 108, idec = 109, isol = 276 без контроля устойчивости и is = 273, idec = 273, isol = 547 с контролем устойчивости. При расчетах с замораживанием матрицы вычислительные затраты следующие: is = 481, idec = 91, isol = 1228 без контроля устойчивости и is = 350, idec = 35, isol = 966 с контролем устойчивости.
При решении задачи (32) методом (4) при <^(у) = 0, т.е. когда в схеме (4) отключена явная часть, для алгоритма ASODE без замораживания матрицы Якоби is = 71, idec = 72, isol = 159, а для алгоритма с замораживанием — is = 151, idec = 22, isol = 341.
Пример 7 ([1, с. 168]).
yi = -Ayi - ву1у3,
у2 = ауЪ
у3 = АУ1 - ВУ1У3 - МСУ2У3 + СУ4, (33)
у4 = ВУ1У3 - СУ4,
t G [0, 360], yi(0) = 1.76 ■ 10-3, У2(0) = У3(0)= У4(0) = 0, ho = 10-6,
где A = 7.89 ■ 10-i0, B =1.1 ■ 107, C = 1.13 ■ 103, M = 106. В данной задаче масштабы переменных имеют очень большой разброс, что приводит к определенным проблемам с
приближенным решением. Переменные у2, у3 и у4 вырастают до величины порядка 10-10, а затем у4 убывает практически до машинного нуля, а у1 и у2 в конце интервала равны между собой, причем у1 = у2 « 8.8612 ■ 10-23. Все это происходит на фоне достаточно большой переменной у1, которая в начале промежутка интегрирования порядка 10-3, а в конце стремится к нулю. Ясно, что в такой ситуации на эффективность алгоритма интегрирования существенное влияние оказывает выбор параметра г в формуле (26). В этом случае значение г имеет смысл задавать по каждой компоненте отдельно. Однако здесь г выбирался одинаковым по всем переменным по следующему правилу. На первом шаге г = 10-3, т.е. порядка у1. Учитывая, что для (33) имеет место закон сохранения у2 - у3 - у4 = 0, на последующих шагах в случае нарушения неравенства |у2 - у3 - у4| < е значение г полагалось равным тт(|у2|, |у3|), где е есть требуемая точность интегрирования. Заметим, что до некоторого момента все три переменные у2, у3 и у4 примерно одного порядка, а затем компонента у4 достаточно быстро убывает к нулю. По этой причине она в выборе параметра г не участвует.
Численное решение (33) осуществлялось методом (4) при ^(у) = 0. Две правильные значащие цифры в решении достигаются при задаваемой точности е = 10-4, при этом is = 26 410, idec = 26 411, iso1 = 65 421 для алгоритма без замораживания матрицы Якоби и is = 41 240, idec = 2073, iso1 = 115 208 для алгоритма с замораживанием. В случае приведения задачи (33) к виду (2) требуются порядка полумиллиона шагов, причем затраты не зависят от того, контролируется устойчивость или нет.
Заключение
Предложенный алгоритм интегрирования разрабатывался по инициативе чл.-корр. РАН В.В. Шайдурова для численного решения задач механики сплошной среды после дискретизации по пространству методом конечных элементов или с помощью конечных разностей. В этом случае в задаче (3) разделение на функции д(у) и ^(у) происходит естественным образом: д(у) есть симметричная часть, описываемая оператором дифференцирования второго порядка, а ^(у) есть несимметричная часть (конвективные слагаемые), описываемая оператором дифференцирования первого порядка. При реализации численной формулы (4) необходимо дважды решать линейную систему алгебраических уравнений. В задачах механики сплошной среды эффективность алгоритма интегрирования может быть достигнута за счет специальных методов решения линейных систем с симметричной матрицей, которая во многих случаях положительно определенная.
Схему (4) можно применять также для решения локально-неустойчивых задач, причем ^(у) в этом случае отвечает за собственные числа матрицы Якоби с положительными вещественными частями. В отличие от А-устойчивых или ¿-устойчивых методов, у которых область неустойчивости обычно небольшая и которые являются А-устойчивыми или ¿-устойчивыми не только в левой, но и в правой полуплоскости плоскости (кЛ), явные методы типа Рунге — Кутты являются неустойчивыми практически во всей правой полуплоскости и поэтому более предпочтительны при определении неустойчивого решения. Для локально-неустойчивых задач в ряде случаев разделение правой части системы обыкновенных дифференциальных уравнений на функции ^(у) и д(у) из физических соображений также не вызывает особых трудностей.
Приведенные результаты численных экспериментов не ориентированы на решение задач механики сплошной среды или локально-неустойчивых задач, а направлены на иссле-
дование возможностей алгоритма интегрирования при решении некоторых общепринятых тестовых примеров. Тестовые примеры выбирались таким образом, чтобы продемонстрировать разные стороны работы алгоритма интегрирования. Если поведение алгоритма на нескольких тестовых задачах совпадало, то выбирался наиболее простой из них. Цель расчетов заключалась в проверке работоспособности алгоритма с переменным шагом и с замораживанием матрицы Якоби, надежности неравенства для контроля точности вычислений, а также в исследовании возможности расчетов с диагональной аппроксимацией матрицы Якоби. В последнем случае вычислительные затраты на шаг в явных методах и построенном алгоритме различаются незначительно. В частности, из анализа результатов расчетов жестких задач следует, что в случае невозможности применения методов с неограниченной областью устойчивости алгоритм (4) существенно эффективнее метода Мерсона — наиболее распространенного среди явных численных схем типа Рунге — Кутты.
Список литературы
[1] Хайрер Э., Ваннер Г. Решение обыкновенных дифференциальных уравнений. Жесткие и дифференциально-алгебраические задачи: Пер. с англ. М.: Мир, 1999. 685 с.
[2] Rosenbrook H.H. General implicit processes for the numerical solution of differential equations // Computer J. 1963. N 5. P. 329-330.
[3] Новиков Е.А. Явные методы для жестких систем. Новосибирск: Наука, 1997. 195 с.
[4] Gear C.W. The automatic integration of stiff ordinary differential equations // Proc. IFIP Congress. 1968. P. 81-85.
[5] Merson R.H. An operational methods for integration processes // Proc. Symp. on Data Proc. Weapons Research Establishment, Salisbury, Australia, 1957.
[6] Kaps P., Rentrop P. Generalized Runge — Kutta methods of order four with stepsize control for stiff ordinary differential equations // Numer. Math. 1979. N 33. P. 55-68.
[7] Новиков Е.А., Шитов Ю.А., Шокин Ю.И. Одношаговые безытерационные методы решения жестких систем // Докл. АН СССР. 1988. Т. 301, № 6. С. 1310-1314.
[8] Cooper G.J., Sayfy A. Additive Runge — Kutta methods for stiff ordinary differential equations // Mathematics of Computation. 1983. Vol. 40, N 161. P. 207-218.
Поступила в редакцию 31 мая 2004 г-
