Алгоритм переменной структуры на основе метода Ческино и l-устойчивой (2,2)-схемы Текст научной статьи по специальности «Математика»

МАТЕМАТИЧЕСКОЕ МОДЕЛИРОВАНИЕ И КОМПЬЮТЕРНЫЕ

ТЕХНОЛОГИИ

УДК 519.622

АЛГОРИТМ ПЕРЕМЕННОЙ СТРУКТУРЫ НА ОСНОВЕ МЕТОДА ЧЕСКИНО И L-УСТОЙЧИВОЙ (2,2)-СХЕМЫ

© 2013 г. Е.А. Новиков

Новиков Евгений Александрович - д-р физ.-мат. наук, профессор, главный научный сотрудник, Институт вычислительного моделирования СО РАН, г. Красноярск. Тел. (391)249-47-24. E-mail: [email protected]

Novikov Evgenii Alexandrovich - Doctor of Physico-Mathe-matical Science, professor, Chief Researcher, Institute of Computational Modeling SB RAS, Krasnoyarsk. Ph. (391)249-47-24. E-mail: [email protected]

Построено неравенство для контроля устойчивости схемы Ческино второго порядка точности. На основе стадий этого метода построена численная формула первого порядка с расширенным до 32 интервалом устойчивости. На основе L-устойчивой (2,2)-схемы и численной формулы Ческино разработан алгоритм переменной структуры, в котором эффективная численная формула выбирается на каждом шаге по критерию устойчивости. Алгоритм предназначен для решения как жестких, так и не жестких задач. Приведены результаты расчетов, подтверждающие эффективность построенного алгоритма.

Ключевые слова: жесткая задача; схема Ческино; (2,2)-метод; контроль точности и устойчивости.

An inequality for stability control of a Ceschino's scheme of second order of accuracy is constructed. A numerical formula of order one is developed that is based on the stages of the this method and its stability interval is extended to 32. On a base of L-stable (2,2)-scheme and a numerical Ceschino's formula, an algorithm of alternating structure, in which an efficient numerical formula is chosen on an every step by a stability criterion, is constructed. The algorithm is intended for solving stiff and non-stiff problems. There are shown results of calculations, confirming efficiency of this algorithm.

Keywords: stiff problem; Ceschino's scheme; (2,2)-method; accuracy and stability control.

Введение

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

Для повышения эффективности расчетов в ряде алгоритмов используется замораживание матрицы Якоби, т. е. применение одной матрицы на нескольких шагах интегрирования [1]. Наиболее успешно этот подход применяется в алгоритмах на основе неявных методов типа Рунге - Кутта и многошаговых численных формул [2]. В алгоритмах интегрирования на основе известных безытерационных методов, к которым относятся схемы типа Розенброка [3] и их различные модификации [4], вопрос о применении одной матрицы на нескольких шагах интегрирования более сложный. В таких алгоритмах матрица Якоби включена в численную схему, и ее аппроксимация может

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

Некоторым аналогом замораживания матрицы Якоби является применение в расчетах алгоритмов интегрирования на основе явных и ¿-устойчивых методов с автоматическим выбором численной схемы [5]. В этом случае эффективность алгоритма может быть повышена за счет расчета переходных участков явным методом. В качестве критерия выбора эффективной численной формулы естественно применять неравенство для контроля устойчивости [6]. Отметим, что применение таких комбинированных алгоритмов полностью не снимает проблему замораживания матрицы Якоби.

Здесь построено неравенство для контроля устойчивости метода Ческино второго порядка точности. На основе стадий данной численной формулы построена схема первого порядка с расширенной до 32 единиц по вещественной оси областью устойчивости. На основе ¿-устойчивой (2,2)-схемы и рассмотренных яв-

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

L-устойчивая (2, 2)-схема

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

У = f (У), y(to) = Уо, to < t < tk

(1)

y«+i = Уп + aki +0,5a k2;

Dnki = hnf (Уп); Dnk2 = hnf (Уп + ak1) " 2ak^

(2)

S = h-

a - 3 ) f *f+ (1" 4 «I ff2

+ 0(h4)

где у и f - вещественные Ж-мерные вектор-функции; t - независимая переменная; рассмотрим (2,2)-метод вида [4]:

Тогда согласно [6] для контроля точности (2) можно использовать оценку ошибки вида

еи = (а-1/3)К2/2/ + O(h3).

Учитывая, что

к2 + (2а - 1)к = (а - 2а2)К2/2/ + O(h3),

величину еп с точностью до членов О(К3) можно оценить по формуле

еп=[(а-1/3)/(а-2а2)][к2+(2а-1)к1].

В результате для контроля точности схемы (2) можно применять неравенство [6]

Кin \k2 +(2a-1)k]||<

