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

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

Двинский Антон Леонидович

Исследование (т, к) -методов с L-устойчивыми промежуточными схемами для решения жестких систем

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

АВТОРЕФЕРАТ

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

Красноярск 2004

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

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

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

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

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

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

Половинкин Владимир Ильич

Защита состоится «28» декабря 2004 года в 15 часов на заседании диссертационного совета К 212.098.03 в Красноярском государственном техническом университете по адресу: 660074, г. Красноярск, ул. Киренского, 26.

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

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

кандидат физико-математических наук, доцент

Левыкин Александр Иванович

Ведущая организация Институт вычислительных

технологий СО РАН. (г. Новосибирск)

Ученый секретарь диссертационного

Сафонов К. В.

ОБЩАЯ ХАРАКТЕРИСТИКА РАБОТЫ

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

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

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

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

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

Известно, что L-устойчивая схема Розенброка на основе к вычислений правой части может иметь порядок точности не выше В то же время, на

основе (т, &)-методов с к вычислениями правой части можно построить Ь-устойчивую схему (/; + 2)-го порядка точности. Исследование (т, к) -методов с L-устойчивыми промежуточными схемами до сих пор не проводилось.

Цель работы - исследовать (ш, &)-методы с Ь-устойчивыми промежуточными схемами.

Цель достигается выполнением следующих задач: доказательство теорем о максимальном порядке (ш, &)-методов, построение методов максимального порядка с Ь-устойчивыми промежуточными схемами, исследование возможности замораживания матрицы Якоби в схемах максимального порядка.

Научная новизна

• Доказаны теоремы о максимальном порядке точности (т, к) -методов с замораживанием и без замораживания матрицы Якоби при к = 1,2 и 3. Построены методы максимального порядка с Ь-устойчивыми промежуточными схемами.

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

• Получена пара (т, к) -методов, в которых для нахождения стадий применяется одна и та же матрица.

Теоретическая значимость

Доказаны теоремы о максимальном порядке точности для (т, к) - методов с аналитическим вычислением матрицы Якоби при к = 1,2,3 и для (т,к)-методов с замораживанием при = 2 и 3. На основе -методов построены алгоритмы интегрирования со второго по пятый порядок точности включительно, имеющие Ь-устойчивые промежуточные схемы. Исследована возможность построения алгоритмов переменного порядка и шага на основе (т,к)-методов. Показано существование (т, &)-методов, в которых для нахождения стадий применяется одна и та же матрица, что делает возможным построение методов интегрирования переменного порядка и шага.

Практическая ценность

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

Методы исследования

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

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

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

• Межвузовская конференция студентов, аспирантов и молодых ученых "Информатика и информационные технологии"(Красноярск. 2000,2001,2002,2003).

• XXXIV краевая научно-техническая студенческая конференция (Красноярск. 2001).

• Региональная научная конференция студентов, аспирантов и молодых ученых "Наука. Техника. Инновации" (Новосибирск, 2001).

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

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

• V Международная конференция молодых ученых по математическому моделированию и информатике (Новосибирск, 2004).

Структура и объем работы. Диссертация состоит из введения, 4 глав основного текста, заключения и одного приложения, содержит 38 таблиц. Библиографический список включает 62 наименования. Объем работы составляет 150 листов.

СОДЕРЖАНИЕДИССЕРТАЦИИ

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

на основе L-устойчивых (то, А;)-методов вида

^ = "<^Мт\] > 1 ,т} е Мьгп> < ¿},~ 2< г< т,

где у и / - вещественные гладкие 1У-мерные вектор-функции, t - независимая переменная, Л - шаг интегрирования. к„ 1 < г < ТО, - стадии метода, - вещественные константы, определяющие свойства точности и устойчивости (2), = д/(уп)/ду- матрица Якоби системы (1). Приведен обзор всех разделов диссертации.

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

Теорема 1. Максимальный порядок точности (то, 1)-методов с аналитическим вычислением матрицы Якоби равен двум.

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

Теорема 3. Максимальный порядок точности (т, 3)-методов с аналитическим вычислением матрицы Якоби равен пяти.

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

