Об одном способе конструирования W-методов для жестких систем ОДУ Текст научной статьи по специальности «Математика»

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

Том 12, № 4, 2007

ОБ ОДНОМ СПОСОБЕ КОНСТРУИРОВАНИЯ W-МЕТОДОВ ДЛЯ ЖЕСТКИХ СИСТЕМ ОДУ*

В.А. Вшивков Институт вычислительной математики и математической геофизики СО РАН, Новосибирск, Россия

e-mail: vsh@ssd.sscc.ru

О.П. СтояновскАя Институт катализа им. Г. К. Борескова СО РАН, Новосибирск Новосибирский государственный университет, Россия e-mail: stop@catalysis.ru

A new approach for construction of W-methods for stiff system of ODE is presented. Using the proposed idea, a W-modification of the one-stage Rozenbrok scheme has been developed. Some computational features of the proposed new method are described.

Введение

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

Г y' = F (x,y), \ y(to) = yo, to < t < tk,

имеет следующие тенденции [1]:

— исследование явления жесткости как такового и создание теоретического аппарата анализа устойчивости методов [2];

— конструирование и усовершенствование методов с учетом специфики решаемых задач (например, в методах могут найти отражения такие особенности, как размерность задачи [3], заполненность соответствующей матрицы Якоби [4], априорные свойства решения — знакоопределенность [5], количество переходных фаз [6] и т. д. );

— конструирование методов, перспективных для распараллеливания [7]. Наиболее полный обзор современного состояния численных методов решения обыкновенных дифференциальных уравнений представлен в монографии Э. Хайрера и Г. Ван-нера [2, 8].

