Внутренняя устойчивость явных методов Рунге-Кутта Текст научной статьи по специальности «Математика»

Литература

1. Фергин, В.Р. Методы оптимизации в лесопильно-деревообрабатывающем производстве / В.Р. Фергин. -М.: Лесн. пром-сть, 1975. - 216 с.

2. Розен, В.В. Цель - оптимальность - решение: математические модели принятия решений / В.В. Розен. -М.: Радио и связь, 1982. - 168 с.

3. Гуслицер, И.И. Основы эффективного функционирования поточных линий по первичной обработке древесного сырья: моногр. / И.И. Гуслицер, В.П. Шмаков, А.П. Меньшиков. - Красноярск: Изд-во КГУ, 1993. -184 с.

4. Каверзин, С.В. Обеспечение работоспособности гидравлического привода при низких температурах / С.В. Каверзин, В.П. Лебедев, Е.А. Сорокин. - Красноярск, 1998. - 240 с.

5. Бутылкин, В.И. Накопители для лесных грузов. Создание и совершенствование / В.И. Бутылкин, В.М. У пит // Проблемы машиностроения и новые материалы: мат-лы Всерос. науч. конф. - Красноярск: ИПЦ КГТУ, 2006. - С. 171-175.

----------♦'------------

УДК 519.622 Е.А. Новиков, Ю.А. Шитов

ВНУТРЕННЯЯ УСТОЙЧИВОСТЬ ЯВНЫХ МЕТОДОВ РУНГЕ-КУТТА1

Для произвольного т получены коэффициенты явных т -стадийных методов типа Рунге-Кутта с первого по третий порядок точности, у которых области устойчивости промежуточных численных формул согласованы с областью устойчивости основной схемы. Построены неравенства для контроля точности и устойчивости.

Введение. Для численного решения задачи Коши для жестких систем обыкновенных дифференциальных уравнений

У = /(*,У\ Ж) = Уо. *о <г < Ч, (1)

в [1] предлагается применять явные методы типа Рунге-Кутта

т г-1

уи+1 = у» + Е Рт&, К = ¥ (г» + аА Уп+Е })> (2)

г=1 }=1

где у и I - гладкие вещественные N -мерные вектор функции; г - независимая переменная;

кг,1 < г < т, - стадии метода; ц, Рт ,1 < г < т, 1 < } < г -1 - коэффициенты, определяющие свой-

ства точности и устойчивости схемы (2). В [1] для произвольного т получены коэффициенты методов с первого по третий порядок, а в [2] - коэффициенты методов четвертого порядка с максимальным интервалом устойчивости. Там же численно показано значительное повышение эффективности алгоритмов интегрирования за счет комбинирования таких численных формул в процессе вычислений на основании критерия

устойчивости. Заметим, что в [1] и [2] не рассмотрен вопрос о выборе коэффициентов , которые влияют

на устойчивость промежуточных численных схем и, в конечном счете, на эффективность алгоритма интегрирования. Авторы ограничились замечанием, что устойчивости промежуточных формул можно добиться за

счет выбора данных коэффициентов достаточно малыми.

Численные схемы. Для упрощения выкладок далее рассмотрим задачу Коши для автономной системы обыкновенных дифференциальных уравнений

У=1 (УХ У(го) = Уо> го <г < гк> (3)

для решения которой применим методы типа Рунге-Кутта, записанные в следующем виде:

г т

Упг = Уп + Е Д+1,А * 1 < 1 < т - 1 Уп+1 = Уп + Е Ртк , (4)

}=1 г=1

1 Работа выполнена при финансовой поддержке РФФИ (грант №05-01-00579-а).

где кг = к/ (уп г-1), 1 < I < т, уп0 = уп. Все полученные ниже результаты можно использовать для неавтономных систем, если в (2) положить

1-1

а = Е вц, 2 < г < т, а1 = 0. (5)

1=1

