CC BY
48Анализ устойчивости явных методов Рунге - Кутты в инструментальной среде ИСМА
Ю. В. Шорников, Е. А. Новиков, Д. Н. Достовалов
Новосибирский государственный технический университет, 630092, Новосибирск, Россия
Разработаны и программно реализованы алгоритмы исследования устойчивости явных методов Рунге - Кутты. Исследованы 13-стадийные методы различных порядков точности на основе стадий Фельберга. Доказана эффективность алгоритмов переменного порядка с контролем точности и устойчивости.
Ключевые слова: явные методы, жесткость, функция устойчивости, интервал устойчивости.
Algorithms of research a stability of explicit methods Runge-Kutta are developed and programmed. Are researched 13-phasic methods of various orders of a precision on the basis of Fehlberg stages. Efficiency of variable-order algorithms with the control of precision and stability is proved.
Keywords: explicit methods, stickiness, function of stability, range of stability.
Введение. В ряде случаев численное решение задачи Коши для жестких систем обыкновенных дифференциальных уравнений целесообразно проводить с использованием явных методов. Как правило, алгоритмы интегрирования на основе неявных или полуявных формул используют обращение матрицы Якоби, что в случае большой размерности системы является отдельной трудоемкой задачей. Однако многие современные алгоритмы на основе явных методов с контролем точности не приспособлены для решения жестких задач, поскольку такой подход приводит к потере эффективности и надежности вследствие обусловленного потерей устойчивости "раскачивания" шага интегрирования на участке установления. Введение в алгоритм интегрирования контроля устойчивости позволяет существенно повысить эффективность расчетов. Для этого необходимо знать характеристики вычислительной устойчивости используемой численной схемы, в частности длину интервала устойчивости. Как правило, процесс исследования устойчивости является достаточно трудоемким, поэтому предложены, программно реализованы и интегрированы в систему ИСМА [1] алгоритмы автоматизации решения данной задачи. Результаты исследования могут быть применены для создания алгоритмов интегрирования переменного порядка, использующих известные численные схемы и построенные на их основе методы с расширенным интервалом устойчивости.
Специфика решения жестких задач. Поведение различных динамических систем удобно описывать с использованием системы обыкновенных дифференциальных уравнений
У = f (t, y), У (to) = У о, (1)
где f : R x Rn ^ Rn - известная непрерывная вектор-функция; RN - N -мерное линейное вещественное пространство; Уо - решение в начальной точке to. Как правило, системы вида (1) являются жесткими. Задача Коши для системы обыкновенных дифференциальных уравнений называется жесткой на некотором интервале I ^ [tо,tk], если для t е I выполняются соотношения Re(Ät) < 0 (1< i < N),
max Re(Ai) / min Re(A,) >> 1 , где Ä, - собственные числа матрицы Якоби df / dy , вычисленной на
1<i<N 1 1<i<N 1
решении y(t) [2]. На практике такие задачи распознаются хорошо. Любая физическая система, имеющая компоненты с существенно различающимися временными константами, приводит к жесткой задаче.
В жестких задачах длина интервала интегрирования связана с медленно меняющимся решением, в то время как существуют быстро затухающие возмущения. Поэтому весь интервал интегрирования можно условно разбить на несколько участков. Для некоторых из них, называемых переходными, характерно наличие больших производных решения при малой длине. Поскольку на переходном участке решение меняется быстро, из условий точности шаг интегрирования выбирается небольшим. На других участках, называемых
участками установления, производные решения малы, а их длина существенно больше длины переходного участка. В этом случае шаг интегрирования по точности будет увеличен, и при отсутствии ограничения на максимальный размер он может стать настолько большим, что рост ошибок вследствие возникающей неустойчивости вызовет переполнение арифметического устройства ЭВМ. Это приводит к необходимости большого количества повторных вычислений решения, в результате шаг выбирается значительно меньше максимально допустимого. Таким образом, применение известных явных схем, обладающих малыми областями устойчивости, для решения жестких задач становится неэффективным или вообще невозможным. Однако введение в алгоритм интегрирования контроля устойчивости позволяет расширить границы применимости явных методов.
Методы Рунге - Кутты. Для решения задачи (1) будем использовать явные методы типа методов Рун-ге - Кутты [2], которые имеют вид
т '—1
Уп+1 = Уп + Е М- , к = У ^п + уп + £ вЧк1), а = 0 (2)
'=1 7=1
(а , в у , Р1, 1 —' — т , 1 — у — ' — 1 - числовые коэффициенты; к - стадии метода; т - число стадий).
Численные формулы (2) хорошо приспособлены для практических расчетов, поскольку не требуют вычисления дополнительных начальных значений и позволяют легко менять шаг интегрирования. Свойства методов определяются набором параметров т , а , в у , р{. Как правило, при каждом заданном т , т. е. при
заданных затратах на шаг интегрирования, коэффициенты а , в у , Р1 выбираются из условия максимального порядка точности численной формулы (2). В ряде случаев часть параметров используется для контроля точности вычислений или расширения области устойчивости метода.
Устойчивость. Поскольку исходная задача обладает малой чувствительностью к малым возмущениям, того же свойства требуют и от метода ее решения. Такие методы называют устойчивыми. Далее будем рассматривать абсолютную устойчивость метода.
Понятие абсолютной устойчивости вводится на линейном скалярном уравнении [2]
у' = Лу, у(0) = у0, г>0, (3)
с комплексным X , Re(X) < 0 .
Точное решение задачи (3) в точке 1п+х записывается в виде у (¿п+1) = е2 у ($п) , г = ИХ , при этом выполняется неравенство | е2 |— 1. В результате применения метода (2) для решения задачи (3) получается уп+1 = Q(г) уп , где Q(2) - функция роста или функция устойчивости метода (2). Функция устойчивости т -стадийного метода (2) представляет собой многочлен степени т [2]:
Q(z) = 1 + Су2 + ... + Ст2т. (4)
Очевидно, что при численном решении задачи (3) также должно быть истинным неравенство | Q(2) |— 1. Метод (2), удовлетворяющий этому условию для данного 2 , называется абсолютно устойчивым. Область Я комплексной плоскости Z называется областью устойчивости метода, если он устойчив при всех 2 е Я . Далее будем использовать понятие интервала устойчивости - пересечение области устойчивости с вещественной осью [2] комплексной плоскости Z .
Поскольку при 2 ^ —со | Q(2) с , области устойчивости явных методов типа методов Рунге - Кут-ты ограничены. Это означает, что шаг интегрирования по устойчивости ограничен неравенством И — П/ | Хтах |, где постоянная П - длина интервала устойчивости; Хтах - максимальное собственное
число матрицы Якоби задачи (1).
Известно [3], что длина интервала устойчивости пропорциональна квадрату числа стадий и обратно пропорциональна порядку точности метода. Так, для т -стадийного метода в случае первого порядка точности
величина П максимальна и равна 2т2.
Коэффициенты функции роста можно вычислить [2] исходя из параметров численного метода. Для этого введем в рассмотрение матрицу В с элементами
Ъ1г = 1, 1 < г < т, Ъы = 0, 2 < к < т , 1 < г < к -1,
1-1
К = ТвуЪк-и , 2 < к < т , к < г < т ,
} =к-1
Т
и
где в у , 2 < г < т (1 < ] < г -1) определены в выражении (2). Обозначим С = (с^---, ст) Р = (рх, —, рт)Т . Тогда имеет место соотношение
ВР = С , (5)
используя которое можно решить обратную задачу - поиск коэффициентов рг (1 < г < т) при заданных ву (2 <г< т , 1 < у <г-1) и сг (1 <г< т).
Таким образом, с помощью выражения (5) можно найти полином устойчивости любого известного метода, затем определить длину интервала устойчивости и расширить границы применимости метода, создав алгоритм интегрирования с контролем устойчивости. С другой стороны, построив область устойчивости нужной конфигурации и получив функцию роста, можно создать численный метод с заданным интервалом устойчивости. В частности, возможно создание методов с расширенными областями устойчивости на основе стадий известных численных схем.
Графическая интерпретация полинома устойчивости. Полином (4) является функцией комплексной переменной, поэтому на плоскости он представляется в виде линий уровня, для которых )| = d при
й = 0,05; 0,30; 0,70; 1,00 (рис. 1). Часть комплексной плоскости, охватываемая линией уровня (г)| = 1,
является областью устойчивости метода.
Область устойчивости, представленная на рис. 1, является "почти" многосвязной. Многосвязная, т. е. разрывная, область устойчивости приводит к сокращению интервала устойчивости [2], который уменьшится
до наиболее близкой к нулю вещественной точки X , где (X)| = 1. В данном случае это точка X = 4,0 .
Сокращение области устойчивости обусловлено следующими причинами. Если при увеличении шага интегрирования происходит переход из ближайшего к нулю участка области, то вероятность оказаться в области неустойчивости высока по сравнению с возможностью попадания в следующую часть области устойчивости. Также в случае, если собственные числа матрицы Якоби расположены достаточно близко, часть их может оказаться в области неустойчивости. Поэтому для получения достаточно большого интервала устойчивости полином устойчивости в экстремальных точках должен принимать значения меньше единицы, например 0,9 (рис. 2).
Рис. 1. Графическая интерпретация полинома устойчивости при т = 3 , к = 2
Рис. 2. Немногосвязная область устойчивости при т = 3 , к = 2
Рис. 3. Многосвязная область устойчивости при m = k = 6
Следует отметить, что при повышении порядка точности область устойчивости становится многосвязной. При m = k = 6 отделяется пара комплексно-сопряженных корней с положительной вещественной частью (рис. 3), а при k = 11 отделяется вторая пара корней. Предполагается, что этот процесс носит регулярный характер.
Расчет коэффициентов функции устойчивости. Для
того чтобы построить m -стадийный метод k -го порядка с расширенной областью устойчивости, необходимо найти соответствующую ему функцию роста, удовлетворяющую
условию Q (z)| < 1 на максимальном вещественном интервале. Для этого одновременно рассмотрим полином Чебышева первого рода T (z') того же порядка m , для которого справедливым является неравенство |T (z')| < 1 при |z '|< 1 [2].
Если выполнить замену переменных z' = 1 — 2z / у , то очевидно, что при Im (z ) = 0 |T (z )|< 1 на интервале R6 (z) e [y,0] . При этом замена переменных тождественна отображению отрезка [—1,1] в отрезок [у, 0]. Приравнивая коэффициенты правой части (4) с одинаковыми степенями к соответствующим коэффициентам полинома Чебышева, получаем систему уравнений для определения коэффициентов функции устойчивости. Учитывая, что при замене переменных T(z' ) = T(z,у) (у = const), всегда можно подо-
Y ,0 будет являться
брать такое значение у = у , при котором |Т (z)| < 1, и вещественный интервал расширенным интервалом устойчивости. Однако, как отмечено ранее, при таком подходе область устойчивости получается "почти" многосвязной, поскольку в вещественных экстремальных точках |Т(z)| = 1. Поэтому будем управлять конфигурацией области устойчивости, задавая значения полинома (4) в экстремальных точках [2]. Тогда, если известно примерное расположение собственных чисел матрицы Якоби решаемой задачи, можно построить область устойчивости соответствующей конфигурации. Например, предполагается, что в окрестности нуля имеется множество комплексных собственных чисел и соответственно корней полинома устойчивости. Кроме того, имеются отрицательные вещественные корни, относительно большие по модулю. Следовательно, область устойчивости должна иметь вид "лампочки", широкая часть которой направлена к нулю (см. рис. 1).
При сравнении (4) с разложением точного решения задачи (3) в числовой ряд видно, что если схема (2) имеет порядок точности к , то первые к коэффициентов полинома Q(z) определяются по формуле
С 1 = 1/1! (1 < I < к). Для отыскания остальных коэффициентов С 1 (к +1 < / < т) обозначим вещественные экстремумы (4) через х^...,Хт_1, причем Х1 > Х2 > ... > хт—1, и потребуем, чтобы в экстремальных точках х1 (к < / < т — 1) многочлен (4) принимал заданные значения ^. Получаем систему нелинейных уравнений
^ (х ) = рг,
[Q'(х.) = 0, к <г< т — 1, ()
где Q' (x ) = 2ic г-
г—1
г=1
Запишем задачу (6) в матричном виде. Обозначим через У , Z , С и Я векторы с координатами
yi = xk+i-1 - Z = ck+i - g, = Fk+i-1 - 1 - Z cjyi - Г = -Z jcj^ j-1 - 1 —i —N - N = m - k -
j=1 j=1
через £i - E2 - диагональные матрицы с элементами e1 = k + i - e2i = VУ, - 1 — i — N - а через A матрицу с элементами aiJ = yk+ j - 1 — i- j — N . Тогда задача (6) принимает вид
Г AZ - G = 0-[E2AE1Z - Я = 0.
В^1разив Z из первого уравнения и подставив во второе- получаем E 2 AE1A G - R = 0 . Для численного
решения данной задачи используем метод установления- состоящий в том- что для стационарной задачи строится некоторый нестационарный процесс- решение которого с течением времени сходится к решению исходной задачи:
У = f (У) - У (0 ) = Го- (7)
где f (у) = E3 (E2AE1 A^1G -R) ; E3 - диагональная матрица с элементами e'3 = (-1)k+i 1 (1 — i — N).
Максимальный размер интервала устойчивости [2] явных m -стадийных методов типа методов Рунге -Кутты равен 2m2 - поэтому в качестве начальных условий данной задачи выберем экстремальные точки многочлена Чебышева- рассматриваемого на отрезке I -2m2- 0 I . С учетом смещения их можно вычислить
[2] по формуле yt = m21 cos—-1 I - 1 — i — m -1. V m J
Задача (7) является жесткой- причем жесткость усиливается с ростом размерности. Поэтому для решения задачи (7) использован алгоритм МК22 [4] на основе L -устойчивого (2-2)-метода. Особенностью алгоритма МК22 является применение одной матрицы Якоби на нескольких шагах интегрирования- что существенно повышает эффективность алгоритма при решении задач большой размерности.
Вычисление длины интервала устойчивости. Поскольку интервал устойчивости явных методов типа методов Рунге - Кутты не превышает 2m2 - величина -D лежит на участке \_-2m2- 0|. Из определения абсолютной устойчивости следует- что Q (x)| — 1 при x е [-D- 0]. Поэтому на интервале
x e_-2m2- 0] функция Q(x)| -1 один раз меняет знак в точке -D . С учетом этого данную задачу будем решать методом дихотомии. Введем обозначения
f (x) = \Q (x)|-1; (8)
a = -2m2- b = 0. (9)
Вычисляя значения функции f (x) в точке a и пробной точке
c = ( a + b У 2- (10)
уточним промежуток существования корня. Процесс уточнения повторяется до тех пор- пока не будет достигнута требуемая точность S или значение функции f (c) не будет малым положительным числом (не больше заданного допуска 5). Данная схема реализуется следующим алгоритмом. Шаг 1. Вычисляются a и b по формулам (9) и f (a) по формуле (8). Шаг 2. Вычисляется c по формуле (10).
Шаг 3. Если (b — a) < 2s , расчет окончен, D — |с|.
Шаг 4. Вычисляется f (с) по формуле (8).
Шаг 5. Если 0 < f (с) < 5 , расчет окончен, D — |с|.
Шаг 6. Если f (a) f (с) < 0, положить b — с . Иначе a — си f (a) = f (с).
Шаг 7. Перейти на шаг 2.
Программа исследования устойчивости. На основе рассмотренного выше математического обеспечения создано программное средство анализа устойчивости явных численных схем, которое позволяет автоматизировать выполнение следующих задач: 1) расчет коэффициентов полинома устойчивости по заданным параметрам численной схемы в соответствии с (5); 2) построение области устойчивости заданной конфигурации, поиск соответствующей функции роста и расчет параметров численного метода с заданной областью устойчивости; 3) графическая интерпретация функции устойчивости (построение области устойчивости); 4) расчет длины интервала устойчивости; 5) создание отчетов об исследовании в формате Microsoft Word или Excel.
Таким образом, использование созданного инструментального средства позволяет существенно сократить временные затраты на создание эффективных алгоритмов численного анализа на основе явных схем.
Новый комплекс программ исследования устойчивости интегрирован в среду моделирования динамических и гибридных систем ИСМА [1] и зарегистрирован в Роспатенте [5]. Правильность работы комплекса доказана конструктивно на всем множестве тестовых примеров из [2]. Наряду с малостадийными методами разного порядка точности успешно исследованы многостадийные схемы Фельберга, в том числе с расширенной областью устойчивости [6, 7], которые будут рассмотрены ниже.
Анализ устойчивости алгоритмов на основе стадий Фельберга. Для иллюстрации практического применения исследования устойчивости рассмотрим схемы на основе стадий Рунге - Кутты - Фельбер-га [8]:
13 1
Уп+i = УП+Z p k >, k = hf (tn+УП +Z j) (11)
'=1 j=1
с коэффициентами
= 0 = _2_ = 1 = 1 = _5_ = 1 = 5 = 1 = 2 = 1 =.
ОСл — 0 , ОС о — , ОС о — , ОС л — , ОС с — , ОСп — , ОС 7 — , ОС о — , ОС(\ — , ОСл п — , ОСл 1 — 1,
1 2 27 3 9 4 6 5 12 6 2 7 6 8 6 9 3 10 3 11
а12 — 0, а13 — 1, ß21 — 2 27 , ß31 — 1 36 , ß32 — 12 , ß« = _ 1 " 24 , ß 42 — 0 , ß43 — 1 — 8 , ß51 — 5 — 12 , ßs2
ßs3 — 25 — 16 , ß — ß54 —16, ß61 — 1 — 20 , ße2 — ßea — 0, ß64 — 1 — 4 , ß65 — 5 ß71 5 — — 25 108 , ß72 — ß73
ß74 — 125 108 , ß —--ß75 — 27 , ß — 125 ß76 — 54 , ß81 — 31 300 , ß82 — ß83 — ß84 !— 0, ß85 — 61 ß — — 225 , ß86 —
ß87 — 13 900 , ß91 — 2 , ß92 — ß93 — 0 , ß94 —* 53 6 , ß — 704 ß P95 - 45 , Pge — 107 9 , ß — - ß ß97 — 90 , 8
ß10,1 = 91 108 , ß10,2 — ß10,3 — 0 , ß10,4 — 23 108 , ß10,5 — -976 ß — _ 1 oc , ^10,6 _ 135 331 ß 54 , ß107 19 ^b^ ß108 —
2 9
17
6 :
1 2383 341 4496 301
До,9 — — ^ ßn,i— ßn,2 — ßu — 0, Au ßn,5 — ßn,6 — — "82,
2133 45 45 18 3
ß117 — 4100 , ß118 — 82 , ß11,9 — 164 , ß1110 — "47, ß121 — 205 , ß12,2 — ß12 3 — ß12,4 — ß12,5 — 0
в — -— в — в — -— в = — в = — в = 0 в --1777
А2,е- 41, А2,7 - 205 , А2,8- 41, - 41, А2,10 - 41, А2Д1-0, Азд - 4100 ,
341 4496 289 2193 51
в13 2 =в13 3 = ^ в13 4 ^^ в13 5 = 1025 , в13 6 =-в13 7 = в138 = 82 ,
33 12
в13,9 = 164 , в13,10 = 41 , в13,11 = 0 , в13,12 = 1 . При значениях коэффициентов
41 34 9
Р1 — ^ , Р2 = Р3 = Р4 = Р5 = 0, Рб , Р7 = Р8 =^7 , 840 105 35
9 41 0
Р9 = Р10 = 280, Рп = 840 , Р12 = Р13 =0
(12)
схема (11) имеет седьмой порядок точности [8], а соответствующий полином устойчивости Q (г) принима-
11
ет вид Q(г) — 1 + ^с^ , где параметры с{ находим из соотношения (5):
с — 1, с2 — 0,5, с7 — 0,19841269841270 • 10-3
г—1
с2 — и> - , с7
с3 — 0,16666666666667, с8 — 0,23165371472663 • 10-4, с4 — 0,41666666666667 •10-1, — 0,23671439526314 • 10-5,
с5 — °,83333333333333•10- , с10 — 0,51829448771964• 10-7,
с6 — 0,13888888888889 •10-2, с11 — -0,43191207309970 • 10-7.
Область устойчивости данного метода приведена на рис. 4, на котором видно, что интервал устойчивости приблизительно равен пяти.
Согласно [8] численная формула (11) с коэффициентами
34 9 9
Р1 — Р 2 — Р3 — Р 4 — Р5 — 0 Р6 —, Р7 — Р8 — , Р9 — Р10 —777 ,
105 35 280 (13)
0 41
Р11 — 0, Р12 — Р13 — —
12
этом
г=1
циенты принимают значения
имеет восьмой порядок точности. В этом случае функция роста имеет вид Q (г) — 1 + ^ с^' , где коэффи-
с — 1, с2 — 0,5, с8 — 0,24801587301587 • 10-4,
с3 — 0,16666666666667, с9 — 0,23490700935724 • 10-5,
с4 — 0,41666666666667 •10-1, с10 — 0,23620053064283 • 10-6,
с5 — 0,83333333333333 •10-2, с11 — -0,25914724385982 • 10-7,
с6 — 0,13888888888889 •10-2, с12 — -0,14397069103323 • 10-7.
с7 — 0,19841269841270 • 10-3,
Область устойчивости приведена на рис. 5, на котором видно, что интервал устойчивости метода восьмого порядка приблизительно равен пяти.
Рис. 4. Область устойчивости метода (11) с коэффициентами (12)
Рис. 5. Область устойчивости метода (11), (13)
На основе стадий численной схемы (11) построим метод первого порядка точности с более широкой областью устойчивости. При т — 13 максимальный интервал устойчивости схемы Фельберга равен 338 [8]. Так как размер области устойчивости по вещественной оси методов седьмого и восьмого порядков равен пяти, то для задач, в которых шаг в основном ограничен по устойчивости, эффективность можно повысить примерно в 67 раз. Однако из результатов расчетов даже простейших задач следует, что шаг интегрирования существенно меньше максимально возможного, поскольку области устойчивости внутренних схем (11) малы, и при больших шагах это приводит к существенным погрешностям в промежуточных вычислениях.
Из результатов численных экспериментов следует [7], что наиболее эффективным является метод первого порядка (11) с коэффициентами при т — 7 . Тогда максимальный интервал устойчивости равен 98,0, а полином Q (г) имеет коэффициенты
с1 — 1, с2 — 0,16326530612245
с5 — 0,43614440711587•10-
-2
с3 — 0,99958350687216 • 10
с4 — 0,29142376293650•10-
сй — 0,32366931882441•10
с7 — 0,94364232893419 • 10
-7
10
Используя соотношение (5), получаем:
Рх — -0,36275196357533, Р2 — +0,19206362739513 , Р3 — +0,84515155604427, Р4 — +0,32207503966499 , р5 — +0,33295308199673 • 10-2, Р6 — +0,13204058892796 • 10-3, р7 — +0,16906205162854 • 10-6.
Область устойчивости данного метода показана на рис. 6. Приведем еще один набор коэффициентов:
(14)
с1 — 1, с2 — 0,17242757067512 с, — 0,12240089243768 • 10-1,
с, — 0,34986371648830•10
-3
с5 — 0,55971678758885• 10
с6 — 0,44431955012931•10-с7 — 0,13862209013567•10
-5
-9
(15)
при которых в вещественных экстремальных точках х{ (1 < / < 6) многочлена устойчивости выполняется условие | Q(х{) | = 0,9 . Подставляя
данные значения в соотношение (5), получаем коэффициенты (15) метода (11) первого порядка точности с интервалом устойчивости, равным приблизительно 90.
Рис. 6. Область устойчивости метода (11), (14)
рх = -0,38057403389775, р2 = +0,11240089243768 • 10-1, р3 = +0,89767190586740,
р4 = +0,38213980246925 , р5 = +0,42473101431975 • 10-2, р6 =+0,18104880747501 -10-3, р7 =+0,24835399856329 -10-6.
Область устойчивости построенного метода первого порядка точности (рис. 7) примерно в 18 раз шире по вещественной оси области устойчивости численной схемы (11) седьмого и восьмого порядков. Кроме того, по числу вычислений правой части задачи (1) метод первого порядка почти в два раза дешевле. Поэтому для задач, в которых шаг ограничен в основном по устойчивости, предполагается теоретическое повышение эффективности в 36 раз.
Алгоритм переменного порядка и шага. Поскольку максимальный интервал устойчивости для явных схем (2) пропорционален квадрату количества т стадий метода и обратно пропорционален порядку точности к, для расширения интервала устойчивости целесообразно использовать многостадийные численные схемы. На основе рассмотренных выше методов (11)—(15) построен и включен в библиотеку численных методов системы ИСМА [1] алгоритм переменного порядка и шага с контролем точности и устойчивости ККР78БТР [6], высокая эффективность которого достигается применением наиболее выгодной численной схемы. На участке установления, где шаг выбирается из условия устойчивости, интегрирование проводится методом первого порядка (11), (15). На переходном участке для выполнения требований точности используется метод седьмого порядка с коэффициентами (12) совместно с методом восьмого порядка (13) для контроля точности.
Для контроля устойчивости применяется неравенство ип < О , где ип - оценка максимального собственного числа матрицы Якоби; О - длина интервала устойчивости используемого на текущем шаге метода. Величина Оп вычисляется с использованием степенного метода по формуле
ип = тах | (12к3 - 18к2 + 6к,), | /1 (к2 - к,).
1<г<#
ное собственное число не обязательно существенно отделено от остальных; 2) в степенном методе применяется мало итераций; 3) дополнительные искажения вносит нелинейность задачи (1). Поэтому здесь контроль устойчивости используется в качестве ограничителя на размер шага интегрирования, откуда шаг по устойчивости И^ выбирается по формуле = г И, где г вычисляется из равенства го = О.
Такая оценка является грубой, поскольку: 1) максималь-
Рис. 7. Область устойчивости метода (11), (15)
Сравнение результатов расчетов с контролем устойчивости
RKF78 RKF78ST RKF78STP
944 229 497 836 18 394
564 907 305 709 59 630
244 898 131 827 12 696
191 313 59 281 19 313
28 859 16 454 4 403
Для конструктивного доказательства эффективности ККР78БТР рассмотрены задачи различной степени жесткости [6]. Результаты расчетов приведены в таблице, в которой в качестве критерия эффективности выбрано суммарное число вычислений правой части задачи (1) на интервале интегрирования. Расчеты проводились с задаваемой точностью 5 = 10 6 с использованием трех алгоритмов: классического метода Рунге - Кутты - Фельберга седьмого порядка с контролем точности (ЫКР78), алгоритма с контролем устойчивости (ККР78БТ) и алгоритма переменного порядка и шага с контролем точности и устойчивости (ККР78БТР).
Из таблицы следует, что контроль устойчивости приводит к повышению эффективности расчетов почти в 2 раза, а поскольку оценка максимального собственного числа матрицы Якоби системы (1) осуществляется через ранее вычисленные стадии и не требует расчета функции / , вычислительные затраты увеличиваются незначительно. Применение на участке установления методов низкого порядка точности с расширенными областями устойчивости позволяет значительно увеличить размер шага интегрирования без увеличения вычислительных затрат. На переходных участках, где определяющую роль играет точность вычислений, целесообразно применять методы более высокого порядка точности, но с небольшой областью устойчивости. Комбинирование методов низкого и высокого порядков с помощью неравенства для контроля устойчивости позволяет существенно повысить эффективность расчетов.
Заключение. Используя разработанную программу, нетрудно проанализировать устойчивость численной схемы или создать собственный метод, имеющий требуемые характеристики устойчивости. Графическая интерпретация функции устойчивости и областей устойчивости не обременяет пользователя дополнительными затратами по обработке данных, что позволяет быстро создавать новые алгоритмы интегрирования. Результаты тестирования подтверждают высокую эффективность алгоритма переменного порядка с контролем устойчивости и многостадийных явных схем для исследования задач повышенной жесткости.
Список литературы
1. Свидетельство об офиц. регистрации программы для ЭВМ № 2005610126. 2005. Инструментальные средства машин-
ного анализа / Ю. В. Шорников, В. С. Дружинин, Н. А. Макаров, К. В. Омельченко, И. Н. Томилов // Программы для ЭВМ. Базы данных. Типологии интегральных микросхем: Офиц. бюл. Рос. агентства по патентам и товар. знакам. 2005.
2. Новиков Е. А. Явные методы для жестких систем. Новосибирск: Наука. Сиб. предпр. РАН, 1997.
3. Новиков Е. А., Новиков В. А. Контроль устойчивости явных одношаговых методов интегрирования обыкновенных
дифференциальных уравнений // Докл. АН СССР. 1984. Т. 277, № 5. С. 1058-1062.
4. Новиков А. Е., Новиков Е. А., Шорников Ю. В. Аппроксимация матрицы Якоби в (2,2)-методе решения жестких
систем // Докл. АН ВШ РФ. 2008. № 1. С. 30-44.
5. Свидетельство о гос. регистрации программы для ЭВМ № 2009610905. Алгоритм численного конструирования облас-
тей устойчивости одношаговых методов / Ю. В. Шорников, Е. А. Новиков, Д. Н. Достовалов // Программы для ЭВМ. Базы данных: Типологии интегральных микросхем: Офиц. бюл. Рос. агентства по патентам и товар. знакам. 2009. № 2.
6. Новиков Е. А., Шорников Ю. В., Никонова О. В. Алгоритм переменного порядка и шага на основе стадий метода
Фельберга седьмого порядка точности // Науч. вестн. Новосиб. гос. техн. ун-та, 2006. № 4. С.105-118.
7. Новиков Е. А., Шорников Ю. В. Контроль устойчивости метода Фельберга седьмого порядка точности // Вычисл.
технологии. 2006. Т. 11. № 4. С. 65-72.
8. Fehlberg E. Classical fifth-, sixth-, seventh- and eighth order Runge - Kutta formulas with step size control // Computing.
1969. V. 4. P. 93-100.
Юрий Владимирович Шорников - канд. техн. наук, доц. кафедры автоматизированных систем управления Новосибирского гос. техн. университета, (383) 346-11-00,
e-mail: shornikov@asu.cs.nstu.ru, shornikov@inbox.ru.