Теорема 5. Максимальный порядок точности (т, 3)-методов с замораживанием матрицы Якоби и возможностью ее численной аппроксимации равен четырем.

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

Исследована (2,2)-схема вида

Оп = Е- аН?п, Опку = н/(уп), Ппк2 = Л/(уп + а^). (3)

Показано, что схема (3) является Ь-устойчивой, имеет второй порядок точности и обладает Ь-устойчивой внутренней схемой. Для контроля точности (3) использовалось неравенство

где - требуемая точность расчетов, -некоторая норма в ЯР.

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

Исследован метод третьего порядка на основе (3, 2)-метода

Уп+1 = Уп + рЛ + Р2^2 + Рзк3, Д, = Е - аЬ,}'п &пкI = Л/Ы, Опк2 ~ ки АЛ = кЦуп + Рцкх + РыЬ) + а32к2.

(4)

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

при которых она имеет третий порядок точности и обладает Ь-устойчивой основной и промежуточной схемами.

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

о = +0.435866521508459, рх = +0.4358665215084, щ = +0.8156217517823, рз = +0.5925925925925, = +0.4358665215084, Д32 = +0.3141334784915,

Неравенство для контроля точности построено на основе ранее вычисленных стадий к\ и с помощью вспомогательной численной схемы второго порядка точности

Для контроля точности и выбора величины шага интегрирования (4) используется неравенство

Применение (5) для контроля точности вычислений не приводит к дополнительным вычислениям правой части и матрицы Якоби, поскольку в данном неравенстве используются ранее вычисленные стадии

Данный метод эффективен при расчетах с требуемой точностью порядка Ю-4, однако при более высокой точности расчетов требуются методы более высокого порядка.

Построен (т, 2)- метод максимального (четвертого) порядка точности. Соответствующая численная схема имеет вид

(6)

где промежуточная численная формула записывается следующим образом

Уп+1 = Уп+031*1+

Найдены аналитические выражения для коэффициентов (б), при которых схема является ¿-устойчивой, имеет четвертый порядок точности и обладает ¿-устойчивой внутренней схемой. Для практических вычислений рекомендуется набор коэффициентов при а = +0.219669914110089, а именно

а = +0.219669914110089, = +0.219669914110089,

Р2 = +0.266835225483347, р3 = +0.401841276140395,

р4 = +0.299682669966500, рь = - 0.108931353514302, (7)

031 = +0.219669914110089, #¡2 = +0.0530330085889911,

а32 = - 2.33854786494379, а42 = +6.85032446594068.

Для контроля точности на основе первых четырех стадий к\, к2, кз, к^ ш-строена схема третьего порядка

г/п+1,3 -Уп + Г А + Г2к2 + Г3к3 + пк4.

Точность (6) контролируется посредством неравенства

На основе (т, к) -методов с тремя вычислениями правой части построен (6,3)-метод максимального порядка точности вида

(8)

Найдены аналитические выражения для коэффициентов (8), при которых схема является ¿-устойчивой и имеет пятый порядок точности. Для практических вычислений рекомендуется набор коэффициентов при а = +0.76677148671184, а именно

а = +0.76677148671184, Р2 = - 2.94136515866556, щ = +1.31497204379270, рй = +0.75454330862107,

Р1 = +1.35502017783518, рз = +1.38436882308358, р5 = - 0.88298387691592, Ра = +1.56879960316762,

(9)

/?42 = - 1.11111156897868, /?43 = +0.36487788897879, Аг = +0.32574077312011, = - 0.26189456959060,

Рез = - 0.41052356997426, ¡Зм = +0.34313548548096,

= +0.21675295014995, а43 = +0.03575255859581.

Неравенство для контроля точности вычислений получено с помощью дополнительной стадии

АА = А'б + а7з^з +

с применением которой построена ¿-устойчивая схема четвертого порядка

Уп+1,4 = Уп + пЬ + Г2к2 + Гзкз + Г4&4 + г5к5 + Гбкб + Г7^7-

Здесь ki: 1 < i< 6 - стадии схемы (8).

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

112/n+i - 2/n+i,4ll < е.

Для контроля точности (8) с коэффициентами (9) использовались следующие коэффициенты вспомогательной схемы

