Численное интегрирование жестких локально-неустойчивых систем обыкновенных дифференциальных уравнений тема автореферата и диссертации по математике, 01.01.07 ВАК РФ

АКАДЕМИЯ НАУК СССР Вычислительный центр

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

РОДИОНОВА ОКСАНА ЕВГЕНЬЕВНА

Численное интегрирование жестких локально-неустойчивых систем обыкновенных дифференциальных уравнений

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

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

Москва 1991

Работа выполнена в Ордена Ленина Институте химической

физики им. Н.Н.Семенова АН СССР

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

старший научный сотрудник ПАВЛОВ Борис Васильевич. Официальные оппоненты: доктор физико-математических наук,

профессор ПОВЗНЕР Александр Яковлевич.

кандидат физико-математических наук старший научный сотрудник ФИЛИППОВ Сергей Сергеевич.

Ведущая организация: Кафедра вычислительной математики и

программирования МАИ

Защита состоится " Щ " ИС&Ъ/ьЛ! 199'г. в Л_часов

на заседании специализированного совета Д 002.32.01 при Вычислительном центре АН СССР по адресу 11796? Москва ГСП 1 ул. Вавилова д. 40

С диссертацией можно ознакомиться в библиотеке Вычислительного центра АН СССР

Автореферат разослан _« 199/г.

Ученый секретарь специализированного совета

д. ф.-м. н. <Zjffa/( Терентьев Е.Д.

%

'. -i . »ш

ТДРЛ

сертаций

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

актуальность темы. Система обыкновенных дифференциальных уравнений /ОДУ/ x=f(x,t); x(tQ )=xQ - локально-неустойчива, если ее решение x(t,zQ) проходит через область неустойчивости вариационной матрицы fx(x,t)=Df(x,t)/Dx. В этой области происходит возрастание погрешности численного интегрирования по мере движения вдоль траектории решения, и это обстаятель-ство требует повышенной точности самого алгоритма интегрирования и аккуратности вычисления его локальной погрешности.

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

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

Целью работы является разработка схем численного решения задачи Коши для жестких, в том числе и локально-неустойчивых систем ОДУ, на основе метода локальной линеаризации / МЯЛ /, с обеспечением надежного контроля за точностью численного решения.

Научную новизну работы представляют :

- формулы для вычисления погрешности численного решения систем ОДУ полученные с помощью рекуррентных соотношений для матриц exp(At) и C(t)=A~1(exp(At-1);

- метод оценки правой границы спектра матрицы At с использованием Tr(exp(Ati)); t

- метод вычисления интегралов V(s) ds основанный на

о

сингулярной аппроксимации ядра exp(A(t-s)) :

AS m , ,

е «ct^HZ^1 1

где ?=t-s; 0(£)-о0общенная функция;

- применение степенной аппроксимации для быстрорастущих под-интегральных функций р. (s);

- алгоритм решения интегрального уравнения МЛЛ для системы

г А(1;-б)

г=СГ(хпгп)+«'е ^ •

где м.(г(з))=Г(хп+2,гп+з)-Аз-Г(2п>гп) , 1=1-^.

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

Апробация работы. Основные результаты диссертационной работы докладывались и обсуждались на V Всесоюзной конференции " Математические методы в химии " (г. Грозный, 1985г.), на I Всесоюзной конференции " Новые подходы к решению дифференциальных уравнений " (г. Дрогобыч, 1987г.), на научных семинарах в ИПМ им. М.В.Келдыша АН СССР, на научных семинарах в математическом отделе ИХФ АН СССР.

Публикации. По теме диссертации опубликовано четыре работы, список которых приведен в конце автореферата.

Структура и объем работы.Диссертация состоит из введения, шести глав, заключения, и списка литературы, включающего 40 наименований. Общий объем диссертационной работы составляет 96 листов машинописного текста.

СОДЕРЖАНИЕ РАБОТЫ

Во введении обосновывается актуальность выбранного направления исследования.

Первая *- -два. "Литературный обзор".Рассмотрены основные метода .используемые для решения задачи Коши для жестких систем ОДУ. Подробно рассмотрены основы метода локальной линеаризации, разработанные до начала данного исследования.

Вторая глава " Интегрирование линейных систем ОДУ с постоянными коэффициентами". В основе МЛЛ лежит решение линейных систем ОДУ, поэтому, они подробно рассмотрены в работе.

Первый параграф " Решение систем вида х=Ах". Решение

такой системы ищется в виде х(1;)=ехр(А1;)х(1;0), где х(г0) -начальное условие. Вычисления матрица ехр(Аг) основано на функциональном уравнении:

А!;.. ЬХ? е =е е ""

Пусть , У(А110)-приближенное значение матрицы

ехр(А110) и вычисляется с помощью т-членов ряда

еАЙ° а У(АЬ0)=Е + АЪд /1!+ ... + /т! при условии, что |А110| < е0< 1. Обозначим через У(АЬП)- приближенное значение матрицы ехр(А1^). Значение У(А1гГ1) вычисляется с помощью рекуррентной формулы

(1) У(АИП)=У(АЬП_1) У(А1гп_^) ,где

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

2п£т+1

(2) АХ(ЬП)= 1й?°ут В^сО^)

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

(3) 14x0^)1 < «М^Н где ет-заданная точность интегрирования.

