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

•' АКАДЕМИЯ НАУК РОССИИ . - СИБИРСКОЕ ОТДЕЛЕНИЕ ВЫЧИСЛИТЕЛЬНЫЙ ПЕНТР

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

УДК 519.622

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

ИССЛЕДОВАНИЕ (ш,к)-МЕТОДОВ РЕШЕНИЯ НЕЯВНЫХ

СИСТЕМ ОДУ

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

АВТОРЕФЕРАТ диссертации на соискание ученой степени кандидата физико - математических наук

НОВОСИБИРСК - 1994

Работа выполнена в Вычислительном центре СО РАН

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

наук Б.А. Новиков

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

наук Ю.М. Лаевский

кандидат физико - математических наук Р. Р. Ахмеров

Ведущая организация - Институт химической кинетики и горения СО РАН

•1-у 1/1с:

Защита состоится 'д!" февраля 1994 г. в 'и. час. на заседании специализированного совета К 002.10.01 по присуждению ученой степени кандидата физико -математических наук в Вычислительном центре СО РАН по адресу: 630090, Новосибирск 90, просп. Академика Лаврентьева, 6 С диссертацией можно ознакомиться в читальном зале библиотеки ВП СО РАН (пр. Академика Лаврентьева, 6).

Автореферат разослан "¿Лг января 1994 г.

Ученый секретарь специализированного совета д. ф. - м. н. J / Ю. И. Кузнецов.

/

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

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

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

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

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

тельно производной.

Апробация работы. Материалы диссертации докладывались и обсуждались на:

- конференции молодых учёных ВЦ СО АН СССР (1987 г.);

- рабочем совещании "Методы решения жёстких систем и их приложения", Красноярск, 1988 г.;

- семинарах "Статистическое моделирование в физике", лаборатории стохастических задач математической физики ВЦ СО РАН.

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

Структура и объём работы. Диссертация состоит из введения, трёх глав, заключения и списка литературы из 62 наименований. Общий объём работы - 153 страницы.

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

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

Р(х,х)=о, х(го) = хо, а « \ (1)

где х, Р- гладкие вещественные N - мерные вектор-функции, X -независимая переменная. Кратко излагается содержание всех разделов.

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

В первом параграфе дана общая форма записи (т,к.)-схем решения задачи Коши (1) и некоторых её конкретных видов записи для неявных систем индекса не выше двух.

Пусть заданы целые числа т, к € Ъ, к ^ т. Рассмотрим множества Мт = {1,..., т),

М1 , =М Л М° .,

т,к т,к т,к

^={га3-11 3>1, т^М® к, т ^ 1}, 1< 1 ^ ш. Тогда (т.к)-схемы численного решения задачи (1) имеют вид:

ш т

приращения к*, к^ вычисляются при и по формулам:

+5А 2 а-Л* + Ь(е-1 )А 2 Тц*?. И .2)

к^ = ^(аЯНк^Ку^Е (3±1к^)-е 2 а 1^-(1-е)Ь 2 7±Х), (1-3)

при 1€ М1 по формулам:

т,к

1^=1/ (аЬ) (1е- 2 (1 + 2 7-,-,^)), (1.5)

х XII ¡ы^3 3 ¿е-т^1-3 3

где р1, р1;), а13> 713 , а - параметры схемы, А1 , А2 - матрицы ,