П - +1.6398200946875, г2 = - 2.0099617353737, г3 = - 0.6411237402454,

П = +0.8346699640250, г3 = +0.13828078911602, гб = 0.0,

г7 = 0.7545433086210, а73 = +1.0725550328766, а75 = - 0.7169404063139.

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

Построен (6, 3)-метод максимального порядка точности, обладающий L-устойчивыми внутренними схемами. Соответствующая численная формула имеет вид

Уп+i = Уп + Pih + P2h + Pzh + Pih + p5fc5 + Рб^б,

D„h = Ц, Dnh = hf(fn+1) + «63^3 +

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

fn+1 = Уп+ /ЗвА + peih + fah + + Ptth.

Найдены аналитические выражения для коэффициентов метода (10), при которых он является L-устойчивым, имеет пятый порядок точности и обладает L-устойчивыми внутренними схемами. Внутренней L-устойчивости удалось добиться за счет изменения вида шестой стадии.

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

а = +0.274551306835187, а именно

а = +0.274551306835187, р2 = - 0.3197890642437, р4 = +0.05866989302597, Рб = +0.6233616729200, 0а = - 3.162478318369, Да = +0.2745513068351, /?бз = +0.8036587906267, /Sss = +0.001094540820854, а« = - 4.421907933163,

pi = +0.2745513068351, Рз = +3.198890217616, К = -0.04502443732750, /?41 = +0.2745513068351, /?43 = +5.165340613251, 062 = +0.2222079136289, /?64 = +0.02732899627387, а« = - 21.22863361480, «65 = - 0.02025913426759.

(П)

Неравенство для контроля точности вычислений получено с помощью дополнительной стадии

Dnkj = к6 + а7зк$ + о^,

с применением которой построена L-устойчивая схема четвертого порядка Уп+1А = Уп + rih + r2k2 + r3k3 + r4fc4 + r5fc5 + r6fce + r7k7.

Здесь ki,l < 1 < 6 - стадии схемы (10).

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

lljfo+i - Утн-мП < е-

Для контроля точности (10) с коэффициентами (И) использовались следующие коэффициенты вспомогательной схемы

г2 = +3.710707537872, г3 = +10.49128588514, rs — - 1.134493098307, r6 = +11.83556197284,

Г1 = +0.2745513068351, г4 = +0.1481385540061, г7 = -11.21220029992, ап = +1.009872457390, а1Ь = - 0.0891885600729.

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

Исследована (2, 2)-схема вида

Dnki = hf(yn), Dnk2 = hf(yn + 0ki)+akh

(12)

где

Бп = Е- аКА„, (13)

а матрица Ап представима в виде.Лп = д/(уп)/ду + АД, + 0(А2). Здесь Вп-некоторая матрица, не зависящая от размера шага интегрирования. Представление Д, в виде (13) обеспечивает возможность замораживания и численного вычисления матрицы Якоби.

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

Для контроля точности (12) использовалось неравенство ПсгО1-^^! + Ъ2к2)\\ < е, 1< jn< 2,

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

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

где Д, представима в виде (13). Схема (14) имеет третий порядок точности, обладает Ь-устойчивой основной и промежуточной схемами при коэффициентах

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

а = +0.435866521508459, д = +0.4358665215084,

Р2 = +0.628266956983082, р3 = +0.75, Дц ~ +0.4358665215084,

Д$2 = +0.230800145158208, а32 = - 1.08551130465539.

Неравенство для контроля точности получено с помощью дополнительной стадии

с применением которой построена численная схема второго порядка точности Уп+1,2 = Уп + + г2к2 + г3Аг3 + г4Лг4,

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

используется неравенство

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

Получена оценка ошибки метода второго порядка при расчетах по схеме (4) и оценка ошибки схемы (15) при расчетах по методу третьего порядка.

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

Метод четвертого порядка имеет вид (6), а метод третьего порядка определяется (4, 2)-схемой

где внутренняя схема имеет вид

Уп+1 ~ Уп + &А +

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

^ _ 3- 57а + №а?- 528а3 + а(4а- 1)2(-162о2 + 1Э6а- 11)

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