В дальнейшем нам потребуется матрица Вт с элементами Ъц вида [1]

Ъ1г = 1, 1 < г < т, Ъкг = 0, 2 < к < т,1 < г < к -1,

1-1

Ъкг = Е вЪк-1,1. 2 < к < т к <г < т (6)

1=к-1

где вц - коэффициенты схемы (2) или (4). Устойчивость одношаговых методов обычно исследуется на линейном скалярном уравнении

у' = Лу, у(0) = Уо, г > 0, (7)

с комплексным Л,Яв(Л) < 0 . Применяя вторую формулу (4) к (7), получим

уп+1 = йт ( ^) уп . 2 = кЛ,

тт йт (2) = 1 + Е Ст£ , Ст1 = Е ЪцРт1 ,1 < 1 < т (8)

г=1 1=

Обозначая Ст = (ст1,•••, стт )Т, Рт = (рт1,..., ртт )Т, соотношения (8) можно переписать в виде

В Р = С , (9)

т т т9 ' >

где элементы матрицы Вт определены соотношениями (6). Для промежуточных численных схем (4) имеем

кк

у,к = а„ (г) у,, (г) = 1+Е V, Сц = Е Ъ„вы,, ,1 < к <т-1. (10)

г =1 ¡=г

Используя обозначения вк = (вк+11,..,вк+1к)Т,ск = (ск1,..,скк)Т, получим, что коэффициенты в у промежуточных схем (4) и коэффициенты соответствующих многочленов устойчивости связаны соотношениями

Вк вк = ск, 1 < к < т -1. (11)

Заметим, что из сравнения (6) и (10) следует, что Ъкг = с-1к-1, то есть элементы (к +1) -го столбца

матрицы Вт совпадают с коэффициентами многочлена устойчивости йк (г). Отсюда следует, что если

заданы коэффициенты многочленов устойчивости основной и промежуточной численных схем, то коэффициенты методов (4) однозначно определяются из линейных систем (9) и (11) с верхними треугольными матрицами Вг, 1 < г < т. Ниже этот факт будет использоваться при определении коэффициентов численных

формул (4) при условии, что области устойчивости промежуточных численных схем согласованы с основной. Разлагая точное и приближенное решения в ряды Тейлора по степеням к, можно записать:

к к

у(*п+1) = у (к)+к/+—//+—[ //2+/ '2/ ]+0(кА),

2 6

тт

уп+1 = уп + (Е Ъ1 ]Рт1 )к/ + (Е Ъ2]Рт1 )к /'п/п + (12)

1=1 1=2

тт

3 г/2 Г , А 7.2 \1 3 г" г2

+(ЕЪ,Рт,)к/+0.5(Ея¡Р,ч)к/’Л+0(к‘),

п

1=3 1=2

где элементарные дифференциалы вычислены на точном у(гп) и приближенном уп решениях соответственно. Из сравнения соотношений (12) видно, что численная формула (4) будет иметь первый порядок точности, если Е т=1 Ъ1 ]Рт] = 1. Требование второго порядка точности (4) означает выполнение дополнитель-

ного условия I т_2 Ь2 ,рт, = 0,5. Наконец схема (4) будет иметь третий порядок, если дополнительно выполняются еще два соотношения

т 1т 1

I Ьз ,р„ = 6 ■ I Х.р,, = 3 ■ (13)

}=3 6 }=2 3

Отсюда следует, что для построения т -стадийных методов первого и второго порядков точности в линейной системе (9) следует положить ст1 = 1 и ст1 = 1, ст2 = 0,5 соответственно. Задача о построении т -стадийных методов третьего порядка сводится к совместному решению линейной системы (9) и второго уравнения (13) при условии ст1 = 1, ст2 = 0,5, ст3 = 1/6. Учитывая, что Ь2 , = а,, 2 < , < т, нетрудно видеть, что совместность (9), (13) эквивалентна совместности соотношений

