Аддитивные алгоритмы решения жестких систем на основе (m,k) - методов тема автореферата и диссертации по математике, 01.01.07 ВАК РФ

На нрапах рукописи

Тузов Антон Олегович

АДДИТИВНЫЕ АЛГОРИТМЫ РЕШЕНИЯ ЖЕСТКИХ СИСТЕМ НА ОСНОВЕ (М,К)~МЕТОДОВ

01.01.07 - вычислительная математика

АВТОРЕФЕРАТ

Диссертации на соискание ученой степени кандидата ф из и ко - математических наук

Красноярск 2007

оозоевзег

003066562

Работа выполнена в Институте вычислительного моделирования СО РАН и Сибирском федеральном университете (г. Красноярск)

Научный руководитель доктор физико - математических

наук, профессор

Новиков Евгений Александрович

Официальные оппоненты доктор физико - математических

наук, профессор Задорин Александр Иванович

доктор физико - математических

наук, профессор Добронец Борис Станиславович

Ве;\ущаи организации Институт вычислительной математики

и математической геофизики СО РАН. (г. Новосибирск)

Защита состоится * ноября 2007 года в 15 часов на заседании диссертационного совета К 212,098.03 в Красноярском государственном техническом университете по адресу: 660074, г.Красноярск, ул. Киренгкого, 26.

С диссертацией можно ознакомиться п библиотеке Политехнического института Сибирского федерального университета

Автореферат разослан сентября 2007 года.

Ученый секретарь диссертационного совета ' . г у,-^' Сафонов К.В.

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

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

Вопросы построения аддитивных методов на основе (т, к) - схем и явных методов типа Рунге - Кутта, ¿-устойчивость внутренних схем и контроль устойчивости явной части метода применительно к аддитивным задачам, а также численное исследование ад дитивных методов с диагональной аппроксимацией матрицы Якоби ранее не исследовались

Цель работы - исследовать методы решения аддитивных жестких задач, построенных на основе (ш, к) - методов и явных методов Рунге - Кутта

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

1 Cooper, G J Additive Runge-Kutta methods for Stiff Ordinary Differential Equations/ G J Cooper, A Sayfy // Mathematics of Computation - vol 40 - №161 - 1983 - P 207-218

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

• Построены новые методы решения аддитивных автономных и неавтономных задач на основе (т, к) - методов и явных методов Рунге - Кутта

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

• Созданы программные реализации алгоритмов переменного шага с диагональной аппроксимацией матрицы Якоби

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

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

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

Апробация работы.

Результаты диссертации докладывались на следующих конференциях

• Международная конференция "Вычислительные и информационные технологии в науке, технике и образовании" (Павлодар, 2006)

• Всероссийская конференция молодых ученых по математическому моделированию и информационным технологиям (Красноярск, 2006)

• Всероссийская студенческая научная конференция "Студенческая наука-взгляд в будущее" (Красноярск, 2006)

• Семинар ИВМ СО РАН (Красноярск, 2007)

• Семинар кафедры прикладной математики Политехнический института Сибирского федерального университета (Красноярск, 2007)

Результаты работы опубликованы в 9 работах, из них три в журналах, включенных в перечень ВАК

Структура и объем работы

Диссертация состоит из введения, 3 глав, заключения и одного приложения, содержит 3 таблицы Библиографический список включает 162 наименования Объем работы составляет 105 листов

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

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

У1 = ¥>(*, У) + у)> о) = Уо, t0<t< tk, (1)

в которых правая часть разбита на нежесткую <p(t, у) и жесткую g(t, у) части В основу предлагаемых алгоритмов интегрирования положены (т, к) - методы и явные схемы типа Рунге - Кутта Класс (т, к) - методов, предложенный в работах Новикова Е А , вводится следующим образом Задаются целые числа т и fc, к < го, и рассматриваются следующие множества

Мт = {1,2, ,т},

Мк = {т, е Мт|1 = mi < т2 < < тк < т},

J, = {т3 - 1 е Mm\j > \,т, € Мк,т} < г} 2 < г < т

Тогда (т, к) - схемы для решения задачи У' = f{t, У), y(ta) = У<з, ta<t< tk, имеют вид

Уп+i -Уп+ Pih +р2к2+ + Рткт, Dn = Е- ahfn,

DnK = hf{yn + ДА) + ]C avk» 1 e Mk' W

j=i

Dnkt = кг-1 + г e Mm \ Mki

J€J,

где E - единичная матрица, f'n = df(yn)/dy - матрица Якоби функции f,h- шаг интегрирования, к,, 1 < г < т, - стадии метода, a,p„0i3,atJ - вещественные константы, определяющие свойства точности и устойчивости (2) Вычислительные затраты на шаг интегрирования в методах (2) следующие - один раз вычисляется матрица Якоби и осуществляется декомпозиция Dn, к раз вычисляется функция /, т раз осуществляется обратный ход метода Гаусса (т, к) - методы так же просты в реализации, как и методы типа Розенброка, однако свойства точности и устойчивости у них лучше

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

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

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

в

Уп+1 =Уп + .

2=1

= МЫ,

Опк2 = Ь<р{уп) + Нд{уп), Опк3 = к2,

з з

к4 = к(р(уп + Д»А) + кз(Уп + «4А)> (з)

А>А% = к4 +

5

ко = +

3=1

Получены коэффициенты схемы (3)

«41 = Аи = /?42 = /?61 = о, «42 = /?62 = 0, а43 = 1 - а,

Р1 = -(1 - 2а)/«, Р2 = (2о2(77-1)-а(77-3)+7)/(6аи), р3 = (—2а2(87 4-1) + а(137 +1) - 2у)/(6аи),

р4 = (-6а2 + 2а(7 + 5) - 1)/(6м), р5 = О 5(2а2 - 4а + 1)/«, (4)

Ре = (1 - 2а)/и, /34з = 1/(6ор5), А = арв/рв.

= 2(1 — /?|3)/(3 — 2/?4з), /?64 = а2/(1-2а), Ди = А - Ди, А» = 0з-а- - (7 + 1)Аи,

при которых она имеет третий порядок точности и обладает Ь - устойчивостью основной и промежуточной численных схем относительно жесткой части системы (1) Всюду далее под Ь - устойчивостью понимается Ь - устойчивость относительно жесткой части В (4) и = 0(7 - 1) + 1, а коэффициенты а и 7 определяются из уравнений

6а3 - 18а2 + 9а - 1 = О,

4а373 - 12а(15а2 - 10а + 1)72 + 3а(374а2 - 228а + 33)7 + 4(813а2 - 486а + 175/3) = О Для практических вычислений рекомендуются коэффициенты при а = +0 43586652150845

«41 = #41 = ßi2 = #61 = О,

а = +0 43586652150845 10°, 7 = -0 45174528144973 • 10\

Pi = +0 91301462909303 Ю-1, р2 = +0 49588787677187 10°,

Рз = +0 75521774748194 10°, р4 = +0 20395977226114 10°,

Ps = +0 12937356107219 10°, рв = -0 91301462909303 Ю-1, (5)

а42 =+0 43586652150845 10°, а4з = +0 56413347849155 10°,

#43 = +0 29556275399511 101, #62 = +0 43586652150845 10°,

#вз = —0 39848721470948 101, #64 =+0 14811267768433 10\ #65 = -0 20987467167967 101,

Контроль точности здесь и для всех построенных ниже методов осуществляется посредством неравенства

где уп,2 - приближение к решению, полученное вложенным методом второго порядка точности, £ - требуемая точность интегрирования, г - положительный параметр Если Ij/JjI < г, то по г-й компоненте решения контролируется относительная ошибка е, в противном случае контролируется абсолютная ошибка г s

Для контроля точности (3) построена вспомогательная численная схема второго порядка точности, не требующая дополнительных вычислений правой части и обращений матрицы Якоби

з

У11+1, 2 = Уп + г'кг + г4&4> »=1

h = h<fi(yn),

Dnk2 = h(v{yn)+g(yn)), (6)

Dnks = k2,

3

¿4 = hvivn + ^ßijkj)

1

Коэффициенты вспомогательной схемы имеют вид ri = -0 b/ßa, r2 = 0 5(4а - l)/a, 7-3 = 0 5(1 - 2а)/а, п = 0 5/#43