Набор коэффициентов схемы (16), соответствующий коэффициентам (7), имеет вид

а = +0.219669914110089, рг = +0.219669914110089, р2 = +0.4126450787451, р3 = +0.5107726296546, щ = +0.08181996293789, #ц = +0.219669914110089, /?32 = 0.53033008588991, а32 = - 9.676674665135,

Данная пара методов может быть использована для построения алгоритмов интегрирования переменного порядка и шага.

В приложении приведены использованные в расчетах тестовые примеры из химической кинетики.

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

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

• Доказаны теоремы о максимальном порядке точности (т, к) -методов с аналитическим вычислением матрицы Якоби при к = 1, 2, 3 и для (т, к)-методов с замораживания матрицы Якоби при = 2 и 3. Построены методы максимального порядка с ¿-устойчивыми промежуточными схемами.

• Построены (т,к)- методы с аналитическим вычислением матрицы Яко-би с второго по пятый порядок точности включительно, обладающие ¿-устойчивыми промежуточными схемами.

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

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

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

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

1. Двинский, А. Л. Об одном алгоритме интегрирования жестких систем ОДУ / А. Л. Двинский // Информатика и информационные технологии.-Красноярск, 1999 .- С. 83-84..

2. Двинский, А. Л. (т,2)-методы решния жестких систем ОДУ / А. Л. Двинский // Информатика и информационные технологии.- Красноярск, 2000.- С. 77-78.

3. Двинский, А. Л. Исследование (3,2) -метода решения жестких систем / А. Л. Двинский // Наука. Тезисы. Инновации. Региональная научная конференция студентов, аспирантов, молодых ученых. Тезисы докладов.- Новосибирск, 2001.- ч. 1 .- С. 57-58.

4. Двинский, А. Л. (4,2)-метод решения жестких систем ОДУ с Ь-устойчивой внутренней схемой / А. Л. Двинский //Информатика и информационные технологии.- Красноярск, 2001.- С. 53-62.

2 5 4 9 9

5. Двинский, А. Л. (4,2)- метод третьего порядка для решения жестких систем / А. Л. Двинский, Е. А. Новиков// труды международной конференции RDAMM-2001, Вычислительные технологии, 2001.- т. б- ч.2.-С. 470-474.

6. Двинский, А. Л. (4,2)-алгоритм интегрирования жестких систем ОДУ / А. Л. Двинский// Материалы XXXIV научной студенческой конференции.- Красноярск, 2001.- С. 53-59.

7. Двинский, А. Л. L-устойчивая (б,3)-схема пятого порядка точности для решения жестких систем/ А. Л. Двинский// IV Всероссийская конференция молодых ученых по математическому моделированию и информационным технологиям. Программа и тезисы докладов.- Красноярск, 2003.- С. 21.

8. Dvinskiy, A.L. ( 4, 2)-Method of Order 3 for Solving Stiff Systems /A. L. Dvinskiy, E. A. Novikov//AMSE, Advances A.- vol. 40.- №3.-2003.-Р. 61-70.

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

10. Двинский, А. Л. Алгоритм интегрирования переменного порядка и шага для решения жестких систем / А. Л. Двинский, Е. А. Новиков// Вестник КГТУ.-выпуск 33.-Красноярск, 2004.- С. 79-87.

11. Двинский, А. Л. (6,3)-схема пятого порядка для решения жестких систем/ А. Л. Двинский// V Всероссийская конференция молодых ученых по математическому моделированию и информационным технологиям. Программа и тезисы докладов.- Новосибирск, 2004.- С. 16-17.

12. Двинский, А.Л. L-устойчивая (6,3)-схема пятого порядка точности для решения жестких систем /А.Л. Двинский, Е.А. Новиков // Вычислительные технологии.-т.9.- Вестник КазНУ,- №3(42).- 2004. - С.228-234.

Усл. печ. л. 1. Тираж 120 экз. № заказа 160

ООО "Печатные технологии" ул. Ленина, 113, офис 413, телефон: (3912) 23-52-73, 29-63-58

Введение

Глава

Исследование (то, к) -методов

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

2. (то, к) -методы.

