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

Вычислительные технологии

Том 16, № 6, 2011

Комбинированный алгоритм третьего порядка для решения жестких задач*

А. Е. Новиков1, Е. А. Новиков2

1 Сибирский федеральный университет, Красноярск, Россия

2

e-mail: aenovikov@bk.ru, novikov@icm.krasn.ru

Построены явная трехстадийная схема типа Рунге — Кутты и Ь-устойчивый (3,2)-метод третьего порядка точности для решения жестких задач. Разработан алгоритм интегрирования переменного шага, в котором выбор эффективной численной схемы осуществляется на каждом шаге с применением неравенства для контроля устойчивости. Приведены результаты расчетов, подтверждающие эффективность предложенного алгоритма.

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

Введение

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

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

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

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

В настоящей работе на основе явного метода типа Рунге —Кутты и ¿-устойчивого (3,2)-метода третьего порядка точности построен комбинированный алгоритм переменного шага, в котором выбор эффективной численной формулы осуществляется по устойчивости, Приведены результаты расчетов, подтверждающие эффективность предложенного алгоритма интегрирования.

1. Класс (т, к)-методов решения жестких задач

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

у' = /(У), У(ьо) = Уo, Ь < Ь < 4, (1)

где у и / — вещественные Ж-мерные вектор-функции, Ь — независимая переменная. Изучение автономной системы (1) не снижает общности, В случае неавтономной задачи введением дополнительной переменной ее всегда можно привести к автономному виду, В методах типа метода Розенброка матрица Якоби введена в численную формулу, что приводит на каждом шаге к необходимости решения линейных систем алгебраических уравнений. Это существенно упрощает программную реализацию алгоритмов интегрирования, но возникают проблемы с использованием одной матрицы Якоби на нескольких шагах интегрирования. Заметим, что в неявных и полуявных методах типа метода Рунге —Кутты необходим итерационный процесс типа ньютоновского, что приводит к более сложному алгоритму интегрирования. Прежде чем описывать класс (т, к)-методов, рассмотрим одностадийный метод типа метода Розенброка

Уп+1 = Уп + рк\, Опк\ = к/(уп), Вп = Е + акЩ, (2)

где Е — единичная матрица, /' = д/(уп)/дУ — матрица Якоби системы (1), к1 — стадия метода, а и р1 — коэффициенты, определяющие свойства точности и устойчивости (2), к1

ставлепие приближенного решения уп+1 в виде ряда Тейлора по степеням к. Разложение точного решения у (Ьп+1) в ряд Тейлора в окрестноети точки Ьп до членов с к4 включительно имеет вид

уЦп+1) = УЦП) + к/ + \к2П + \к\г2! + /72)+