Полагаться на то, что вся жесткость сосредоточена в функции g(t,y), a <p(t, у) есть нежесткая часть, вообще говоря, нельзя Поэтому все построенные в данной работе методы оснащены неравенством для контроля устойчивости явной части численной схемы

1 14 - d'il ^ о

!-г max —Г77 < 2,

|а32| !<•<* К - Ml -

где di = hip(yn + a2iki), d2 = htp(y„ + a3ih + a32di), a2i - a3i + a32, а числом 2 приблизительно ограничен интервал устойчивости явного метода типа Рунге-Кутта третьего

порядка точности В конкретных расчетах были выбраны «21 = 2/3, а31 = а32 = 1/3 Неравенство для контроля устойчивости грубое, поэтому оно используется как ограничитель на рост величины шага интегрирования Прогнозируемый шаг кп+1 по точности с ограничением по устойчивости вычисляется по формуле

/1п+1 = тах{Л„, тт{/гасс, Л^}},

где касс ~ шаг, выбранный исходя из требования точности, /га( - шаг, выбранный исходя из требования устойчивости

Для улучшения свойств устойчивости схемы (3) рассмотрена четвертая стадия в виде

з

з

АЛ = Му» + X) + +а4Л)

В результате новая схема с улучшенными свойствами устойчивости имеет вид

б

кх = Ыр(уп), Бпк2 = Ь,<р(уп) + кд(уп), А^з = кг,

з

з

АЛ = 1ир(уп + р^к,) + 1гд(уп + ^ а4зк,),

(7)

ДА = к4 + укз,

5

Получены коэффициенты схемы (7) третьего порядка точности, Ь - устойчивой и с Ь -

а41 = /?41 = /?в1 = /?ю = О,

а42 = Аг = о, «43 = 1 - а, рг = а, 7 = 2а(о + 1)/(6а3 - 18а2 + 9а - 1), Рз = (а2 — 4а/3 + 1)/(1 — а), р4 = (6а3 - 20а2 -1- 11а - 1)/(6а - 6а2), р5 = (6а3 - 18а2 + 9а - 1)/(6а2 - 6а), А = (а- 1)/(6а3 - 16а2 + 7а - 1),

Аз = А - а, (8)

А = (1-Д?)/(15-А), Рб = (0 5 - А/3)/А,Р1 = -Рв, А = (6АРв)-1,

А = [1/6 - а(2А - А2)/3]/рв, 0№ = [в(А - 2А) + А - А]/(°7 + а), Аз = А - А - 7А5, А4 = А — А51

где коэффициент а определяется из уравнения

24а4 - 96а3 + 72а2 - 16а + 1 = О

Для практических вычислений рекомендуется комплект коэффициентов при значении а равном 0 57281606248213, то есть

«41 = Ра = А1 = Аг = О, а = +0 57281606248213 10°, рх = -0 48695861160292 10°, Рз = 0 13211252622010 101, Ре = 0 42438423735834 10°, а42 = 0 57281606248213 10°, 042 ~ 0 57281606248213 10°, Аз = 0 25149936861896 101, Аб = 0 91371881359681 10°

7 = -0 28918950092395 101,

р2 = +о 57281606248213 10°,

р4 = -0 91050904025002 • Ю-1,

рв = +0 48695861160292 10°, (9)

а4з = +0 42718393751787 10°,

Аз = -0 18882050162851 10°,

0и = -0 22405291307056 Ю-1,

4 _

Уп+1, 2 = Уп + 12 + Г5*5 ,

1=1

к1 = Ьф(уп)>

АЛ = МйО + Ьд(уп),

АЛ = Аа, (10)

3 , 3

АЛ = МУ» + АЛ) + + X)

3=1 3=1

АЛ =

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

П = 0, г2 = о, гз = 1 —в —г4 = 2 - а+ (г/ — 0 5)/а, г5 = ь-г4, (11)

где V = 0 5/04 На каждом шаге интегрирования метод (10) требует дополнительно один обратный ход в методе Гаусса и не требует дополнительных вычислений правой части Поскольку на каждом шаге интегрирования в основном методе (7) производится одна декомпозиция матрицы Якоби и четыре обратных хода в методе Гаусса, то для задач большой размерности увеличение вычислительных затрат, связанных с использованием вложенного метода, будет незначительным

В методах (3),(7) в четвертой стадии функции <р и д вычисляются в произвольных, не связанных друг с другом точках (уп + рА]к3) и (уп + оц}к}), то есть (р(уп + /Мч) и з{Уп + а4)к}) Обозначив

з з

04 = ОЦ = У^Оаз.

3=1 3=1

эти точки можно интерпретировать следующим образом

з з

у(и + /?4Л) « Уп + + &Уп + ¿>А>

3=1 3=1

где у{1„ + /?4/г) и у{гп 4- а4й) есть точное решение в точках (¿п + /?4Л) и {1п + оцк), соответственно Точное решение удовлетворяет соотношению у1 = ц>(у) 4- д{у), в котором обе функции вычисляются в одной и той же точке Естественно потребовать выполнения того же свойства и от приближенного решения, то есть положить £'=1 Р*з — ац, В работах Соорег'а это условие называется условием согласованности

Рассмотрим численную схему вида (7), но при определении числовых коэффициентов учтем условие согласованности Коэффициенты, при которых данная схема имеет третий

порядок точности и являться £ - устойчивой вместе с промежуточными схемами, определяются так

а« = А1 = /?б1 = Аг = 0, А = 2/3, »42 = Аг = Р2 = О, а43 = Аз = 2/3 - а, 7 = (4а2 — 2а — 1)/(1 — За), « = (7 + 1)/(3(1 - а)7), (12)

р4 = (6а - 1)/(4а), р5 = 3/4-р4 рз = 1/4 - а - 7р5, р6 = 1/(4«), рх = -р6, А5=-1/т, Аз = 1-й, А4 = Щ-Аб,

где « = (7 + 1)/(3(1 - 0)7), а коэффициент а определяется из уравнения 4а2 - 9а + 3 = 0 Для практических вычислений рекомендуются коэффициенты при а = (9 — \/33)/8, то есть

ощ = А1 = А1 = Аг = 0, а = +0 40692966918275 10°, рх = -0 37323757000745 10°, р3 = +0 55049743857359 10°, р5 = -0 13564322306092 10°, а42 = +0 40692966918275 10°, /342 = +0 40692966918275 10°, Аз = +0 33018532942702 10°, Аб = -0 19174162478890 10°

7 = +0 52153516540863 101,

р2 = +0 4069296691827 10°,

р4 = +0 88564322306092 10°,

р6 =+0 37323757000745 10°, (13)

«43 = +0 25973699748392 10°,

Аз = +0 25973699748392 10°,

/0в4 = +0 86155629536189 10°,

Для контроля точности воспользуемся вложенной схемой вида (10) Коэффициенты, при которых вложенная схема имеет второй порядок точности и являться Ь - устойчивой вместе с промежуточными схемами, задаются в (11)

Из анализа результатов расчетов с диагональной аппроксимацией матрицы Якоби следует, что алгоритм, построенный с учетом условия согласованности, является наиболее эффективным среди всех рассматриваемых методов Поэтому он рекомендуется для решения жестких задач с диагональным преобладанием в матрице Якоби, а во всех построенных ниже методах учитывается условие согласованности В случае диагональной аппроксимации вычислительные затраты на шаг интегрирования построенных алгоритмов практически не отличаются от соответствующих затрат явных методов Поэтому сравнение эффективности проводилось на задачах средней жесткости с известными явными методами Мерсона, Фельберга и Дорманда-Принса Из результатов расчетов следует преимущество построенных алгоритмов по вычислительным затратам более чем на порядок Третья глава посвящена исследованию двух семейств численных методов решения задачи Коши для жестких аддитивных неавтономных систем При построении методов учтены условия согласованности Построены алгоритмы интегрирования переменного шаг га с контролем точности вычислений Приведены результаты численных экспериментов

Исследована шестистадийная численная схема вида

Уп+i = Уп + ^2Piku !=1

кг = hip(tn,yn), Dk2 = h(p(tn, уп) + hg(t„, уп), Dk3 = hi,

з з

ki = h<p{tn + c4/i, Уп + ^2 043k}) + hg(tn + cbh, yn + a^k,),

j=i

?=i