3. Максимальный порядок точности (то, к) -методов с аналитическим вычислением матрицы Якоби

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

5. Построение неравенства для контроля точности вычислений

6. Алгоритмы интегрирования.

Глава

Алгоритмы интегрирования с аналитическим вычислением матрицы Якоби

1. (2,1)-метод второго порядка точности.

2. (2,2)-метод второго порядка с внутренней L -усточивостыо

3. (3,2)-метод третьего порядка с внутренней L -устойчивостью

4. (4,2)-метод третьего порядка с внутренней L- устойчивостью

5. (5,2)-метод четвертого порядка с внутренней L -устойчивостью

6. (6,3)-метод пятого порядка.

7. (6,3)-метод пятого порядка с внутренней L -устойчивостью

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

Глава

Алгоритмы с замораживанием матрицы Якоби

1. (2,2)-метод второго порядка

2. (3,2)-метод третьего порядка.

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

Глава

Алгоритмы переменного порядка

1. Метод переменного порядка с переключением по оценке локальной ошибки

2. (т, к) -методы с одной матрицей для нахождения стадий

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

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

Стремление ко все более точному описанию физических процессов приводит к постоянному росту размерности и жесткости соответствующей системы дифференциальных уравнений и предъявляет все возрастающие требования к методам интегрирования. Применение многошаговых схем в некоторых случаях нежелательно из-за эффекта "срезания экстремумов". Реализация неявных методов Рунге-Кутты сложна и приводит к итерационному процессу для нахождения стадий. В работах [1]-[4] предложен способ интегрирования обыкновенных дифференциальных уравнений, основанный на локальном многочленном приближении, где также используется итерационный процесс. Применение явных методов типа Рунге-Кутты приводит к обременительному ограничению на размер шага интегрирования ввиду ограниченного размера области устойчивости. Построение явных методов с расширенными и согласованными областями устойчивости [5, 6] не способно полностью решить проблему.

При решении задачи Коши для жестких систем обыкновенных дифференциальных уравнений широкое распространение получили методы типа Розенброка [7]. Данные численные формулы получены из полуявных методов типа Рунге-Кутты [8, 9], в которых при решении нелинейной системы алгебраических уравнений используется одна итерация метода Ньютона (см., например, [10]). В результате при вычислении каждой стадии вместо нелинейной системы алгебраических уравнений нужно решать линейную систему, а требуемая точность вычислений достигается за счет подходящего выбора величины шага интегрирования. Для снижения вычислительных затрат при решении линейной системы обычно используют LU -разложение. Другие способы решения линейных систем можно найти в работах [11]-[21].

Максимальный порядок точности методов типа Розенброка не превышает (к + 1) , где число к задает количество вычислений правой части системы дифференциальных уравнений на шаге. Несмотря на более сложную (по сравнению с явными методами) реализацию и необходимость вычисления и обращения матрицы Якоби, методы Розенброка обладают неограниченной областью устойчивости, их применение позволяет снять ограничения по устойчивости на размер шага интегрирования.

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

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

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

Известно, что при замораживании матрицы Якоби в методах Розенброка максимальный порядок точности не превышает двух [24].

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

В настоящее время на базе (га, к) -методов с замораживанием матрицы Якоби удалось построить алгоритмы интегрирования до четвертого порядка порядка точности [26]. Показано, что максимальный порядок точности (га, к) -методов равен (к + 2) в случае аналитического вычисления матрицы Якоби и (& + 1) при ее замораживании. Здесь к - количество вычислений правой части на шаг интегрирования. Построены (т, к) -схемы максимального порядка точности, обладающие свойством L -устойчивости.

При построение методов интегрирования жестких систем необходимо учитывать, что классический порядок точности метода, основанный на сравнении разложений приближенного и точного решений в ряды Тейлора, в случае жестких систем оказывается завышенным по сравнению с фактическим. На это явление впервые было указано в работе [27], а затем в [28] было введено понятие В -сходимости, способное предоставить оценку глобальной ошибки метода, не зависящую от жесткости решаемой задачи. Как показано в [10, с.267], методы Розенброка не могут быть В -сходящимися.