тт

1а,Рт! = 2 • ^Рт, = 3 • (14)

}=2 2 ,=2 3

где а,, 2 < , < т определены формулами (5).

Коэффициенты методов первого и второго порядков точности. Будем предполагать, что заданы коэффициенты многочленов устойчивости

Qk (г) = 1 +1 сыг , 1 < к < т. (15)

¡=1

При выборе коэффициентов многочленов (15) можно руководствоваться различными соображениями. Если, например, основная численная схема (4) применяется с максимальным интервалом устойчивости, то имеет смысл выбрать полиномы устойчивости промежуточных схем с аналогичными свойствами. Если же область устойчивости основной схемы "раздута" по мнимой оси, то и области устойчивости промежуточных схем должны быть устроены подобным образом. Это позволяет избежать неустойчивости промежуточных численных формул, когда шаг интегрирования по устойчивости основной схемы близок к максимальному. Коэффициенты нужных многочленов можно вычислить предложенным в [1] алгоритмом. В [3] приведены коэффициенты многочленов Qk (г), 2 < к < 13 с максимальным интервалом устойчивости.

Для каждого к,1 < к < т обозначим через | ук | длину такого максимального интервала [ук, 0],

что для всех г е [ук, 0] имеет место неравенство Qk (г) < 1. Учитывая, что г = кЛ, в (15) для каждого

Qk (г),1 < к < т проведем замену к на (кук 1ут). В результате вместо (15) получим:

й\(г) = 1 +1 сиг1, с'и =(Ук/ГтУси> 1 < к < т (16)

=1

Далее данные полиномы будут использоваться в качестве многочленов устойчивости методов (4).

Замена к на (кук/ут) означает, что приближенное решение по промежуточным схемам (4) вместо точек

(гп + ск1к),1 < к < т _ 1 будет вычисляться в точках (гп + с к1к),1 < к < т _ 1. В этом случае максимальный шаг из условия устойчивости основной схемы будет максимальным и для промежуточных численных формул.

Определение коэффициентов методов (4) будем осуществлять по следующему алгоритму. С использованием [1] вычислим коэффициенты полиномов (15), удовлетворяющими некоторым заданным свойствам. Затем вычислим коэффициенты многочленов (16) с применением соответствующей замены переменных. Учитывая, что коэффициенты (к +1) -го столбца матрицы Вт совпадают с коэффициентами многочлена

устойчивости @к(г), сформируем матрицу Вт, которая имеет вид:

Г 111 ... 1 ^

в =

0

с 11 с 21 ♦♦♦ с т_\,\

\

0 0 с 22 . с т _1 , 2

0 0 0 ••• с т_1 т_1

(17)

Используя в (11) вектор с\ = (с'къ• ••,с\к)Т, вместо ск все коэффициенты метода (4) однозначно

определяются из линейных систем (9) и (11).

Коэффициенты методов третьего порядка точности. При построении методов третьего порядка точности с согласованными областями устойчивости отличие от предыдущего случая заключается в том, что необходимо выполнение совместности уравнений (14). Поэтому будем полагать, что число с'и в матрице (17) есть пока неопределенный коэффициент полинома @ (г) = 1 + с 11г. Из линейной системы (9) последовательно определим ртт, ртт_1,.., рт3. Учитывая, что ц = с ¡-и, 2 < I < т, перепишем соотношения (14) в виде

т т

О 11рт2 = 12 - Е с']~1,1рт] , с Прт2 = 1/3 - Е ^)~1,1Рт} •

}=3 }=3

Разделив второе уравнение на первое, получим:

тт

с \1 = (113 - Е с ' 2-ир,„) / (1/2 - Е с',-1.А,).

}=3 }=3

Данное с 11 обеспечивает совместность (14). Значения рт2 и рт1 последовательно определим из второго и первого уравнений (9). Оставшиеся коэффициенты в,, 3 < I < т,1 < , < I -1 вычислим последовательным решением систем линейных алгебраических уравнений Вк@к = с\, 2 < к < т -1, где