Dk 5 = fc4 + 7&з,

= hip(tn + c6/i, j/n +

j=i

Получены коэффициенты схемы (14)

(14)

7 = 4Й-3, «41 = 0, /?4i = /?6i = 0, р2 = (За-1)/(2а), «42 = о, /?42 = 0, /?б2 = -а(3а - 1) г, р3 = -(9а - 5)/(4а), »43 = 2/3 -а, Ди = 2/3, /?ез = (15а2 - 10а + 1)/2 г, р4 = (За — 1)/(4а), с4 = сб= 2/3, /364 = 3а2/2г, р5 = 1/(4а), /Зб5 = -(6а - 1)/4 г, Pt = -Pl=T, с6 = 0

(15)

при которых она имеет третий порядок точности и обладает Ь - устойчивостью основной и промежуточной схем относительно жесткой части В (15) г = 0 25/(/?б4 +/?ея). г = 1/ [(6а2-6а 4- 1)т], а коэффициент а определяется из уравнения 6а3 — 18а2 + 9а — 1 = 0 Для практических вычислений рекомендуется коэффициенты при а равном 0 43586652150845, то есть

«41 = All = Р42 = Pei = се = 0,

а = +0 43586652150845 10°,

Pi = -0 10674632795531 101,

Рз = +0 61785045034886 10°,

р5 = +0 57357009006977 10°,

а42 = +0 43586652150845 10°,

/?43 = +0 66666666666667 10°,

Да = +0 50155967944070 10°,

/?б5 = +0 79584003971969 10°,

с4 = с5 = 2/3,

7 = -0 12565339139662 • 10\ р2 = +0 35285981986046 10°, р4 = +0 17642990993023 10°, р6 = +0 10674632795531 101, а4з = +0 23080014515822 10°, /?62 = +0 26424022694610 10°, /364 = - + 56163994610649 10°,

4

Уп+1, 2 = Уп + r'fc' + г=1

fa = htp(tn, уп), Dnk2 = htp(tn, уп) + hg(tn, уп),

Dnh = fe, (I7)

з з

fa = h<p(tn + c4h, Уп + ^2 04]k3) + hg{t„ + c6h, + <*4зк3), 3=1 1=1

Dnkb = k4 + 7&з, = fa,

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

3 [4 а2 — 3 о (7 + 4) Ч- 7 4- 4]

п =0, г7 = - , , At—> i-4 = 0 75 - г7,

4 [о (3 7 - 1) + 7 + 4]

П = -1 - (37+2)7-7, г2 = 0 25 — гз — уг7

На каждом шаге интегрирования метод (17) требует дополнительно один обратный ход в методе Гаусса и не требует дополнительных вычислений правой части

Для улучшения свойств устойчивости схемы (14) четвертая стадия изменена следующим образом

з з

Dnfa = hifi{tn + c4h, уп + Pijkj) + hg(tn + c5h, + a4Зк3)

з=1 3=1

В результате новая схема записывается в виде б

Уп+1 = Уп + У^ргкг,

г=1

h = h(p(tn,yn), Dfa = h<p(tn, yn) 4- hg(tn, yn),

Dk3 = fa, (18)

3 3

Dfa - hp(tn + c4h, Уп + ^2 Pijk3) + hg(tn + c5h, yn + ^2 a^k,), 3=1 3=1

Dfa — fa + 7 fa,

5

fee = hp{tn + c$h,yn + ^2(k]k}), 3=1

Коэффициенты, при которых схема (18) имеет третий порядок точности, является Ь -устойчивой и имеет Ь - устойчивые внутренние схемы, вычисляются по формулам

«41 = /?41 = /?61 = Сб = О,

с4 = с5 = 2/3, /?42 = а, /З43 = 2/3 - о, 7 = (12о3 - 48а2 + 28а - 4)/(За - I)2, «42 = (~а( 37 + 4) + 7 + 4)/3, «43 = (а(37 + 4)-7-2)/3, Р! = (7(3а2 - 4а + а42))/(4(о7 + а42)),

р2 = (а(-37 + 2) + 7 + 1)/(4о), (19)

Рз = (а(б7-1)-27-1)/(4а), Р4 = 3/2 — 1/(4а), Рб = -3/4 + 1/(4»), р6 = -Рь

062 = -1 — 1/7 _ (3(а— 1))/(4р6),

063 = 2 + 1/7 + (За — 4)/(4рб),

064 = 1/7+ 1/(4рб), 0б5 = -1/7>

а коэффициент а определяется из уравнения а4 - 4а3 + За2 - ^а + ^ = О

Для практических вычислений рекомендуются коэффициенты при а = 0 57281606248213, то есть

«41 = 041 = 061 = Сб = 0, а = +0 57281606248213 10°, Р1 = -0 11816662081845 10°, р3 = -0 24546934012082 10\ р5 = -0 31355972471041 • 10°, а42 = +0 12448344453974 101, 042 = +0 57281606248213 10°, 062 = +0 20659784402781 101, 064 = +0 17610018794415 101,

с4 = с5 = 2/3,

7 = -0 28196432554463 Ю1,

р2 = +0 18205668382489 101,

р4 = +0 10635597247104 10\

Ре = +0 11816662081845 10°, (20)

а43 = -0 57816777873078 10°,

043 = +0 93850604184537 Ю-1,

063 = -0 31816351144595 101,

0в5 = +0 35465479473988 10°,

4

2/71+1, 2 = Уп + +

1=1

Опк2 = Ь.1р(1п, уп) + кд(гп, уп),

АЛ = Ь, (21)

3 3

Дгк4 = + с4Л, у„ + ^2 04}к3) + К д{гп + с5/г, уп + ацА),

7=1 7=1

ДА = £4 + 7/1:3, =

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

3 [4 о2 — 3 а (7 + 4) + 7 + 4] п __

Г1 = 0, г7 = 1 — А^ Л ■ т-4 = 0 75 - гг, 4 [а (3 7 - 1) + 7 + 4]

г3 =-1 - (37 + 2)Г7, г2 = 0 25 - Г3 - ~/г7

На каждом шаге интегрирования метод (21) требует дополнительно один обратный ход в методе Гаусса и не требует дополнительных вычислений правой части В заключении сформулированы основные результаты В приложении приведены использованные в расчетах тестовые примеры

ОСНОВНЫЕ РЕЗУЛЬТАТЫ ДИССЕРТАЦИИ

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

• Исследованы три семейства методов решения жестких автономных и два семейства методов решения жестких неавтономных задач

• Получены коэффициенты пяти Ь - устойчивых шестистадийных численных схем третьего порядка точности

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

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

Работа инициирована член-корреспондентом РАН Шайдуровым В В и выполнена при финансовой поддержке грантов РФФИ №06-08-920, №05-01-00579, НШ-3428 2006 9

Публикации по теме диссертации

По теме диссертации опубликовано 9 работ Основные результаты диссертации опубликованы в следующих статьях

1 Тузов, А О Шестистадийный метод третьего порядка для решения аддитивных жестких систем / Б А Новиков, А О Тузов // Сиб журн вычисл математики / РАН Сиб отд-ние - Новосибирск, 2007 - Т 10, № 3 -С 307-316

2 Тузов, А О Неоднородный метод третьего порядка для аддитивных жестких систем / Б А Новиков, А О Тузов // Математическое моделирование - Москва, 2007 - Т 19, № 6 - С 61-70

3 Тузов, А О Численный метод третьего порядка точности для решения автономных аддитивных жестких систем /АО Тузов // Вестник КрасГАУ -выпуск 14 - Красноярск, 2006 - С 467-473

4 Тузов, А О Метод третьего порядка для автономных аддитивных жестких систем / Е А Новиков, А О Тузов // Вычислительные и информационные технологии в науке, технике и образовании Труды международной конференции - Павлодар 2006 -С 77-84

5 Тузов, А О Аддитивная шестистадийная численная схема третьего порядка точности для жестких систем /АО Тузов // VII Всероссийская конференция молодых ученых по математическому моделированию и информационным технологиям Программа и тезисы докладов - Красноярск, 2006 - С 30

6 Тузов, А О Метод третьего порядка для решения жестких аддитивных систем / А О Тузов / / Конференция молодых ученых 2006 материалы конференции молодых ученых Института вычислительного моделирования СО РАН - Красноярск, ИВМ СО РАН, 2006 - С 92-98