a - 2a

a-1/3

(3)

где а = 1- \/2 / 2 ; hn - шаг интегрирования; к1 и к2 -стадии метода; Dn=E-ah„An; Е - единичная матрица; Ап - матрица, представимая в виде Ап = fn + hnBn + + O(h2n); /п = д/(у,)/ду - матрица Якоби системы (1); Вп - не зависящая от шага интегрирования матрица. Использование матрицы Ап, представимой в виде Ап = f'n + hnBn + O(h2n), позволяет применять метод (2) с замораживанием как аналитической, так и численной матрицы Якоби [7]. При использования матрицы Якоби/п-к, вычисленной к шагов назад, имеем:

Ап = Гп-ккп/:/п + (К2);

Вп =-кГ^„, 0 < к < ^ ,

где ik - максимальное число шагов с замороженной матрицей Якоби; /"п /п = [д2/(уп) / ду2\/(уп). Если матрица Якоби вычисляется численно с шагом г = сД где с, 1 < ] < Ж, - постоянные числа, то получим

Ап = /п + hnGn + О(К2), Вп = Gn, элементы g„,iJ■ матрицы Gn определяются по формулам

gn,y■ = 0,5^/^)/ду2, 1 < /■, ] < Ж.

В случае замораживания численной матрицы Якоби можно записать

Ап = /п - К ^п - к//) + К), Вп = Gn - к/п/п.

Так как при записи схемы (2) матрица Вп произвольная, то вопрос о замораживании и численной аппроксимации матрицы Якоби можно рассматривать одновременно.

Локальную ошибку 5п метода (2) можно записать в виде

где 1 < }п < 2; е - требуемая точность интегрирования; ||-|| - некоторая норма в ЯЖ, а целочисленный параметр }п выбирается наименьшим, при котором выполняется неравенство (3).

Оценка максимального собственного числа ^п0 = = КХп,тах матрицы Якоби системы (1), необходимая для перехода на явную численную схему, вычисляется по формуле

wn, 0 = h

N

= h max У

1<i<N j=1

¿У,

Метод Ческино

Для решения системы (1) рассмотрим явную формулу типа Рунге - Кутта вида

Уп+1 = Уп + Рт1к1 + ... + Рт4к 4;

к1 = К/ (tn, Уп); к2 = К/(tn + 0,25К, Уп + 0,25^); (4)

кз = К/ (^ + 0,5К, уп + 0,5к2);

к4 = К/^п + К,Уп + к -2к2 + 2кз),

где К - шаг интегрирования; к - стадии метода; рт -числовые коэффициенты; т - порядок точности метода. При коэффициентах [8]

Р21 = ^ Р22 = -2, Р23 = 2 , Р24 = 0

(5)

схемы (4), (5) являются ¿-устойчивыми и имеют второй порядок точности. Схема (4) с коэффициентами р41 = р44 = 1/6, р42 = 0 и р43 = 2/3 имеет четвертый порядок. Тогда для контроля точности схемы второго порядка можно использовать оценку ошибки 5п2 вида

Ж

8п,2 =е (рм - рл )к, .

ь

В результате для контроля точности вычислений применяется неравенство ||5п,2 ||< е, где ||-|| - некоторая норма в Еы, е - требуемая точность расчетов. Учитывая, что имеет место соотношение 5п2 = О(И3), шаг Иас по точности выбирается по формуле Иас = qh, где q находится из уравнения q3||Sn,2|| = е. Если q <1, то происходит повторное вычисление решения (возврат) с шагом И, равным qh. В противном случае вычисляется приближенное решение, а прогнозируемый шаг Ип+1 вычисляется по формуле Ип+1 = qh. Неравенство 5и 2 < е хорошо зарекомендовало себя при решении многих практических задач, и будет использоваться ниже.

Контроль устойчивости метода Ческино

Построим неравенство для контроля устойчивости схемы (4). Для этого применим (4) для решения линейной задачи у' =Ау с постоянной матрицей А. Первые три стадии ^ и ^ для у'=Ау имеют вид:

к, = Хуп , k2 = [X + 0,25X2]у„ , ^ = [X + 0,5X 2 + 0,125X3]уп ,

где X = hA. Нетрудно видеть, что имеют место соотношения:

- 2^ + ^ = 0,125X3уп , 0,5(^ - k1) = 0,125X2уп .

Оценку максимального собственного числа матрицы Якоби можно вычислить степенным методом [6]. Введем обозначение

hn+l = max

w„, = 2 max

1<i< N

(k -2k2 + k3)i

(k2 k1 )i

(6)

hn, min (hac, hst)], (7)

где hn есть последний успешный шаг интегрирования.

Формула (7) применяется для прогноза шага интегрирования Ип+1 после успешного вычисления решения с предыдущим шагом hn и поэтому фактически не приводит к увеличению вычислительных затрат. Если шаг по устойчивости меньше последнего успешного, то он уменьшен не будет, потому что причиной этого может быть грубость оценки максимального собственного числа. Однако шаг не будет и увеличен, потому что не исключена возможность неустойчивости численной схемы. Формула (7) позволяет стабилизировать поведение шага на участке установления решения, где определяющую роль играет устойчивость.

Из результатов расчетов алгоритмом интегрирования с контролем точности и дополнительным контролем устойчивости следует, что фактическая точность вычисления решения на интервале установления значительно выше задаваемой. Это естественно, потому что старые ошибки подавляются за счет контроля устойчивости, а новые ошибки невелики за счет малости производных решения. В такой ситуации эффективнее проводить вычисления методом более низкого порядка точности с широкой областью устойчивости.

Метод первого порядка точности

Для численного решения задачи (1) рассмотрим схему вида

Уп+1 = Уп + r1k1 + ... + r4k4

(8)

Тогда для контроля устойчивости метода Ческино можно применять неравенство м>„1 < Б, где число Б ограничивает интервал устойчивости.

Устойчивость методов типа Рунге - Кутта обычно исследуется на скалярном тестовом уравнении у' = Ху, где Х есть произвольное комплексное число, Re(Х) < 0. Смысл Х - некоторое собственное число матрицы Якоби задачи (1). Применяя (4), (5) для решения у=Ху, получим, что функция устойчивости 02(х) метода второго порядка точности имеет вид Q2(x) = 1 + х + + 0,5х2 + 0,125х3, х = hХ. Интервал устойчивости метода второго порядка равен двум, а метода четвертого порядка приблизительно равен 2,8. Поэтому в неравенстве м>„1 < Б положим Б = 2. Учитывая, что ^пд=О(И), шаг к1' по устойчивости можно выбирать по формуле И—И, где г вычисляется из равенства г^пд= 2.

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

где стадии ^ вычислены по формулам (4), а коэффициенты Г1 подлежат определению. Заметим, что при П=1, г2= -2, г3=2 и г4=0 численная формула (8) имеет второй порядок точности и совпадает с (4) с коэффициентами (5). Построим менее точную схему с максимальным интервалом устойчивости. Для этого применим (8) для решения скалярного тестового уравнения у' = Ху. Получим уп+1= Q(x)yn, где функция устойчивости Q(x) имеет вид

Q( х) = 1 + (г1 +... + г4 ) х + (0, 25г2 + 0,5г3 + +г4) х 2 +(0,125г3 + 0, 5г4 ) х3 + 0,25г4 х4.

Требование первого порядка приводит к соотношению Г1+...+г4=1, которое ниже будем считать выполненным. Оставшиеся коэффициенты Г1 выберем таким образом, чтобы метод (8) имел максимальный интервал устойчивости. Для этого рассмотрим многочлен Чебышева 7,4(г)=824-822+1 на промежутке [-1,1]. Проведем замену переменных т=1-2х/у, получим

32 160

256

128

T4(x) = 1--x + —- x--— x + —- x

У

У

У

У

при этом отрезок [у, 0] отображается на [-1,1]. Нетрудно показать, что среди всех многочленов Р4(х)