Наряду с В-сходимостью было введено понятие BSI -устойчивости [29], которое играет существенную роль для В-сходимости и связано с поведением внутренних схем. В [22, с. 194] установлено, что однократно неявные методы Рунге-Кутты, к которым относятся и методы Розенброка, BSI -устойчивы при некотором ограничении на шаг интегрирования. Данное ограничение можно обойти, потребовав L -устойчивость промежуточных численных формул.

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

Следует заметить, что свойства В -сходимости и BSI -устойчивости зависят от класса решаемых задач, и проверка этих условий для достаточно общих классов задач трудоемка, в то время как во многих случаях применение L -устойчивых схем с L -устойчивыми внутренними схемами, построенных на основе классического порядка точности, вполне себя оправдывает.

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

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

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

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

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

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

Основные результаты диссертации опубликованы в работах [30]-[41].

Заключение

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

• Доказаны теоремы о максимальном порядке точности (ш, к) -методов при к = 1,2,3 с замораживанием и без замораживания матрицы Якоби, а также с возможностью ее численной аппроксимации. Показано, что методы максимального порядка могут обладать свойством внутренней L -устойчивости.

• На основе (ш, к) -методов с аналитическим вычислением матрицы Якоби построены алгоритмы второго, третьего, четвертого и пятого порядков точности, обладающие свойством внутренней L -устойчивости.

• На основе (т, к) -методов с замораживанием матрицы Якоби и ее численной аппроксимацией построены алгоритмы второго и третьего порядков точности, обладающие свойством внутренней L -устойчивости.

• На десяти тестовых примерах из области химической кинетики показано почти полуторократное повышение эффективности методов за счет L -устойчивости промежуточных численных схем.

• Исследована возможность построения методов переменного порядка и шага на основе (ш, к) -методов. Проведены численные расчеты, подтверждающие повышение эффективности алгоритмов интегрирования за счет внутренней L -устойчивости.

1. Залеткин, С.Ф. О построении многочленных приближений при численном решении обыкновенных дифференциальных уравнений /С.Ф. Залеткин, Татевян С.К., Сорокин Н.А.// Вычислительные методы и программирование.-т.2.-2001.-С.56-64. (http://num-meth.srcc.msu.su)

2. Залеткин, С.Ф. Формула численного интегрирования Маркова и ее применение в ортогональных разложениях /С.Ф. Залеткин, Татевян С.К., Сорокин Н.А.// Вычислительные методы и программирование.-Т.2.-2001.-С.131-158. (http://num-meth.srcc.msu.su)

3. Залеткин, С.Ф.Численное интегрирование обыкновенных дифференциальных уравнений с использованием рядов Чебышева /С.Ф. Залеткин, Татевян С.К., Сорокин Н.А.// Вычислительные методы и программирование.-т.З.-2002.-С.52-81. (http://num-meth.srcc.msu.su)

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

5. Новиков, Е. А. Явные методы второго порядка с согласованными областями устойчивости/ Е. А. Новиков, JL Н. Контарева// Вычислительные технологии, 2001.- т. 6.- JYH.- С. 40-50.

6. Rosenbrock, Н.Н. Some general implicit processes for the numerical solution of differential equations /Н. H. Rosenbrock// Computer, №5,1963. p. 329-330.

7. Butcher, J. C. Implicit Runge-Kutta Processes /J. C. Butcher// Math. Сотр. 1964.-18.- С. 50-64.

8. N0rsett, S. P. Semi-Explicit Runge-Kutta Methods//S. P. N0rsett .-Tech.Rept.-Math and Cornp.-№6/74.-1974.-Univ. of Trondheim.

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

10. Axelsson, О. A general incomplete block-factorization method/O.Axelsson //Linear Alg. Appl., 74 P.179-190, 1986.

11. Axelsson, O. Iterative Solution Methods/O. Axelsson// Cambridge University Press.- New York.- 1994.

12. Bank, R. E. An analysis of the composite step biconjugate gradient method/ R. E. Bank, T. F. Chan// Numer. Math. 66: P.295-319.- 1993.

13. Brezinski, C. Projection Methods for Systems of Equations/ C. Brezinski//North-Holland.- Amsterdam.- 1997.