* Работа выполнена при финансовой поддержке проекта "Развитие научного потенциала высшей школы" Р.Н.П.2.1.1.1969, Российского фонда фундаментальных исследований (грант № 05-01-00665). ( Институт вычислительных технологий Сибирского отделения Российской академии наук, 2007.

Мы ограничимся рассмотрением группы одношаговых в-стадийных методов вида

Уп+1 = уп + П^Ьгк-

г=1

Здесь Ьг — константы, удовлетворяющие соотношению ^г Ьг = 1, которые наряду с константами, используемыми для вычисления кг, определяют порядок аппроксимации и свойства устойчивости метода. Способ вычисления стадий кг определяет один из возможных классов метода.

Неявные методы Рунге—Кутты имеют вид

кг = ^ уп + к аИ к

\ 3 = 1

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

кг = ^ |V + к == аг3к^ , 1 < г < в,

требуют меньших затрат на матричные вычисления. В случае в = г для нахождения кг необходимо решение нелинейной системы, размерность которой совпадает с размерностью исходной системы ОДУ.

Полунеявные методы типа Розенброка [9] с фиксированным числом обращений матриц на каждом шаге интегрирования имеют вид

— Мцкг = ^ ^уга + к == агзк¿^ + к^у == дуку, 1 < г < в,

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

Тем не менее обращение матриц на каждом шаге по времени представляется слишком трудоемким для систем большой размерности, поэтому в работе [10] был впервые предложен еще один класс методов.

Методы с неточной матрицей Якоби, или Ш-методы:

(г-1 \ г—1

УП + к агзкЛ + кА ^ дцку, 1 < г < в,

3=1 3=1

где W(к, дгг, А) = I — кдггА; А — такая произвольная матрица, что W обратима. Подобный подход позволяет, избегая пересчета частных производных и обращения матриц

на каждом шаге, получить методы второго и более высоких порядков аппроксимации. Недостатки этих методов, получивших развитие в работах [3, 7] и др., связаны со сложностью исследования устойчивости методов, поскольку аппарат функций устойчивости оказывается неприменим (см. разд. 3.2 статьи, где приведены дополнительные теоретические сведения). Это обусловливается невозможностью одновременной диагонали-зации матриц А и С//Су, поэтому тестовое уравнение Далквиста перестает быть приемлемой моделью поведения решения.

Из приведенного обзора ясно, что движение в сторону минимизации затрат на матричные вычисления на каждом шаге по времени сопровождается потерей вычислительной устойчивости методов и приводит к необходимости уменьшения шага интегрирования. Таким образом, для эффективного использования численных методов для ОДУ важен вопрос управления длиной шага интегрирования, этот вопрос разрабатывается вычислителями. Большинство алгоритмов предсказания оптимальной величины шага основано на построении оценок локальной погрешности решения. Кроме традиционно применяемого правила Рунге используется вложенный метод, впервые предложенный в работе Р. Мерсона в 1957 году [8]. Явные методы с расширенными областями устойчивости, например, Рунге—Кутты—Чебышева, которые также применяются для решения жестких задач [11], предполагают дополнительную оценку устойчивости метода и ее учет в алгоритме управления шагом. В работе Л. Шампайна и А. Витта в 1995 году [12] также был подробно рассмотрен простой универсальный метод выбора величины шага, не требующий оценки локальной погрешности решения. Этот метод представляет собой интегрирование по дуге и может применяться в тех случаях, когда нет простого способа оценки локальной погрешности решения. Более того, метод интегрирования по дуге использовался значительно раньше, в работе Ю.А. Березина, В.А. Вшивкова в

К основным особенностям Ш-методов относится их "безматричность". Это означает, что при использовании этих методов интегрирования нет необходимости вычислять, хранить и выполнять преобразования матрицы Якоби системы на каждом шаге по времени. Исторически первый способ реализации Ш-методов состоял в использовании одной и той же обратной матрицы на протяжении нескольких шагов, после чего вычисля-

Тем не менее ясно, что приближение матрицы А к точной матрице Якоби С//Су на каждом шаге по времени улучшает свойства устойчивости методов этого класса и делает их более приемлемыми для решения жестких систем ОДУ. Поэтому идея аппроксимировать матрицу Якоби путем численного дифференцирования правых частей показала свою высокую эффективность. При этом построение аппроксимации матрицы С//Су размерности п х п в виде матрицы А = (Бт ((,Б С Япхк) меньшей размерности к х к, являющейся проектором на соответствующее подпространство Крылова, позволяет существенно экономить вычислительные ресурсы [7, 14].

Но для систем средней и большой размерности, возникающих в химической кинетике, матрица Якоби, как правило, очень разреженная [4]. Поэтому использование в расчетах точной матрицы частных производных не оказывается слишком дорогостоящим в тех случаях, когда нет необходимости ее обращать. В настоящей работе предлагается способ конструирования Ш-методов, который сочетает в себе неявный подход с неполным обращением матрицы. Частичная неявность метода дает возможность применять его для жестких систем. Сократить вычислительные затраты на каждом шаге интегри-

1976 году [13].

лась точная матрица Якоби и находилась точная обратная матрица

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

Идея построения приближения к обратной матрице на каждом шаге по времени и способ ее использования для конструирования W-методов представлен в разд. 1. В разд. 2 описан пример построения численной схемы второго порядка аппроксимации на основе правила средней точки. Подходы к анализу устойчивости метода представлены в разд. 3. В разд. 3.4 приведено сравнение вычислительных характеристик двух методов — метода Розенброка второго порядка аппроксимации и нового W-метода с неполным обращением матрицы, сконструированного на его основе. Важно отметить, что в настоящей работе авторы не ставят целью представить самый экономичный метод интегрирования жестких систем ОДУ. Цель работы — показать, что предложенный способ конструирования W-методов позволяет получить выигрыш по сравнению с базовыми методами Розенброка за счет замены полного обращения матрицы двумя матричными умножениями. Актуальность работы состоит в том, что эта методика обладает потенциалом для дальнейшего повышения эффективности расчетов, в первую очередь, за счет использования ускоренных алгоритмов умножения матриц [15].

1. Новый способ реализации №-методов

Основная цель разработки W-методов — повышение вычислительной устойчивости за

( д/\ —1

счет максимального приближения матрицы W—1 = (I — кдА)-1 к (I — кд— ) при

V дУ)

минимальной сложности вычислений на каждом шаге по времени. В предлагаемом методе построение приближения к обратной матрице опирается на идею, изложенную в работе [16].

Определение 1. Пусть В — (и х и)-матрица, тогда определим Я(В):

Я(В) = I — В А; В называют достаточно хорошей оценкой для А-1, если

||Д(В)|| = и(В) < 1. Теорема 1. Пусть В — достаточно хорошая оценка для А-1, тогда

(I + Е(В) + Я2(В) + ...)В ^ А-1,

т. е. если

Бк = (I + Е(В) + Я2(В) + ... + Як (В))В, то справедливо при к ^ ж

||А—1 9 и < 11В №(В) ^ 0

||А — Бк1 — и{В) ^ 0 Такой способ построения обратной матрицы в виде суммы ряда дает возможность

/ дfп+1 V11

строить приближение к (I — кдгг—-— ) , используя в качестве достаточно хорошего

V ^ )

( дfn^ —1 начального приближения приближение к I — кдц——

V ^

Этот метод оказывается также принадлежащим к классу Ш-методов:

уп+1 = уп + ^Ьгкг,

г=1

(г-1 \ г-1

Уп + о*!кА + кА ^ Сг3к3, 1 < г < в,

3 = 1 ) 3 = 1

Ш(к,Си,А) = I — НСцА,

-1

( С/г' 1

если положить, что приближенное значение ( I — кСц —— ) является (I — кСггАп) .

V СУ )

Опишем предлагаемый способ реализации Ш-методов.

Будем считать, что Сгг = С, и положим, что Вп = (I — кСАп)-1, тогда

(г-1 \ г-1

уп + к ^ Огз кЛ + кАп^ Сгз к3

3=1 3=1

может быть преобразовано к виду

(г-1 \ г-1

уп + к ^ ОгзкЛ + кВпАп ^ Сгзку.

3=1 3=1

При этом нет необходимости вычислять и хранить Ап, поскольку в силу

В— кСАп) = I, кСВ пАп = Вп — I,

получим

кг = ВП (^ (уп + к ^ Огз к3Л + С ^ Сгз к3Л — С ^ С*з ку.

3=1 С 3=1 С 3=1

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

_ ( С/п+1 \-1

(I — кСАп) 1 представляет собой хорошее приближение для I — кС—-— , т. е. име-

Су ет место

\\Яп+1\\ = III — Вп ^I — м/р^ || < 1.

Тогда

-1

— кС*/р =(I + Кп+1 + Я2п+1 + ...)Вп,

положим

Вп+1 = (I + Яп+1)Вп = (^Н — Вп ^I — кС/р ^ В тем самым определив Ап+1: Вп+1 = (I — кСАп+1)-1.

2. Реализация метода второго порядка аппроксимации на основе правила средней точки

Будем исходить из метода Розенброка второго порядка аппроксимации:

у"+1 = у"+^—К Ь)

который получен из модифицированного метода Эйлера (правила средней точки):

уп+1 — уп 1

у-- = 1(/(уп+1) + / (уп))-

т 2

1 этап. Получение схемы с произвольным разбиением матрицы Якоби и проверка ее аппроксимации

Этот этап подразумевает конструирование W-метода на основе метода Розенброка. Стоит отметить, что для применения техники неполного обращения возможно использование уже найденных методов более высокого порядка аппроксимации, например, [7]. Полагаем

/ = яп = яп + яп,

где яп « яп;

к \ уп+1 — уп

I — 2(яп + яп))у кУ = / (уп),

п\ —1

у п+1 — уп 1 1 1

у ку = 2 яп уп+1 — 2 япуп + 2 яп (уп+1 — уп)

2 яп уп+1 — 2 яп уп + 2 яп (уп+1 — уп) + / (уп).

Выполнив разложение

1 яп (уп+1 — уп) = 1 яп (уп + + о(к2) — уп) = 2 яп (к/(уп) + 0(к2)),

получаем

I — к яп) (уп+1 — уп) = [I + к яъ) к/(уп).

Для проверки аппроксимации полученной схемы заметим, что

% = У (у) = = я/(у)'

дЬ2 дЬ ду дЬ

тогда, учитывая

.,п+^„,п, ъ (ду\ , „п , г, ^У) ,П/73ч

у"+ = у"+Ч*; [ж*) +0(к ) = у"+к/> + т {д¥) +0(к ^

получим

I — к яп) (к(ди)П + к* (ду)П + 0(к3)) — к/(уп) — к* яп / (уп) =

к* (д*у\П к* ^пгг, пч ,

т дф - т{у") + 0(к ) = 0(к

2 этап. Более эффективная организация вычислений по схеме с разбиением

( ь \-1

Положим Вп = ( I — —ЯП) — известная матрица, тогда

в пяч = -(вп — I).

Полученная схема

I — ЬЯп) (уп+1 — уп) = (I + ЬЯ2) ! (Уп)

в этом случае перепишется в следующем виде:

уп+ = уп + вп(I + ЬЯ2) !(уп).

Исключим Впд^:

2

в пяп = в п(яп — Яп) = впяп — -(вп — I),

Ь впдп2 = гп к~— ' ' к

таким образом получим формулу

вп + - впяп = вп + - впяп — вп +1 = I + - впяг

уп+ = уп + ^ + Ьвп!(уп). 3 этап. Рекуррентное вычисление матрицы вп

п = 0, в0 = ^ — -Я^ 1 п > 1, вп = (-I — вп-1^I — Ьо ) ) вп~х

3. Результаты численных экспериментов 3.1. Тестовые задачи

В связи с большим интересом к разработке эффективных численных методов для жестких систем к настоящему моменту создано несколько наборов задач, включающих системы ОДУ разной размерности и жесткости, которые могут быть использованы в качестве тестовых для изучения вычислительных особенностей разрабатываемых методов. На страницах http://pitagora.dm.uniba.it/testset/ и http://www.unige.ch/hairer/testset /testset.html, а также в [2] даны описания серий задач-тестов для жестких систем. Результаты, представленные в настоящей статье для иллюстрации особенностей метода, получены на решении двух тестовых задач.

3.1.1. Система Робертсона

Классическим примером жесткой задачи является реакция Робертсона [2]:

А ^ в,

-в ^ с + в, в + с ^ А + с,

с константами скоростей соответственно: к1 = 0.04, к2 = 3 • 107, к3 = 104. Заметим, что задача является жесткой за счет большого разброса порядков в константах скоростей реакции.

Изменение концентраций веществ в ходе реакции Робертсона описывается следующей системой ОДУ:

= —к1у1 + кзу2уз,

^у2 = к1у1 — кзу2уз — Ьу1,

Луз , 2 ИЬ = к2у22

= 1, 0, 0.

Матрица Якоби для этой системы имеет вид

/ -0.04

104уз

V 0

6 ■ 107У2

104У2 \

0.04 -104y3 - 6 ■ 107y2 -104y2

0

Ее собственные числа легко могут быть найдены аналитически.

Можно заметить, что Ai = 0, а два следующих собственных значения — это корни квадратного уравнения:

A2 + (0.04 + 6 • 107y2 + l04y3)A + 6 • 107у2(0.04 + 104у2) = 0.

Данная задача не имеет аналитического решения. "Точное" решение, которое можно использовать в процессе тестирования, было получено методом RADAU 5 с относительной и абсолютной локальной погрешностью Ю-18 (взято с http://www.unige.ch/hairer/ testset / testset. html).

3.1.2. Система уравнений HIRES

Система ОДУ HIRES [2] также взята из химической кинетики, она описывает химические превращения восьми реагентов:

dy f( , dt = f (У), У(0) = Уо,

у е Я8, 0 < г < 321.8122,

причем

I (у) =

/ -1.71у1 +0.43у2 +8.32у3 +0.0007 \

1.71у1 -8.75у2

-10.03уз +0.43у4 +0.035у5

8.32у2 1.71уз -1.12у4

-1.745уь +0.43у6 +0.43у7

-280у6у8 +0.69у4 +1.71у5 -0.43у6 +0.69у7

280у6у8 -1.81у7

V -280у6у8 + 1.81У7 /

Уо = (1, 0, 0, 0, 0,0, 0, 0.0057)т.

3.2. Особенности вычислительной устойчивости

Приведем ряд определений устойчивости, которые используются при исследовании линейной устойчивости численных методов для систем ОДУ [2].

Они вводятся на линейном скалярном уравнении, называемом тестовым уравнением Далквиста:

У' = Ау, , у (го) = 1,

где А — комплексное, И,е(А) < 0.

Определение 2. Функцией устойчивости метода Я(г) называется численное решение после одного шага тестового уравнения Далквиста:

уп+1 = Я(г)уп, г = АН.

Например, для явного метода Эйлера

уп+1 — уп

= АУП,

уп+1 = (1 + АН)уп функция устойчивости будет иметь вид

Я(г) = 1 + г.

Определение 3. Метод называется абсолютно устойчивым для г = АН, если Я(г) < 1. Область Я комплексной плоскости называется областью устойчивости метода, если он устойчив при всех г е Я.

Пересечение области устойчивости с вещественной осью называется интервалом устойчивости.

Определение 4. Метод называется А-устойчивым, если его область абсолютной устойчивости включает всю полуплоскость Ке(АН) < 0.

Определение 5. Одношаговый метод называется Ь-устойчивым, если он А-ус-тойчив и если Б,е(АЬ) ^ 0 при |АН| ^ то.

Однако, как уже было отмечено во введении, аппарат функций устойчивости, разработанный для исследования методов Рунге—Кутты, неприменим к исследованию линейной устойчивости Ш-методов. Поэтому, как правило, ограничиваются рассмотрением предельного свойства устойчивости, когда матрица А или Q1 представляет собой точную матрицу Якоби системы [7]. Рассматриваемый метод имеет предельную А-устойчивость [2], но не обладает предельной ¿-устойчивостью.

Более того, при использовании метода необходимо иметь представление о его "внутренней устойчивости": вп = (I — кд,Ап)-1 представляет собой хорошее приближение

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

Схема проведения вычислительных экспериментов

1. Выбирается достаточно мелкий шаг кьавю (порядка 10-6 ... 10-4).

2. Методом с неполным обращением матрицы рассчитываются решение и матрица в с указанным шагом. При этом на каждом шаге проводятся следующие оценки.

3. Вычисленное значение фиксируется как начальное.

4. Выбирается более мелкий (внутренний) шаг к интегрирования (порядка 10-6), запускается цикл кпе„ = гк.

4.1. Из начального значения с внутренним шагом к рассчитывается значение у1, по у1 строится Q1 и вычисляется (I — кQ1)-1.

4.2. По матрице в вычисляется матрица в1 с внутренним шагом к:

4.3. Вычисляются спектральная норма матрицы и максимум суммы абсолютных значений элементов матрицы по столбцам

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

На рис. 1 и 2 слева представлена зависимость неспектральной нормы матрицы В от внутреннего шага к для ряда внешних шагов кьавю = 10-6 при решении системы уравнений Робертсона. Видно, что существует участок решения, на котором немонотонность величины ||О| наиболее ярко выражена. На этом участке применение традиционного алгоритма предсказания длины шага по правилу Рунге окажется неэффективным, поскольку уменьшение шага интегрирования может сопровождаться ослаблением внутренней устойчивости метода. Стоит отметить, что этот "аномальный" участок решения представляет собой переходную фазу, где резко изменяется вторая компонента решения системы.

В

в1 — (I — кQl)-1

3.3. Замечания об алгоритме управления длиной шага

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

1. Задаем начальный шаг h.

2. Вычисляем решение yl с шагом h.

2.1. По найденному y\ вычисляем Qn+l и делаем оценку:

stabi = ||I - Bn ^I — hQn+l

2.2. Вычисляем решение y2 из того же начального значения, но с шагом 0.5h.

2.3. Аналогично 2.1 вычисляем stab2l, stab22.

2.4. Если stab = max(stabl, stab2l, stab22) > 1, то переходим к 1, полагая, что

hnew = 0.7h,

иначе переходим к 2.5.

2.5. Оцениваем локальную ошибку

Ь2 — yi\

err =

3

2.6. Если err < tol, где tol — значение допустимой локальной погрешности, то начальный шаг h считается принятым, переходим к вычислению следующего шага из решения yl, считая, что h = hnew, предварительно полагая facmax = min(1.1,1 + (1 — stab)a), в противном случае повторяем вычисления из того же начального значения с шагом h = hnew, при этом производя переприсваивание facmax = 1.

2.7. Вычисляем оптимальный шаг

hnew = h min ( facmax, max ( facmin, 0.7\ —

err

Значения facmin и коэффициента запаса (0.7) выбираются произвольно, таким образом, чтобы каждое из них было меньше единицы. В программе использовались соответственно 0.3 и 0.7. При этом a — параметр сглаживания, выбирается большим единицы. Для системы Робертсона использовался a = 1.3, для системы HIRES — a = 1.8.

3.4. Сравнительные оценки вычислительной эффективности

Чтобы показать работоспособность сконструированного метода и выявить его вычислительные особенности, было проведено сравнение времени счета исходного метода типа Розенброка с предложенной модификацией метода (неполным обращением матрицы). Авторы ограничились состязанием двух методов "в легком весе", поскольку предложенный способ конструирования Ж-методов позволяет получить выигрыш по сравнению с базовыми методами Розенброка за счет замены полного обращения матрицы двумя матричными умножениями.

В табл. 1 представлены результаты работы обоих методов с постоянным шагом интегрирования. Время счета метода с неполным обращением матрицы оказалось на 10 % меньше, чем время счета исходного метода.

Рис. 3. Поведение принятых по внутренней устойчивости шагов Н на участке [0,1] при решении системы Робертсона: графики приведены для одной (а), трех (б) и пяти (в) итераций; отброшенные по точности шаги представлены отрицательными

Таблица 1. Время счета, с, на Pentium 4 CPU 2.6 GHz, 512 MB, методов неполного обращения матрицы и правила средней точки системы HIRES с постоянным шагом интегрирования, 1 000 000 шагов

Метод с неполным Метод Розенброка

обращением матриц (правило средней точки)

5.56 6.15

100% - 111%

Таблица 2. Время счета, с, на Pentium 4 CPU 2.6 GHz, 512 MB, метода с неполным обращением матрицы при решении задачи Робертсона на отрезке [1,10] c автоматическим выбором шага

Число Время, Принято Отброшено шагов

итераций с шагов Внутренняя неустойчивость Точность

1 0.98 508 176 295

2 0.55 231 172 103

3 0.40 126 38 11

4 0.30 119 12 0

Отметим, что при использовании большего числа итераций

Вп,1 = (21 - Вп-Ц I - h Qn

В

п- 1

k > В

n,k+1

2I - В

n,k

i - hQ

В

n,k

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

Вместе с тем увеличение числа итераций приводит к возрастанию вычислительной работы, связанной с преобразованием матриц, но позволяет сократить число шагов. Это проиллюстрировано рис. 3 и табл. 2. В табл. 2 представлены время счета и число шагов, выполненных предложенным методом для решения задачи Робертсона на отрезке [1,10]. На рис. 3 представлено поведение принятых по внутренней устойчивости шагов к на участке [0,1] при решении системы Робертсона для одной, трех и пяти итераций. Видно, что менее затратным с вычислительной точки зрения оказывается применение на некоторых участках решения большего числа итераций.

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

Заключение

Предложен новый способ реализации Ш-методов с использованием рекуррентного построения приближения к обратной матрице.

На основе правила средней точки сконструирован одностадийный метод второго порядка аппроксимации для решения жестких систем ОДУ. Исследованы свойства устойчивости предложенного метода.

Для реализации управления величиной шага интегрирования модифицирована формула предсказания длины шага, в которую включена оценка внутренней устойчивости метода.

Авторы выражают благодарность В.Н. Снытникову за поддержку и постоянное внимание к работе.

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

[1] Butcher J.C. Numerical methods for ordinary differential equations in the 20th century // J. of Comput. and Appl. Math. 2000. N 125. P. 1-29.

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

[3] Weiner R., Schmitt B.A., Podhaisky H. ROWMAP — a ROW-code with Krylov techniques for large stiff ODEs // Appl. Numerical Math. 1997. N 25. P. 303-319.

[4] Nejad L.A.M. A Comparison of stiff ode solvers for astrochemical kinetics problems // Astrophys. and Space Sci. 2005. N 299. P. 1-29.

[5] Bertolazzi E. Positive and conservative schemes for mass action kinetics // Computers Math. Applic. 1996. Vol. 32, N 6. P. 29-43.

[6] Novati P. An explicit one-step method for stiff problems // Computing. 2003. N 71. P. 133-151.

[7] Podhaisky H., Schmitt B.A., Weiner R. Design, analysis and testing of some parallel two-step W-methods for stiff systems // Appl. Numerical Math. 2002. N 42. P. 381-395.

[8] Хайрер Э., Нёрсетт С., Ваннер Г. Решение обыкновенных дифференциальных уравнений, нежесткие задачи. М.: Мир, 1990.

[9] КАлиткин Н.И. Численные методы решения жестких систем // Математическое моделирование. 1995. Т. 7, № 5. С. 8-11.

[10] Steihaug T., Wolfbrandt A. An attempt to avoid exact jacobian and nonlinear equations in the numerical solution of stiff differrential equaion // Math. of Comput. 1979. Vol. 33, N 146. P. 521-535.

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

[12] Shampine L.F., Witt A. A simple step size selection algorithm for ODE codes // J. of Comput. and Appl. Math. 1995. N 58. P. 345-354.

[13] Березин Ю.А., Вшивков В.А. О критических параметрах ударных волн в плазме // Прикл. математика и техн. физика. 1976. № 2. C. 27-36.

[14] Schmitt B.A., Weiner R. Matrix-free W-method using a multiple Arnoldi iteration // Appl. Numerical Math. 1995. N 18. P. 307-320.

[15] Smith J.R. The Design and Analysis of Parallel Algorithms. N.Y.; Oxford: Oxford Univ. Press, 1993.

[16] ВЕржвицкий В.М. Численные методы, линейная алгебра и нелинейные уравнения. М.: Высшая школа, 2000.

Поступила в редакцию 2 февраля 2007 г., в переработанном виде —18 апреля 2007 г.