7 Новиков, Е А Численный метод решения задачи Коши для автономных аддитивных жестких систем / Е А Новиков, А О Тузов // Ресурсосберегающие технологии механизации сельского хозяйства - Приложение к вестнику КрасГАУ - выпуск 4 -Красноярск, 2007 - С 156-160

8 Тузов, АО L- устойчивый метод третьего порядка для жестких систем /АО Тузов // Студенческая наука - взгляд в будущее материалы Всероссийской студенческой научной конференции ч 2/Краснояр гос аграр ун-т-Красноярск, 2005 - С 141143

9 Тузов, А О Метод третьего порядка решения жестких автономных систем /АО Тузов // Студенческая наука - взгляд в будущее материалы Всероссийской студенческой научной конференции ч 2 /Краснояр гос аграр ун-т- Красноярск, 2005-С 140-141

ЛП №040943 от 02 03 99 Подписано в печать 5 04 07 Формат бумаги 60 х 841/]6 Уел печ л 1 Тираж 100 экз Заказ 12

Отпечатано на ризографе ИВМ СО РАН 660036, Красноярск, Академгородок

Введение

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

1.1 Основные понятия и определения.

1.2 Численная схема второго порядка доя автономных систем.

1.2.1 Условия второго порядка точности.

1.2.2 Исследование устойчивости

1.2.3 Связь схемы второго порядка с явным методом Рунге - Кугта и ш, к) - методом.

1.2.4 Контроль точности вычислений.

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

1.3 О применимости к неавтономным задачам схем. построенных для автономных задач.

1.4 Численная схема второго порядка для неавтономных систем.

1.5 Анализ результатов расчетов.

2 Методы третьего порядка точности для автономных систем

2.1 Численная схема 1.

2.1.1 Условия третьего порядка точности.

2.1.2 Исследование устойчивости

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

2.1.4 Связь схемы 1 с явным методом Рунге - Кутта и (т, к) - методом

2.2 Численная схема 2.

2.2.1 Условия третьего порядка точности.

2.2.2 Исследование устойчивости

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

2.3 Численная схема 3.

2.3.1 Исследование устойчивости

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

2.4 Анализ результатов расчетов.

3 Методы третьего порядка точности для неавтономных систем

3.1 Численная схема 1.G

3.1.1 Условия третьего порядка точности.

3.1.2 Исследование устойчивости

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

3.2 Численная схема 2.

3.2.1 Условия третьего порядка точности.

3.2.2 Исследование устойчивости

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

3.3 Анализ результатов расчетов.

При решении ряда задач, таких как проектирование радиоэлектронных схем, моделирование кинетики химических реакций, расчет динамики механических систем и других возникает необходимость численного решения задачи Коши для жестких систем обыкновенных дифференциальных уравнений y' = f(t,y), y(to) = Vo, t0<t<tk, (0.0.1) где у и / - вещественные N - мерные вектор - функции, у0 - начальное условие, t -независимая переменная, которая меняется на заданном конечном интервале. Учет все большего числа факторов при построении математических моделей физических процессов, имеющих компоненты с сильно различающимися временными константами, приводит к жестким системам обыкновенных дифференциальных уравнений все большей размерности. Кроме того, для некоторых задач правая часть может быть разрывной.

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

Для построения эффективных алгоритмов интегрирования необходимо

1) выбрать методы, соответствующие классу решаемых задач,

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

Каждый из предложенных в данной работе методов основывается на двух численных формулах - явной и неявной. Из них строится метод, называемый аддитивным [1-4,50-55]. Для применения аддитивного метода правая часть / исходной дифференциальной задачи (0.0.1) разбивается на две части f(t,y) = ip(t,y) + g{t,y), которые будем называть "нежесткой" и "жесткой", соответственно. Аддитивный метод строится так, что явный метод используется для решения задачи с "нежесткой" частью <p(t,y), а неявный - с "жесткой" g(t,y). Поэтому, говоря о методе, "жесткую" часть g(t,y) будем называть неявный частью, а "нежесткую" часть ip(t, у) - явной. В качестве неявной составляющей аддитивного метода используются (т, к) - методы [5-20], а в качестве явной - явные методы типа Рунге - Кутта. Таким образом, в основу алгоритмов интегрирования положены одноша-говые безытерационные методы вида

Уп+i = Уп + hipf(tn, уп, h), п = 0,1,2,., где начальное условие у0 задано, п - текущая точка интегрирования, h - шаг интегрирования, <р/ - заданная вектор - функция, зависящая от правой части исходной задачи. В данной форме можно представить не только безытерационные методы, но и неявные итерационные методы типа Рунге - Кутта с фиксированным итерационным процессом, в котором число итераций не зависит от номера шага интегрирования. Одношаговые безытерационные методы имеют следующие преимущества перед многошаговыми методами:

• многошаговые методы усредняют решение ("срезают экстремумы"), что для некоторых задач неприемлемо,

• многошаговые методы малоэффективны для задач с разрывной правой частью,

• безытерационность метода позволяет оценить вычислительные затраты на шаг интегрирования до проведения расчетов,

• отсутствие итерационного процесса значительно упрощает программную реализацию алгоритма интегрирования.

Достаточно полный обзор работ по решению (0.0.1) многошаговыми методами содержится в [25,26, 29, 30,148-162]. Поэтому далее на многошаговых методах останавливаться не будем.

Многие современные способы управления величиной шага основаны на контроле точности численной схемы. Обзор способов оценки ошибки приведен в [21,25-27,32-48]. Здесь для оценки глобальной ошибки метода р - го порядка точности используется вложенный метод (р — 1) - го порядка, а неравенство для контроля точности имеет вид l?Ai,p Уп.р—il max —г1-.—;-< е.

У<i<N \Уп,р\ + Г где уп,р и уп,р~ 1 ~ соответственно, приближения к решению в точке tn методами р - го и (р — 1) - го порядков, г - положительный параметр. Если \у1п\ > г, то контролируется относительная ошибка е, в противном случае контролируется абсолютная ошибка е-r. Указанный способ хорошо зарекомендовал себя в [29,30,32-37,44-48] и будет использоваться ниже.

Поскольку в аддитивные методы входят явные численные формулы, возникают определенные проблемы с устойчивостью численной схемы. При решении жестких систем явными методами на участках установления, которые, как правило, составляют большую часть интервала интегрирования, шаг ограничен не требованием точности, а требованием устойчивости численной схемы. При выборе величины шага исходя только из требования точности на этих участках возникает неустойчивость, приводящая к раскачиванию шага, что в лучшем случае снижает эффективность алгоритма интегрирования. Этот недостаток можно устранить дополнительным контролем устойчивости численной формулы [21]. Поскольку неравенство для контроля устойчивости грубое, то оно используется как некоторый ограничитель на рост шага. В результате прогнозируемый шаг Л„+1 по точности с ограниченном но устойчивости вычисляется по формуле hn+i = max{hn, min{hacc,hst}}, где hacc - шаг, выбранный исходя из требования точности, hst - шаг, выбранный исходя из требования устойчивости, hn - предыдущий шаг.

Обычно для решения умеренно жестких задач применяются явные методы с расширенными областями устойчивости [81-135]. Однако явные методы имеют ограниченную область устойчивости, следовательно, при решении задач высокой степени жесткости требование устойчивости накладывает слишком обременительное ограничение на размер шага интегрирования [25]. Поэтому, явные методы применяются для решения задач средней жесткости [21].

Для решения жестких задач обычно используют L - устойчивые методы [5-20,56-76]. В случае большой размерности системы дифференциальных уравнений в L - устойчивых методах общие вычислительные затраты фактически полностью определяются временем вычисления и обращения матрицы Якоби. Затраты можно значительно уменьшить, если использовать одну и ту же матрицу Якоби на нескольких шагах интегрирования, то есть применять алгоритмы с замораживанием матрицы Якоби. В итерационных методах замораживание матрицы не влияет на порядок точности численной схемы, а определяет скорость сходимости итерационного процесса. Поэтому данный подход широко используется при реализации методов такого типа [26]. В безытерационных методах, к которым относятся методы типа Розенброка [77] и их различные модификации (например, [9], [32]), матрица Якоби включена непосредственно в численную формулу. Поэтому ее аппроксимация может приводить к потере точности численной схемы, что приводит к определенным проблемам с ее замораживанием. Методы типа Розенброка для автономной системы