14. Dongarra, J. J. Numerical Linear Algebra for High-Performance Computers/J. J. Dongarra, I. S. Dulf, D. C. Sorensen, H. A. van der Vorst// SIAM.- Philadelphia, PA.- 1998.

15. Elman, H. С. A stability analysis of incomplete LU factorizations/H. C. Elman// Math. Сотр.- 47.-P. 191-217.- 1986.

16. Hackbusch, W. Iterative Solution of Large Linear Systems of Equations/ W. Hackbusch// Springer Verlag.- New York.- 1994.

17. Hestenes, M. R. Methods of conjugate gradients for solving linear systems/ M. R. Hestenes,E. L. Stiefel//J. Res. Nat. Bur. Stand., Section В.- 49.-P.409-436.-1952.

18. Lanczos, C. Solution of systems of linear equations by minimized iterations/ C. Lanczos//J. Res. Nat. Bur. Stand.-49.-P.33-53.- 1952.

19. Meurant,G. Computer solution of large linear systems/ G. Meurant//North-Holland.- Amsterdam.- 1999.

20. Greenbaum, A. Iterative Methods for Solving Linear Systems/ A. Greenbaum//SIAM.- Philadelpha, PA.-1997.

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

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

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

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

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

26. Prothero, A. On the Stability and Accuracy of One-Step Methods for Solving Stiff Systems of Ordinary Differential Equations/A. Prothero, A. Robinson//Math, of Coput.-vol.28.-1974.-P.145-162.

27. Frank, R. The Concept of Б-convergence /R. Frank, J. Schneid, C. W. Ueberhuber//SIAM J. Numer. Anal.- vol. 18.-1981.-P.753-780.

28. Frank, R. Stability Properties of Implicit Runge-Kutta Methods/R. Frank, J. Schneid, C. W. Ueberhuber//Bericht Nr 52/82.-Institut fur Numerische Mathematik, TU-Wien.-1982.

29. Двинский, А. Л. Об одном алгоритме интегрирования жестких систем ОДУ / A. JI. Двинский // Информатика и информационные технологии.- Красноярск, 1999 .- С. 83-84.

30. Двинский, A. J1. (т,2)-методы решния жестких систем ОДУ / A. JT. Двинский // Информатика и информационные технологии.- Красноярск, 2000.- С. 77-78.

31. Двинский, A. JI. Исследование (3,2) -метода решения жестких систем / A. JI. Двинский // Наука. Тезисы. Инновации. Региональная научная конференция студентов, аспирантов, молодых ученых. Тезисы докладов.- Новосибирск, 2001 .- ч. 1 С. 57-58.

32. Двинский, A. JI. (4,2)-метод решения жестких систем ОДУ с L-устойчивой внутренней схемой / A. JI. Двинский //Информатика и информационные технологии.- Красноярск, 2001.- С. 53-62.

33. Двинский, А. Л. (4,2)- метод третьего порядка для решения жестких систем / A. J1. Двинский, Е. А. Новиков// труды международной конференции RDAMM-2001, Вычислительные технологии, 2001.- т. 6.- ч. 2.- С. 470-474.

34. Двинский, A. JI. (4,2)-алгоритм интегрирования жестких систем ОДУ / A. JI. Двинский// Материалы XXXIV научной студенческой конференции.- Красноярск, 2001.- С. 53-59.

35. Dvinskiy, A.L. ( 4, 2)-Method of Order 3 for Solving Stiff Systems /А. L. Dvinskiy, E. A. Novikov//AMSE, Advances A.- vol. 40.- №3.-2003.- P. 61-70.

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

37. Двинский, A. JI. Алгоритм интегрирования переменного порядка и шага для решения жестких систем / A. JI. Двинский, Е. А. Новиков// Вестник КГТУ.-выпуск ЗЗ.-Красноярск, 2004.- С. 79-87.

38. Двинский, АЛ. L-устойчивая (6,3)-схема пятого порядка точности для решения жестких систем /АЛ. Двинский, Е.А. Новиков // Вычислительные технологии.-т.9.- Вестник КазНУ.- №3(42).- 2004. С.228-234.