+^4(//3/ + /772 + з/772 + Г/3) + о(к5), (з)

где элементарные дифференциалы вычислены на точном решении у(£п). При условии Уп = У(£п), сравнивая полученное представление приближенного решения Уп+1 с рядом Тейлора для точного решения у(£п+1), получим, что схема (2) будет иметь второй порядок точности, если выполнены соотношения р1 = 1 и ар1 = 0.5 ил и р1 = 1 и а = 0.5, Применяя (2) для решения тестовой задачи у' = Ау, у(0) = у0, £ > 0, получим уп+1 = ^(г)уп, где г = ЛА, а функция устойчивости имеет в ид = [1 + (р1 — а)г]/(1 — аг). Отсюда с учетом условий порядка следует, что одностадийный метод типа метода Розенброка второго порядка точности является А-устойчивым, Требование ¿-устойчивости приводит к соотношению р1 = а, которое противоречит второму порядку. Обычно от второго порядка отказываются в пользу ¿-устойчивости.

Теперь рассмотрим численную формулу следующего вида:

Уп+1 = Уп + ЬЛ + Ы2, ад = Л/(Уп), = кь (4)

Разлагая ^ и к2 в ряды Тейлора и подставляя в первую формулу (4), получим Уп+1 = Уп + (Ь1 + Ь2)Л/п + а(Ь1 + 2^2) Л2 / /п + 0(Л3).

Нетрудно видеть, что схема (4) будет иметь второй порядок точности, если Ь1 + Ь2 = 1 и а(Ь1 + 2Ь2) = 0.5, где параметр а удовлетворяет условию ¿-устойчивости Ь1 = а, которое с учетом условий порядка переписывается в виде а2 — 2а + 0.5 = 0. В результате имеем коэффициенты Ъ\ = а и Ь2 = 1 — а, при которых схема (4) имеет второй порядок точности и является ¿-устойчивой. Уравнение а2 — 2а + 0.5 = 0 имеет два корня: а1 = 1 —0.5л/2 и а2 = 1 + 0.5-\/2. В расчетах обычно применяется параметр а = а\, так как в этом случае величина коэффициента в локальной ошибке меньше.

При большой размерности исходной задачи (1) основные вычислительные затраты приходятся на декомпозицию матрицы Оп (порядка N3 арифметических операций). На

фоне декомпозиции обратный ход в методе Гаусса (порядка N арифметических опе-

/

вычиелительным затратам отличается от (2) на один дополнительный обратный ход

¿

ленной формулы (2) имеет второй порядок. По этой причине она эффективнее (2), Более того, для метода (4) легко построить неравенство для контроля точности вычислений, которое имеет вид ||к2 — А1|| < е [9], где е — требуемая точность интегрирования, || • || — некоторая норма в Ям.

Рассмотрение стадий вида Дпк2 = А1 привело к идее (т, к)-методов [10], которые определяются следующим образом. Пусть Z есть множество целых чисел и заданы т, к € Z, к < т. Обозначим через Мт множество чисел г € Z, 1 < г < т, а через Мк и ¿¿, 1 < г < т — подмножества из Мт гада Мк = (тг € Мт|1 = т1 < т2 < ... < т^ < т} и ¿г = (т^-1 € Мт|^' > 1, ту € Мк, ту < г}. Рассмотрим следующие численные схемы:

т

Уп+1 = Уп + ^ рА ,

г=1

(г—1 \ г—1

Уп + ^ вгу М + Ьу + Л/4 ^ сц Ьу, г € Мк,

3 = 1 ' М=1

г—1

ДА = Ьг—1 + ^ ау Ьу + / ^ сцЬу, г € Мт—к, (5)

Ме^ 3=1

где Мт-к = Мт\Мк, Ьг, 1 < г < т, — стадии метода; а, Рг, в%з, а^ и с^ — постоянные коэффициенты; к — шаг интегрирования; Е — единичная матрица; / = д/(уп)/ду — матрица Якоби системы (1); к — количество вычислений функции / на шаге; т — число стадий или количество обратных ходов в методе Гаусса, На каждом шаге интегрирования осуществляются одно вычисление матрицы Якоби и одна декомпозиция матрицы Оп. Допускается аппроксимация матрицы Якоби / матрицей Ап гад а Ап = /+кВп+О(к2), где матрица Вп не зависит от величины шага интегриров ания. Так как кит полностью определяют затраты на шаг, а набор чисел т1,..., из множества Мк только распределяет их внутри шага, то методы типа (5) названы (т, к)-методами. Отметим, что при к = т и а^ = с^ = 0 методы (5) совпадают со схемами типа схемы Розенброка [4], а при к = т и а^ = 0 — с КОШ-методами [1]. В отличие от КОШ-методов в схемах (5) более точно определены затраты на шаг интегрирования и более правильно описана область определения коэффициентов численных формул, что упрощает их исследование и делает предпочтительнее,

к=т

когда число стадий и количество вычислений правой части задачи (1) совпадают, В случае схем (5) при т > к можно построить численные схемы более высокого порядка точности [11]. Основное отличие приведенных здесь методов (5) от существующих численных формул заключается в том, что в них стадия метода не связывается с обязательным вычислением правой части исходной задачи (1), за счет чего свойства данных методов могут быть улучшены,

2. ¿-устойчивый (3,2)-метод

Рассмотрим (5) при т = 3, т.е. численную формулу вида

Уп+1 = Уп + Р1к1 + Р2к2 + Рзкз, Бпк1 = к/(уп), В пк2 = к1, Вкз = к/(уп + вз1к1 + вз2 к2) + аз2к2, (6)

где матрица Оп определена в (5), Разложим стадии кг, 1 < г < 3, в ряды Тейлора к

соотношение:

Уп+1 = Уп + [Р1 + Р2 + (1 + аз2)рз]к2к/+ + [ар1 + 2ар2 + (а + вз1 + вз2 + 3ааз2)рз]к2 //+

+ [а2р1 + 3а2р2 + (а2 + 2авз1 + 3авз2 + 6а2аз2)рз]кз /2/+

+ ¿(031 + 032 ?р^гпй + + /зз2)3рз]^4/г/п + +а(вз1 + вз2)(вз1 + 2вз2)Рзк4/'/'„,/+

+ [азр1 + 4азр2 + (аз + 3а2вз1 + 6а2вз2 + 10азаз2)рз]к/з/п+ 1

+ -а(/3 31 + + 0(к5). (7)

Ряд Тейлора для точного решения y(tn+1) в окрестности точки tn имеет вид (3), Полагая yn = y(tn) и сравнивая ряды для точного и приближенного решений до членов с h3 включительно, найдем условия третьего порядка точности схемы (6):

Р1 + Р2 + (1 + «32)Р3 = 1,

api + 2ар2 + (а + /531 + /З32 + Заск32)р3 =

a2pi + За2р2 + (а2 + 2а/331 + 3 а/332 + 6а2а32)р3 = 7,

6

(/Зз1+/ЗЗ2)2РЗ = ^. (8)

Положим в системе (8) коэффициенты а, в31 и в32 свободными, и исследуем ее на совместность, Значение p3 выразим из последнего уравнения (8), Умножим первое соотношение (8) на —а и сложим со вторым. Затем умножим первое уравнение на -а2 и сложим с третьим. В результате имеем

ар2 + 1аа32р3 = а - (Зр3,

2а2р2 + Ъа2а32р3 = \ - а2 - а(2/331 + 3/332)р3, 6

где в = в31 + в32- Разрешая данную линейную систему отноеительно p2 и а32, значение p1 определим из первого соотношения (8). В результате можно записать

_ /32(18а2 - 9а + 1) + 2a(/?3i - а)

Pl

p2 =

6а2в2

в2(—18а2 + 15а — 2) + 2а(в32 — £31)

6а2в

2

1 в2 (6а2 — 6а +1) — 2авз2 = «з2 =-—2-• (9)

Исследуем устойчивость схемы (6) па линейном скалярном уравнении у' = Ау. Применяя (6) для решения данного уравнения, получим уп+1 = ^(г)уп, где г = ЛА, а функция устойчивости записывается в виде

_ (-а3 + а2р\ + а2р3 - а(33гр3)г3 (1 - аг)3

[За2 - 2ар\ - ар2 + (/З31 + /З32 - 2а)р3]г2 (1 - аг)3 [-За + рх + р2 + (1 + о;32)р3]г + 1 (1 - аг)3

Из вида следует, что для ¿-устойчивости схемы (6) необходимо выполнение соотношения а2 — а(р1 + р3) + в31р3 = 0. Подставляя сюда коэффициенты (9), получим

а

6а3 — 18а2 + 9а — 1 = 0.

(10)

Сравнивая представление приближенного и точного решений до членов с к4 включительно, видим, что локальная ошибка будет минимальной, если

(Рз1 + /Зз2)3Рз = р а(р 31 + /З32)(/З31 + 2/332)рз =

Теперь отсюда и (9) получим набор коэффициентов

130а2 - 33а + 6 -54а2 + 21а - 4 16

Р1 =-54а2-' Р2 =-18^-' = 27'

48а - 3 3 - 24а 54а2 - 30а + 6

~~ 32а ' ~~ 32а ' ^ ~ 32^ '

при которых локальная ошибка 6п,з схемы (6) имеет вид

^ = 1"12а+^°2"24а3т+тг+о(Л»), (п)

где значение а определяется из условия ¿-устойчивости (10), Уравнение (10) имеет три вещественных корня:

а1 = 2.40514957850286, а2 = 0.158983899988677, аз = 0.435866521508459.

Аа 1/3 < а < 1.0685790, поэтому выбираем корень а = аз, В результате имеем следующий комплект коэффициентов:

а = +0.43586652150846, рх = 0.15902052285216 ■ 101, р2 = -0.14930556622438 ■ 101, рз = 0.59259259259259, вз1 = 0.12849112162238 ■ 101, вз2 = -0.53491121622384, аз2 = 0.52356010690630.

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

В жестких задачах поведение ошибки определяется элементарным дифференциалом /'з/. В силу этого согласно [13] при построении оценки аналога глобальной ошибки будем учитывать только первое слагаемое в (11). Для контроля точности вычислений и автоматического выбора величины шага интегрирования используем идею вложенных методов, для чего рассмотрим двухетадийный метод (4) в виде уп+1,1 = уп + Ь1к1 + Ь2к2, где приближение уп вычислено по формуле (6). Отметим, что в данной численной формуле применяются стадии метода (6), и поэтому ее применение практически не приводит к увеличению вычислительных затрат. С использованием (3) и (7) нетрудно видеть, что при значениях коэффициентов Ь1 = 0.5(4а - 1)/а и Ь2 = 0.5(1 - 2а)/а данная схема имеет второй порядок точности, а ее локальная ошибка 6п,2 имеет вид

6а2 - 6а + 1 з 2 4

5п,2 =-«-Ь:! ! + 0{Н ).

6

еп

1 — 12а + 36а2 — 24а3

4(6а2 - 6а + 1)

"(Уп+1 — Уп+1,1).

Подчеркнем особенность построенной оценки ошибки. Из вида функции устойчивости схемы (6) следует, что при г — —то выполняется соотношение — 0. Для

точного решения у(£п+1) = ехр(г)у(£п) задачи у' = Ау выполняется аналогичное евой-

еп

г — —то. Однако для построенной оценки имеет место еп = 0(1), так как оценивается

еп

рассмотрим оценки еп(^п), 1 < < 2, вида [8]

1 - 12а + 36а2 - 24а3 , ■ £пЫ)- 4(6а2 — 6а + 1)

При = 1 оценка еп(^п) совпадает с еп и будет А-устойчивой, а при = 2 — ¿-устойчивой, Теперь неравенство для контроля точности имеет вид

11^—"(Уп+1 — Уп+1,1)|| < Се, (12)

где с = 4-|6а2 — 6а+1|/|1 —12а+36а2 — 24а31 ~ 3, е — требуемая точность интегрирования, а значение параметра выбирается наименьшим, при котором выполняется неравенство (12), Заметим, что в смысле главного члена оценки еп(1) и еп(2) совпадают. Кроме того, неравенство (12) при = 2 проверяется редко, в основном при резком увеличении шага интегрирования. Поэтому применение еп(^п) вместо еп не приводит к значительному росту вычислительных затрат. Более того, за счет применения еп(^п) вместо еп эффективность расчетов повышается примерно на 10 ^ 15 %, Это связано с устранением некоторых повторных вычислений решения, возникающих за счет неправильных свойств оценки еп. Норма ||£|| в левой части неравенства (12) вычисляется по формуле ||£|| = тах |£г|/(|уп| + V). При выполнении неравенства |уп| < V по г-й компоненте

1<г<М

решения контролируется абсолютная ошибка vе, в противном случае контролируется

е.

Сделаем некоторые замечания по выбору шага численного дифференцирования. При численном вычислении матрицы Якоби данный шаг гг, 1 < г < N, выбирается по формуле Гц = тах[гтщ, |уц|], 1 < ] < N, где гт;п — минимальный шаг численного дифференцирования, зависит от разрядности ЭВМ, При расчетах с двойной точностью величина гт;п полагалась равной 10—14, Тогда ^'-й столбец ап численной матрицы вычисляется по формуле ап = /(у1, • • • , Уу + Гц, • • • , у^) — /(у1, • • • , у3-, • • • , у^)/г3-, т. е, для задания матрицы требуется N вычислений правой части системы уравнений (1),

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

7 11 11 ,

1»п,о = п т- = II тах > —- .

ду 1<г<М ^ ду]

Ниже данная оценка будет применяться для автоматического выбора численной схемы.

3. Явный метод третьего порядка

Для численного решения задачи (1) рассмотрим явный трехетадийный метод типа метода Рунге — Кутты вида

Уп+1 = Уп + РЛ + Р2к2 + Рзкз, к1 = к/(уп), к2 = к/(уп + в21к1), кз = к/(уп + взЛ + вз2к2),

(13)

где к1; к2 и кз — стадии метода, р1,р2,рз, в21, вз1, вз2 _ числовые коэффициенты, определяющие свойства точности и устойчивости (13), В случае неавтономной системы

У' = /(^У), У(*о) = Усъ ¿о < * < ¿к, схема (13) записывается в виде

Уп+1 = Уп + Р1к1 + рк + Рзкз, к1 = к/(¿п, Уп), к2 = к/(¿п + в21к,уп + в21к1), кз = к/(¿п + [вз1 + вз2]к, Уп + взЛ + вз2к2).

Ниже для сокращения выкладок будем рассматривать задачу (1), Однако построенные методы можно применять для решения неавтономных задач. Получим соотношения на

к1

к2 кз к4

пия в первую формулу (13), В результате получим

Уп+1 = Уп + (Р1 + Р2 + Рз)к/п + [в21Р2 + (вз1 + вз1)рз]к2/п/п+ +^3[/321/332р3/^/га + ¿(/32>2 + (&1 + Рз2)2РЗ)Ш2] +

+к4

1

в21вз2рз/п/п' У)2 + в21(вз1 + вз2)вз2рз/п ^^

0(к5).

Здесь элементарные дифференциалы вычислены на приближенном решении уп. Точное решение у(*п+1) в окрестности точки ¿п имеет вид (3), Сравнивая полученные ряды для точного и приближенного решений до членов с кз включительно при условии уп = у(*п), запишем условия третьего порядка точности схемы (13):

Р1 + Р2 + Рз = 1, в21Р2 + (вз1 + вз2)Рз :

в221Р2 + (вз1 + вз2)2Рз

/?21 /З32р3 = 76

1

2

1

Локальную ошибку численной задачи (13) можно вычислить по формуле = у(^п+1) — Уп+1- Учитывая представления у(£га+1) и уга+1 в виде рядов Тейлора, получим

= ¿//3/ +

- ^21/Зз2РЗ

/772 +

^ - - ^(/Зз1 + /З32)3рз

о - /521 (/5з1 + /З32)/З32р3

/7 4 + о(^5).

/772+

(15)

В нелинейной системе алгебраических уравнений (14) имеются два свободных коэффициента, Исследуем три варианта.

Вариант 1. Положим в21 = в31 + в32 и в31 = в32- Это означает, что приращения к2 и к3 будут вычислены в одной и той же точке ¿п + в21^, причем вклад и к2 при определении к3 учитывается одинаково. Тогда нелинейную систему (14) можно переписать в виде

1 2'

1) Р1 + Р2 + Р3 = 1, 2) в21(Р2 + Р3) =

з) /322!(р2 +Рз) = 4) /321/332Рз = ^ 3 6