вк = (вк+1,1’-"’вк+1,к ) ’ с к = (с к1,•••> с кк) .

Контроль точности вычислений. Для контроля точности методов первого порядка будем использовать оценку локальной ошибки 8п1. С применением (12) видим, что для т -стадийного метода она имеет вид

8„ = (0,5 - Е % Ь ,рт)к‘ /}' + 0(Н1) = (0.5 - ст1)Н2 // + 0(к1), где ст2 - коэффициент при г2 многочлена устойчивости (8). Оценку £п1 данной ошибки можно вычислить многими способами. Разлагая к{,1 < I < т в ряды Тейлора нетрудно видеть, что к{ - к, = (а{ - а, )к2/'п/п + 0(к3), I Ф Отсюда следует, что еп1 можно вычислить по приближенной формуле

£п,1 = [(0,5 - ст2)/(аI - а,)] (к - к,), 1 < I,} < т, I ф (18)

Для жестких задач характерным является быстрое изменение решения на небольших промежутках. Для того чтобы избежать неудовлетворительной точности, естественно в оценке ошибки использовать стадии, вычисленные в крайних точках отрезка [гп, гп+1]. Стадия к1 вычисляется в точке гп и поэтому в (18) положим } = 1,1 > }. Для методов с согласованными областями устойчивости имеет место неравенство а1 <а2 < ■■■ <ат < 1. Отсюда следует, что ни одна стадия в точке гп+1 не вычисляется. С другой стороны, применение в (18) стадий с достаточно большими номерами будет приводить к значительному увеличению вычислительных затрат в случае повторных вычислений решения вследствие невыполнения требуемой точности. С целью повышения эффективности расчетов будем поступать следующим образом. В качестве предварительной оценки будем применять величину £ п1 = [(0,5 -ст2)/а2](к2 -к1). Так как к1 зависит от размера шага линейно, то нарушение неравенства II £п111< £ будет приводить всего лишь к одному дополнительному вычислению правой части задачи (3). Здесь £ - требуемая точность расчетов, II ■ II - некоторая норма в Ям. Учитывая, что к/(уп+1) - к1 = к2/'п/п + 0(к3), окончательное решение по точности будем принимать на основе неравенства II £"nЛII<£, где £"пЛ = (0,5 - ст2)(к/ (уп+1) - к1). Оценку ошибки методов второго порядка будем вычислять с помощью схем первого порядка точности. Нетрудно видеть, что в этом случае нужно подправить соответствующим образом оценку ошибки схемы первого порядка.

Контроль точности методов третьего порядка будем осуществлять с помощью численных формул второго порядка точности. Для схем третьего порядка шаг интегрирования по устойчивости значительно меньше, чем для формул первого порядка. Поэтому при выборе шага ограничимся одной оценкой ошибки. С

использованием разложения стадий к1 и к2 в ряды Тейлора нетрудно видеть, что численная формула Уп+\ 2 = Уп + [(а2 -0,5)/а2]к1 + [0,$/а2]к2 имеет второй порядок точности. Тогда для контроля точности вычислений можно применять оценку ошибки £п3 = уп+1 - уп+12, где приближение уп+1 к точному решению y(tn+1) вычислено с третьим порядком.

Литература

1. Новиков, Е.А. Явные методы для жестких систем / Е.А. Новиков. - Новосибирск: Наука, 1997. - 197 с.

2. Golushko, M.I. Explicit fourth-order methods for stiff system / M.I. Golushko, E.A. Novikov // Russ. J. Numer. Anal. Math. Modelling. - 1999. - P. 71-85

3. Новиков, В.А. О построении явных методов типа Рунге-Кутта с расширенными областями устойчивости / В.А. Новиков, Е.А. Новиков. - Красноярск, 1998.

♦