Например, при ет=е0=10~3,т=3 , число возможных шагов по формуле (1) п =< 40.

Существует возможность уточнить значение матричной экспоненты:

(4) е * - У(АЙП){Е + В0 +[ -^ттг В0]2 +...}

Добавление даже одного члена ряда в формуле (4) позволяет значительно увеличить число вычислений по рекуррентной формуле (1), и , например, при ет=е0=1СГ3, т=3 , это позволяет решать очень плохо обусловленные системы, такие что

где Алах и \rnin - максимальное и минимальное ( по модулю)

собственные значения матрицы А.

Второй параграф "Решение системы вида х=Ах+Ь".Решение данной системы ищется в виде Хд^Хд+г, где г(т:)=С(т:)Ъ,

1 кЫ а, -1 Ат С(т)= /е 1 'йз = А (е - Е) о

Для дискретных значений 1^= t-tn, С^А!^) - приближенное значение матрицы С (А!^), такое что

2П£пн-1

С(АЬП)=С1 (Айп)+А"1 [ечВ - Е1 , -щитт-Рекуррентное соотношение (5) для вычисления С1(АЬд)

(5) (Айп)=С1 (АЬ^ )+СЕ+С1 (А11п_1 ШО, (Айп_1)

эквивалентно рекуррентному соотношению (1) и , следовательно, формулы для определения погрешности интегрирования в параграфе 1 применимы и в данном случае.

Третий параграф. Асимптотические свойства матриц У(А11П) и С1(А!^). Показано, что для устойчивой матрицы А :

|е -V ) | =»0 и |С(А11П)-01(А11Г1 )| • О 11 П =» со П =» 00

Поэтому, на практике, для интегрирования жестких линейных систем с постоянными коэффициентами можно ограничиться использованием матриц УСА!^) и С^Айд).

Четвертый параграф "Дробно-рациональное приближение матрицы еА1:". Если известно, что 11е(?^1;)<М , где М ~ 1, за начальное приближение ехрСА!^) принимается

(6) 6АЛ°=(Е-аА110)(Е-7АЬ0)~Р > 1^=1/2".