Из второго и третьего уравнений данной системы имеем в21 = 2/3, Из соотношений в21 = в31 + в32 и в31 = в32 заппшем в31 = в32 = 1/3. Из четвертого уравнения системы получим р3 = 3/4. Из равенства р2 + р3 = 3/4 имеем р2 = 0. Наконец, из первого соотношения системы получим р1 = 1/4. В результате коэффициенты схемы (13) определяются однозначно:

21

021 = д, Рз! = /З32 =

Р1

Р2 = 0, Р3 =

1

4' " 4"

При данных соотношениях локальная ошибка схемы (13) имеет вид

1

(16)

=

216

^4(9/'3/ + 773 — 3/772 + 3/772) + о(^5)

Вариант 2. Минимизируем локальную ошибку (15). Для этого, учитывая вид (15), вместо (14) рассмотрим расширенную нелинейную систему:

1) Р1 + Р2 + Р3 =1, 2) в21Р2 + (в31 + в32)р3

1

з) /321р2 + (/?31 + /З32)2р3 = 4) /321/332р3 =

1

1

5) /ЗДаРз = у^ 6) /321(/331 + /З32)/З32р3 = -.

(17)

При 1.5в21 = в31 + в32 Два последних уравнения (17) совпадают. Из четвертого и пятого соотношений (17) имеем в21 = 0.5. Из второго и третьего равенств получим р2 = 1/3 и р3 = 4/9. Из первого уравнения (17) запишем р1 = 2/9, а из четвертого имеем в32 = 3/4. Наконец, из соотношения в31 + в32 = 3/4 запишем в31 = 0. В результате коэффициенты метода (13) с минимальной локальной ошибкой имеют вид

