CC BY
45
УДК 519.622
Е. А. Новиков
АЛГОРИТМ ИНТЕГРИРОВАНИЯ ПЕРЕМЕННОЙ СТРУКТУРЫ ДЛЯ РЕШЕНИЯ ЖЕСТКИХ ЗАДАЧ НА ОСНОВЕ ЯВНОГО И ¿-УСТОЙЧИВОГО МЕТОДОВ1
Построены явная двухстадийная схема типа Рунге-Кутта первого порядка с расширенной областью устойчивости и L-устойчивый (т, к)-метод второго порядка точности. Разработан алгоритм переменного порядка и шага, в котором выбор наиболее эффективной численной схемы осуществляется на каждом шаге с применением неравенства для контроля устойчивости. Приведены результаты расчетов, подтверждающие эффективность алгоритма при расчетах с точностью 1% и ниже.
Современные методы решения жестких задач, как правило, используют вычисление и обращение матрицы Якоби системы дифференциальных уравнений. В случае достаточно большой размерности системы эффективность этих методов практически полностью определяется временем обращения этой матрицы. Для повышения эффективности расчетов одна матрица обычно используется на нескольких шагах интегрирования [1]. Наиболее успешен этот подход в алгоритмах на основе многошаговых методов [2]. Эта проблема не вызывает особых трудностей при построении алгоритмов интегрирования на основе других численных схем, если вычисление стадий в них проводится с участием матрицы Якоби в итерационном процессе. В этом случае матрица Якоби не влияет на порядок точности метода, а определяет только сходимость итераций. Необходимость в ее пересчете возникает лишь при значительном замедлении скорости сходимости итерационного процесса.
В известных безытерационных численных схемах, к которым относятся, например, методы типа Розенбро-ка [3], замораживание матрицы Якоби является более сложной задачей. Следует отметить, что безытерацион-ные методы более просты с точки зрения реализации на ЭВМ, чем алгоритмы с применением итераций, однако в методах вида [3] матрица Якоби влияет на порядок точности численной схемы и поэтому возникают трудности с ее замораживанием. Если вопрос об использовании одной матрицы на нескольких шагах интегрирования оставить не решенным, то нужно заведомо ограничиться задачами небольшой размерности. В [4; 5] эта проблема изучалась применительно к методам типа Розенброка и было доказано, что максимальный порядок точности данных методов с замораживанием матрицы Якоби равен двум.
Некоторым аналогом замораживания матрицы Якоби является применение в расчетах алгоритмов на основе явных и L-устойчивых методов с автоматическим выбором численной схемы. В этом случае эффективность может быть повышена за счет расчета переходного участка, соответствующего максимальному собственному числу матрицы Якоби, явным методом. В работе [6] в качестве критерия выбора эффективной численной формулы применялось неравенство для контроля устойчивости и на основе явного метода типа Рунге-Кутта и
L-устойчивого (2, 1)-метода был построен алгоритм переменной структуры, в котором допускалось замораживание как численной, так и аналитической матрицы Якоби.
Для численного решения задачи Коши для жестких систем обыкновенных дифференциальных уравнений
У = /(У), У(*о) = У» *о ^ ^ *к, (1)
где у и f - вещественные ^-мерные векторы-функции; I - независимая переменная, был предложен класс (т, к)-методов [7]. С точки зрения реализации (т, к)-ме-тоды столь же просты, как и схемы [3], однако в данном классе более просто решается проблема замораживания матрицы Якоби, а также ее численной аппроксимации. Кроме того, (т, к)-методы обладают лучшими свойствами точности и устойчивости при незначительном увеличении вычислительных затрат.
Для решения (1) рассмотрим (2, 1)-схему вида
Уп+1 = Уп+РА + Р2К АА = V(Уп I АА = (2)
где к1 и к2 - стадии метода; Dn = Е - акАп; Е - единичная матрица; h - шаг интегрирования; Ап - некоторая матрица, имеющая вид Ап = fn + hBn + О(к2); fn = д/(уп) / ду -матрица Якоби задачи (1); Вп - не зависящая от шага интегрирования произвольная матрица; а, р1 и р2 - числовые коэффициенты. Использование матрицы Ап, представленной в виде А =/ + кВ + О(к2), позволяет применять (2) с замораживанием как аналитической, так и численной матрицы Якоби [8]. В случае применения матрицы Якоби / п к, вычисленной к шагов назад, имеем Вп = -к/X / Г/п = д 2/(у) / ду/(уп). Если матрица Якоби вычисляется численно с шагом г = ск, то элементы Ь ..
1 1 ’ п, ч
матрицы Вп имеют вид Ьп ч = 0,5с. 2/(уп) / ду2.
Получим коэффициенты L-устойчивой численной схемы (2) второго порядка и неравенство для контроля точности вычислений. Разложение точного решения у(гп + 1) в ряд Тейлора в окрестности точки 1п до членов с к3 имеет вид
У(*п+1) = У(*п) + / + 0,5к2// + ^ (/Л/ + //2) + 0(к4), (3)
6
где элементарные дифференциалы/,///2/и/'/2 вычислены на точном решенииу(г). Для нахождения коэффициентов а,р ир2 схемы (2) запишем разложения стадий к1 и к2 в ряды Тейлора в окрестности точки уп до членов с к3 включительно и подставим их в (2):
1 Работа поддержана грантами РФФИ № 08-01-00621 и Президента НШ-3428.2006.9.
Уп+1 = Уп + (Р1 + Р 2 )Н/п + а(Р1 + 2 Р 2 )Л 2 ГЛ. +
+ а2( Р1 + 3 Р2)Л3 /'^/п + а( Р1 + 2 Р2)^3 вп/п + 0(Л4Х ^ где элементарные дифференциалы вычислены на приближенном решенииуп. Полагаяуп = у(£п) и сравнивая (3) и (4) до членов с к2 включительно, получим условия второго порядка точности схемы (2):
Р1 + Р2 = 1, аР1 + 2ар2 = 0,5. (5)
Исследуем устойчивость численной формулы (2). Применяя ее к задаче у' = лу, Яв(Х) < 0, получим уп+1 = Q(x)yn, х = кХ, где функция устойчивости Q(x) имеет вид Q(x) = = [1 + (р1 + Р2 -2 а)х + а(а - _р1)х2] / (1 - ах)2. Схема (2) будет L-устойчивой, еслир1 = а. Подставляя это соотношение в (5), получим набор коэффициентов
Р1 = а, р2 = 1 - а, (6)
где а определяется по условию L-устойчивости:
а2 -2а + 0,5 = 0. (7)
Сравнивая (3) и (4) до членов с к3 включительно, получим, что локальная ошибка 5 численной схемы (2) с коэффициентами (6) будет
5п = Л3 [(а - V/ +1 //2 -1 Вп/] + 0(И4). (8)
3 6 2
Уравнение (7) имеет два корня: а1 = 1 - 0, ^л/2 и а2 = 1 +0, 5л/2. Выберем а = а1, так как в этом случае коэффициент в главном члене (а -1/3) к3/г2/'меньше ошибки (8).
Контроль точности вычислений построим по аналогии с [9]. Для этого введем обозначение
К Л) = А?' (к2 - *1), (9)
где к1 и к2 вычисляются по формулам (2). Тогда, согласно [9], для контроля точности вычислений на каждом шаге нужно проверять неравенство
||уОп ) || < е, 1 < ]п < 2, (10)
где £ - требуемая точность расчетов; ||-|| - некоторая норма в Я"; целочисленная переменная 1 выбирается наименьшей для выполнения неравенства (10).
Отметим важную особенность построенной оценки ошибки (9). Схема (2) L-устойчива, т. е. для ее функции устойчивости Q(x) при х ^ -го имеет место соотношение Q(x) ^ 0. Так как для точного решения у(^п+1) = ехр(х)у(£п) задачи у' = Ху, Яв(Х) < 0, выполняется аналогичное свойство, то естественным будет требование стремления к нулю оценки ошибки при х ^ -го. Однако для у(1) это свойство не выполняется: данная оценка ведет себя А-устойчивым образом. Для исправления асимптотического поведения ошибки вместо у(1) введена оценка v(/n), 1 < 1п < 2. В этом случае поведение оценки ошибки при 1п = 2 будет согласовано с поведением точного решения тестовой задачи при х ^ -го. Подчеркнем, что в смысле главного члена оценки v(1) и v(2) совпадают и использование v(1) фактически не приводит к увеличению вычислительных затрат. Это связано с тем, что v(jn) при1п = 2 проверяется только в том случае, если оно нарушено при
1 = 1. Не такая ситуация встречается достаточно редко, в основном при быстром росте величины шага интегрирования.
Оценку максимального собственного числа „ = кХ
' п,0 п, тах
матрицы Якоби системы (1), необходимую для перехода на явную формулу, оценим через ее норму тмп0 = к ||Э/0у||.
Данная оценка будет применяться ниже для автоматического выбора численной схемы.
Для численного решения задачи (1) рассмотрим схему вида [10]
Уп+1 = Уп + Г1к1 + r2k2, к1 = Л/ (Уп ), к2 = Л/ (Уп + К)- (11)
Отметим, что при г1 = г2 = 0,5 численная формула (11) имеет второй порядок точности. Построим менее точную схему с максимальным интервалом устойчивости. Для этого применим схему (11) для решения скалярного тестового уравнения у'= Ху, Яв(Х) < 0. Получим уп+1 = = Q(x)y , где функция устойчивости Q(x) имеет вид Q(x) = 1 + (г1 + г2)х + г2 х2, х = кХ. Требование первого порядка точности приводит к соотношению г1 + г2 = 1, которое ниже будем считать выполненным. Теперь выберем г2 таким образом, чтобы метод (11) имел максимальный интервал устойчивости. Для этого рассмотрим многочлен Чебышева Т2^) = ^2 - 1) на промежутке [-1, 1]. Проведем замену переменных, полагая z = 1 - 2x1 у Получим Т2(х) = 1 - 8х / у + 8х2 / у2, при этом отрезок [у, 0] отображается на отрезок [-1, 1]. Нетрудно показать, что среди всех многочленов вида Р2(х) = 1 + х + с2 х2 для Т2(х) неравенство |Т2(х)| < 1 выполняется на максимальном интервале [у, 0], у = -8. Потребуем совпадения коэффициентов Q(x) и Т2(х) при у= -8. Это приведет к соотношениям г1 + г2 = 1 и г2 = 1/8. В результате имеем коэффициенты г1 = 7/8 и г2 = 1/8 метода первого порядка точности с максимальным интервалом устойчивости, локальная ошибка 5п которого имеет вид
Дп = 3 Л2// + 0(къ).
Для контроля точности численной формулы первого порядка будем использовать оценку локальной ошибки. Учитывая, что к2 - к1 = к2//п + О(к3) и вид локальной ошибки, неравенство для контроля точности запишем в виде ||к2 - к1|| < 8£ / 3, где ||-|| - некоторая норма в Я"; £ - требуемая точность расчетов.
Построим неравенство для контроля устойчивости метода первого порядка. Для этого рассмотрим вспомогательную стадию к3 = к/у ). Заметим, что к3 совпадает со стадией к1, которая применяется на следующем шаге интегрирования и поэтому ее использование не приводит к дополнительным вычислениям правой части системы (1). Запишем к1, к2и к3 для задачи у' = Ау, где А - матрица с постоянными коэффициентами:
к = ХУп, к2 = (X + X 2)Уп, кз = (X + X2 + 0,125Х3) Уп, где X = кА. Легко видеть, что
к2 - к1 = X 2Уп , 8(к3 - к2 ) = X3Уп .
Тогда, согласно [6], оценку максимального собственного числа w = кХ матрицы Якоби системы (1) можно вычислить по формуле
^п1 = 8тах(|к\ -к\ |/|к[ -к'1 |). (12)
1</<Л^
Интервал устойчивости схемы (11) первого порядка точности равен восьми. Поэтому для контроля ее устойчивости можно применять неравенство ^п2 < 8. В случае использования данного неравенства для выбора шага следует учитывать грубость оценки (12), потому что вовсе не обязательно, чтобы максимальное собственное
число было сильно отделено от остальных. Кроме того, в степенном методе применяется мало итераций и дополнительные искажения вносит нелинейность задачи (1). Поэтому контроль устойчивости используется как ограничитель на размер шага интегрирования. В результате прогнозируемый шаг кп+1 будем вычислять следующим образом. Новый шаг кас по точности определим по формуле кас = дкп, где кп - последний успешный шаг интегрирования, а д, учитывая соотношение к2 - к1 = 0(кп2), задается уравнением д2||к2 - к1| | = 8е / 3. Шаг к“1 по устойчивости зададим формулой к= dк , где d определяется по равенству dwn2 = 2. Тогда прогнозируемый шаг кп+1 вычисляется по формуле
кп+1 = тах[кп ,тт(кас, к51)]. (13)
Заметим, что формула (13) применяется для прогноза величины шага интегрирования кп+1 после успешного вычисления решения с предыдущим шагом кп и поэтому фактически не приводит к увеличению вычислительных затрат. Если шаг по устойчивости меньше последнего успешного, то он уменьшен не будет, поскольку причиной этого может быть грубость оценки максимального собственного числа. Однако шаг не будет и увеличен, так как не исключена возможность неустойчивости численной схемы. Если шаг по устойчивости должен быть уменьшен, то в качестве следующего шага будет применяться последний успешный шаг кп, в связи с чем для выбора шага и предлагается формула (13). Данная формула позволяет стабилизировать поведение шага на участке установления решения, где определяющую роль играет устойчивость.
На основе ¿-устойчивого метода (2) и численной формулы (11) первого порядка точности легко сформулировать алгоритм интегрирования с автоматическим выбором численной схемы. При расчетах по явному методу нарушение неравенства < 8 вызывает переход на схе-
му (2). А при вычислениях по ¿-устойчивой схеме передача управления явному методу происходит в случае выполнения неравенства ^п 0 < 8. Численную формулу (2) без потери порядка точности можно применять с замораживанием матрицы Якоби. Попытка замораживания матрицы осуществляется после каждого успешного шага. Размораживание матрицы происходит в следующих случаях:
- если нарушена точность расчетов;
- если число шагов с замороженной матрицей достигло заданного максимального числа iqк^,
- если прогнозируемый шаг больше последнего успешного в дк раз.
Числа iqк и дк могут влиять на перераспределение вычислительных затрат: при iqк = 0и дк = 0 замораживания не происходит, при увеличении iqк и дк число вычислений правой части возрастает, а количество обращений матрицы Якоби убывает. Построенный алгоритм с автоматическим выбором явной или ¿-устойчивой численной схемы будем называть ИКМК2.
Расчеты проводились на 1ВМ РС АШоп(1т)ХР 2000 с двойной точностью. В конкретных расчетах левая часть неравенства для контроля точности вычислялась по формуле
| к - к |
|| к2 - к, || тах—-——,
2 1 1*<N | Уп | +Г
где г - положительный параметр. Если по г-й компоненте решения выполняется неравенство 1уп'1 < г, то контролируется абсолютная ошибка ге, в противном случае - относительная ошибка £. В расчетах параметр г выбирался таким образом, чтобы по всем компонентам решения фактическая точность была не хуже задаваемой. Расчеты проводились с точностью £ = 10-2. Это связано с тем, что в построенном алгоритме применяются схемы низкого порядка точности и осуществлять расчеты с более высокой точностью данным методом нецелесообразно.
Сравним эффективность алгоритма ИКМК2 с эффективностью известного метода Гира в реализации А. Хин-дмарша DLSODE из коллекции ODEPACK1.
Ниже через г/и и гуа будут обозначены соответственно суммарное число вычислений правой части и количество обращений матрицы Якоби задачи (1), которые позволяют объективно оценить эффективность алгоритма интегрирования. В качестве тестового примера выбрана простейшая модель реакции Белоусова-Жаботинского [11]
У' = 77,27(У2 - уЛ + у -8,37510-6у2),
У2 = ( У2 - У1У2 + Уз ) / 77,27,
У3З = 0,161(У! - У3Х
1 е [0,300], У1(0) = 4, У2(0) = 1,1, У3(0) = 4, к0 = 10-3. (14)
Расчеты проводились с численной матрицей Якоби. Решение данной задачи алгоритмом ИКМК2 вычислено с затратами г/и = 1 154и гуа = 83. При расчетах только по ¿-устойчивой схеме (2) затраты г/и = 926 и гуа = 109. Фактическая точность расчетов в конце интервала интегрирования не хуже задаваемой. Решение (14) удалось вычислить явным методом с затратами г/и = 2 121 221. Конечно, данная задача слишком жесткая для явных методов, однако результаты расчетов приведены здесь, чтобы продемонстрировать принципиальную возможность применения явных методов с контролем устойчивости для решения достаточно жестких примеров. Данные алгоритмы на некоторых задачах большой размерности могут быть эффективнее ¿-устойчивых методов.
При расчетах программой DLSODE требуемая точность 10-2 достигается при задаваемой точности 10-4с затратами г/и = 1 129 и гуа = 107. При более высокой точности расчетов DLSODE эффективнее алгоритма ККМК2, потому что последний алгоритм основан на схемах низкого порядка точности и проводить с его помощью высокоточные расчеты нецелесообразно. Но при задаваемой точности 10-2 алгоритм ККМК2 почти в 1,5 раз эффективнее метода DLSODE по числу обращений матрицы Якоби, в то время как количество вычислений правой части задачи (14) для RKMK2 и DLSODE различается незначительно. А в случае большой размерности задачи (1) построенный алгоритм интегрирования по времени счета может быть существенно эффективнее DLSODE. Такая же тенденция сохраняется при интегрировании 10 тестовых примеров из [11].
Алгоритм ЯКМК2 предназначен для расчетов с небольшой точностью (порядка 1% и ниже). В этом случае
1 Подробная информация об этой коллекции представлена на сайте http://www.netlib.org/odepack/index.html.
достигается его максимальная эффективность. В этом алгоритме с помощью признака можно задавать различные режимы расчета:
- явным методом первого порядка точности с контролем или без контроля устойчивости;
- ¿-устойчивым методом с замораживанием или без замораживания как аналитической, так и численной матрицы Якоби;
- с автоматическим выбором численной схемы.
Все это позволяет применять данный алгоритм для
решения как жестких, так и нежестких задач. При расчетах с автоматическим выбором численной схемы вопрос о том, является ли задача жесткой или нет, перекладывается на алгоритм интегрирования.
Использование неравенства для контроля устойчивости фактически не приводит к увеличению вычислительных затрат, потому что оценка максимального собственного числа матрицы Якоби системы (1) осуществляется через ранее вычисленные стадии и не приводит к росту числа вычислений функции f Такая оценка получается грубой, однако применение контроля устойчивости в качестве ограничителя на рост шага позволяет избежать негативных последствий грубости оценки. Более того, в некоторых случаях это даже приводит к нестандартно высокому повышению эффективности алгоритма, например когда на участке установления за счет контроля устойчивости старые ошибки стремятся к нулю, а новые невелики за счет малости производных решения. В некоторых случаях вместо оценки максимального собственного числа оценивается число, следующее по порядку. Шаг интегрирования становится больше максимально допустимого и интегрирование осуществляется с таким шагом до тех пор, пока не нарушится неравенство для контроля точности. Как правило, число таких шагов невелико. Но величина шага может на порядок превышать максимальный шаг по устойчивости. После нарушения неравенства для контроля точности шаг уменьшается до максимально возможного. Такой эффект может повторятся многократно в зависимости от длины участка установления. В результате средний шаг интегрирования может превышать максимально допустимый.
Библиографический список
1. Хайрер, Э. Решение обыкновенных дифференциальных уравнений. Жесткие и дифференциально-алгебраические задачи / Э. Хайрер, Г. Ваннер. М. : Мир, 1999.
2. Byrne, G. D. ODE solvers: a review of current and coming attractions / G. D. Byrne, A. C. Hindmarsh // J. of Comput. Physics. 1987. № 70. P. 1-62.
3. Rosenbrock, H. H. Some general implicit processes for the numerical solution of differential equations / H. H. Rosenbrock // Computer. 1963. № 5. P. 329-330.
4. Новиков, E. А. Замораживание матрицы Якоби в методе типа Розенброка второго порядка точности / E. А. Новиков, В. А. Новиков, Л. А. Юматова // Журн. вы-числ. математики и мат. физики. 1987. Т. 27, № 3. С. 385-390.
5. Новиков, E. А. Замораживание матрицы Якоби в методах типа Розенброка / E. А. Новиков, А. Л. Двинский // Вычисл. технологии. 2005. Т. 10. С. 108-114.
6. Новиков, E. А. Явные методы для жестких систем / E. А. Новиков. Новосибирск : Наука. Сиб. изд. фирма Сиб. отд-ния Рос. акад. наук, 1997.
7. Новиков, E. А. Одношаговые безытерационные методы решения жестких систем / E. А. Новиков, Ю. А. Шитов, Ю. И. Шокин //Докл. АН СССР. 1988. Т. 301, № 6. С. 1310-1314.
8. Новиков, E. А. Алгоритм интегрирования жестких систем на основе (m, к)-метода второго порядка точности с численным вычислением матрицы Якоби : препр. / E. А. Новиков, Ю. А. Шитов ; ВЦ Сиб. отд-ния АН СССР. Красноярск, 1988.
9. Демидов, Г. В. Оценка ошибки одношаговых методов интегрирования обыкновенных дифференциальных уравнений / Г. В. Демидов, E. А. Новиков // Числ. методы механики сплош. среды. 1985. Т. 16, № 1. С. 27-42.
10. Кнауб, Л. В. Алгоритм интегрирования переменного порядка и шага на основе явного двухстадийного метода Рунге-Кутта / Л. В. Кнауб, Ю. М. Лаевский, E. А. Новиков // Журн. вычисл. математики. 2007. Т. 10, №2. С. 177-185.
11. Enright, W. H. Comparing numerical methods for the solutions of systems of ODE’s / W. H. Enright, T. E. Hull // BIT. 1975. Vol. 15. P. 10-48.
E. A. Novikov
AN ALGORITHM OF VARIABLE STRUCTURE FOR SOLVING STIFF SYSTEMS BASED ON AN EXPLICIT AND L-STABLE METHOD
An explicit Runge-Kutta type scheme with extended stability region is developed and an L-stable (m, k)-method of order two’s built. An algorithm of variable step and order based on these methods is obtained where the most effective numerical scheme is chosen for each step basing on stability and accuracy control inequality. The results are given that confirm the effectiveness of this algorithm for accuracy of 1 % and below.