У' = /Ы, У {to) = Уо, t0<t<tk, имеют вид т i-l

Уп+i =Уп + J^Pifci, Dni — Е - Oihf'(yn + i= 1 j=1 i-l

Dnih = hf{yn + j=l где E - единичная матрица, /' = df /ду - матрица Якоби правой части, Pi,ca, Pij^ij ~ параметры. Наиболее эффективные реализации методов такого типа возникают тогда, когда все аг равны между собой и = 0. В этом случае на каждом шаге требуется вычисление и обращение только одной матрицы Dn = Е — ahf'(yn).

Проблема замораживания матрицы Якоби значительно проще решается в рамках (т, к) - методов [5-20]. Заметим, что (т, к) - методы так же просты в реализации, как и методы типа Розенброка, а свойства точности и устойчивости у них лучше [16]. Класс т, к) - методов вводится следующим образом. Зададимся целыми числами m и к, к < т, и рассмотрим следующие множества

Мт = {1,2, .,т},

Mk = {mi е Мт|1 = т,\ < т2 < • • • < тк < т},

Ji = {trij — 1 6 Mm\j > 1, nij e Mk,mj <i} 2 < i < m.

Тогда (m, A;) - схемы имеют вид

Уп+1 = Уп +Piki+ р2к2 + . + Pmkmi Dn = E- ahf'n, t-i

Dnki = Д/(г/„ + £ + XI» e Л/t, (0.0.2) i=i jeJi

Dnki = k^i + X ctijfcj, i G il/m \ АД. jeJi

Здесь E1 - единичная матрица, f'n = Of(yn)/dy - матрица Якоби функции f,h- шаг интегрирования, 1 < i < т, - стадии метода, a,pi,/3ij,aij ~ вещественные константы, определяющие свойства точности и устойчивости (0.0.2). Для описания традиционных одношаговых методов достаточно одной константы т - числа стадий. В данных схемах для описания затрат на шаг необходимо введение двух чисел тп и к. Вычислительные затраты на шаг интегрирования в методах (0.0.2) следующие - один раз вычисляется матрица Якоби и осуществляется декомпозиция Dn, к раз вычисляется функция /, тп раз осуществляется обратный ход метода Гаусса.

Отметим, что при к = m и а^- = 0 данные методы совпадают со схемами типа Розен-брока. В этом случае (к, к) - схему (0.0.2) можно поставить в соответствие к - стадийной полуявной численной формуле типа Рунге - Кутта, при реализации которой на каждом шаге используется одна матрица размерности N. Относительно таких формул известно, что нельзя построить к - стадийную схему выше (к + 1) - го порядка точности. Отметим, что схема максимального порядка может быть только А - устойчивой, что недостаточно для построения эффективных алгоритмов интегрирования. Поэтому в настоящее время применяются методы типа Розенброка, в которых порядок точности совпадает с числом вычислений правой части дифференциальной задачи. Однако если рассматривать схемы (0.0.2) при m > к, то при к равном 1 можно построить L - устойчивый (2, 1) - метод второго порядка, а при к равном 2 и 3 можно построить L - устойчивые методы (к + 2) - го порядка точности, соответственно.

Независимость порядка точности (т, к) - метода от числа шагов с замороженной матрицей достигается за счет того, что при получении условий порядка предполагается, что при реализации метода вместо матрицы Якоби правой части системы дифференциальных уравнений применяется матрица Л„, представимая в виде Лп = f„ + hRn, где Rn -некоторая матрица, не зависящая от шага интегрирования. Это естественное требование в случае применения одной и той же матрицы Якоби на нескольких шагах интегрирования. Отметим, что при замораживании матрицы Якоби максимальный порядок точности L -устойчивой численной схемы с к вычислениями правой части дифференциальной задачи равен (к + 1).

Попытка повысить эффективность алгоритма интегрирования отделением жесткой и нежесткой компонент решения предпринята в проекционных методах. Идея проекционных методов состоит в разделении жесткой системы у' — f(y) на жесткую систему небольшой размерности у'а = /а{Уа,Уь) и нежесткую систему большой размерности у'ь = /ь(уа,уь) с целью применения к жесткой системе неявного метода, а к нежесткой - явного. Трудность этого подхода заключается в проблеме разделения на жесткую и нежесткую системы.

Построенные в данной работе аддитивные методы предусматривают разбиение правой части f{t,y) па жесткую g(t,y) и нежесткую <p(t,y) части

Предполагается, что большая часть жесткости сосредоточена в функции g{t,y). Методы конструируются так, чтобы неявный метод использовался для решения задачи с жесткой частью, а явный - с нежесткой. В качестве неявной составляющей аддитивного метода используется (ш, к) - методы [5,13-16,18], а в качестве явной - явные методы типа Рунге -Кутта. При разбиении правой части задачи (0.0.1) на жесткую и нежесткую части можно выделить два случая

1) Естественное разбиение. При численном решении ряда задач, описываемых дифференциальными уравнениями в частных производных, после дискретизации по пространственным переменным (методом конечных разностей или конечных элементов) возникает необходимость решения задачи Коши для системы обыкновенных дифференциальных уравнений вида (0.0.3), в которой правая часть естественным образом разбивается на жесткую и нежесткую составляющие. При этом разбиении выполняется предположение о том, что большая часть жесткости сосредоточена в функции g(t, у). Например, в некоторых задачах механики сплошной среды жесткая часть появляется после дискретизации дифференциального оператора второго порядка, а нежесткая - после дискретизации дифференциального оператора первого порядка. Кроме того, в некоторых задачах вида (0.0.3) матрица Якоби функции g(t, у) имеет специальный вид (например, в задачах механики сплошной среды она симметрична).

2) Искусственное разбиение. Полагая g(t,у) = By, <p(t,y) = f(t,y) — By, где В = df(t,y)/dy, можно считать, что вся жесткость сосредоточена в функции g(t,y) [1-4, 50-53]. Тогда система (0.0.1) примет вид

В данной работе и в [4,50-53] в качестве В может быть использована произвольная аппроксимация матрицы Якоби df(t,y)/dy. При "разумной" аппроксимации этой матрицы можно считать, что основная жесткость задачи сосредоточена в функции g(t,y). Однако считать, что вся жесткость сосредоточена в функции g(t,y), а ip(t,y) есть нежесткая часть, вообще говоря, нельзя. Поэтому все построенные в данной работе методы оснащены неравенством для контроля устойчивости явной части численной схемы, которое в ряде случаев повышает эффективность алгоритма интегрирования. Контроль устойчивости явной части может оказаться полезным и для случая естественного разбиения правой части. Если есть уверенность, что вся жесткость сосредоточена в функции g(t, у), то контроль устойчивости следует отключить, у' = <p{t,y) + g{t,y), y{to) = yo, t0<t<tk.

0.0.3) y' = [f(t,y)-By} + By,

0.0.4) потому что в некоторых случаях он может приводить к понижению эффективности алгоритма интегрирования.

Условия порядка для всех предложенных здесь схем выводились без каких-либо предположений относительно вида функций <p(t, у) и g(t, у). В результате, для построенных аддитивных методов несущественно, естественным образом или искусственным была разбита правая часть, а имеет значение лишь то, как распределена жесткость задачи относительно функций <p(t,y) и g(t,y). В идеале вся жесткость должна быть сосредоточена в функции g(t,y). Возможность произвольной аппроксимации матрицы Якоби исходной задачи, без снижения порядка точности, в аддитивных методах достигается за счет произвольности представления задачи (0.0.1) в виде (0.0.4). В частности, метод допускает замораживание матрицы Якоби.

Хотя совершенно произвольная аппроксимация матрицы Якоби не снижает порядка точности метода, она может понизить его эффективность.

1. Например, для задачи (30) из [26, с. 167] матрица В = diag(df(y)dy) слишком грубо показывает направления изменения решения из-за особенности этой задачи - когда компонента решения у2 случайно становится отрицательной, она стремится к -ос. В результате, для данной задачи алгоритм интегрирования с диагональной аппроксимацией матрицы Якоби менее эффективен, чем с полной матрицей Якоби, несмотря на экономию вычислительных затрат при ее обращении.