в:

21

32 031 = 0, /З32 = -, р 1 = -,

14

Р2 = -, Р3 = --

(18)

1

При данных соотношениях локальную ошибку схемы (13) можно записать следующим образом:

^ = ^4(12//3/-//73) + 0(/г5).

При использовании (13) с наборами коэффициентов (16) и (18) ни одна стадия не вычисляется в точке ¿п+1. При быстром изменении решения это приводит к понижению точности расчетов.

Вариант 3. Положим в21 = 0.5 и вз1 + вз2 = 1- Тогда на каждом шаге приращения к1; к2 и кз вычисляются в точках ¿п, ¿п + 0.5к и ¿п + к соответственно, В этом случае условия третьего порядка записываются в виде

11

1) Р\ +Р2 + Рз = 1, 2)-р2+рз = -, 3) ^Р2+Рз = 4) /332Рз =

Из второго и третьего уравнений этой системы имеем р2 = 2/3 и рз = 1/6, из первого и последнего — р1 = 1/6 и вз2 = 2. Из равенства вз1 + вз2 = 1 следует вз1 = — 1В результате коэффициенты метода (13) имеют вид

1 12 1

/321 =/Зз1 = ~1, Рз2 = 2, = Р2=з> = (19)

При данных соотношениях локальную ошибку схемы (13) можно записать следующим образом:

Л ^ 24

//3/ - /772 -

0(к5).

4. Контроль точности и устойчивости

Неравенство для контроля точности вычислений метода третьего порядка построим с использованием идеи вложенных методов. Для этого рассмотрим вспомогательную схему

Уп+1,1 = Уп + Пк! + Г2к2,

к1 к2

точности. Разложение приближенного решения уп+1,1 в виде ряда Тейлора по степеням к к2

Уп+1,1 = Уп + (п + Г2)к/п + в21 г 2 к2 /п /п + 0(кз).

Сравнивая ряды для у(*п+1^ уп+1,1, видим, что требование второго порядка точности будет выполнено, если г1+г2 = 1 и в21г2 = 0.5. Отсюда получим г2 = 0.5/в21 и г1 = 1 — г2, где значение в21 определено в (16), (18) или (19),

С помощью идеи вложенных методов ошибку еп,з метода третьего порядка можно оценить по формуле [13]

£п,з = Уп+1 — Уп+1,1 = (Р1 — Г1)к1 + (р2 — Г2)к2 + рзкз.

Тогда неравенство для контроля точности вычислений запишется в виде

Ь - + (Р2 - Г2)&2 + РзМ < е,

где || • || — некоторая норма в е — требуемая точность интегрирования, В конкретных расчетах как наиболее надежный применялся метод (13) с коэффициентами (19), В этом случае неравенство для контроля точности имеет вид

||^1 - 2^2 + М < 6е.

Неравенство для контроля устойчивости численной формулы (13) построим предложенным в [8] способом. Запишем стадии к1; к2 и к3 применительно к задаче у' = Ау, где А есть матрица с постоянными коэффициентами, В результате получим

¿1 = Xy„, k2 = (X + 2)y„, ¿3 = [X + (вэ! + вз2)Х2 + в21вз2Х где X = hA Найдем коэффициенты d1; d2 и d3 из условия выполнения равенства

diki + d2^2 + d3^3 = X V Легко видеть, что данное требование выполняется при

d = /Зз1 + /З32 ~ /?2i ^ = /З31+/З32 ^ = 1

в21в32 ' 021032 ' 021032

Также видно, что

- fci) = Х2уп.

в21

Тогда согласно [8] оценку максимального собственного числа wn,3 = hAn max матрицы Якоби системы (13) можно вычислить по формуле

w„,3 = в21 max (|d1ki + d2&2 + d3&3l/|k2 - ¿11). (20)

Интервал устойчивости численной схемы (13) приблизительно равен 2,5, Поэтому для контроля устойчивости можно применять неравенство wn,3 < 2.5,

Полученная оценка (20) является грубой, так как: 1) вовсе не обязательно максимальное собственное число сильно отделено от остальных, 2) в степенном методе применяется мало итераций, 3) дополнительные искажения вносит нелинейность задачи (1), Поэтому здесь контроль устойчивости используется как ограничитель на размер шага интегрирования. Прогнозируемый шаг hn+1 будем вычислять следующим образом. Новый шаг hac то точности определим по формуле hac = q1hn где hn есть последний успешный шаг интегрирования, a учитывая соотношение ега,3 = О(Л-П), задается уравнением q3||£n,3|| = £■ Шаг hst то устойчивости зададим формулой hst = q2hn, где q2, учитывая равенство wn,3 = O(hn), определяется из уравнения q2wn,3 = 2.5, В результате прогнозируемый шаг hn+1 вычисляется по формуле