вида Р4(х) = 1+х+с2х2+с3х3+с4х4 для Т4(х) неравенство |Т4(х)|<1 выполняется на максимальном интервале [у, 0], у = -32. Потребуем совпадения коэффициентов Q(x) и Т4(х) при у = -32. Это приводит к соотношениям Г!+...+г4=1, 0,25г2 + 0,5г3 + г4= 5/32, 0,25г4 = 1/8192, 0,125г3 + 0,5г4 = 1/128. В результате имеем коэффициенты Г\ = 895/2048, г2 = 257/512, г3 = 31/512 и г4 = = 1/2048 метода первого порядка точности с максимальным интервалом устойчивости, локальная ошибка 5п которого имеет вид 5п = 9К2/'/32 + О(К3). С учетом соотношения к2 - к\ = 0,25К2/ '/п + О(К3) и вида локальной ошибки неравенство для контроля точности записывается в виде || к2-ка||<е, где ||-|| - некоторая норма в ЯЖ, е - требуемая точность расчетов. Интервал устойчивости численной схемы (8) равен 32. Поэтому для ее контроля устойчивости можно применять неравенство < 32.

Алгоритм переменной структуры

На основе построенных явных методов первого и второго порядков точности легко сформулировать алгоритм переменного порядка и шага. Расчеты всегда начинаются методом второго порядка как более точным. Переход на схему первого порядка осуществляется при нарушении неравенства < 2. Обратный переход на метод второго порядка происходит в случае выполнения неравенства < 2. При расчетах по методу первого порядка наряду с точностью контролируется устойчивость, а выбор прогнозируемого шага производится по аналогии с методом второго порядка точности по формуле типа (7).

В случае использования схемы (2) формулировка алгоритма интегрирования также не вызывает трудностей. Нарушение неравенства < 32 вызывает переход на ¿-устойчивую схему (2). Передача управления явным методам происходит в случае выполнения неравенства ^п0 < 32, где оценка ^п0 вычисляется через норму матрицы Якоби.

Численную формулу (2) без потери порядка точности можно применять с замораживанием матрицы Dn. При замораживании матрицы Якоби величина шага интегрирования остается постоянной с целью сохранения свойства ¿-устойчивости метода (2). Попытка замораживания матрицы Dn осуществляется после каждого успешного шага. Размораживание матрицы происходит в следующих случаях: 1) нарушение точности расчетов; 2) если число шагов с замороженной матрицей достигло заданного максимального числа ih; 3) если прогнозируемый шаг больше последнего успешного в qh раз. Числами ih и qh можно влиять на перераспределение вычислительных затрат. При ih=0 и qh=0 замораживания не происходит, при увеличении ih и qh число вычислений правой части возрастает, а количество обращений матрицы Якоби убывает. В расчетах использовались значения ih = 10 и qh = 2. Норма ||ф|| в неравенствах для контроля точности вычисляются по формуле

||ф||= max {|фг|/(I y'n\ + г)} ,

1<i<N I vi i ))

где i - номер компоненты; r - положительный параметр. Если по i-й компоненте решения выполняется неравенство \y„'\<r, то контролируется абсолютная ошибка re, в противном случае - относительная ошибка е. Ниже построенный алгоритм переменного порядка и шага, а также с автоматическим выбором явной или ¿-устойчивой численной схемы будем называть MK_CES2.

Результаты расчетов

Расчеты осуществлялись с задаваемой точностью е = 10-2 на PC l„tel(R) Core i7-3770S [email protected] с двойной точностью. Невысокая точность расчетов связана с тем, что в построенном алгоритме применяются схемы низкого порядка точности, и поэтому данным методом осуществлять расчеты с более высокой точностью нецелесообразно. Сравнение эффективности проводилось с известным методом Гира в реализации А. Хиндмарша D¿SODE из коллекции ODEPACK [9] и алгоритмом RKMK2 [7] на основе ¿-устойчивого и явных двухстадийных методов. В качестве тестовой задачи выбрана модель реакции Белоусова - Жаботинского, для которой на промежутке интегрирования характерно наличие трех переходных процессов. Тестовый пример имеет вид:

У = 77,27 (y2 - yiУ2 + y - 8,375 • 10-6y2 ) ;

У2 =(уз -У2 -У1У2)/77,27 ;

У3 = 0,161(yi - Уз) , t е [0,300], у(0) = 4;

У2(0) = 1,1, У3(0) = 4, h0 = 2 •Ю-3.