При р=2, а=1+У2~ , 7=1+У0.Ь ,при р=3, а=(1+УЗ)/2 , 7=(3+уЗ)/2 . В случае интегрирования жестких систем это позволяет значительно увеличить начальный шаг Для р=2 показано, что в случае устойчивой матрицы А при п=10, матрица ехр(АЪ02к) будет удовлетворительно ( с точностью не хуже чем 1.8-10"3) приближать точное значение при к ? 5. При наличии положительных собственных значений в спектре матрицы.

А , приближение с заданной относительной точностью е может быть получено на интервале <- (0,1.5-е~122п).Таким обра-ом ,во второй главе обосновано применение рекуррентных формул (1) и (5) для получения решения систем ОДУ с постоянными коэффициентами , в том числе жестких и неустойчивых.

Третья глава " Интегрирование линейных систем ОДУ с

переменными коэффициентами" Первый параграф." Решение системы вида х=Ах+<р(1;)" Решение данной системы ищется в виде хп=х0+г , где

2=С(т)[ср(1;0)+Ах0]+<11 , т=1;-1;0 , х(^)=х0

х

г А(Т-з)

<1.,= £ е Ф(э)с1з , ф(э)=ф(1;0+8)-<р(1;0).

А(Т-з)

При вычислении «I., приближается ядро е на некотором фиксированном интервале т=1г, s40.li] , в предположении о том, что спектр матрицы А расположен вблизи действительной оси спектральной плоскости . Тогда , для устойчивой матрицы А , при больших 1г ядро сосредоточено в окрестности точки з^й , так что

11

г А (11-3)

(7) ^ ф(в) с1з а С(П)ф(1х) О

Для сохранения этого важного асимптотического свойства в аппроксимацию ядра введена обобщенная функция ©(1п-з). На зЧ0,11] приближение имеет вид:

(8) е = Нт(Аз) = 0^(8)+ I р^ООз^^ОНР^з)

1=1

причем коэффициенты % и р1 определяются в узлах э:{11,11/2,

1г/4,... ,1г/2ш> системой уравнений Ьр Ьр

/Нт(Аз) йз = /еАз йз , р=0,1,...ш , 11р=2 11 0.0

Данная аппроксимация обладает следующими свойствами:

1.Для устойчивой матрицы А коэффициенты р1(¡1) экспоненциально стремятся к нулю при 1г-*оо , таким образом .приближение

(8) обладает требуемым асимптотическим свойством (7)

2.Точность аппроксимации (8) проверяется по точности вычисления интегралов с приближенным ядром на векторе вида

к1

ф±(8)=(в/Й) Хф(Ь) и зависит от показателя к^.

3.Качество предложенных формул можно оценить , проводя сравнение собственных значений моментов матричной экспоненты, вычисленных по точным и приближенным формулам. Непосредственные вычисления показали , что при т=4 с помощью Н4(Лз) моменты до 50 порядка включительно вычисляются с погрешностью менее 1% на интервале Лт <- (-со,1].

4.Аппроксимация (8) удовлетворительно приближает ядро лишь при условии

(9) Мк1

В работе показано , что для проверки условия (9) достаточно проверить условие

(10) СТге2А11-ТгеА11][Тге2А11-ТгеА11+0.25] < 25 Второй параграф. "Решение систем вида х=А(1;)х+ф(1;)" На

конкретном примере показан "механизм" получения качественно неверных решений при интегрировании жестких локально-неустойчивых систем ОДУ конечно-разностными методами. Рассмотрена схема решения линейной системы ОДУ вида

х=А(1)х+ф(г) , х(г0)=х0 с использованием аппроксимации (8).

Четвертая глава "Схемы метода локальной линеаризации" Рассматривается задача Коши для системы ОДУ вида

(11) х=Г(хД) , х(г0)=х0

Существует две основные трудности при решении жестких локально-неустойчивых систем ОДУ широко распространенными конечно-разностными методами:

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