2. Несмотря на контроль устойчивости явной части численной схемы, алгоритм интегрирования будет достаточно эффективен, когда основная часть жесткости сосредоточена в функции g(t,y), соответствующей неявной части схемы. Например, теоретически возможным, но практически бесполезным, было бы положить 5 = 0. Тогда g{t,y) = 0 и вся нагрузка легла бы на явный метод Рунге - Кутта, который, по указанным выше причинам, малоэффективен для решения жестких задач.

3. Диагональную аппроксимацию целесообразно использовать для матрицы Якоби с диагональным преобладанием.

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

Цель настоящей работы состоит в разработке эффективных алгоритмов численного решения задачи Коши для жестких систем обыкновенных дифференциальных уравнений большой размерности. Построенные алгоритмы призваны расширить область применения, занимаемую явными методам (возможно, обладающими расширенной областью устойчивости и оснащенных неравенством для контроля устойчивости [21]), в сторону задач большей степени жесткости, а не конкурировать с методами, использующими полную аппроксимацию матрицы Якоби. При этом их вычислительные затраты на шаг интегрирования должны быть сопоставимыми с явными методами.

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

В исходной задаче Коши правая часть может быть представлена как в виде (0.0.1), •гак и (0.0.3). В первом случае задача приводится к виду (0.0.4), то есть производится искусственное разбиение правой части. Во втором случае правая часть уже разбита исходя из некоторых соображений предметной области, согласно которым большая часть жесткости сосредоточена в функции g(t,y), то есть разбиение было произведено естественным образом.