Расчеты проводились с численной матрицей Яко-би. Сравнение эффективности осуществлялось по числу вычислений правой части ifu и количеству декомпозиций матрицы Якоби ija на интервале интегрирования. Решение этой задачи алгоритмом MK_CES2 вычислено с затратами ifu = 1 077 и ija = 58. При вычислениях по ¿-устойчивой схеме (2) затраты ifu = 1072 и ija = 82. Фактическая точность расчетов в конце интервала интегрирования не хуже задаваемой. Решение задачи удалось вычислить явными методами переменного порядка и шага с затратами ifu = = 978 524. Данная задача слишком жесткая для явных методов. Однако результаты расчетов приведены здесь с целью демонстрации принципиальной возможности применения явных методов с контролем устойчивости и переменным порядком для решения достаточно жестких примеров. Явные методы на некоторых жестких задачах большой размерности могут быть эффективнее ¿-устойчивых методов. Решение данной задачи алгоритмом RKMK2 [7] вычислено с затратами ifu = 1 214 и ija = 65. Решение задачи явными методами переменного порядка и шага из алгоритма RKMK2 вычислено с затратами ifu = 2 112 678.

При расчетах программой DLSODE требуемая точность е=10-2 достигается при задаваемой точности е = 10-4 с затратами = 1 129 и уа = 107. При более высокой точности расчетов DLSODE эффективнее построенного алгоритма. Это является следствием низкого порядка точности построенных численных формул.

При точности е = 10-2 алгоритм MK_CES2 более чем в 2 раза эффективнее алгоритма DLSODE по числу декомпозиций матрицы Якоби, в то время как количество вычислений правой части для MK_CES2 и DLSODE различается незначительно. В случае большой размерности задачи (1) построенный алгоритм интегрирования по времени счета может быть существенно эффективнее программы DLSODE.

Заключение

Построенный алгоритм MK_CES2 предназначен для расчетов с небольшой точностью - порядка 1 % и ниже. В этом случае достигается его максимальная эффективность. В MK_CES2 с помощью признака можно задавать различные режимы: 1) явными методами первого или второго порядков точности с контролем или без контроля устойчивости; 2) явными методами с переменным порядком и шагом; 3) L-устойчивым методом с замораживанием или без замораживания аналитической или численной матрицы Якоби; 4) с автоматическим выбором численной схемы. Все это позволяет применять данный алгоритм для решения как жестких, так и нежестких задач. При расчетах с автоматическим выбором численной схемы вопрос о том, является ли задача жесткой или нет, перекладывается на алгоритм интегрирования.

Применение на участке установления явного метода первого порядка точности с расширенной областью устойчивости позволяет в 16 раз увеличить размер шага интегрирования по сравнению с явным методом второго порядка без увеличения вычислительных затрат. На переходных участках, где определяющую роль играет точность вычислений, более

Поступила в редакцию

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

Работа поддержана РФФИ (проекты 11-01-00106 и 11-01-00224).

Литература

1. Hairer E., Wanner G. Solving ordinary differential equations. Stiff and differential-algebraic problems. Springer-Verlag. 1996. 601с. [Хайрер Э., Ваннер Г. Решение обыкновенных дифференциальных уравнений. Жесткие и дифференциально-алгебраические задачи: пер. с англ. / под ред. С.С. Филиппова. М., 1999. 685 с.]

2. Byrne G.D., Hindmarsh A.C. ODE solvers: a review of current and coming attractions // J. of Comput. Physics. 1987. № 70. P. 1 - 62.

3. Rosenbrock H.H. Some general implicit processes for the numerical solution of differential equations // Computer. 1963. №5. P. 329 - 330.

4. Новиков Е.А., Шитов Ю.А., Шокин Ю.И. Одношаговые безытерационные методы решения жестких систем // Докл. АН СССР. 1988. Т. 301, № 6. С. 1310 - 1314.

5. Новиков Е.А. Построение алгоритма интегрирования жестких систем дифференциальных уравнений на неоднородных схемах // Докл. АН СССР. 1984. Т. 278, № 2. С. 272 - 275.

6. Новиков Е.А. Явные методы для жестких систем. Новосибирск, 1997. 197 с.

7. Новиков A.E., Новиков E.A. Численное решение жестких задач с небольшой точностью // Мат. моделирование. 2010. Т. 22, № 1. С. 46 - 56.

8. Ceschino F., Kuntzman J. Numerical solution of initial value problems. New Jersey: Prentice-Hall, Englewood Clis. 1966. 287 p.

9. Hindmarsh A.C. ODEPACK, a systematized collection of ODE solvers // Lawrence Livermore National Laboratory. 1982. Preprint UCRL-88007.

21 августа 2013 г.