2.схема т-го порядка плохо приближает компоненты х(т), которые на данном интервале т [0,1г1 растут быстрее чем 1го, а это приводит к потере контроля за точностью численного решения в области неустойчивости матрицы {Б1/Вх>.

Первый параграф."Интегральное уравнение МЛЛ". Для

системы (11) запишем уравнение в приращениях

¿(г)=Г (х^^НАг+ц^.т) , и(0)=0

где введены обозначения:т=1;-1;п ,х(1;п)=хп, г(т)=х(1;)-х(1;п) ,

А=Г'х(х0,г0) , |1(2,т;)=Г(хп+2,гп+1;)-Г(хпД11)-А2. Система

приводится к эквивалентному интегральному уравнению

г

(12) 2(т)=С(т)Г(хпДп)+ХеА(т-£3)м.(2(з)) ¿д,

которое служит основой для построения всех схем ММ.

Второй параграф."Схемы метода локальной линеаризации" Набор неявные схемы МЛЛ можно получить, используя различные аппроксимации функции ц(г(з)). В работе использовуется степенное приближение. Показатель степенного роста ^ определяется по-компонентно в процессе решения неявной схемы. Пусть т=11 ,

(13) ц1(а(а),з)=а1з 1 Коэффициенты а. и к . определяются в точке з=Ь.

к1

р-гШ)^)]! к,-=-

1 ц.(з(И),И)

Таким образом, получаем неявную схему МЛЛ вида: 2(Ь)=0(11)Г(хГ111^)+Г(Ь.к)ц(г1Ь)

г А(Ь-в) к. к.

О

к. к.

где в(з /Ь )-диагональная матрица с элементами вида к- к.

(0)^=6^.(3 Х/Ь. х) ^ - символ Кронекера. Для вычисления Г(1г,к) используется приближение ядра (8).

Третий параграф."Разностный аналог схем МЛЛ". Если интегрирование ведется с шагом Ь. , таким что |А1г|<1 ,то вычисление ь

г А(Ь-а)

Je (1(г(з)) йз о

существенно упрощается , т.к. ехр(АЬ) можно вычислить

разложением в ряд

АЪ. т т

е =Е+АЬ+...+А 11 /т!

с небольшим числом членов (обычно т равно двум или трем). Применение разностного аналога схем МЛЛ существенно сокращает число арифметических операций и возможно тогда , когда система ОДУ стала нежесткой.

Четвертый параграф."Выбор шага интегрирования". Это особенно важно для численного решения задачи Коши для жестких локально-неустойчивых систем. Рассмотрим некоторую неявную схему МЛЛ. (14) г№)-0(Ь)Г(хп,1п)+Р(11,г(Ь))

IV ГУ

где г (11)-приближенное значение г(11). О свойствах г(Л) можно судить по его реакции бг на малое возмущение правой части

(14)-Д:

-1

©и=[Е-Р*(г.Ь)] А

и

Если спектральный радиус матрицы р^ меньше единицы,

то матрица [Е-Р^и.й) Г1 заведомо ограничена, и реакция е2

будет малой. В соответствии с этим, решение а(Ь) считается "правильным", если для всех т с [О.Ь] выполняется условие

(15) Н(Р£(а<г)рг)) < 1

При вычислениях, условие (15) проверяется только при 1=1^, где ^-пробные значения 1г. Поэтому, условие (15) можно записать

(16) Н(Р^(а(Ь),Ь)) < 1

Это условие автоматически проверяется решением (14) методом прямых итераций:

~ (пн-1) ~ (т) (0)

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

Пятый параграф."Оценка точности интегрирования". Для

исходного интегрального уравнения

ъ.