39. Новиков, Е. А. М,3-метод третьего порядка для жестких неавтономных систем ОДУ/ Е. А. Новиков, М. И. Голушко// Вычислительные технологии, 1998.- т. 3.- №3.- С. 48-54.

40. Новиков, Е. А. Одношаговые безитерационные методы решения жестких систем /Е. А. Новиков.- Новосибирск: ВЦ СО РАН, дисс. д.ф.м.н. .- 1991.- 327с.

41. Новиков, Е. А. Построение (т, к )-методов решения линейных систем обыкновенных дифференциальных уравнений/ Е. А. Новиков// Вычислительные технологии, сборник научных трудов ИВТ СО РАН: Новосибирск .- 1993.- т. 2.-№7.- С. 140-155.

42. Richardson, L. F. The Deferred Approach to the Limit/L. F. Richardson// Phil. Trans. A .1927.- vol. 226 .-P.299-349.

43. Richardson, L. F. The Approximate Arithmetical Solution by Finite Differences od Physical Problems Including Differential Equations, with an Application to the Stresses in a Masonry Dam/L. F. Richardson//Phil. Trans. A .1910.- vol. 210 .-P.307-357.

44. Cescino, F. Numerical Solution of Initial Value Problems /F. Cescino, J. Kuntzman//Prentice-Hall, Englewood Cliffs, A, New Jersey, 1966.

45. Fehlberg, E. Low-Order Classical Runge-Kutta Formulas with Step Size Control and their Application to Some Heat Transfer Problems/E. Feglberg//NASA Technical Report 315.-1969.

46. England, R. Error Estimates for Runge-Kutta Type Solutions to Systems of Ordinary Differential Equations/ R. England//The Computer J. .- 1969.-vol.12.-P. 166-170.

47. Sarafyan, D. Error Estimation for Runge-Kutta Methods through Pseudo-Iterative Formulas/D. Sarafyan//Techn. Rep. No.l4.-1966.-Lousiana State Univ.-New Orleans.

48. Dormand, J. R. A Family of Embedded Runge-Kutta Formulae/J. R. Dormand, P. J. Prince//J. Сотр. Appl. Math. .- vol.6 .-1980.-P.19-26.

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

50. Shampine, L. F. The Art of Writing a Runge- Kutta Code/ L. F. Shampine, H. A. Watts//II, Appl. Math. Comput.,1979, vol.5, P. 93-121.

51. Gear, C. W. The Automatic Integration of Stiff Ordinary Differential Equations/C. W. Gear//Inform. Proc.-1969.-P.187-193.

52. Byrne, G. D. Stiff ODE Solvers: a Review of Current and Coming Attractions/G. D. Byrne, A. C. Hindmarsh//J. Сотр. Phys.-1986.- №70.-P.l-62.

53. Артемьев, С. С. Алгоритм переменного порядка и шага для численного решения жестких систем обыкновенных дифференциальных уравнений /С. С. Артемьев, Г. В. Демидов//ДАН СССР.-1978.-Т.238,№3.-С. 517-520.

54. Демидович, Б. П. Численные методы анализа/Б.П. Демидович, И. А. Марон, Э. 3. Шувалова.-М.:Наука.-1967

55. Демидов, Г. В. Оценка ошибки одношаговых методов интегрирования обыкновенных дифференциальных уравнений. /Г. В. Демидов, Е. А. Новиков//Числ. мет. мех. сплош. среды.-Новосибирск.-1985.-т. 16.-т.-С. 27-39.

56. Новиков, Е. А. Апроксимация матрицы Якоби в (М. К.)-методах третьего порядка точности // Е. А. Новиков, М. И. Голушко, Ю. А. Шитов .- Препринт Ml : Красноярск, ВЦ СО РАН .- 1991 .- 36 с.

57. Swart, J. J. В. On the Construction of Error Estimators for Implicit Runge-Kutta Methods/ J. J. B. de Swart, G. Soderlind //MAS-R9704.-1997.

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

59. Демидов, Г. В. Исследование точности неявных одношаговых методов // Г. В. Демидов, JI. А. Юматова .- Препринт №11 : Новосибирск, ВЦ СО АН СССР.- 1976 .- 22 с.