hn+1 = max[hn, min(hac, hst)]. (21)

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

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

На основе построенных L-устойчивого и явного методов третьего порядков точности легко сформулировать алгоритм переменного шага и структуры. Нарушение неравенства wn,3 < 2.5 вызывает переход с явной численной формулы на схему (6), Передача управления с L-устойчивой схемы явному методу происходит в случае выполнения неравенства wn,0 < 2.5, Нижепоетроенный алгоритм переменного шага с автоматиче-

L

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

Компьютерные вычисления проводились на Intel Core 2 Quad CPU с двойной точностью, Расчеты проводились с точностью е = 10-4, при которой наиболее эффективны методы третьего порядка. Ниже через ifu и ija обозначены соответственно суммарное число вычислений правой части и количество обращений (декомпозиций) матрицы Якоби задачи (1), позволяющие объективно оценить эффективность алгоритма интегрирования, В качестве первого тестового примера выбрана простейшая модель реакции Белоусова — Жаботинского [14]

у' = 77.27(y2 - У1У2 + yi - 8.375 ■ 10-6у2), у2 = (-У2- У1У2 + уз)/77.27, у3 = 0.161(yi - уз),

t G [0,300], yi(0) = 4, У2(0) = 1.1, уз(0) = 4, ho = 2 ■ 10-3. (22)