х(гп+1)=еАЬх(гп)+-ГеА(Ь_3[1(х(з),з) йэ

после аппроксимации ц(х(з),з)), в общем случае

получим схему :

АЬ.

тогда уравнение для глобальной ошибки в т. 1; +1 будет иметь ввд: ч ^

ехп+1 = [Е~\+1] [е

Локальная ошибка 8лок на шаге интегрирования Ь. складывается из ошибки, связанной с приближением ядра ехр(А(1г-з)) и ошибки , связанной с приближением ц(2(з),з), з с [0,Ы. Если решение г(11) на шаге 11 определяется как

Ь

то 0

11 Ч

^=4+1 №-в)в((а/Ь) )ц(2(1г),1г) йз + я(Ю

где 0^.,- ошибка в приближении ядра ехр(АШ-з)) функцией порядка т, q(h)-oшибкa, связанная с приближением функции г(1г).

Шестой параграф."Основные свойства схем МЛЛ". Приведены основные свойства набора неявных одношаговых схем МЛЛ второго порядка точности.

Пятая глава."Алгоритм". Разработан алгоритм МЛЛ для решения автономных систем ОДУ химической кинетики. Для этого класса жестких систем вычисление правых частей и якобиана системы возможно производить в стандартном виде. Формулы алгоритма с небольшими изменениями применимы и для неавтономных систем. В работе предлагается неоднородный алгоритм, состоящий из неявного метода Эйлера, для выбора начального шага интегрирования, разностного двухшагового метода, для интегрирования нежестких ситсем и набора неявных схем МЛЛ различного порядка для жестких систем ОДУ.

Первый параграф."Выбор начального шага". Для исходной системы ОДУ

(1Т) х=Г(х) , ха0)=х0

выбор начального шага осуществляется с помощью неявного

метода Эйлера

(18) ' х^+ЬоКх,)

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

Шаг выбирается так, что й £ е«1. После выбора шага , решение х^ уточняется по двухшаговой неявной схеме вида:

-^(Хд)]]^/»}11-1^) ,где (п)- • (п)-, г • г (п)

^[г^-Г-и., ^/[^(^-[ЗД ЫЛ*]]

Второй параграф."Неявная разностная схема". Счет по двухшаговой схеме ведется в предсказывагаце-исправляющем режиме. Для вычисления предсказывающего значения используется

локальная экстрополяция на отрезке [0,1гп+]ап_1 ]

хЭ1(^+з)=хп_1 +(ь+3)Г.(хп_1)+У1.(^1) ,1

.где

VI• ^пИ^п • У1 ^-^п-Т "Г(хп-1

Имея значение х®+1 можно оценить величину |Л.тахЬ.|, так как согласно параграфу 3 главы 4, данная разностная схема надежна при условии

(19) ^ < 1

Обозначим ). Выполнение одновременно условий

2 2

1г|ТгА., | < 1 и 11 |ТгА., I < 1 гарантирует выполнение условия (19).

Если в исходное интегральное уравнение

подставить полученное экстраполяционное .значение , то получим первое исправленное значение х

где ^-вычисляется с использованием квадратурной формулы Гаусса с четырьмя узлами. Если локальная ошибка

| велика, то для вычисления хп+1 используется неявная схема г,

(1+1) 1 г (1) Э

хп+1 =хп+1+;^2 (гп+з))-Г(х (гп+з)) аз

и для вычисления х^1 ^ также используется квадратурная формула Гаусса.

После того как закончены итерации, локальную ошибку на 1г0 определяем как разность между предсказывающим и исправляющим значениями. Если ошибка велика, то необходимо делить шаг. Если нарушено условие (19), то надо переходить на счет по МЛЛ.

Третий параграф." Вычисление С(1г0)". При требуемой точности интегрирования е, вычисляется из условия

Ь0=е|1^(х ¿Г1, тогда

СК^Н^Е + АЪф / 2! / 3!]

Четвертый параграф."Вычисление С(1г)". Если 1г-это шаг,с которым получено последнее решение по разностной схеме, для вычисления С№) применяется рекуррентная формула (5), т.е.

С(1г)=С(11/2)+[Е+С(]1/2)А]С(1г/2)

Несколько последних матриц С(Ь/21) запоминаются и используются далее для приближения ядра исходного интегрального уравнения.

Пятый параграф."Простая неявная схема МЛЛ". Эта схема имеет вид:

- 12 -

