CC BY
67
дельной цепи А совпадают, поэтому способ вычисления стационарного распределения цепи А является также и одним из способов определения стационарного распределения сети N. Преимуществом этого способа является возможность его применения для вычисления стационарного распределения достаточно широкого класса сетей обслуживания, а недостатком — необходимость выполнения большого объема вычислений даже для сетей обслуживания средней размерности. Для сетей обслуживания, обладающих свойством локального равновесия, к которым относится сеть N, эффективным в отношении объема вычислений является способ вычисления стационарного распределения с использованием мультипликативной формы стационарного распределения. При решении практических задач анализа сетей обслуживания не всегда требуется определение стационарного распределения. Поэтому в статье предложен также способ вычисления стационарных характеристик сети N с использованием рекурсивного метода анализа сетей обслуживания.
Библиографический список
1. Kelly F. P. Reversibility and stochastic networks. London: Wiley, 1979. 230 p.
2. Уолренд Дж. Введение в теорию сетей массового обслуживания. М.: Мир, 1993. 335 с.
3. Митрофанов Ю. И. Анализ сетей массового обслуживания. Саратов: Научная книга, 2005. 177 с.
4. Henderson W, Pearce C.E.M., Taylor P. G., Dijk N.M. Closed queueing networks with batch services // Queueing Systems. 1990. Vol. 6. P. 59-70.
5. Henderson W, Taylor P. G. Product form in networks of queues with batch arrivals and batch services // Queueing Systems. 1990. Vol. 6. P. 71-88.
6. Boucherie R. J., Dijk N. M. Product forms for queueing networks with state-dependent multiple job transitions // Advances in Applied Probability. 1991. Vol. 23, № 1. P. 152-187.
7. Serfozo R. F. Queueing networks with dependent nodes and concurrent movements // Queueing Systems. 1993. Vol. 13. P. 143-182.
8. Miyazawa M. Structure-reversibility and departure functions of queueing networks with batch movements and state dependent routing // Queueing Systems. 1997. Vol. 25. P. 45-75.
9. Woodward M. E. Towards the accurate modelling of high-speed communication networks with product-form discrete-time networks of queues // Computer Communications. 1998. Vol. 21. P. 1530-1543.
10. Гурьянов А. И., Митрофанов Ю.И. Определение параметров замкнутых линейных сетей систем массового обслуживания // Системное моделирование. Новосибирск, 1970. Вып. 1. C. 39-49.
УДК 519.622
АЛГОРИТМ ПЕРЕМЕННОГО ПОРЯДКА И ШАГА НА ОСНОВЕ ЯВНОГО ТРЕХСТАДИЙНОГО МЕТОДА ТИПА РУНГЕ - КУТТА
Е. А. Новиков
Институт вычислительного моделирования СО РАН, Красноярск,
отдел вычислительной математики E-mail: novikov@icm.krasn.ru
Получено неравенство для контроля устойчивости трехстадий-ного метода Рунге - Кутта третьего порядка точности. Построен метод первого порядка с расширенной областью устойчивости. Сформулирован алгоритм интегрирования переменного порядка. Приведены результаты расчетов жестких задач, подтверждающие повышение эффективности алгоритма с переменным порядкам по сравнению с расчетами по фиксированной схеме.
Ключевые слова: жесткие задачи, явный метод, контроль точности и устойчивости, переменный порядок.
Variable Order and Step Algorithm Based on a Stages of Runge - Kutta Method of Third Order of Accuracy
E. A. Novikov
Institute of Computational Modeling SB RAS, Krasnoyarsk, Department of Calculus Mathematics E-mail: novikov@icm.krasn.ru
An inequality for the stability control of 3-stage Runge - Kutta method of 3th order of accuracy is obtained. Method of first order with expanded stability domain is constructed. Algorithm of variable order is formulated. The results of stiff system computations are provided, which confirm an increase in efficiency for the variable order method as compared to a calculation with fixed scheme.
Key words: stiff problems, explicit method, stability and accuracy control, variable order.
ВВЕДЕНИЕ
При решении задачи Коши для жестких систем обыкновенных дифференциальных уравнений большой размерности в ряде случаев возникает необходимость применения алгоритмов на основе явных методов. Алгоритмы интегрирования на основе неявных или полуявных формул, как правило, используют обращение (декомпозицию с выбором главного элемента по строке или столбцу, а иногда и по всей матрице) матрицы Якоби, что в данном случае есть отдельная трудоемкая задача. Затраты на декомпозицию порядка Ж3 арифметических операций, где N есть размерность исходной системы. Кроме того, получение элементов матрицы Якоби и составление подпрограммы ее нахождения требуют от вычислителя больших затрат личного времени. В такой ситуации предпочтительнее применять алгоритмы на основе явных формул, если жесткость задачи позволяет за разумное время получить приближенное решение.
Современные алгоритмы на основе явных методов в большинстве своем не приспособлены для решения жестких задач по следующей причине. Обычно алгоритм управления шагом интегрирования строится на контроле точности численной схемы. Это естественно, потому что основным критерием является точность нахождения решения. Однако при применении алгоритмов интегрирования на основе явных формул для решения жестких задач этот подход приводит к потере эффективности и надежности, потому что на участке установления вследствие противоречивости требований точности и устойчивости шаг интегрирования раскачивается. В лучшем случае это приводит к большому количеству возвратов (повторных вычислений решения), а шаг выбирается значительно меньше допустимого. Этого можно избежать, если наряду с точностью контролировать и устойчивость численной схемы.
В настоящее время можно выделить два подхода к контролю устойчивости [1,2]. Первый связан с оценкой максимального собственного числа матрицы Якоби /у через ее норму с последующим контролем (наряду с точностью) неравенства ||Н/У|| < О [1], где положительная постоянная Э зависит от размера области устойчивости. Ясно, что для явных методов, где матрица Якоби не участвует в вычислительном процессе, это приводит дополнительно к ее нахождению и, следовательно, к значительному увеличению вычислительных затрат. Второй подход основан на оценке максимального собственного числа Лтах матрицы Якоби степенным методом через приращения правой части системы дифференциальных уравнений с последующим контролем неравенства |НЛтах| < О [2]. Во всех рассмотренных ситуациях такая оценка фактически не приводит к увеличению вычислительных затрат [3].
Целью данной работы является разработка алгоритма интегрирования переменного порядка и шага на основе трехстадийной схемы типа Рунге - Кутта.
Для численного решения задачи Коши для системы обыкновенных дифференциальных уравнений
где у и / — вещественные Ж-мерные вектор-функции, £ — независимая переменная, Н — шаг интегрирования, , и — стадии метода, рьр2,р3, 021, 031, в32 — числовые коэффициенты, определяющие свойства точности и устойчивости (2). В случае неавтономной системы
1. МЕТОД ТРЕТЬЕГО ПОРЯДКА ТОЧНОСТИ
у' = /(у), у(£о) = уо, £о < £ < ^,
(1)
рассмотрим явный трехстадийный метод типа Рунге - Кутта вида
уп+1 = Уп + Р1 + Р2 &2 + Р3&3, &1 = Н/(уп ), ^2 = Н/(уп + 021&1), = Н/(уп + 031 + 032&2 ),
(2)
у' = /(£,у), у(£о)= уо, £о < t < ^,
схема (2) записывается в следующем виде:
уп+1 = Уп + Р1 + Р2 &2 + Р3&3,
= hf (tn, Уп), &2 = hf (tn + 021 h, Уп + 021 kl), k3 = hf (tn + [031 + 032]h, Уп + 031 ki + 032^2)-
Ниже для сокращения выкладок будем рассматривать задачу (1). Однако построенные методы можно применять для решения неавтономных задач. Получим соотношения на коэффициенты метода (2) третьего порядка точности. Для этого разложим стадии k1, k2 и k3 в ряды Тейлора до членов с h4 включительно, т.е.
k1 = hfn, k2 = hfn + 021 h2 fnfn +1021 h3fn fn +1031 h4fn'fn + O(h5 ),
2 b
k3 = hfn + (031 + 032 )h2fn fn + 021032 h3 fn2fn + 2(031 + 032 fh3 f f2n + 0.5021032 h4 fnfn f2 +
+021(031 + 032 )032h4fn fn fn + 1(031 + 032 )h4fn'f3n + O(h5). Подставим полученные разложения в первую формулу (2). В результате получим
yn+1 = yn + (Р1 + Р2 + P3)hfn + [021Р2 + (031 + 031 )P3 ]h2 Щ fn + h3[021032P3 fnfn + 0.5(0^2 +
+ (031+032 )2p3)fn fn]+h4[O.50221032P3fn fn f +021(031 +032 Ш fn + 1(031 +032)3P3 f' f3 ]+O(h5).
Здесь элементарные дифференциалы вычислены на приближенном решении yn, fn = df (yn)/dy,
fn = d2f (yn)/dy2, fn = d3f (yn)/dy3. Точное решение y(tn+1) в окрестности точки tn имеет вид
h3 h4
y(tn+1) = y(tn) + hf + 0.5h2 ff + hr[f '2 f + ff2] + ^^^[f3 f + fff2 + 3fff2 + f'f3] + O(h4),
где элементарные дифференциалы вычислены на точном решении y(tn). Сравнивая полученные ряды для точного и приближенного решений до членов с h3 включительно при условии yn = y(tn), запишем условия третьего порядка точности схемы (2), которые имеют вид
P1 + P2 + P3 = 1, 021P2 + (031 + 032 )P3 = 1, 021P2 + (031 + 032 )2P3 = 1, 021032P3 = 1 - (3)
2 3 b
Локальную ошибку ¿п численной формулы (2) можно вычислить по формуле ¿п = у(Ьп+{) — уп+1. Учитывая представления у{Ьп+1) и уп+1 в виде рядов Тейлора, получим
&п = Ь4{ 24 /'3/ + — 2в221 /'/'//2 + [1 — 021 (031 + 032)032Р3] /'//'/2 +
+ [24 — 1 в31Р2 — 1(031 + 032^] /'''/3} + 0(Н5). (4)
В нелинейной системе алгебраических уравнений (3) присутствуют два свободных коэффициента. Исследуем три варианта.
Вариант 1. Положим 021 = 031 + 032 и 031 = 032. Это означает, что приращения к2 и к3 будут вычислены в одной и той же точке + 021Н, причем вклад к1 и к2 при определении к3 учитывается одинаково. Тогда нелинейную систему (3) можно переписать в виде
1) Р1 + Р2 + Р3 = 1, 2) 021 (р2 + Р3) = 1, 3) 021 (Р2 + Р3) = 1, 4) 021032Р3 = 1.
2 3 Ь
Из второго и третьего уравнений имеем 021 = 2/3. Из соотношений 021 = 031 + 032 и 031 = 032 запишем 031 = 032 = 1/3. Из четвертого уравнения системы получим Р3 = 3/4. Из равенства Р2 + Р3 = 3/4 имеем р2 = 0. Наконец, из первого соотношения системы получим р1 = 1/4. В результате коэффициенты схемы (2) определяются однозначно и имеют вид
2 11 3
021 = 3, 031 = 032 = 3, Р1 = 4, Р2 =0, Р3 = 4. (5)
При данных соотношениях локальную ошибку ¿п схемы (2) можно записать следующим образом:
¿п = 24 н4 / '3/ + ^ Н4 /''//3 — -2 Н4 Г г''Г 2 + Н4/'/'/2 + О(Н5).
Вариант 2. Минимизируем локальную ошибку (4). Для этого, учитывая вид (4), вместо (3) рассмотрим следующую расширенную нелинейную систему:
1) Р1 + Р2 + Рз = 1, 2) 021Р2 + (вз1 + вз2)рз = 1, 3) 021Р2 + (031 + вз2)2рз = 1, 1 1 2 1 3 4) 0210з2Рз = 6, 5) вз2Рз = 12, 6) в21 (вз1 + вз2)вз2Рз = £ •
(6)
При 1.5021 = вз1 + вз2 два последних уравнения (6) совпадают. Из четвертого и пятого соотношений (6) имеем 021 = 0.5. Из второго и третьего равенств получим р2 = 1/3 и рз = 4/9. Из первого уравнения (6) запишем р1 = 2/9, а из четвертого имеем 0з2 = 3/4. Наконец, из соотношения 0з1 + 0з2 = 3/4 запишем 0з1 = 0. В результате коэффициенты метода (2) с минимальной локальной ошибкой можно записать в виде
1 3 2 14
в21 = ^' вз1 = вз2 = 4, Р1 = 9, Р2 = 3, Рз = 9 • (7)
При данных соотношениях локальную ошибку ¿п схемы (2) можно записать следующим образом:
Л4 Л4
¿п = f'зf - ;;з +
При использовании (2) с наборами коэффициентов (5) и (7) ни одна стадия не вычисляется в точке £п+1. При быстром изменении решения это может приводить к понижению точности расчетов.
Вариант 3. Положим 021 = 0.5 и 0з1 + 0з2 = 1. Тогда на каждом шаге приращения , и вычисляются в точках £п, £п + Л/2 и £п + Н соответственно. В этом случае условия третьего порядка записываются в виде
1 1 . 1 1 1
1) Р1 + Р2 + Рз = 1, 2) 2^2 + Рз = 2 > 3) 4Р2 + Рз = 3, 4) 0з2Рз = 3 •
Из второго и третьего уравнения данной системы имеем р2 = 2/3 и Рз = 1/6. Из первого и последнего уравнений имеем р1 = 1/6 и 0з2 = 2. Из равенства 0з1 + 0з2 = 1 следует 0з1 = -1. В результате коэффициенты метода (2) можно записать в следующем виде:
1 12 1
^21 = ~, вз1 = -1, вз2 =2, Р1 = т, Р2 = ^, Рз = т • (8)
2 6 3 6
При данных соотношениях локальную ошибку ¿п схемы (2) можно записать следующим образом:
Л4 г 1-,
¿п = НН^ ^/3f - f ^' f2 - 3 f''' f 1 + °(Н5).
2. КОНТРОЛЬ ТОЧНОСТИ ВЫЧИСЛЕНИИ
Неравенство для контроля точности вычислений метода третьего порядка построим с использованием идеи вложенных методов. Для этого рассмотрим вспомогательную схему
Уп+1,1 = Уп + Г1&1 + Г2 ,
где и определены в (2). Потребуем, чтобы данный метод имел второй порядок точности. Разложение приближенного решения уп+1,1 в виде ряда Тейлора по степеням Н до членов с Н2 включительно имеет вид
Уп+1,1 = Уп + (г1 + г 2 )НД + 021^2 Н2 fn ^ + 0(Нз). Сравнивая ряды для у(^п+1) и уп+1;1, видим, что данное требование будет выполнено, если г1 + г2 = 1, в21г2 = 0.5. Отсюда получим г2 = О.5021, г1 = 1 - г2, где значение 021 определено в (5), (7) или (8).
С помощью идеи вложенных методов оценку ошибки еп,з метода третьего порядка можно оценить по формуле £п,з = Уп+1 - Уп+1,1 = (Р1 - ^1)^1 + (Р2 - ^2)^2 + Рз[3]. Тогда неравенство для контроля точности вычислений имеет вид
II(Р1 - Г1 + (Р2 - Г2)^2 + Рз&з II < е,
где || - | — некоторая норма в Ян, е — требуемая точность интегрирования. В конкретных расчетах применялся метод (2) с коэффициентами (8), как наиболее надежный. Тогда неравенство для контроля точности имеет вид
11|&1 - 2^2 + II < е. 6
3. КОНТРОЛЬ УСТОЙЧИВОСТИ ЧИСЛЕННОЙ СХЕМЫ
Неравенство для контроля устойчивости численной формулы (2) построим предложенным в [3] способом. Запишем стадии кг, к2 и кз применительно к задаче у' = Ау, где А есть матрица с постоянными коэффициентами. В результате получим
кг = Хуп, к2 = (X + в21Х2 )уп, кз = [X + (вз1 + вз2)Х2 + в21вз2 X 3]уп,
где X = НА. Найдем коэффициенты , ¿2 и ¿з из условия выполнения равенства
¿г кг + ¿2 к2 + ¿з кз = X 3уп.
Легко видеть, что данное требование выполняется при следующих значениях:
¿г = взг + вз2 - в2г ¿2 = _ взг + вз2 ¿з = 1
в21 вз2 в21 вз2 ' в21 вз2
Нетрудно видеть также, что
-^(к2 _ кг ) = X2 уп.
в21
Тогда согласно [3] оценку максимального собственного числа г>п,з = НАптах матрицы Якоби системы (1) можно вычислить по формуле
Уп,з = в21 тах (|¿1 к1 + ¿2к2 + ¿зк3|/|к2 _ к11). (9)
1<г<^
Область устойчивости метода третьего порядка приведена на рис. 1.
Интервал устойчивости численной схемы (2) приблизительно равен 2.5. Поэтому для ее контроля устойчивости можно применять неравенство г>п,з < 2.5.
Полученная оценка (9) является грубой, потому что: 1) вовсе не обязательно максимальное собственное число сильно отделено от остальных, 2) в степенном методе применяется мало итераций, 3) дополнительные искажения вносит нелинейность задачи (1). Поэтому здесь контроль устойчивости используется как ограничитель на размер шага интегрирования. Прогнозируемый шаг Нп+г будем вычислять следующим образом. Новый шаг Нас по точности определим по формуле Нас = q1hn, где Нп есть последний успешный шаг интегрирования, а q1, учитывая соотношение еп,з = О(Нп), задается уравнением q3||£n,3|| = е. Шаг Наг по устойчивости зададим формулой = q2Нп, где q2, учитывая равенство г>п,з = 0(Нп), определяется из уравнения q2^п,з = 2.5. В результате прогнозируемый шаг Нп+г вычисляется по следующей формуле:
Нп+г = тах[Нп, тт(Нас, )]. (10)
где Нп есть последний успешный шаг интегрирования. Отметим, что формула (10) применяется для прогноза величины шага интегрирования Нп+г после успешного вычисления решения с предыдущем шагом Нп, и поэтому фактически не приводит к увеличению вычислительных затрат.
Если шаг по устойчивости меньше последнего успешного, то он уменьшен не будет, потому что причиной этого может быть грубость оценки максимального собственного числа. Однако шаг не будет и увеличен, потому что не исключена возможность неустойчивости численной схемы. Если шаг по устойчивости должен быть уменьшен, то в качестве следующего шага будет применяться последний успешный шаг Нп. В результате для выбора шага и предлагается формула (10). Данная формула позволяет стабилизировать поведение шага на участке установления решения, где определяющую роль играет устойчивость. Собственно говоря, именно наличие данного участка существенно ограничивает возможности применения явных методов для решения жестких задач.
Рис. 1. Область устойчивости метода третьего порядка точности
4. МЕТОД ПЕРВОГО ПОРЯДКА ТОЧНОСТИ
Для численного решения задачи (1) рассмотрим схему вида
Уп+1 = Уп + Г1 + Г2 к2 + гзкз, = /(Уп ), = /(Уп + 021&1), кз = /(Уп + вз1 + вз2&2 ),
(11)
где коэффициенты в21, вз1 и вз2 заданы при описании метода третьего порядка точности, а гь г2 и гз подлежат определению. Построим менее точную схему с максимальным интервалом устойчивости. Для этого применим (11) для решения скалярного тестового уравнения у' = Ау,Ле(А) < 0, где Л интерпретируется как некоторое собственное число матрицы Якоби задачи (1). Получим уп+1 = ^(х)уп, где функция устойчивости ф(х) имеет вид
д(х) = 1 + (Г1 + Г 2 + Гз)х + [в21Г2 + (вз1 + вз2 )гз ]х2 + 021 вз2 ГзХз, X = ЛА.
Требование первого порядка точности приводит к соотношению г1 + г2 + гз = 1, которое будем считать выполненным. Теперь выберем г2 и гз таким образом, чтобы метод (11) имел максимальный интервал устойчивости. Для этого рассмотрим многочлен Чебышева Тз(г) = (4гз — 3г) на промежутке [-1,1]. Проведем замену переменных, полагая г = 1 — 2х/7. Получим Тз(х) = 1 — 18х/т + 48х2/72 — 32хз/7з, при этом отрезок [7,0] отображается на отрезок [—1,1]. Нетрудно показать, что среди всех многочленов вида Рз(х) = 1+ х + с2х2 + сзхз для Тз(х) неравенство |Тз(х)| < 1 выполняется на максимальном интервале [7,0], 7 = —18 [3]. Потребуем совпадения коэффициентов ф(х) и Тз(х) при 7 = —18. Это приводит к соотношениям г1 + г2 + гз = 1, 021г2 + (вз1 + вз2)гз = 4/27, 021 вз2гз = 4/729. В результате имеем коэффициенты
4
гз =
729в21вз2'
г2 =
4/27 — (вз1 + вз2 )гз 021 '
г1 = 1 — г2 — гз
метода первого порядка точности с максимальным интервалом устойчивости, локальная ошибка <5п>1 которо-
19
го имеет вид 6п 1 = —- Л2/' / + 0(Лз). Для контро-
54
ля точности численной формулы первого порядка будем использовать оценку локальной ошибки. Учитывая, что к2 — к1 = в21Л2/п/п + 0(Лз) и вид локальной ошибки, неравенство для контроля точности записывается в виде
19
54|в21|
||к2 — к11| < е,
Рис. 2. Область устойчивости метода первого порядка точности
где || ■ || — некоторая норма в Л7 , е — требуемая точность расчетов.
Построим неравенство для контроля устойчивости метода первого порядка. Область устойчивости метода первого порядка приведена на рис. 2.
Интервал устойчивости численной схемы (11) первого порядка точности равен 18 [3]. Поэтому для ее контроля устойчивости можно применять неравенство г>п,з < 18, где г>п,з определяется по формуле (9).
5. АЛГОРИТМ ПЕРЕМЕННОГО ПОРЯДКА И ШАГА
На основе построенных методов первого и третьего порядков точности легко сформулировать алгоритм переменного порядка и шага. Сначала для расчетов используют метод третьего порядка как дающий более точные результаты. Переход на схему первого порядка осуществляется при нарушении неравенства г>п,з < 2.5. Обратный переход на метод третьего порядка происходит в случае выполнения неравенства г>п,з < 2.5. При расчетах по методу первого порядка наряду с точностью контролируется устойчивость, а выбор прогнозируемого шага производится по аналогии с методом третьего порядка точности по формуле (10).
Построенный алгоритм интегрирования будем называть RK3PP(Ac,St). При значениях параметра Ac, равном 1 или 3, расчеты проводятся соответственно методом первого или третьего порядка точности, при Ac =0 используют расчеты с переменным порядком. Если St = 0, то устойчивость контролируется, при St=1 — не контролируется. При расчетах с переменным порядком устойчивость всегда контролируется.
6. РЕЗУЛЬТАТЫ РАСЧЕТОВ
Расчеты проводились на Intel(R) Core(TM)2 Quad CPU с двойной точностью. В конкретных расчетах левая часть неравенства для контроля точности вычисляется по формуле
их и l^iI
||5n|| = max 1 nl
i<i<N |yi | + r'
где i — номер компоненты, r — положительный параметр. Если по i-й компоненте решения выполняется неравенство |уП| < r, то контролируется абсолютная ошибка re, в противном случае — относительная ошибка е. В расчетах параметр r выбирался таким образом, чтобы по всем компонентам решения фактическая точность соответствовала задаваемой. Значение е полагалось равным 10-3. При такой точности наиболее эффективны методы третьего порядка.
Сравнение эффективности алгоритма переменного порядка и шага RK3PP(0,0) проводилось с трех-стадийным методом Рунге - Кутта третьего порядка точности с контролем (алгоритм RK3PP(3,0)) и без контроля (алгоритм RK3PP(3,1)) устойчивости. Далее через isa, iwo и ifu обозначены соответственно суммарное число шагов интегрирования, число повторных вычислений решения (возвратов) вследствие нарушения требуемой точности расчетов и число вычислений правой части задачи (1). Пример 1 [4].
У1 = -0.04yi + 0.01у2Уз, у2 = 400yi - 100у2уз - 3000у2, y3 = 30у2,
(12)
t е [0,40], yi(0) = 1, У2(0) = уз(0) = 0, ho = 10-5.
Алгоритмом переменного порядка и шага RK3PP(0,0) решение задачи (12) удалось вычислить с затратами isa = 6 889, iwo = 124, ifu = 20 792. Для метода третьего порядка RK3PP(3,1) без контроля устойчивости затраты равны isa = 44 441, iwo = 11 758, ifu = 156 839, а для алгоритма RK3PP(3,0) с контролем устойчивости — isa = 44 951, iwo = 655, ifu = 136 163. Пример 2 [4].
у1 = Уз - 100yiУ2, У2 = Уз + 2y4 - 100У1У2 - 20000y2, y3 = -уз + 100yiy2,
(13)
t е [0, 20], yi(0)=0, у2(0)= уз(0)=0, ho = 2.5 ■ 10-5.
Алгоритмом RK3PP(0,0) решение задачи (13) вычислено с затратами isa = 366, iwo = 6, ifu = 1 105. Для метода третьего порядка RK3PP(3,1) без контроля устойчивости затраты равны isa = 2 240, iwo = 555, ifu = 7 830, а для алгоритма RK3PP(3,0) с контролем устойчивости — isa = 1 032, iwo = 20, ifu = 3 136. Пример 3 [4].
yi = -0.013yi - 1000y2уз, У2 = -2500у2уз, уз = -0.013yi - 1000yiуз - 2500у2уз,
5 (14)
t е [0, 50], yi (0) = 0, у2(0) = уз(0) =0, ho = 2.9 ■ 10-5.
Алгоритмом RK3PP(0,0) решение задачи (14) вычислено с затратами isa = 12 689, iwo = 106, ifu = 38 173. Для метода третьего порядка RK3PP(3,1) без контроля устойчивости затраты равны isa = 74 281, iwo = 19 555, ifu = 261 953, а для алгоритма RK3PP(3,0) с контролем устойчивости — isa = 61 047, iwo = 1 686, ifu = 186 513.
Пример 4 [4].
y1 = 77.27(y2 - ViV2 + Vi - 8.375 ■ 10-6у?, у2 = (-y2 - yiy2 + Уз)/77.27, у3 = 0.161(yi - уз),
3 (15)
t е [0,300], yi(0) = 4, y2(0) = 1.1, Уз(0) = 4, ho = 10-3.
Задача (15) является простейшей математической моделью с периодическим решением для описания реакции Белоусова - Жаботинского (орегонатор). Она является «слишком» жесткой для явных методов. Тем не менее пример (15) приведен здесь с целью подчеркнуть преимущество методов с контролем устойчивости, а также алгоритма переменного порядка и шага.
Алгоритмом RK3PP(0,0) решение задачи (15) вычислено с затратами isa = 438 941, iwo = 965, ifu = 1 317 819. Для метода третьего порядка RK3PP(3,1) без контроля устойчивости затраты isa = 2 904 014, iwo = 768 860, ifu = 10 249 762, а для алгоритма RK3PP(3,0) с контролем устойчивости — isa = 2 871 743, iwo = 11 653, ifu = 8 638 535.
Из сравнения результатов расчетов жестких задач следует, что контроль устойчивости всегда приводит к повышению эффективности. Это является следствием устранения некоторых возвратов (повторных вычислений решения), возникающих из-за неустойчивости численной формулы. В конце интервала интегрирования фактическая точность вычислений для всех алгоритмов лучше задаваемой. Такая же тенденция сохраняется при интегрировании всех рассмотренных задач [4-5].
ЗАКЛЮЧЕНИЕ
Использование неравенства для контроля устойчивости фактически не приводит к увеличению вычислительных затрат, потому что оценка максимального собственного числа матрицы Якоби системы (1) осуществляется через ранее вычисленные стадии и не приводит к росту числа вычислений функции f. Такая оценка получается грубой. Однако применение контроля устойчивости в качестве ограничителя на рост шага позволяет избежать негативных последствий грубости оценки. Более того, в некоторых случаях это приводит к нестандартно высокому повышению эффективности алгоритма. На участке установления за счет контроля устойчивости старые ошибки стремятся к нулю, а новые невелики за счет малости производных решения. В некоторых случаях вместо оценки максимального собственного числа оценивается следующее по порядку. Шаг интегрирования становится больше максимально допустимого и с таким шагом осуществляется интегрирование до тех пор, пока не нарушается неравенство для контроля точности. Как правило, число таких шагов невелико. Однако величина шага может на порядок превышать максимальный шаг по устойчивости. После нарушения неравенства для контроля точности шаг уменьшается до максимально возможного. Такой эффект может повторяться многократно в зависимости от длины участка установления. В результате средний шаг интегрирования может превышать максимально допустимый.
Работа выполнена при финансовой поддержке РФФИ (проекты 11-01-00106 и 11-01-00224).
Библиографический список
1. Shampine L.M. Implementation of Rosenbrock method // ACM Transaction on Mathematical Software. 1982. Vol. 8, № 5. P. 93-113.
2. Новиков Е.А., Новиков В. А. Контроль устойчивости явных одношаговых методов интегрирования обыкновенных дифференциальных уравнений // Докл. АН СССР. 1984. Т. 277, № 5. С. 1058-1062.
3. Новиков Е.А. Явные методы для жестких систем. Новосибирск: Наука, 1997. 197 с.
4. Eright W. H., Hull T. E. Comparing numerical methods for the solutions of systems of ODE's // BIT. 1975. № 15. P. 10-48.
5. Хайрер Э., Ваннер Г. Решение обыкновенных дифференциальных уравнений. Жесткие и дифференциально-алгебраические задачи. М.: Мир, 1990. 685 с.