Для данной задачи характерным является наличие на промежутке интегрирования трех участков, на которых решение изменяется быстро. По этой причине шаг выбирается достаточно малым, что позволяет проводить расчеты на данных участках по явной схеме. Вычисления проводились с численной матрицей Якоби, Решение данной задачи алгоритмом MKRK3 вычислено с затратами ifu = 2518 и ija = 411. При расчетах только по L-устойчивой схеме (6) ifu = 2501 и ija = 701. Фактическая точность расчетов в конце интервала интегрирования не хуже задаваемой. Из результатов расчетов следует, что по числу декомпозиций матрицы Якоби построенный алгоритм эффективнее L

части задачи (22) для обоих алгоритмов примерно одинаково.

Решение задачи (22) было вычислено явным методом переменного шага (13) с контролем точности вычислений и устойчивости численной схемы с затратами ifu = 10 497424. При расчетах без контроля устойчивости ifu = 13 250 508. Данная задача слишком

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

В качестве второго примера выбрано известное уравнение Ван дер Поля

у' = У2, у2 = 102[(1 — у2)У2 — ш], I е [0,11],

У1(0) = 2, У2(0) = 0, ^ = 10-6. (23)

Для данной задачи также характерным является наличие участков с быстрым изменением решения, которые можно просчитывать по явной схеме. Расчеты проводились с численной матрицей Якоби, Решение задачи (23) алгоритмом МКНКЗ вычислено с затратами г/и = 19 432 и г^'а = 5010. При расчетах только по ¿-устойчивой схеме (6) г/и = 18 670 и г^'а = 5671. Фактическая точность расчетов в конце интервала интегрирования не хуже задаваемой. Из результатов расчетов следует, что по числу декомпозиций матрицы Якоби построенный алгоритм снова эффективнее ¿-устойчивого метода (6), хотя преимущество не столь значительно. Число вычислений правой части задачи (23) для обоих алгоритмов примерно одинаковое. Решение задачи (23) было вычислено явным методом переменного шага (13) с контролем точности вычислений и устойчивости численной схемы с затратами г/и = 22 030 302. При расчетах без контроля устойчивости г/и = 27 350 638. Задача (23) также слишком жесткая для явных методов.