(20) а(т)=С(Ь)Г(хп)+С(11)и.(а(ш)) , 2(0)=С(ЮГ(Хд)

решается прямыми итерациями на кавдом шаге интегрирования.

Показатель сходимости схемы (20) Н(С(Ь)ц^(г)) определяется на

последовательности векторов ае® и вычисляется как й=|аега|, где

(т) (т-1) , . ,

ае =С(Ь)^(г)ае /|эеда1;|

Требование

(21) |ае(т)| < 1

ограничивает шаг интегрирования, гарантируя получение качественно правильного решения. Если вместо условия (21) выполняется более сильное условие

(22) |ае(т)||££

то неявная схема дает достаточно точное решение исходного уравнения хп+1=хп+а, и можно переходить к следующему шагу интегрирования. В качестве локальной ошибки 6лок можно принять величину Если условие (22) не выполнено, то необхо-

димо применение более сложной схемы ММ. Если нарушено условие (21), то необходимо делить шаг.

Шестой параграф."Схемы ММ". Рассматривается набор схем МЛЛ второго порядка точности. Счет, также как и в случае разностной схемы, ведется в предсказывающе-исправляющем режиме. Для вычисления предсказывающего значения рассмотрим исходную систему уравнений в приращениях

¿=Г(хп)+А2+ц(г)

где ц(г)=ц(хп,хп+1)=Г(хп+г)-Аг-£(хп) , 2=хп+1-хп , и строим степенную экстраполяцию це(2(з))

(23) ц®(г(а))=а1(в/11п) 1 ,з <- £0,1^]

II—2 71-1

¥

г _к1<У2+У1+У

К'= <У2+У1 >£ ^УА^Уг+У 1

к.=1п

Н^п-г^

13 -

1 /ш Г

/ 1 \-2+ЬП-1

Подставляя (2э) в исходное интегральное уравнение (12) получим

и

(24) 2э(Ю=С(ЮГ(Хп)+Г(П>к)=С(Ь)Г(хп)+£еА(11 3 (з/Ьп)К)(1з

Для вычисления интегралов типа Г(1г,к) используется приближение ядра ехр(А(11-з)) вида (8). Подставляя (8) в (24) получаем предсказывающее значение :

(25) ^

О

Неявная схема ММ , аналогично (25) будет иметь вид

(1) (1-1) (1-1)

(2б) > Е (1-1)

Решая неявную схему (26) одновременно определяем порядок ядра , т.е. величину т=т^. Итерационный процесс (26) осуществляется до сходимости итераций с точностью е1. После каждого вычисления х^;] происходит перевычисление ц(хп,х^^;]) по формуле

Седьмой параграф."Уточнение решение" Значение вектора

фактически является решением исходной системы ОДУ на

данном шаге йд. Однако, схема (26) не является консервативной

, так как для аппроксимации функции м.(г(з)), для каждой

компоненты 1 вычисляется свой показатель степенного роста

Для получения консервативного решения , приближаем х(1;п+з) на

з 10,Ьд] функцией вида'

(27) х(гп+з)=хп+С(з)Г(хп)+у(з).

к2, (1-1) где у±(в>-уО±(в/Ьп) } уО±=Гп 0)^+1 ) , .

- 14 -

кг.^А.уо+цс^,^ )]./У0.(Ьп) Подставляя аппроксимацию (27) в исходное интегральное уравнение (12) получим .

+у0.о((Б/11п) г)) Ла

Вектор хп+1 вычисляется с помощью квадратурной формулы Гаусса , а величину б=|хп+1-х^]| можно принять за локальную ошибку.

Восьмой параграф."Выбор нового шага интегрирования". Если локальная ошибка 6лок на данном шаге слишком мала,т.е елок<0.18, где е-заданная точность интегрирования, то шаг надо удваивать. Удвоение шага возможно при соблюдении двух условий :

(28) Не^Ьд^) < 1

(29) С(ЬГ1)А=АС(1пп)

Проверка условия (28) производится при помощи условия (10). Проверка условия (29) осуществляется непосредственно.Если локальная ошибка 6лок>£, то надо повторить интегрирование с более мелким шагом 1^=1^/2. Если удвоение шага невозможно, то необходимо оценить эффективность счета с данной матрицей А0. Возможно, новая линеаризация позволит значительно увеличить шаг.