аппроксимирующие производные Fr¡x=дí (хп,уп)/дх, 3? =дР(хп1уп)/(9у, к - число вычислений функции Р, т - число обратных ходов в методе Гаусса. Так как шик полностью определяют затраты на шаге интегрирования, а набор чисел т1,...»т^ только распределяют их внутри шага, то схемы типа (1.1)-(1.5) называются (ш.к)-схемами. Параметр е равен 0 или 1. Схемы при е = 1 предпочтительней для вычислений, поскольку требуют при реализации меньше умножений вектора на матрицу, при 8 =0 использовались при исследовании условий на параметры, обеспечивающих заданный порядок точности метода . Число параметров схемы равно тк + 1.

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

класса общих линейных безитерационных схем решения (1), включающего ЖЖ-методы и метода Розенброка. Доказано утверждение об эквивалентности рассматриваемого класса схем классу (т.к)-схем, удовлетворяющих условию <^=0. Получены формулы для перехода от (т,к)-схем с параметром е =0 к случаю е =1.

Во втором параграфе для неявной системы индекса 1

х = Г(х,у) (1.6)

О = е(х,у) ,

где х(г0)=х0?Нн, у(г0)=у0с Н.м, г ^ - невырождена,с помощью помеченных деревьев и соответствующих элементарных дифференциалов задачи (1.6), записано точное решение в точке г ^ через решение и его производные в точке 1; .

В третьем параграфе получено разложение в ряд Тейлора матрицы производных

А,

А

11

"21

12

22

УхГуе>

х?=х(1;п+111),у$= у(гп+:ш),1<0 в окрестности точки решения (*П»УП) относительно параметра схемы шага й. С помощью этого разложения получена рекуррентная по порядку элементарных дифференциалов процедура вычисления производных численного решения (т.к)-схем задачи Коши (1.6), использующих "замороженные" матрицы производных. Приведены условия сходимости численного решения задачи (1.6), даваемого (т.к)-схемами, к точному, условия согласованности схем до четвёртого порядка точности.

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

ц,=1, цг=1-а/2,"а=1±У2/2. • (1.7)

При ш=3 получена Ь-устойчивая (3,2)-схема третьего порядка точности с двумя свободными параметрами а, рг), определяемая мно-

жествами М° 2={1,2}, М^ 2={3}.

Для решения жёсткой задачи Коши

х=Г(х), х(го)=хо, Ха $ % $ (1.8)

при ш=4, М°2={1,3), М^ г={2,4) построена А-устойчивая, а при т=5, 2={Г,3), М^ г=£2,'4,5} Ь-устойчивая (га,2)-схемы четвертого порядка точности, удовлетворяющие дополнительному условию согласованности, полученному в третьем параграфе, которое может быть полезным при очень высокой жесткости задачи (1.8).

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

Доказано утверждение, что при любом ш нельзя построить (ш.З)-метод четвёртого порядка точности.