Заключение

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

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

матрицы Якоби системы (1) осуществляется через ранее вычисленные стадии и не при/

применение контроля устойчивости в качестве ограничителя на рост шага позволяет избежать негативных последствий грубости оценки. Более того, в некоторых случаях это приводит к нестандартно высокому повышению эффективности алгоритма. На участке установления за счет контроля устойчивости старые ошибки стремятся к нулю, а новые невелики за счет малости производных решения, В некоторых случаях вместо оценки максимального собственного числа оценивается следующее по порядку. Шаг интегрирования становится больше максимально допустимого и с таким шагом осуществляется интегрирование до тех пор, пока не нарушается неравенство для контроля точности. Как правило, число таких шагов невелико. Однако величина шага может на порядок превышать максимальный шаг по устойчивости. После нарушения неравенства для контроля точности шаг уменьшается до максимально возможного. Такой эффект может

повторяться многократно в зависимости от длины участка установления, В результате

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

Список литературы

[1] Хайрер Э., Ваннер Г. Решение обыкновенных дифференциальных уравнений. Жесткие и дифференциально-алгебраические задачи. М.: Мир, 1999. 685 с.

[2] Byrne G.D., Hindmarsh А.С. ODE solvers: a review of current and coming attractions // J. Comput. Phvs. 1987. No. 70. P. 1-62.

[3] Современные численные методы решения обыкновенных дифференциальных уравнений / Под ред Дж. Холла и Дж. Уатт. М.: Мир, 1979. 312 с.

[4] Rosenbrock Н.Н. Some général implicit processes for the numerical solution of differential équations // Computer. 1963. No. 5. P. 329-330.

[5] Новиков E.A., Новиков В.A., Юматова Л.A. О повышении эффективности алгоритма интегрирования на основе формулы типа метода Розенброка второго порядка точности за счет замораживания матрицы Якоби. Новосибирск, 1985. (Препр. АН СССР. Сиб. отд-ние. ВЦ. № 592).

[6] Новиков Е.А., Двинский А.Л. Замораживание матрицы Якоби в методах типа метода Розенброка // Вычисл. технологии. 2005. Т. 10. С. 108-114.

[7] Новиков В.А., Новиков Е.А. Контроль устойчивости явных одношаговых методов интегрирования обыкновенных дифференциальных уравнений // Докл. АН СССР. 1984. Т. 277, № 5. С. 1058-1062.

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

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

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

[11] Новиков Е.А. Исследование (m, 2)-методов решения жестких систем // Вычисл. технологии. 2007. Т. 12, № 5. С. 103-115.

[12] Демидов Г.В., Юматова Л.А. Исследование некоторых аппроксимаций в связи с А-ус-тойчивостью полуявных методов // Численные методы механики сплошной среды. 1977. Т. 8, № 3. С. 68-79.

[13] Демидов Г.В., Новиков Е.А. Оценка ошибки одношаговых методов интегрирования обыкновенных дифференциальных уравнений // Там же. 1985. Т. 16, № 1. С. 27-42.

[14] Enright W.H., Hull Т.Е. Comparing numerical methods for the solutions of svstems of ODE's // BIT. 1975. Vol. 15. P. 10-48.

Поступила в редакцию 1 марта 2011 г., с доработки — 26 апреля 2011 г.