Девятый параграф."Оценка эффективности счета с данной матрицей А0".Носит эмпирический характер, и подбирается так, чтобы новая линеаризация проводилась не очень часто.

Десятый параграф."Новая линеаризация".Так как в процессе счета по неявным схемам МЛЛ все время ограничивается правая граница спектра йе(А.+11п) < 1,то для вычисления начального значения матрицы С(АНЬ^)), где Ан=(БГ/1)х)х=хп , можно использовать дробно-рациональное приближение (6) матрицы ехр(Ан1г ),

С(АН110)=(Е-^7Ь0)"1 <ЬЕ -А^МЕ-А^)"1 Для вычисления матрицы (Е-А^!^)-1 используется модифицированный метод Ершова с выбором главного элемента.

Шестая глава."Тестовые примеры". Рассматриваются результаты рассчетов тестовых примеров по программе "БОЬУ 3",составленной на основе алгоритма, изложенного в предыдущей главе.

Приведенные тестовые примеры можно разделить на два класса:

1. Небольшие системы ОДУ, имеющие точное решение. Они позволяют оценить точность и эффективность работы отдельных блоков программы. Примеры подобраны так, чтобы они частично моделировали реальные процессы, возникающие в химической кинетике и связанные с этими процессами математические особенности ОДУ.

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

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

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

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

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

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

приближения ядра интеграле /eA^h"3^(s) ds до пятого

О

порядка включительно; для этих формул указан способ оценки точности приближения.

-Разработан способ определения правой границы спектра матрицы Ah по Tr(eAhi).

-Разработан алгоритм степенной аппроксимации сильно растущих компонент решения с определением коэффициента степенного роста в процессе интегрирования систем ОДУ.

-Предложен надежный алгоритм для решения жестких нелинейных систем ОДУ , как устойчивых так и. неустойчивых, вида:

x=í(t,x) , x(t0)=xQ с автоматическим выбором шага интегрирования , порядка схемы интегрирования , с автоматическим переключением со схем МЛЛ

на разностные схемы и обратно; в процессе инегрирования производится оценка локальной ошибки на каждом шаге h.

-Записаны уравнения накопления глобальной ошибки для используемых схем МЛЛ.

-Разработана экспериментальная программа "S0LV 3" и проведено испытание алгоритма на ряде жестких неустойчивых систем ОДУ.

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

ПУбЛИКАЦИИ ПО ТЕМЕ ДИССЕРТАЦИИ

1. Родионова O.E..Павлов Б.В. О численном интегрировании локально неустойчивых систем обыкновенных дифференциальных уравнений // Математические метода в химии / Тез.докл.-Грозный. 1985.С.31.

2. Павлов Б.В..Родионова O.E. Метод локальной линеаризации при численном решении систем обыкновенных дифференциальных уравнений // Ж.вычисл.матем.и матем.физ.1987.т.27..№ 5.С.688-699.

3. Павлов Б.В..Родионова O.E. Численное интегрирование жестких локально-неустойчивых систем ОДУ методом локальной линеаризации // Новые подходы к решению дифференциальных уравнений / Тез.докл. M.I987.С.83-84.

4. Павлов Б.В..Родионова O.E. О численном решении жестких локально-неустойчивых систем обыкновенных дифференциальных уравнений // В сб.:Численное решение обыкновенных дифференциальных уравнений. М.:Ин-т прикл. матем. им. М.В. Келдыша АН СССР. 1988.С.84-95.

О.Е.Родионова

Численное интегрирование жестких локально-неустойчивых систем обыкновенных дифференциальных уравнений

Подписано в печать 24.04.91. Заказ 29 Тирал 100 экз. Формат бумаги 60x90 I/I6 Уч.-изд. л. 0,72. Усл.-печ. л. 1,3. Бесплатно

Отпечатано на ротапринтах в ВЦ АН СССР I17333, Москва, ул. Вавилова, 40