При ш=5, М°3=£1,3,4), Мд 3={2,5> получена Ь-устойчивая (5,3)-схема четвёртого порядка точности решения задачи (1) для случая невырожденной матрицы Р . При ш=4, д={1,3,4,5}, М^- А =(2) получена Ь-устойчивая (5,4)-схема с семью свободными параметрами четвертого порядка точности решения задачи (1.6).

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

Доказано утверждение, что при любом т нельзя построить (т,2)-метод третьего порядка точности решения задачи (1.6) с параметрами не зависящими от числа 1 - шагов, на которых матрица Вп при приращениях схемы (1.1)-(1.5) не меняется.

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

При т=5, 3={1,3,4), М^ ={2,5} получена Ь-устойчивая (5,3)-схема четвертого порядка точности решения задачи (1), использующая "замороженные" матрицы производных, для случая невырожденной матрицы Р . Показано, что при т=5, к=4 не существует (т,Ю-схемы четвйртого порядка точности решения задачи (1.6) с замораживанием матриц.

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

В первом параграфе для неявной системы индекса 2 х = Г(х,У)

(2.1)

О = g(x) ,

где x{t0)=x0€RN, y(t0)=yQe Нм, t <stk, g^fy - невырождена , с помощью помеченных деревьев и соответствующих элементарных дифференциалов задачи (2.1) , записано точное решение в точке t через решение и его производные в точке t .

Во втором параграфе получено разложение в ряд Тейлора матрицы производных

х =х(1;п+Ш,у£= у("1п+!Ш),КО в окрестности точки решения (*П»УП) относительно параметра схемы шага 1г. С помощью этого разложения получена рекурентная по порядку элементарных дифференциалов процедура вычисления производных численного решения (т.к)-схем задачи Коши (2.1), использующих "замороженные" матрицы производных. Приведены условия сходимости численного решения задачи (2.1), даваемого (т.к)-схемами, к точному, условия согласованности схем до третьего порядка точности.

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

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

Третья глава посвящена построению алгоритмов переменного шага решения задачи Коши (1) и (1.6) индекса 1, в том числе с замораживанием матриц производных.

В первом параграфе рассмотрены два способа получения неравенств для контроля точности вычислений (т.к)-методами решения

и

21

12

g^VV

wv

о

о

задачи Коши (1.6) индекса 1.

Пусть решение УП(Р_,. хп,р. Уп,р задачи (1.6) в

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

^_1ip(t> = cg(t),

где с - вычисляемый параметр, ?(t) - некоторая функция, не зависящая от шага интегрирования. Тогда контроль точности вычислений осуществляется проверкой неравенств

l|enip|| S £> (3.2)

-1

ilßn

lahg^.

« E (3.3)

в случае задачи (1.6), неравенства (3.2) и Ц0~1Рп+1!| < е, в случае задачи (1 )■ индекса 1, где еп_р=1/с(хп_р-хп>р_1),

- значения соответствующих функций в точке хп+1 , уп+1 р, в - задаваемая точность вычислений, Ц-1| - некоторая норма.

Во втором параграфе на основе исследованных в первой главе (т-Ю-схем и предложенного контроля точности построены три алгоритма решения задачи Коши для систем индекса 1:

А) На основе (2,1)-схемы (1.7) алгоритм БОРИ второго порядка точности с параметрами р^а =1-у2/2, рг=1-а. Для контроля точности вычислений используются неравенства

^г"^ >И * 7£ • 1 Чг* 2' (3•4}

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

Б) На основе Ъ-устойчивой (5,2)-схемы (1.1)-(1.5) алгоритм БОРТ третьего порядка точности, использующий точные матрицы производных. Параметры схемы выбирались из условий внутренней Ь-ус-тойчивости и минимизации уклонения

Н= тах ¡ЩгЬехрСг)'! (3.5)

оператора перехода схемы й(г) от экспоненты на вещественной полуоси (г < 0) и имеют следующие значения:.

а = 3(1~У2/2)/4 = .2196699141101,

р, = .21889220743866,

рг = .73291371161961, р3 = .39052429175127, рл = .32231663874474, р5 = -.12024833790342, р31 = а, рзг = 3/4 - а = .5303300858899, азг = 1.643990879855,

а., = 5.0524042692361, а„ = 21.119178189642.

АС. Ьс.

Для контроля точности вычислений используются неравенства

Д^-Р^к^! < 6. 1^2, (з.б)

^ е- (3.7)

где параметры Р1 вложенной схемы выбирались так, чтобы константа с в (3.1) равнялась 1 и имеют следующие значения:.

Р, = .19907547727838, Р2 = .85281143732257, Р., = .245105262331, Р„ = .50726356613138,

3 4

Р5 = -.1597762358698.

В) На основе ь-устойчивой (5,3)-схемы (1.1)-(1.5) алгоритм БОРТХ третьего порядка точности, использующий "замораживание" матриц производных. Параметры выбирались из условий (3.1) на локальную ошибку и (3.5) на вид оператора перехода схемы и имеют следующие значения:.

а = 1/4 ,

р1 = 1, р£ = 4.4246100341335,

Р3 = -25> Ра = Р5 = 2/3' Р31 = Рзг = 1 '

(341 = 1 .54719608240 , ¡3Д2 = -.9693921648,

Раз = -1875'

азг = -^149542272,

адг = 3.971960824, адз = -.875 = -11.50513731840, = .75.

Ос Ьо

Для контроля точности вычислений используются неравенства (3.6), (3.7), где параметры Р вложенной схемы имеют следующие значения:.

Р1 = 2.2687470078894, Р2 = .90367715360981,

Р„ = .39681610116676, Р, = .19575480155567,

о 4

Р5 = .470911865111.

Попытка замораживания матриц производных осуществляется после каздого шага. Следующие причины могут приводить к их размораживанию [Новиков Е.А. Одношаговые Сезитерационные методы решения жёстких систем. Автореферат дисс. докт. физ.-мат. наук.-Новосибирск. -ВЦ СО, РАН, 1991. ]:

1) в случае повторного вычисления решения при невыполнении точности вычислений;

2) если прогнозируемый шаг больше произведения некоторой постоянной Нр, Ер > 1, на величину предыдущего шага, то есть если шаг достаточно быстро возрастает;

3) если количество шагов с замороженной матрицей достигает своего максимального значения;

4) если р;1Рп+,|| > е.

