CC BY
56
УДК 519.622
E. A. Новиков
АЛГОРИТМ ИНТЕГРИРОВАНИЯ С ПРИМЕНЕНИЕМ L-УСТОЙЧИВОГО И ЯВНЫХ МЕТОДОВ1
Аннотация. Актуальность и цели. При моделировании кинетики химических реакций, расчете электронных схем и электрических сетей и других важных приложений возникает необходимость решения задачи Коши для жестких систем обыкновенных дифференциальных уравнений. Материалы и методы. Для решения таких задач применяются Z-устойчивые численные схемы. В таких методах при большой размерности системы дифференциальных уравнений основные вычислительные затраты приходятся на декомпозицию матрицы Якоби. Сокращения затрат достигают замораживанием матрицы Якоби, т.е. применением одной матрицы на нескольких шагах интегрирования. Дополнительного сокращения затрат добиваются за счет применения алгоритмов интегрирования на неоднородных схемах. В состав таких алгоритмов включаются явные и L-устойчивые методы. Эти алгоритмы сами распознают, является задача жесткой или нет. Эффективная численная схема выбирается на каждом шаге по критерию устойчивости. Здесь разработан неоднородный алгоритм интегрирования на основе Z-устойчивого и явных двухстадийных методов. Построено неравенство для контроля устойчивости схемы Рунге - Кутта второго порядка точности. На основе стадий этого метода предложена численная формула первого порядка с расширенным до 8 интервалом устойчивости. На основе Z-устойчивой (2,2)-схемы и численных формул типа Рунге - Кутта первого и второго порядков точности разработан алгоритм переменной структуры, в котором эффективный метод выбирается на каждом шаге по критерию устойчивости. При расчетах по Z-устойчивому методу допускается замораживание матрицы Якоби, которая может вычисляться как аналитически, так и численно. Алгоритм предназначен для решения как жестких, так и нежестких задач. Результаты. Приведены результаты расчетов, подтверждающие эффективность построенного алгоритма.
Ключевые слова: жесткая задача, (да,&)-схемы, методы Рунге - Кутта, контроль точности и устойчивости
E. A. Novikov
THE INTEGRATION ALGORITHM USING THE L-STABLE AND EXPLICIT METHODS
Abstract. Background. When modeling the kinetics of chemical reactions, calculating electronic circuits and electrical networks and other critical applications there is a need to solve the Cauchy problem for stiff systems of ordinary differential equations. Materials and methods. To solve these problems L-stable numerical schemes are applied. In such methods basic computing costs fall on the decomposition of the Jacobi matrix due to high dimensionality of the differential equations system. Cost reduction is achieved by freezing the Jacobian matrix, i.e. using the same matrix in several integration steps. Additional cost savings are achieved by applying heterogeneous integration algorithms. The structure of such algorithms includes explicit and L-stable methods. These algorithms automatically recognize whether the problem is stiff or not. An efficient numerical scheme is chosen at each step according to
1 Работа выполнена при поддержке Российского фонда фундаментальных исследований (проект № 11-01-000106).
a criterion of stability. The heterogeneous integration algorithm based on the L-stable and the explicit two-stage methods is constructed. The inequality for stability control of the second order Runge-Kutta scheme is constructed. The numerical formula of the first order with the extended up to 8 interval of stability based on stages of this method has been proposed. Based on the L-stable (2,2)-scheme and numerical formulas of the first and second orders of accuracy Runge - Kutta method, the algorithm of variable structure in which the most efficient method is chosen at each step according to the criterion of stability has been constructed. The Jacobi matrix freezing is allowed while using the L-stable method. Besides, the matrix can be calculated both analytically and numerically. The algorithm is designed for the solution of both stiff and non-stiff problems. Results. Numerical results confirm the efficiency of the algorithm.
Key words: stiff problems, (m, k)-schemes, the Runge - Kutta method, the accuracy and stability control.
Введение
Во многих приложениях возникает проблема численного решения задачи Коши для обыкновенных дифференциальных уравнений. Для решения жестких задач в основном применяются неявные методы, в которых основные затраты приходятся на декомпозицию матрицы Якоби. Для повышения эффективности расчетов в ряде алгоритмов используется замораживание матрицы Якоби, т.е. применение одной матрицы на нескольких шагах интегрирования [1]. Наиболее успешно этот подход применяется в алгоритмах на основе многошаговых численных формул [2] и для методов, в которых стадии вычисляются с участием матрицы Якоби в некотором итерационном процессе. В алгоритмах интегрирования на основе известных безытерацион-ных численных схем, к которым относятся методы типа Розенброка [3] и их различные модификации [1, 4-5], матрица Якоби включена непосредственно в численную формулу, и ее аппроксимация может приводить к понижению порядка точности. В работе [6] для методов типа Розенброка доказано, что максимальный порядок точности данных численных схем равен двум, если в алгоритме интегрирования одна матрица Якоби применяется на нескольких шагах интегрирования. Это означает, что применение методов типа Розенброка эффективно при решении задач небольшой размерности или при небольшой точности расчетов.
Некоторым аналогом замораживания матрицы Якоби является применение в расчетах алгоритмов интегрирования на основе явных и Z-устойчивых методов с автоматическим выбором численной схемы [7-8]. В этом случае эффективность алгоритма может быть повышена за счет расчета некоторых переходных участков по явному методу. В качестве критерия выбора эффективной численной формулы естественно применять неравенство для контроля устойчивости [9-10]. Применение комбинированных алгоритмов полностью не снимает проблему замораживания матрицы Якоби [11].
Здесь построено неравенство для контроля устойчивости метода типа Рунге - Кутта второго порядка точности. На основе стадий данной численной формулы построена схема первого порядка с расширенной до восьми единиц по вещественной оси областью устойчивости. На основе Z-устойчивой (2,2)-схемы и рассмотренных явных численных формул разработан алгоритм переменной конфигурации, в котором эффективная численная формула вы-
бирается на каждом шаге по критерию устойчивости. Алгоритм предназначен для решения жестких и нежестких задач. Приведены результаты вычислений, подтверждающие эффективность построенного алгоритма.
1. Исследование (2,2)-схемы
Для численного решения задачи Коши для системы обыкновенных дифференциальных уравнений
у = f(y), У = У0, (0 < 1 < 1к, (1)
где у и f - вещественные ^-мерные вектор-функции; 1 - независимая переменная, рассмотрим (2,2)-метод вида [4]
Уп+1 = Уп + Р1к1 + Р2к2 , °пк1 = кпГ(Уп ) , Впк2 = кпГ(Уп + Рк1) + ак1 (2)
с коэффициентами
Р1 = в = а = 1 -72/2, Р2 = 0,5 / а , а = -2а. (3)
Коэффициенты (3) приведены для наглядности, а ниже их выбор обоснован. Здесь кп - шаг интегрирования; к\ и к2 - стадии метода; Вп = Е - акпЛп; Е - единичная матрица; Лп - матрица, представимая в виде
Лп = Гп + кпВп + 0(кп),
f 'п = д/(уп)/ду - матрица Якоби системы (1); Вп - не зависящая от шага интегрирования матрица; а, а, в, р1 ир2 - числовые коэффициенты.
Использование матрицы Лп, представимой в виде Лп = Гп + кпВп + 0(к2п), позволяет применять (2) с замораживанием как аналитической, так и численной матрицы Якоби [5]. В случае использования матрицы Якоби Г’„_к, вычисленной к шагов назад, имеем
Лп = ЭГ (ду-к) = Г - ккпГГ + (Лп2), Вп =-тп, 0 < к < к, дУ
где - максимальное число шагов с замороженной матрицей Якоби, ГГ = [дГ(уп)/ду2] • Г(уп). Если матрица Якоби вычисляется численно с шагом Г = с,к, где с, 1 <, < N - постоянные числа, то получим
Лп = + Квп + 0(кп ) , Вп = вп ,
где элементы gr^ij матрицы Сп определяются по формулам gnЛj = 0,5суЗГуп)/ду/,
1 < i,, < N. В случае замораживания численной матрицы Якоби можно записать
Так как при записи схемы (2) матрица Вп произвольная, то вопрос о замораживании и численной аппроксимации матрицы Якоби можно рассматривать одновременно.
Теперь перейдем к исследованию схемы (2). Ниже для сокращения записи будем полагать к = кп, если это не будет приводить к недоразумениям.
№ 3 (27), 2013 Физико-математические науки. Математика
Разлагая к! и к2 в ряды Тейлора в окрестности точки уп до членов с к3 включительно и подставляя в первую формулу (2), получим
Уп+1 = Уп + [ Р1 + (1 + а) Р2 ] к/п + [аР1 + (а + Р + 2аа) Р2 ] к2 Л/,п +
+
а2р1 + (2+2аР+3а2а) )Гп ^ + 2р2Р2^ /п/'п +
+а [/1 +(1 + 2а) Р2 ]к>Вп/п + 0(к4^
где элементарные дифференциалы вычислены на приближенном решении Уп. Разложение точного решения У(1п+1) в ряд Тейлора в окрестности точки 1п до членов с к3 включительно имеет вид
У('п+1) = У(1п) + к/ + 2к2// + 6к3 [/'2/ + //2 ] + 0(к4),
2 6 ь J
где элементарные дифференциалы вычислены на точном решении У(1п).
Запишем условия второго порядка точности схемы (2). Сравнивая полученные ряды при условии Уп = у(1п) до членов с к2 включительно, имеем
Р1 + (1 + а) Р2 = 1, аР1 + (а + р + 2аа)Р2 = 2.
В приближенное решение входит элементарный дифференциал В/, а в точном решении он отсутствует. Данный дифференциал появился за счет применения «испорченной» матрицы Якоби. Потребуем, чтобы коэффициент при В/п в ряде Тейлора дляУп+1 был равен нулю, т.е. положим
Р +(1 + 2а)Р2 = 0 .
Теперь исследуем устойчивость схемы (2). Для этого применим ее для решения скалярной тестовой задачи У '= Ху, у(0) = У0, Яе(Х) < 0. Получим Уп+1 = б(х)Уп, где х = кХ, а функция устойчивости Q(x) с применением условий порядка имеет вид
1 + (1 - 2a )х + a2 - a (p + Р2 ) + вР2
Q( x) =----------------------------------------------------^-2-:
х2
(1 - ах)
Отсюда условие а2 - а(р1 + р2) + вр2 = 0 обеспечивает Х-устойчивость метода (2).
В процессе вычислений по схеме (2) приближение к решению Уп+в = Уп + Ркг вычисляется в промежуточной точке 1п+р. Основная схема (2) Х-устойчивая, поэтому естественно того же свойства потребовать от промежуточной численной формулы. В противном случае ошибки в промежуточных вычислениях могут существенно влиять на итоговую ошибку. Применяя схему Уп+р = Уп + вк для решения уравнения У' = Ху, получим Уп+р = [1 + (в -- а)х)]Уп/(1 - ах). Отсюда следует, что промежуточная схема будет Х-устой-чивой, если в = а.
Исследуем совместность условий порядка и устойчивости. Для этого перепишем их в следующем виде:
1) Pl +(1 + a)p2 =1;
2) a _ pi + (1 + a) P2 ] + aap2 + PP2 = 0,5; 3) a _pi + (1 + a)p2 ] + aap2 = 0 ;
4) a2 - a _ pi +(l + a)p2 ] + aap2 + PP2 = 0; 5) P = a .
(4)
Вычитая из второго равенства (4) третье, получим Рр2 = 0,5. Отсюда запишем р2 = 0,5/а. Из первого и третьего равенств (4) получим ар2 = -1. Тогда имеем а = -2а. Из первого соотношения выразим р1 = 2 - 0,5/а, а из четвертого соотношения (4) получим уравнение а2 - 2а + 0,5 = 0, которое обеспечивает Х-устойчивость схемы (2). Теперь локальную ошибку 5„ метода (2) можно записать в виде
8n =>
a -Ї )■Г2/ + [1 -~4 a I-//2
+ O(h4).
Тогда согласно [9] для контроля точности (2) можно использовать оценку ошибки вида 8п = (а - 1/3)к2/2/ + 0(к3). Учитывая, что к2 + (2а - 1)кх =
2 2 2 3 3
= (а - 2а2)к// + 0(к ), величину еп с точностью до членов 0(к ) можно оценить по формуле еп = [(а - 1/3)/(а - 2а2)][к2 + (2а - 1)к1 ]. В результате для контроля точности схемы (2) можно применять неравенство [9]
D1-jn _k2 +(2a-1)]
<
a - 2a
a -1/3
e, 1 < jn < 2,
(5)
где 8 - требуемая точность интегрирования; ||-|| - некоторая норма в ^, а значение целочисленного параметра ,п выбирается наименьшим, при котором выполняется неравенство (5). Уравнение а - 2а + 0,5 = 0 имеет два корня
а! = 1 - 0,5 л/2 и а2 = 1 + 0,5 л/2 . Выберем а = аь потому что в этом случае меньше коэффициент в главном члене (а - 1/3)к3/'2/’локальной ошибки 5п. Теперь коэффициенты схемы (2) определяются однозначно и совпадают с (3)
Оценку максимального собственного числа ^п0 = кХп,тах матрицы Якоби системы (1), необходимую для перехода на явную схему, оценим по формуле
wn, 0 =>
д/
дУ
= h • max
1<i< N
N
z
j=1
д/-
дУ,
2. Явный метод второго порядка
Для численного решения задачи (1) рассмотрим явный метод вида
Уп+1 = Уп + Р1к1 + Р2к2 , к1 = к/(Уп ), к2 = к/(Уп +вк1) . (6)
Получим соотношения на коэффициенты метода (6) второго порядка точности. Для этого разложим стадии к1 и к2 в ряды Тейлора и подставим в первую формулу (6). Получим
включительно, запишем условия второго порядка точности схемы (6), т.е. Р1 + р2 = 1 и вр2 = 0,5. Тогда локальная ошибка 5п,2 схемы (6) будет иметь вид
Для контроля точности (6) применим схему первого порядка Уп+1,1 = Уп + к\. С помощью идеи вложенных методов оценку ошибки £п 2 метода второго порядка можно вычислить по формуле £п,2 = Уп+1 - Уп+1,1 = = р2(к2 - к1). Для повышения надежности данной оценки выберем в = 1. Тогда стадия к\ вычисляется в точке п а к2 - в точке ^п+1. Как показывают расчеты, использование информации в крайних точках шага приводит к более надежному неравенству для контроля точности вычислений. При в = 1 коэффициенты метода второго порядка определяются однозначно Р1 = Р2 = 0,5, а локальная ошибка 5п,2 и неравенство для контроля точности имеют соответственно вид
Теперь построим неравенство для контроля устойчивости схемы (6) предложенным в [9] способом. Для этого рассмотрим вспомогательную стадию к3 = к/(Уп+1). Заметим, что к3 совпадает со стадией кь которая применяется на следующем шаге интегрирования, и поэтому ее использование не приводит к дополнительным вычислениям правой части системы (1). Запишем стадии кь к2 и к3 применительно к задаче У' = Ау, где А есть матрица с постоянными коэффициентами. Получим
Тогда согласно [8] оценку максимального собственного числа ^п,2 = кХптах матрицы Якоби системы (1) можно вычислить степенным методом по формуле
5n 2 = (2 -3P)h3// + бh3f2/ + O(h4).
12 б
5n 2 = h3// + бh3/f2/ + O(h4), 0,51| k2 - kl ||<e,
12 б
где || • || - некоторая норма в ; е - требуемая точность интегрирования.
kl = Xyn , k2 = (X + X2)Уп , k3 = (X + X2 + 2XЗ)Уп ,
гдеX = НА. Легко видеть, что
k2 - kl = X2 Уп , 2 (з - k2 ) = X3 Уп .
Интервал устойчивости схемы (6) второго порядка точности приблизительно равен двум. Поэтому для ее контроля устойчивости можно применять неравенство ^п2 < 2. Полученная оценка ^п2 является грубой, потому что не обязательно максимальное собственное число сильно отделено от остальных, в степенном методе применяется мало итераций и дополнительные искажения вносит нелинейность задачи (1). Поэтому здесь контроль устойчивости используется как ограничитель на размер шага интегрирования. В результате прогнозируемый шаг Ип+\ будем вычислять следующим образом. Новый шаг Иас по точности определим по формуле Иас = д\Ип, где Ип есть последний успешный шаг интегрирования, а д\, учитывая соотношение к2 - к! = 0(Ип2), задается уравнением ^12||к2 - к^| = 2е. Шаг И по устойчивости зададим формулой Ий = д2Ип, где д2 определяется из равенства д2^п,2 = 2. Тогда прогнозируемый шаг кп+1 вычисляется по формуле Нп+\ = тах{Ип, тт[Иас, Ий]}.
3. Метод первого порядка точности
Для численного решения задачи (1) рассмотрим схему вида
Уп+1 = Уп + г1к1 + г2к2 , к1 = И/(Уп ), к2 = И/(Уп + к1) . (7)
Заметим, что при Т\ = г2 = 0,5 численная формула (7) имеет второй порядок точности. Построим менее точную схему с максимальным интервалом устойчивости. Для этого применим (7) для решения тестового уравнения У '= Ху. Получим Уп+1 = Q2(x)yn, где функция устойчивости 02(х) имеет вид 02(х) = 1 + (г1 + г2)х + г2х2, х = ИХ. Требование первого порядка точности приводит к соотношению Т\ + г2 = 1, которое ниже будем считать выполненным.
Теперь выберем г2 таким образом, чтобы метод (7) имел максимальный интервал устойчивости. Для этого рассмотрим многочлен Чебышева Т2(г) = = 2г2 - 1 [12] на промежутке [-1, 1]. Проведем замену переменных, полагая г = 1 - 2х/у. Получим Т2(х) = 1 - 8х/у + 8х2/у2, при этом отрезок [у,0] отображается на [-1, 1]. Нетрудно показать, что среди всех многочленов вида Р2(х) = 1 + х + с2х2 для Т2(х) неравенство |Т2(х)| < 1 выполняется на максимальном интервале [у, 0] при у = -8. Потребуем совпадения коэффициентов Q2(x) и Т2(х) при у = -8. Это приводит к соотношениям Г1с + г2 = 1 и г2 = 1/8. В результате имеем коэффициенты Т\ = 7/8 и г2 = 1/8 метода порядка точности с максимальным интервалом устойчивости, локальная ошибка 5п,1 которого имеет вид 5п,1 = 3И2///8 + 0(И3). Для контроля точности численной формулы первого порядка будем использовать оценку локальной ошибки. С учетом к2 - к1 = И2/'п/п + 0(И3) и вида локальной ошибки неравенство для контроля точности записывается в виде 3||к2-к1||/8 < е, где е есть требуемая точность расчетов.
Построим неравенство для контроля устойчивости метода первого порядка точности. Для этого рассмотрим вспомогательную стадию к3 = И/у^). Заметим, что снова к3 совпадает со стадией кь которая применяется на следующем шаге интегрирования, и поэтому ее использование не приводит к дополнительным вычислениям правой части системы (1). Запишем стадии кь к2 и к3 применительно к задаче у' = Ау, где А есть матрица с постоянными коэффициентами. В результате получим
kl = Xy, , k2 = (X + X2)У, , k3 = (X + X2 + 0.125X3)y, ,
гдеX = НА. Легко видеть, что
k2 - k1 = X2yn , S(k3 - k2 ) = X3yn .
Тогда согласно [9] оценку максимального собственного числа ^п,1 = ИХптах матрицы Якоби системы (1) можно вычислить по формуле
Интервал устойчивости схемы (7) первого порядка равен восьми. Поэтому для ее контроля устойчивости можно применять неравенство ^п,1 < 2. Выбор прогнозируемого шага выбирается по аналогии с методом второго порядка.
На основе построенных явных методов первого и второго порядков точности легко сформулировать алгоритм переменного порядка и шага. Расчеты всегда начинаются методом второго порядка как более точным. Переход на схему первого порядка осуществляется при нарушении неравенства ^п2 < 2. Обратный переход на метод второго порядка происходит в случае выполнения неравенства ^пд < 2. При расчетах по методу первого порядка наряду с точностью контролируется устойчивость ^пд < 8, а выбор прогнозируемого шага производится по аналогии с методом второго порядка точности по формуле вида Ип+ = тах{Ип, тт[Иас, ИА]}.
В случае использования схемы (2) формулировка алгоритма интегрирования также не вызывает трудностей. Нарушение неравенства ^пд < 8 вызывает переход на Х-устойчивую схему (2). Передача управления явным методам происходит в случае выполнения неравенства ^п0 < 8, где оценка ^п0 вычисляется через норму матрицы Якоби.
В расчетах шаг численного дифференцирования г}- по }-й компоненте выбирается по формуле Г} = тах(1014, У}]-10~7). Расчеты предполагается проводить с двойной точностью, т.е. при представлении чисел мантисса содержит 14 значащих цифр. Постоянное число 10-7 введено для того, чтобы выдвинуть шаг численного дифференцирования на середину разрядной сетки.
Численную формулу (2) без потери порядка точности можно применять с замораживанием матрицы Вп. Отметим, что при замораживании матрицы Якоби величина шага интегрирования остается постоянной с целью сохранения свойства Х-устойчивости метода (2). Попытка замораживания матрицы Вп осуществляется после каждого успешного шага. Размораживание матрицы происходит в следующих случаях: 1) если нарушена точность расчетов; 2) если число шагов с замороженной матрицей достигло заданного максимального числа /И; 3) если прогнозируемый шаг больше последнего успешного в дИ раз. Числами гИ и дИ можно влиять на перераспределение вычислительных затрат. При гИ = 0 и дИ = 0 замораживания не происходит; при увеличении гИ и дИ число вычислений правой части возрастает, а количество обращений матрицы Яко-
4. Алгоритм интегрирования переменной структуры
би убывает. В случае большой размерности системы дифференциальных уравнений (1) имеет смысл выбирать 7ь и достаточно большими числами. В расчетах использовались значения 7ь = 10 и = 2.
Норма ||ф|| в неравенствах для контроля точности вычисляются по формуле ||ф|| = шах1<г<^{|фг|/[[у„г| + V]}, где V - положительный параметр, 7 - номер компоненты. Если по 7-й компоненте решения выполняется неравенство У„7| < V, то контролируется абсолютная ошибка V • £, в противном случае - относительная ошибка £. Для опытных вычислителей рекомендуется вместо константы V выбирать постоянные v7, 1 < 7 < Ы, по каждой компоненте решения. Построенный ниже алгоритм переменного порядка и шага, а также с автоматическим выбором явной или Ь -устойчивой численной схемы будем называть шк2гк21.
5. Результаты расчетов
В качестве тестового примера рассмотрим уравнение Ван дер Поля для имитации осциллирующих физических процессов. Здесь это уравнение изучается в виде системы двух уравнений первого порядка
У1 = У2, у2 =[(1 - У12)У2 - У1]/^
с начальными условиями у1 (0) = 2 и у2(0) = 0, ^ е [0, 11]. Изменением параметра ц варьируется жесткость модели. В табл. 1 приведены результаты расчетов при различных значениях ц. Вычислительные затраты приведены в форме 1А(у), где через if и у обозначены соответственно число вычислений правой части и декомпозиций матрицы Якоби на интервале интегрирования. Сравнение эффективности построенного алгоритма шк2гк21 проводилось с методом Гира в реализации А. Хиндмарша dlsode из коллекции ОББРЛСК [13]. Расчеты проводились таким образом, чтобы в точном и приближенном решениях совпадали две значащие цифры, где под точным понимается решение при точности вычислений е = 10-11 различными методами интегрирования.
Таблица 1
Вычислительные затраты при различных значениях ц
м 10-1 10-2 10-3 10-4 10-5 10-6
mk2rk21 2412(0) 5745(0) 8279(182) 9701(265) 11718(358) 13041(451)
dlsode 1756(118) 4634(310) 5213(461) 6503(600) 8021(733) 9357(841)
Из анализа результатов расчетов следует, что при небольшой жесткости задачи в построенном алгоритме работают только явные методы. Начиная с ц = 10-3 в зависимости от условия устойчивости в процессе вычислений комбинируются явные и Ь-устойчивые схемы. Из изучения пошаговых вычислений видно, что на переходных участках расчеты осуществляются по явным формулам, а на участках установления - по Ь-устойчивой схеме. По количеству декомпозиций матрицы Якоби построенный алгоритм шк2гк21 при всех значениях ц эффективнее алгоритма dlsode, но по числу вычислений правой части if метод шк2гк21 уступает. Это означает, что на задачах большой размерности шк2гк21 может быть эффективнее.
Заключение
Максимальная эффективность построенного алгоритма достигается при небольшой точности расчетов - порядка 1 % и ниже. В mk2rk21 с помощью признака можно задавать различные режимы расчета: 1) явными методами первого или второго порядков точности с контролем или без контроля устойчивости; 2) явными методами с переменным порядком и шагом; 3) Z-устойчивым методом с замораживанием или без замораживания численной или аналитической матрицы Якоби; 4) с автоматическим выбором численной схемы. Все это позволяет применять данный алгоритм для решения как жестких, так и нежестких задач. При расчетах с автоматическим выбором численной схемы вопрос о том, является задача жесткой или нет, перекладывается на алгоритм интегрирования.
Использование неравенства для контроля устойчивости фактически не приводит к увеличению вычислительных затрат, потому что оценка максимального собственного числа матрицы Якоби системы (1) осуществляется через ранее вычисленные стадии и не приводит к росту числа вычислений функции f. Такая оценка получается грубой. Однако применение контроля устойчивости в качестве ограничителя на рост шага или для выбора эффективной численной формулы позволяет избежать негативных последствий грубости оценки. В некоторых случаях это приводит к нестандартно высокому повышению эффективности явных методов. На участке установления за счет контроля устойчивости старые ошибки стремятся к нулю, а новые невелики за счет малости производных решения. В некоторых случаях вместо оценки максимального собственного числа оценивается следующее по порядку. Шаг интегрирования становится больше максимально допустимого, и с таким шагом осуществляется интегрирование до тех пор, пока не нарушается неравенство для контроля точности. Как правило, число таких шагов невелико. Однако величина шага может на порядок превышать максимальный шаг по устойчивости. После нарушения неравенства для контроля точности шаг уменьшается до максимально возможного. Такой эффект может повторяться многократно в зависимости от длины участка установления. В результате средний шаг интегрирования может превышать максимально допустимый. Применение на участке установления явного метода первого порядка точности с расширенной областью устойчивости позволяет в 4 раза увеличить размер шага интегрирования по сравнению с явным методом второго порядка без увеличения вычислительных затрат. На переходных участках, где определяющую роль играет точность вычислений, более эффективным является метод второго порядка точности, хотя и с небольшой областью устойчивости. Комбинирование методов низкого и высокого порядков с помощью неравенства для контроля устойчивости позволяет повысить эффективность расчетов.
Список литературы
1. Хайрер, Э. Решение обыкновенных дифференциальных уравнений. Жесткие и дифференциально-алгебраические задачи / Э. Хайрер, Г. Ваннер. - М. : Мир, 1999. - 685 с.
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. Н овиков, Е. А. Одношаговые безытерационные методы решения жестких систем / Е. А. Новиков, Ю. А. Шитов, Ю. И. Шокин // ДАН СССР. - 1988. -Т. 301, № 6. - С. 1310-1314.
5. Н овиков, Е. А. Аппроксимация матрицы Якоби в (да,2)-методах решения жестких задач / Е. А. Новиков // Журнал вычислительной математики и математической физики. - 2011. - Т. 51, № 12. - С. 2194-2208.
6. Новиков, В. А. Замораживание матрицы Якоби в методе типа Розенброка второго порядка точности / В. А. Новиков, Е. А. Новиков, Л. А. Юматова // Журнал вычислительной математики и математической физики. - 1987. - Т. 27, № 3. -С. 385-390.
7. Н овиков, Е. А. Построение алгоритма интегрирования жестких систем дифференциальных уравнений на неоднородных схемах / Е. А. Новиков // ДАН СССР. - 1984. - Т. 278, № 2. - С. 272-275.
8. Novikov, A. E. Numerical Integration of Stiff Systems with Low Accuracy / A. E. Novikov, E. A. Novikov // Mathematical Models and Computer Simulations. -2010. - Vol. 2, № 4. - P. 443-452.
9. Н овиков, Е. А. Явные методы для жестких систем / Е. А. Новиков. - Новосибирск : Наука, 1997. - 197 с.
10. Новиков, Е. А. Численное моделирование пиролиза этана явным методом третьего порядка точности / Е. А. Новиков // Известия высших учебных заведений. Поволжский регион. Физико-математические науки. - 2010. - № 4. -С. 64-72.
11. Новиков, Е. А. Компьютерное моделирование жестких гибридных систем / Е. А. Новиков, Ю. В. Шорников. - Новосибирск : Изд-во НГТУ, 2012. - 450 с.
12. Бахвалов, Н. С. Численные методы / Н. С. Бахвалов, Н. П. Жидков, Г. М. Кобельков. - М. : Наука, 1987. - 598 с.
13. Hindmarsh, A. C. ODEPACK, a systematized collection of ODE solvers / A. C. Hindmarsh. - Lawrence Livermore National Laboratory, 1982. - Preprint UCRL-88007.
References
1. Khayrer E., Vanner G. Reshenie obyknovennykh differentsial'nykh uravneniy. Zhestkie i differentsial'no-algebraicheskie zadachi [Solution of regular differential equations. Stiff and differential algebraic problems]. Moscow: Mir, 1999, 685 p.
2. Byrne G. D., Hindmarsh A. C. J. of Comput. Physics. 1987, no. 70, pp. 1-62.
3. Rosenbrock H. H. Computer. 1963, no. 5, pp. 329-330.
4. Novikov E. A., Shitov Yu. A., Shokin Yu. I. DAN SSSR [Reports of the Academy of Sciences of the USSR]. 1988, vol. 301, no. 6, pp. 1310-1314.
5. Novikov E. A. Zhurnal vychislitel'noy matematiki i matematicheskoy fiziki [Journal of calculus mathematics and mathematical physics]. 2011, vol. 51, no. 12, pp. 2194-2208.
6. Novikov V. A., Novikov E. A., Yumatova L. A. Zhurnal vychislitel'noy matematiki i matematicheskoy fiziki [Journal of calculus mathematics and mathematical physics]. 1987, vol. 27, no. 3, pp. 385-390.
7. Novikov E. A. DAN SSSR [Reports of the Academy of Sciences of the USSR]. 1984, vol. 278, no. 2, pp. 272-275.
8. Novikov A. E., Novikov E. A. Mathematical Models and Computer Simulations. 2010, vol. 2, no. 4, pp. 443-452.
9. Novikov E. A. Yavnye metody dlya zhestkikh sistem [Explicit methods for stiff systems]. Novosibirsk: Nauka, 1997, 197 p.
10. Novikov E. A. Izvestiya vysshikh uchebnykh zavedeniy. Povolzhskiy region. Fiziko-matematicheskie nauki [University proceedings. Volga region. Physical and mathematical sciences]. 2010, no. 4, pp. 64-72.
11. Novikov E. A., Shornikov Yu. V. Komp’yuternoe modelirovanie zhestkikh gibridnykh system [Computer modeling of stiff hybrid systems]. Novosibirsk: Izd-vo NGTU, 2012, 450 p.
12. Bakhvalov N. S., Zhidkov N. P., Kobel'kov G. M. Chislennye metody [Numerical methods]. Moscow: Nauka, 1987, 598 p.
13. Hindmarsh A. C. ODEPACK, a systematized collection of ODE solvers. Lawrence Livermore National Laboratory, 1982. Preprint UCRL-88007.
Новиков Евгений Александрович
доктор физико-математических наук, главный научный сотрудник, Институт вычислительного моделирования Сибирского отделения Российской академии наук (Россия, г. Красноярск, Академгородок, д. 50, стр. 44)
E-mail: [email protected]
Novikov Evgeniy Aleksandrovich Doctor of physical and mathematical sciences, senior scientific associate, Institute of Computational Modeling of Siberian Branch of Russian Academy of Sciences (44 Building,
50 Akademgorodok, Krasnoyarsk, Russia)
УДК 519.622 Новиков, Е. А.
Алгоритм интегрирования с применением ^-устойчивого и явных методов / Е. Л. Новиков // Известия высших учебных заведений. Поволжский регион. Физико-математические науки. - 2013. - № 3 (27). - С. 58-69.