В случае естественного разбиения правой части эффективность алгоритма интегрирования достигается за счет того, что отделение более жесткой части g(t,y) от менее жесткой (f(t,y) позволяет использовать неявную схему, входящую в аддитивный метод, преимущественно для решения задачи у' = g(t,y). Следовательно, неявная схема требует вычисление и обращение матрицы Якоби только функции g(t,y), а не всей правой части исходной задачи. Во многих прикладных задачах матрица dg(t,y)/dy, в отличие от d((p(t,y) + g(t,y))/dy, имеет специальный вид (например, в некоторых задачах механики сплошной среды она симметрична). В этом случае эффективность можно повысить за счет выбора специального метода решения линейных систем алгебраических уравнений, учитывающего информацию о виде матрицы Якоби (например, метода квадратного корня для симметричных матриц [31]), и, возможно, использующего параллельный алгоритм вычислений [136-147]).

В случае искусственного разбиения правой части затраты на обращение матрицы В можно значительно уменьшить за счет использования в качестве В не полной матрицы df(t,y)/dy, а некоторой ее аппроксимации специального вида и применения соответствующего этой аппроксимации метода решения линейных систем алгебраических уравнений. Применение аппроксимации матрицы Якоби специального вида не только уменьшит затраты на ее вычисление, но и позволит применить соответствующий метод решения линейных систем алгебраических уравнений, что может весьма значительно уменьшить затраты на ее обращение. В данной работе использовалась диагональная аппроксимация матрицы Якоби, поскольку именно в этом случае затраты на ее обращение наиболее низки и аддитивный метод по затратам на шаг интегрирования становится сопоставимым с явными методами, а по свойствам устойчивости, относительно неявной части, значительно превосходит их.

Возможно сочетание описанных способов, например, использование трехдиагональ-ной аппроксимации матрицы Якоби, применение метода прогонки для решения систем с этой матрицей и ее замораживание. При решении жестких задач на эффективность алгоритма интегрирования существенно влияют свойства устойчивости не только основной, но и промежуточных численных формул [21,26,77]. В данной работе этому также уделено немало внимания. Для многих задач эффективность алгоритма интегрирования можно повысить за счет контроля устойчивости явной части численной схемы.

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

Заключение

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

• Исследованы три семейства методов решения жестких автономных и два семейства методов решения жестких неавтономных задач.

• Получены коэффициенты L - устойчивых шестистадийных численных схем третьего порядка точности.

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

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

1. Cooper, G.J. Additive methods for the numerical solution of ordinary differential equations / G.J. Cooper, A. Sayfy // Mathematics of Computation.- vol. 35.- №152.1980.- P. 1159-1172.

2. Cooper, G.J. Additive Runge Kutta methods for Stiff Ordinary Differential Equations / G.J. Cooper, A. Sayfy // Mathematics of Computation.- vol. 40.- №161.- 1983.- P. 207-218.

3. Cooper, G.J. Semiexplicit A stable Runge - Kutta methods / G.J. Cooper, A. Sayfy // Mathematics of Computation.- vol. 33.- №146.- 1979.- P. 541-556.

4. Новиков, Е.А. Некоторые методы решения жестких систем, индуцированные одним и двумя вычислениями правой части / Е. А. Новиков, Ю.А. Шитов // Математические модели и методы решения задач механики сплошной среды.- Красноярск, 1986.- С. 11-19.

5. Новиков, Е.А. Исследование (т, к) методов решения жестких систем с одним и двумя вычислениями правой части / Е. А. Новиков, Ю.А. Шитов // Препринт №15: Красноярск, ВЦ СО АН СССР,- 1987,- 41 с.

6. Новиков, Е.А. Об одном классе одношаговых безытерационных методов решения жестких систем / Е. А. Новиков // Актуальные проблемы вычислительной и прикладной математики.- Новосибирск, 1987.- С. 138-139.

7. Новиков, Е.А. О классе (т, к) методов решения жестких систем / Е. А. Новиков, Ю.А. Шитов, Ю.И. Шокин // ЖВМ и МФ,- 1989.- т.29.- №2,- С. 194-201.

8. Новиков, Е.А. Численные методы решения дифференциальных уравнений химической кинетики / Е. А. Новиков // Математические методы в химической кинетике.-Новосибирск: Наука.- 1990.- С. 53-68.

9. Новиков, Е.А. Алгоритм интегрирования жестких систем на основе (тп, к) метода второго порядка точности с численным вычислением матрицы Якоби / Е. А. Новиков, Ю.А. Шитов // Препринт №20: Красноярск, ВЦ СО АН СССР.- 1988.- 27с.

10. Двинский, АЛ. Замораживание матрицы Якоби в (3, 2) методе решения жестких систем / A. J1. Двинский, Е. А. Новиков // Вычислительные технологии. 2003. т.8. Региональный вестник Востока.- 2003.- №3.- (Совм. выпуск.- ч. 2.)- С. 272-278.

11. Новиков, Е.А. L устойчивая (5, 2) - схема четвертого порядка / Е. А. Новиков, A. JI. Двинский // Вопросы математического анализа.- выпуск 8, Красноярск, ИПЦ КГТУ, 2004.- С. 134-142.

12. Новиков, Е.А. Максимальный порядок точности (т, 3) методов с замораживанием матрицы Якоби / Е. А. Новиков, A. JI. Двинский // Вопросы математического анализа.- выпуск 9, Красноярск, 2006.- С. 54-64.

13. Новиков, Е.А. (2, 2) метод с замораживанием матрицы Якоби для жестких систем / Е. А. Новиков, А. Л. Двинский, Ю.А. Шитов // Вестник КрасГАУ: Ресурсосберегающие технологии,- выпуск 3.- Красноярск, 2005.- С. 100-104.

14. Новиков, Е.А. Замораживание матрицы Якоби в методах типа Розенброка второго порядка точности / Е. А. Новиков, В. А. Новиков, Л. А. Юматова // ЖВМ и МФ,-1987.- т.27,- JV3.- С. 385-390.

15. Хайрер, Э. Решение обыкновенных дифференциальных уравнений. Нежесткие задами. / Э. Хайрер, С. Нерсетт, Г. Ваннер // М.: Мир,- 1990.- с. 512.

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

17. Shampine, L. Numerical Solution of Ordinary Differential Equations / L. Shampine // New York: Chapman к Hall.- 1994.

18. Деккер, К. Устойчивость методов Рунге Кутты для жестких нелинейных дифференциальных уравнений / К. Деккер, Я. Вервер // М.: Мир,- 1988.

19. Штеттер, X. Анализ методов дискретизации для обыкновенных дифференциальных уравнений , X. Штеттер // М.: Мир.- 1978.

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

21. Бахвалов, Н.С. Численные методы / Н.С. Бахвалов, II.П. Жидков, Г.М. Кобельков // М.: Наука,- 1987.- с. 598.

22. Kaps, P. Generalized Runge Kutta methods of order four with stepsize control for stiff ordinary differential equations / P. Kaps, P. Rentrop // Numer. Math.- 1979.- JV«33.- P. 55-68.

23. Gear, C.W. The automatic integration of stiff ordinary differential equations / C.W. Gear // Proc. IFIP Congress.- 1968.- P. 81-85.

24. Merson, R.H. An operational methods for integration processes /' R.H. Merson // Proc. Svinp. on Data Proc.- Weapons Research Establishment.- Salisbury.- Australia.- 1957.

25. Fehlberg, E. Classical fifth-, sixth-, seventh- and eighth order Runge Kutta formulas with step size control / E. Fehlberg // Computing.- 1969.- vol.4.- P. 93-106.

26. Enright, W.H. Comparing numerical methods for the solutions of stiff systems of ODE's / W.H. Enright, Т.Е. Hull // BIT.- 1975.- №15.- P. 10-48.

27. Estep, D. Global error control for the continuous Galerkin finite element method for ordinary differential equations / D. Estep and D. French // Model. Math. Anal. Numer.-1994,- j\*e28.- P. 815-852.

28. Estep, D. Estimating the error of numerical solutions of systems of reaction diffusion equations / D. Estep, M.G. Larson, R. Williams // Mem. Amer. Math. Soc.- 2000.- №146.-P. 1-109.

29. Gustafsson, K. A PI stepsize control for the numerical solution of ordinary differential equations / K. Gustafsson, M. Lundh and G. Soderlind // BIT.- 1988.- №28.- P. 270-287.

30. Shintani, H. On a one-step method of order 4 / H. Shintani / J. Sei. Hiroshima Univ.-1966.- №30,- P. 91-107.

31. Shintani, H. Two step processes by one - step methods of order 3 and of order 4 / H. Shintani // J. Sei. Hiroshima Univ.- 1966.- №30,- P. 183-195.

32. Rosenbrock H.H. Some general implicit processes for the numerical solution of differential equations // Computer J., 1963, No 5.- P. 329-330.

33. Тузов, А.О. Численный метод третьего порядка точности для решения автономных аддитивных жестких систем / А.О. Тузов // Вестник КрасГАУ.-выпуск 14,- Красноярск, 2006.- С. 467-473.

34. Новиков,Е.А. Метод третьего порядка для автономных аддитивных жестких систем / Е.А. Новиков, А.О. Тузов // Вычислительные и информационные технологии в науке, технике и.образовании: Труды международной конференции,- Павлодар, 2006.-С. 77-84.

35. GO. Lambert, J.D. Computational Methods in Ordinary Differential Equations / J.D. Lambert ,// Wiley.- New York.- 1972.

36. Liniger, W. Efficient, numerical integration of stiff systems of ordinary differential equations / W. Liniger, R.A. Willoughby // Technical Report RC-1970.- Thomas J. Watson Research Center.- New York.- 1976.

37. Butcher, J.C. Diagonally implicit multi - stage integration methods / J.C. Butcher // Appl. Numer. Math.- 1993.- №11.- P. 347-363.

38. Butcher, J.C. Experiments with a variable order type 1 DIMSIM code / J.C. Butcher, P. Chartier, Z. Jackiewicz // Numer. Algorithms.- 1999.- №22.- P. 237-261.

39. Butcher, J.С. Automatic selection of the initial step size for an ODE solver / I. Glachvell, L.F. Shampine, RAV. Brankin // J. Comput. Appl. Math.- 1987,- №18.- P. 175-192.

40. G8. Shampine, L.F. Numerical Solution of Ordinary Differential Equations / L.F. Shampine // Chapman к Hall.- London.- 1994.

41. Shampine, L.F. The Matlab ODE suite / L.F. Shampine, MAY. Rcichelt // SIAM J. Sci. Comput.- 1997.- №18.- P. 1-22.

42. Артемьев, С.С. Минимизация овражных функций численным методом для решения жестких систем уравнений / С.С. Артемьев, Г.В. Демидов, Е.А. Новиков // Препринт №74: Новосибирск, ВЦ СО АН СССР.- 1980,- 13 с.

43. Альшин, А.Б. Числнное решение сверхжестких дифференциально алгебраических систем / А.Б. Альшин, Е.А. Алыпина, Н.Н. Калиткин, А.Б. Корягипа // ДАН.-2006,- т.408.- .№4,- С. 1-5.

44. Альшин, А.Б. Схемы Розенброка с комплексными коэффициентами для жестких дифференциально алгебраических систем / А.Б. Альшин, Е.А. Алыпина, Н.Н. Калиткин, А.Б. Корягина // ЖВМ и МФ,- 2006,- т.46.- №8.- С. 1410-1431.

45. Alshina, Е.А. Integration of differential algebraic stiff system / E.A. Alshina, N.N. Kalitkin, А.В Koryagina // Methematical Modelling and Analisis.- Inter. Conf.-MMA2005&CMAM2.- 2005.- P. 301-307.

46. Alshina, E.A. New schemes for differential algebraic stiff system / E.A. Alshina, N.N. Kalitkin, А.В Koryagina // Springer Verlag.- Progress in Industrial Mathematics at ECMI.- 2005.

47. Dormand, J.R. High-Order Embedded Runge Kutta - Nystrom Formulae / J.R. Dormaml, E.A. El - Mikkawy, P.J. Prince // IMA J. Nuiner Anal.- 1987.- №7.- P. 423-430.

48. S31 Новиков, В.А. О построении явных методов тина Рунге Кутта с расширенными областями устойчивости / В.А. Новиков, Е.А. Новиков // Препринт №9: Красноярск, ВЦ СО АН СССР,- 1988.

49. Новиков, В.А. Численное конструирование областей усустойчивости явных методов / В.А. Новиков, Е.А. Новиков // Препринт №15: Красноярск, ВЦ СО АН СССР,-1988.

50. Лебедев, В.И. Явные разностные схемы с переменными шагами по времени для решения жестких систем уравнений / В.И. Лебедев // Препринт №177, М.: ОВМ АН СССР,- 1987.

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

52. Новиков, В.А. Два эффективных алгоритма численного решения задачи Коши для систем обыкновенных дифференциальных уравнений / В.А. Новиков, Е.А. Новиков // Препринт №5: Новосибирск, ИТПМ СО АН СССР,- 1984.

53. Новиков, В.А. Об алгоритме переменной структуры на основе явных формул тина Рунге Кутта первого и второго порядков точности / В.А. Новиков, Е.А. Новиков // Препринт №112: Новосибирск, ВЦ СО АН СССР,- 1985.

54. Новиков, В.А. Явные методы для решения жестких систем обыкновенных дифференциальных уравнений / В.А. Новиков, Е.А. Новиков // Препринт №629: Новосибирск, ВЦ СО АН СССР,- 1985.

55. Новиков, Е.А. Численные методы решения дифференциальных уравнений химической кинетики / Е.А. Новиков // Математические методы в химической кинетике, Новосибирск: Наука.- 1990, С. 53-68.

56. Novikov, E.A. The Explicit Methods: Algorithms with stability Control / E.A. Novikov, Kontareva L.N. // M.: Dynamics of non-homogeneous systems.- 2001.- vol.5.- P. 107-118.

57. Novikov, E.A. The Explicit Methods: Integration Algorithms with Accuracy Control / E.A. Novikov, Kontareva L.N. // M.: Dynamics of non-homogeneous systems.- 2001.-vol.5.- P. 119-138.

58. Novikov, E.A. The program NODE for solution of ODE stiff systems / E.A. Novikov // Numerical Analysis.- Novosibirsk: NCC Publisher, Bulletin of the Novosibirsk computing center.- 2002,- №11,- P. 95-101.

59. Новиков, E.A. Численные методы для жестких систем со специальными свойствами устойчивости / E.A. Novikov // Вопросы математического анализа.- выпуск 6, Красноярск, КГТУ, 2003,- С. 187-198.

60. Новиков, Е.А. Согласование областей устойчивости явных методов / E.A. Novikov // Информационные системы и технологии.- Новосибирск, НГТУ, 2003.- С. 137-143.

61. Новиков, Е.А. Внутренняя устойчивость явных методов решения жестких систем / Е.А. Novikov // Международная конференция но вычислительной математике.-Новосибирск, 2004,- ч.2.- 567-573.

62. Новиков, Е.А. Алгоритм переменного порядка и шага на основе стадий метода Фель-берга седьмого порядка точности / Е.А. Новиков, Ю.В. Шорников, О.В. Никонова // Научный вестник НГТУ,- Новосибирск,- 2006.- vol.25.- №2,- С. 105-119.

63. Verwer, J.G. Explicit Runge Kutta methods for parabolic partial differential equations / J.G. Verwer // Appl. Numer. Math.- 1996,- №22.- P. 359-379.

64. Ansorge, R. Zur stabilitat des Nystromschen verfahren / R. Ansorge, W. Tornig // Z. Angew. Math. Mech.- I960,- №40.- P. 568-570.

65. Chawla, M.M. Absolute stability of explicit Runge Kutta - Nystrom methods for y" = f{x,y,y') / M.M. Chawla, S.R. Sharma // J. Сотр. Appl. Math.- 1984,- №10.-P. 163-168.

66. Chawla, M.M. Nurnerov made explicit has better stability / M.M. Chawla // BIT.- 1984.-№24,- P. 117-118.

67. Chawla, M.M. An explicit sixth order method with phase - lag of order eight for y" = f(t,y) / M.M. Chawla, P.S. Rao // J. Сотр. Appl. Math.- 1987.-№17.- P. 365-368.

68. Jeltsch, R. Largest disk of stability of explicit Runge Kutta methods / R. Jeltsch, O. Nevanlinna // BIT.- 1978.- №18,- P. 500-502.

69. Jeltsch, R. Stability and accuracy of time discretizations for initial value problems / R. Jeltsch, 0. Nevanlinna // Report HTKK MAT - А187,- Helsinki, Univ. of Technology.-1981.

70. Jeltsch, R. Stability of explicit time discretizations for solving initial value problems / R. .Jeltsch, O. Nevanlinna // Num. Math.- 1981.- №37,- P. 61-91.

71. Riha, W. Optimal stability polynomials / Riha W. // Computing.- 1972.- №9.- P. 37-43.

72. Sommeijer, B.P. On the economization of stabilized Runge Kutta methods with applications to parabolic initial value problems / B.P. Sommeijer, P.J. Van der Houwen // ZAMM.- 1981,- №61.- P. 105-114.

73. Van der Houwen, P.J. Explicit Runge Kutta methods with increased stability boundaries / P.J. Van der Houwen // Numer. Math.- 1972,- №20,- P. 149-164.

74. Van der Houwen, P.J. Stabilized Runge Kutta methods for second - order differential equations without first derivatives / P.J. Van der Houwen // SIAM J. Numer. Anal.-1979,- №16.- P. 523-537.

75. Van der Houwen, P.J. A special class of multistep Runge Kutta methods with extended real stability interval / P.J. Van der Houwen, B.P. Sommeijer // IMA J. Numer. Anal.-1982.- №2.- P. 183-209.

76. Van der Houwen, P.J. Predictor corrector methods with improved absolute stability regions / P.J. Van der Houwen, B.P. Sommeijer // IMA J. Numer. Anal.- 1983.- .№3.- P. 417-437.

77. Verwer, J.D. A class of stabilized three step Runge - Kutta methods for the numerical integration of parabolic equations / J.D. Verwer // J. Comput. Appl. Math.-1977.-№3.-P. 155-166.

78. Eriksson, K. Explicit time stepping for stiff ODEs / K. Eriksson, C. Johnson, A. Logg // SIAM J. Sci. Comput.- 2003.

79. Verwer, J.D. An implementation of a class of stabilized explicit methods for the integration of parabolic equations / J.D. Verwer // ACM Trans. Math. Software.- 1980.- №3.- P. 188205.

80. Verwer, J.D. A note on a Runge Kutta - Chebyshev method /' J.D. Verwer // ZAiMM.-1982,- №62.- P. 561-563.

81. Hoffmann, W. Approximating Runge Kutta matrices by triangular matrices / W. Hoffmann, J.J.B. Swart // Preprint NM-R9517, CWI.- Amsterdam.- 1995.

82. Houwen, P.I. Embedded diagonally implicit Runge Kutta algorithms on parallel computers / P.J. Houwen, B.P. Somineijer, W. Couzy // Math. Сотр.- 1992.- №58.-P. 135-159.

83. Houwen, P.J. Triangularly implicit iteration methods for ODE IVP solvers / P.J. Houwen, J.J.B. Swart // Preprint NM-R9510, CWI.- Amsterdam.- 1995.

84. Jackson, K.R. The potential for parallelism in Runge Kutta methods. Part 1: RK formulas in standard form / K.R. Jackson, S.P. Norsett //' SIAM J. Numer. Anal.- 1995.-№32,- P. 49-82.

85. Estep, D. Accurate parallel integration of large sparse systems of differential equations / D. Estep, R. Williams // Math. Models Methods Appl. Sci.- 1996.- №6,- P. 535-568.

86. Butcher, J.C. Order and stability of parallel methods for stiff problems / J.C. Butcher // Adv. Comput. Math.- 1997,- №7.- P. 79-96.

87. Burrage, K. Parallel methods for initial value problems / K. Burrage // Appl. Numer. Math.- 1993,- №11.- P. 5-25.

88. Burrage, K. Parallel and sequential methods for ordinary differential equations / K. Burrage // Clarendon Press, Oxford.- 1995.

89. Gear, C.W. Parallelism across time in ODEs / C.W. Gear, Xu Xuhai // Appl. Numer. Math.- 1993.- №11,- P. 45-68.

90. Serban, R. CVODES, the sensitivity enabled ode solver in SUNDIALS / R. Serban, A. C. Hindmarsh // Nonlinear Dynamics and Control.- In Proceedings of the 5th International Conference on Multibody Systems.- Long Beach.- CA.- 2005.

91. Cohen, S.D. Stiff/Nonstiff ODE Solver in С / S.D. Cohen, A.C. Hindmarsh // Computers in Physics.- 1996.- vol.10.- №2.- P. 138-143.

92. Brown, P.N. Consistent Initial Condition Calculation for Differential Algebraic Systems / P.N. Brown, A.C. Hindmarsh, L.R. Petzold // SIAM J. Sci. Сотр.- 1998.- vol.19.- P. 1495-1512.

93. Hindmarsh, A.C. Simulation of Chemical Kinetics Transport in the Stratosphere in Stiff Differential Systems / A.C. Hindmarsh, J.S. Chang, N.K. Madscn // R. A. Willoughby.-New York.- 1973,- P. 51-65.

94. Hindmarsh, A.C. Numerical Solution of Stiff Ordinary Differential Equations / A.C. Hindmarsh // AIChE Today Series.- American Institute of Chemical Engineers.- New York.- 1977.

95. Hindmarsh, A.C. GEARS: A Package for the Solution of Sparse Stiff Ordinary Differential Equations / A.C. Hindmarsh, A.H. Sherman //in Electrical Power Problems: The Mathematical Challenge.- SIAM.- Philadelphia.- 1980.- P. 190-200.

96. Hindmarsh, A.C. The ODEPACK Solvers, in Stiff Computation / A.C. Hindmarsh // Oxford U. Press.- 1985,- P. 167-174.

97. Hindmarsh, A.C. ODE Solvers for Use with the Method of Lines / A.C. Hindmarsh // in Advances in Computer Methods for Partial Differential Equations IV, R. Vichnevetsky and R. S. Stepleman, eds.- IMACS.- New Brunswick.- NJ.- 1981.- P. 312-316.

98. Hindmarsh, A.C. ODE Solvers for Time Dependent PDE Software / A.C. Hindmarsh // in PDE Software: Modules, Interfaces, and Systems, B. Engquist and T. Smedsaas, eds.- North-Holland.- Amsterdam.- 1984,- P. 325-337.

99. Hindmarsh, A.C. Stiff System Problems and Solutions at LLNL / A.C. Hindmarsh // in Stiff Computation, ed.- Oxford U. Press.- 1985.- P. 24-29.

100. Hindmarsh, A.C. Matrix Free Methods for Stiff Systems of ODE's / A.C. Hindmarsh, P. N. Brown // SIAM J. Num. Anal.- 1986,- №23,- P. 610-638.