5) если ¡¡V(1)|| > ||У(2)|| , где V(Л) вычисляется по формуле

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

Л

грирования. Однако, если используются матрицы, вычичисленные несколько шагов назад, то оценка ошибки V(2) = D~'v(1) может оказаться хуже за счёт применения "испорченной" матрицы . Поэтому, если ||V(1 )|| > ||V(2)|| и ||V(2)|| < s, то полученное приближение к решению принимается , но следующий шаг выполняется с перевычисленными матрицами производных.

Во третьем параграфе приведены расчбты

а) восьми задач из радиоэлектроники;

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

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

г) на модельной задаче химической кинетики, приведенной в качестве тестовой в описании пакета LSODI.

Все рассмотренные задачи можно записать в виде

A(x)x=i(x,t), x(tQ)=x0, tQ «с t с tk (3.8)

и относятся к неявным системам (1) индекса 1. Три задачи а) соответствуют нелинейным транзисторным схемам, в их записи (3.8) матрица А постоянная вырожденная, а функция f разрывна. Остальные примеры линейные с плохо обусловленной матрицей А и соответствуют RC-схемам. В задачах б) и г) матрица А постоянная вырожденная, в задаче в) матрица А - зависит от решения и плохо обусловлена.

Сравнение построенных алгоритмов проводилось с алгоритмом ADISE2, основанном на схеме типа Розенброка второго порядка точности и программой LSODI переменного порядка и шага, основанной на неявных многошаговых схемах. Из результатов расчётов следует, что фактическая точность решения, полученного построенными алгоритмами, как правило , не ниже задаваемой. Все построенные методы выигрывают по числу вычислений функций и их производных у алгоритма ADISE2 и при точности s=I0-2 у алгоритма LSODI, при точности 10алгоритм SOFT выигрывает у LSODI по числу' вычислений функций, a S0PTZ по числу шагов и обращений матрицы Dn.

Таким образом, построенные алгоритмы являются конкурентно-способными с известными методами.

В заключении сформулированы основные результаты, которые

сводятся к следующему.

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

2. Для систем индекса не выше двух получены условия согласованности (m,к)-схем, использующих как аналитические так и "замороженные" матрицы производных.

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

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

Автор выражает глубокую благодарность научному руководителю д.ф.-м.н. Новикову Е.А. за постоянное внимание к данной работе.

По результатам, представленным в диссертации опубликованы работы:

1. Левыкин A.M., Новиков Е.А. -О (т.к)-методе второго порядка точности для решения неявных систем обыкновенных дифференциальных уравнений.- Новосибирск, 1987, 16 с . ( Препринт / АН СССР. Сиб. отд-ние. ВЦ; 768)

2. Левыкин A.M., Новиков Е.А. Одношаговый метод третьего порядка точности для решения неявных систем обыкновенных дифференциальных уравнений.// Моделирование в механике. 1989 Т.3(20). №4. с. 90-101.

3. Левыкин А.И., Новиков Е.А. Исследование (т.к)-методов третьего порядка точности для решения неявных систем обыкновенных дифференциальных уравнений.- Новосибирск, 1990 , 24 с. ( Препринт / АН СССР. Сиб. отд-ние. ВЦ; 882)