автореферат диссертации по информатике, вычислительной технике и управлению, 05.13.18, диссертация на тему:Метод вариационных сплайнов для сингулярных дифференциальных уравнений

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

Мартыненко Юлия Вячеславовна

МЕТОД ВАРИАЦИОННЫХ СПЛАЙНОВ ДЛЯ СИНГУЛЯРНЫХ ДИФФЕРЕНЦИАЛЬНЫХ УРАВНЕНИЙ

Специальность 05.13.18 - математическое моделирование, численные методы

и комплексы программ

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

15

Ульяновск - 2008

003456915

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

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

профессор

Горбунов Владимир Константинович

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

профессор

Леонтьев Виктор Леонтьевич

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

Коноплева Ирина Викторовна

Ведущая организация: ГОУ ВПО Московский авиационный

институт (государственный технический университет, МАИ)

Защита состоится "26" декабря 2008 г. в на заседании

диссертационного совета Д 212.278.02 при Ульяновском государственном университете по адресу: ул. Набережная р. Свияги, 106, корп. 1, ауд. 703.

С диссертацией можно ознакомиться в библиотеке Ульяновского государственного университета, авторефератом - на сайте вуза http://www.uni.ulsu.ru

Отзывы по данной работе просим направлять по адресу: 432000, г. Ульяновск, ул. Л. Толстого, д.42, УлГУ, Управление научных исследований.

Автореферат разослан " ^ 3" ноября 2008 г.

Ученый секретарь диссертационного совета // г М. А. Волков

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

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

Проблема численного решения сингулярных НОДУ стала активно разрабатываться в нашей стране и за рубежом в последние три десятилетия. Большой вклад в проблемы численного решения сингулярных НОДУ внесли работы H.H. Лузина, Ф.Р. Гантмахера, Ю.Е. Бояринцева и его школы (М.В. Булатов, В.Ф. Чистяков)1, В.И. Шалашилина, Е.Б. Кузнецова2, C.W. Gear, S. Campbell, L.R. Petzold, E. Hairer, G. Wanner' и др.

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

1 Бояринцев Ю.Е.. Данилов В. А.. Логинов A.A.. Чистяков В.Ф. Численные методы решения сингулярных систем. - Новосибирск: Наука, 1989.

" Шалашилин В.И.. Кузнецов Е.Б. Метод продолжения решения по параметру и наилучшая параметризация. -М: Эдиториал УРСС. 1999

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

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

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

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

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

4 Горбунов В.К. Метод нормальной сплайн-коллокации // ЖВМ и МФ. 1989. Т. 29. №2. С.212-224.

5 Горбунов В.К. Метод параметризации задач оптимального управления // ЖВМ и МФ. 1979. Т. 19. №2. С.292-303.

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

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

Предметом исследования являются численные методы решения задач математического моделирования.

Объектом исследования являются сингулярные обыкновенные дифференциальные уравнения.

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

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

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

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

Основные положения, выносимые на защиту:

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

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

3. Сравнительный анализ методов вариационных сплайнов и параметризации.

4. Комплекс программ, реализующий метод вариационных сплайнов.

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

Апробация работы. Результаты диссертации докладывались и обсуждались на следующих конференциях: XIII Байкальской международной школе-семинаре "Методы оптимизации и их приложения" (Иркутск, Байкал, 2005); VI международной конференции "Математическое моделирование физических, технических, экономических, социальных систем и процессов" (Ульяновск, 2005); VII международной конференции "Дифференциальные уравнения и их приложения" (Саранск, 2006); IX Р>сероссийском симпозиуме по прикладной и промышленной математике (Кисловодск, 2008); The Fourth International Conference 'Inverse problems: Modelling and Simulations" (Fethiye-Turkey, 2008).

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

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

Структура и объем диссертации. Диссертация объемом 108 страниц состоит из введения, четырех глав, приложения и списка литературы из 68 наименований. Работа содержит 11 рисунков и 21 таблицу.

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

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

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

Вторая глава посвящена изложению теории и вычислительной схемы метода вариационных сплайнов (далее ВС). В первом параграфе рассматривается начальная задача для системы НОДУ:

/¡\х,х,0 = 0, (0<1<Т, (1)

где х е К", с условием

*('„) = (2)

Функция Р(х,х,1) считается непрерывной по всем переменным и дважды непрерывно дифференцируемой по х и х. Матрица Якоби системы (1)

дР(х,хА)

у. (3)

ох

может произвольно вырождаться в области, содержащей решение (1), (2).

Метод ВС основан на переходе к задаче минимизации интегральной невязки системы

/

J{x{■))=\\F{x{t)Лt),t)tdt (4)

при условии (2) на некотором параметрическом классе гладких (непрерывно дифференцируемых) или кусочно-гладких функций х(0.

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

В общем случае, ввиду возможного вырождения матрицы Якоби (3), задача (2), (4) должна рассматриваться как некорректная, т.е. она должна быть регуляризирована некоторым образом.

Если решение (1), (2) не существует, то элементы множества М, = Лг§ттЩх(-)):(2)} (5)

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

||х(Г)-г(/)|2Л, (6)

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

поиска г-нормального решения может быть приближенно решена с помощью минимизации Тихоновского сглаживающего функционала

т

Л;« (*(■)) = |{|И*(0, *(0, 'Г + «1И0 - г(0||2 }Л, а > 0 (7)

'о

с условием (2). Из теоремы 3, доказанной в работе6, следует, что при а 1О множество решений (2), (7) аппроксимирует множество х-нормальиых решений и таким образом обеспечивается сильная сходимость метода ВС, т.е. сходимость по решению х.а —> х..

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

['о. Л

/0 </, нГ (8)

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

(*) = х,* (*, V *) = £ у кш{1~1ку, I=г,п, 1к < I < , (9)

а=О

—к ( к к к 1 1 где V = {1>(0,у115...,у , г = 1,п),с условиями непрерывности

= (10) и, возможно, непрерывности первых производных

хк+\(к+],ум) = хк(Гк+]-0У). (11)

Таким образом, сплайн (7) определяется набором параметров -векторов V* и скаляров 1к. Если условия (11) присутствуют, этот сплайн имеет дефект р -1, в противном случае - дефект р.

Определение 1. Назовем вариационным сплайном степени р функцию (9), доставляющую минимум функционалу (4) при условиях (2), (10) и, возможно, (11).

Горбунов В.К. Регуляризация нелинейных некорректных задач с параметризованными данными // Нелинейный анализ и нелинейные дифференциальные уравнения / Под. ред. В.А. Треногина и А.Ф. Филиппова. - М.: Физматлит. 2003. С.418-447.

В третьем параграфе приводится обоснование разрешимости задачи минимизации (4) при условиях (2), (9), (10), и возможно, (11). Оно основано на условиях сходимости метода параметризации. Поставленная задача о ВС легко сводится к задаче оптимального управления заменой

х, (/) = «,(/)»' = и.,и;

II 1,2 0 2)

Если исходная задача (1), (2) разрешима, то минимум функционала невязки на классе гладких, а следовательно, и кусочно-гладких функций равен нулю и достигается на этом решении. При этом непрерывная производная х(/) решения (1), (2) может быть сколь угодно точно равномерно приближена сплайном (8), непрерывным или кусочно-непрерывным. Это означает, что выполняется основное условие сходимости МП7, которое обеспечивает сходимость по функционалу при стремлении диаметра сетки разбиения (8) к нулю. Таким образом, задача (4) при условиях (2), (9), (10), и возможно, (11), разрешима.

В четвертом параграфе приведена основная вычислительная схема метода ВС. Интегральный функционал (4) аддитивен относительно отрезка интегрирования [/0>^]> поэтому его минимизацию можно проводить последовательно. Для каждого промежутка ) вводится функция частичной пошаговой невязки

уравнения

= (13)

Чтобы обеспечить на всем отрезке ,7"] требуемый уровень

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

Горбунов В.К., Лутошкин И.В. Развитие и опыт применения метода параметризации в вырожденных задачах динамической оптимизации // Известия РАН: теория и системы управления. 2004. № 5. С.67-84.

<Pk(ykAj<s. (14)

Функция (13) при сделанных предположениях относительно функции F(x,x,() имеет вторые непрерывные производные по всем переменным, и ее минимизация может выполняться градиентным методом или методом Ньютона. Задача упрощается, если применять

эти методы относительно переменных групп vk при фиксированном значении tk+v Если при этом точность (14) не достигается, то следует

уменьшить промежуток за счет приближения правого конца

tk+] к левому tk.

В пятом параграфе рассматривается вопрос вычисления частичной пошаговой невязки и ее производных. Начальное условие (2) и условие непрерывности (10) позволяют последовательно определить компоненты {v*0} и перейти к задаче безусловной минимизации (13) по параметрам

Если на сплайн (9) наложены условия непрерывности производной (11), то можно последовательно определить компоненты

( fe л

{v,,} и сократить размерность задачи безусловной минимизации (13) до р(п-1).

Пусть аналитическое представление минимизируемой на каждом шаге k = 0,N-\ функции частичной пошаговой невязки (13) отсутствует. В этом случае вычисление ее значений, а также значений ее производных может выполняться на основе аппроксимации интеграла, зависящего от неизвестных параметров vk, приближенным значением, полученным с помощью квадратурной формулы Симпсона. Тогда минимизируемая функция (13) заменяется интегральной суммой

Т — / ~п'

¥к <Ук,) s -т-^ IX Х7^ )У (г,- Ш, (15)

6m 7?о tf где коэффициенты квадратурной формулы

к. =

1, г = 0,2т;

2, г = 2,4,...,2»; -2; 4, г = 1,3,...,2т-1.

Далее обозначим V = V *. Доказана следующая теорема.

Теорема 1. Если функция непрерывно

дифференцируема по х и х, то первые производные функции (15) по переменным V имеют вид: дЩ (Т-(0)-2"'

ар, Й/^ Зс,

(16)

и если ^(х,х,0 дважды непрерывно дифференцируема по х и х, то вторые производные имеют вид:

Зт

ы о

£

Э/7, ¿Зх, ¿У7, дх, + -

дху д\\р дх, ду,р

дР, дх,. дГ, дхг

дхг д\'„ дхг

+ Р

дх. дх„ д2/7 дх, дх„ Э2/1,

(17)

Эу.о Эу„„ дх.дх.

ду,я дуга дх,дхг

дх, дхг д2Р, дх, дх„ +-------— + -

г

о г,

дуга ду,р дх,дхг ду,р дут Зс,дкг

Здесь производные по параметрам дхг / дуга и дхг / в точке (определяются для сплайна степени р по формулам:

дх. , ^ дх.. . —

= a{t-tk) , г = 1,ЛГ; а = \,р.

Таким образом, полученные формулы (16), (17) позволяют применить для минимизации (15) методы эффективный метод Ньютона.

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

А(0х(0 + В(0х(0 =/({), с условием

х(^) = х°,(й<(<Т.

Тогда функция(15) принимает вид

2т п ( п 1

(18) (19)

^ = ХЫт.^+ХВ/Т.^Т.МГ,) I . (20)

ЬШ 1-0 1=1 ^ ,=1 /=] )

Доказана теорема 2.

Теорема 2. Формулы (16) и (17) для линейного уравнения (18) принимают вид:

дц/к 74

т п \ п

3 т

^[А./т^УТ^Щ/Т,^^)]-/^,)

1=0 ,=} {

Зс, дх

(21)

га у

2т

3?п

/=0 ;=|

(22)

<7У ¿/V

^ ("Г/ ^ К 1-/-У

Из формулы (21) следует, что система дц/к/ду = 0 линейна

относительно параметров V, следовательно, оптимальное значение параметров может быть найдено как решение системы линейных уравнений

Су = с1. (23)

Матрица С - {с } этой системы совпадает с симметричной матрицей вторых производных и находится по формуле (22), а правая часть системы с/ = (с/,} вычисляется по формуле (24) в случае непрерывного сплайна

/у » \ 2т п ( п ^

/=0 /=1 ^ у=1

\

га га J

и по формуле (25) в случае непрерывно дифференцируемого:

2т п ( п ^

(Т / п _п

(25)

где г'(/\а) = рп{г -1) + а.

В седьмом параграфе исследуются вопросы сходимости и

устойчивости для уравнения, разрешенного относительно производной

х = Дх(/), /), х(Г0) = , < / < Т. (26)

Пусть решение этой задачи х (?) существует и единственно и для некоторого заданного уровня невязки е удалось построить вариационный сплайн хЛ (?), определяемый коэффициентами V, на сетке /0,...,/Л>.

Обозначим хд=х^(/^). Тогда можно интерпретировать конечный набор значений х = {х0,х,,...,хЛ.} как результат решения

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

Пусть /гтах - максимальная длина промежутков разбиения (8). Введем величину

Д, =тах{|к/)-х*(0|Ь'*

Определение 2. Назовем величину ЛЛ' =тах погрешностью

к

сплайна, а А^ = max х^. - х (tk) погрешностью сплайна на сетке

к " "

Очевидно, что Д* <ДЛ. Имеет место следующая оценка этих погрешностей.

Теорема 3. Пусть для уравнения (26) построен сплайн x'v (/) с уровнем невязки s и функция f(x,t) удовлетворяет условию Липшица по переменной х, т.е. j|/(x,f) - - Ц\х _>|> с

постоянной L такой, что L < 1//г,тах. Тогда для погрешности сплайна Дч верна оценка

д. idtB- ,р-А11ю. (27,

Следствие 3. Пусть начальная задача для линейной системы (18) обладает следующими свойствами:

1) Решение задачи существует и единственно:-

2) Матрица A(t) обратима на [?0,Т] и |д^(/)|<ог;

3) \B„{t)\<p.

Тогда, если ар <1 //гтач, то для погрешности сплайна Xх (t) с уровнем невязки G выполняется соотношение

у N-ap{N-\)hmm ГГ, Г8)

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

Теорема 4. Метод параметризации для кусочно-гладкого сплайна первой степени на равномерной сетке с шагом h является I-устойчивым, А-устойчивым и не является L-устойчивым.

Третья глава посвящена практическим и теоретическим аспектам, связанным с предлагаемым методом.

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

Если при этом не достигается условие (14), то длина промежутка уменьшается в два раза и снова решается задача минимизации (13). При этом не учитывается возможное изменение "жесткости" искомого решения, и после уменьшения шага на некотором участке невозможно его увеличить. Это может привести к излишнему измельчению шага, например при решении жесткой задачи с пограничным слоем на левом конце отрезка решения. Поэтому добавляется контроль за величиной частичной пошаговой невязки на предыдущих шагах. Если на нескольких предыдущих промежутках это значение убывало, то делается попытка увеличить длину текущего промежутка в два раза. Это позволяет восстанавливать шаг в процессе решения. Численные эксперименты показывают, что такой контроль значительно уменьшает объем вычислений и время решения.

Во втором параграфе оценивается трудоемкость метода ВС по двум характеристикам: необходимому объему вычислений функции F(д:,x,/) и ее производных и числу арифметических операций, необходимых для решения задачи минимизации частичной пошаговой невязки. Для сравнения аналогичные характеристики были оценены для неявного 5-стадийного метода Рунге-Кутгы.

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

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

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

Оба метода имеют свои вычислительные характеристики и требуют углубленного исследования как теоретически, так и экспериментально. Была проведена сравнительная оценка их трудоемкости.

Предположим, что оба метода используют метод Ньютона для минимизации функционала невязки. Ключевым становится вопрос о вычислении первых и вторых производных по параметрам. В МП для этого предварительно нужно решать задачи Коши для нахождения сопряженных систем и матричных импульсов. Пусть для этого используется 8-стадийный метод Рунге-Кутты. После этого необходимо вычислять интегралы, используя некоторую квадратурную формулу. В методе ВС производные находятся непосредственным дифференцированием функционала невязки и вычисляются на основе приближенного вычисления интегралов по квадратурной формуле. Если оба метода используют одну и ту же квадратурную формулу, то в МП требуется в 5 " раз больше вычислений функции .Р(х,х,/) = 0 и ее производных. Численные эксперименты показывают, что при одинаковом порядке невязки и погрешности метод ВС требует гораздо меньше времени.

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

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

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

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

Для уравнений в нормальной форме получена оценка точности метода и проведено исследование устойчивости.

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

Все полученные теоретические и алгоритмические результаты реализованы в разработанном комплексе программ.

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

Публикации в изданиях, входящих в список ВАК

1. Горбунов В.К., Мартыненко Ю.В. Метод вариационных сплайнов для сингулярных дифференциальных уравнений // Обозрение прикладной и промышленной математики. 2008. Т.15. Вып. 3. С. 461-462.

Публикации в прочих изданиях

2. Горбунов В.К., Лутошкин И.В., Мартыненко Ю.В. Метод параметризации для решения неявных дифференциальных уравнений // Обратные и некорректные задачи прикладной математики: труды XIII Байкальской межд. школы-семинара <Методы оптимизации и их приложения> (Северобайкальск, 2-8 июля, 2005 г.). Иркутск: Изд-во ИСЭМ СО РАН. 2005. Том 3.С. 100-105.

3. Горбунов В.К., Лутошкин И.В., Мартыненко Ю.В. Метод параметризации для сингулярных обыкновенных дифференциальных уравнений // Труды Средневолжского математического общества. Саранск: СВМО. 2006. Т.8. Вып.1. С.36-50.

4. Горбунов В.К., Мартыненко Ю.В. Метод вариационных сплайнов для неявных дифференциальных уравнений // Вестник Самарского государственного технич. ун-та. Сер. Математическая. 2007. №2(6). С. 16-28.

5. Мартыненко Ю.В. Вариационно-параметрические сплайны в вырожденных системах дифференциальных уравнений // Труды шестой международной конференции "Математическое моделирование физических, технических, экономических, социальных систем и процессов". Ульяновск: УлГУ, 2005. С.90-91.

6. Мартыненко Ю.В. Метод вариационных сплайнов для дифференциальных уравнений в нормальной форме // Труды Средневолжского математического общества. Саранск: СВМО. 2007. Т.9. Вып.2.С.110-120.

7. Мартыненко Ю.В. Программный комплекс для решения дифференциально-алгебраических уравнений [Текст] / Ю.В. Мартыненко // Инновации в науке и образовании (Телеграф отраслевого фонда алгоритмов и программ). - №9(44). - С. 50.

8. Gorbunov V.K., Lutoshkin I.V., Martynenko Y.V. The parametrization method for numerical solution of singular differential equations // International Conference "Tikhonov and Contemporary Mathematics", Moscow, June 19-25, 2006. Sect. №4. C. 64-65.

9. Gorbunov V.K., Martynenko Y.V. The variational splain method for singular differential equations // The Fourth International Conference "Inverse problems: Modelling and Simulations", Fethiye-Turkey, May 26-30, 2008. P. 53.

Подписано в печать 18.11 08. Формат 60x84/16. Гарнитура Times New Roman Уел печ. л. 1,0 Тираж 100 экз Заказ № 118

Отпечатано с оригинал-макета в Издательском центре Ульяновского государственного университета 432000, г. Ульяновск, ул Л Толстого, 42

Введение

1 Сингулярные дифференциальные уравнения и методы их решения

1.1 Уравнения в нормальной форме и методы их решения

1.1.1 Разностные методы.

1.1.2 Проекционные методы.

1.1.3 Метод сплайн-коллокации.

1.1.4 Жесткие задачи

1.2 Неявные и дифференциально-алгебраические уравнения

1.2.1 Особенности НОДУ.

1.2.2 Основные методы решения НОДУ.

1.3 Метод параметризации

2 Метод вариационных сплайнов

2.1 Общая схема метода.

2.2 Вариационный сплайн.

2.3 Сходимость.

2.4 Основная вычислительная схема.

2.5 Вычисление пошаговой частичной невязки и ее производных

2.6 Производные по параметрам для линейного уравнения

2.7 Метод вариационных сплайнов для уравнения, разрешенного относительно производной.

2.7.1 Сходимость.

2.7.2 Устойчивость.

3 Алгоритмические вопросы метода ВС

3.1 Алгоритмы построения вариационного сплайна.

3.2 Оценка сложности метода ВС.

3.3 Сравнение метода ВС и метода параметризации.

4 Вычислительный эксперимент

4.1 Жесткие задачи

4.2 Линейные дифференциально-алгебраические уравнения

4.3 Нелинейные ДАУ.

А Описание комплекса программ

Диссертация посвящена разработке численных методов решения обыкновенных дифференциальных уравнений в неявной форме (НОДУ)

F(x,x,t) = 0, xeRn, (1) несводимых к нормальной форме Коши x = f{x,t). (2)

Эта проблема стала активно разрабатываться как в СССР, так и за рубежом в последние три десятилетия. Прикладные задачи, моделируемые в классе НОДУ, возникают в различных областях науки, техники, экономики, в частности в электродинамике и механике (кинематические уравнения) [45, 49, 63].

Неразрешенность уравнения (1) относительно производной х может вызываться как сложностью соответствующих преобразований, так и вырожденностью матрицы Якоби dF (x,x,t) , д± ' {6) В последнем случае система (1) называется сингулярной. Постановка начальной задачи для такого уравнения требует учета структуры вырождения.

В литературе для систем вида (1) используются различные термины: дифференциально-алгебраические уравнения (ДАУ), неструктурированные ДАУ, алгебро-дифференциальные уравнения, сингулярные ДУ и др. Далее будем называть уравнение (1) НОДУ в общем случае и сингулярным в случае вырождения матрицы Якоби (3).

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

A(t)x(t) + B{t)x(t) = f(t), (4) и его матрица Якоби (3) равна А (£). Она не зависит от решения х (£).

Наиболее важным классом НОДУ, который часто встречается в приложениях, являются так называемые структурированные ДАУ, или просто ДАУ. В этом случае динамическая система представляется двумя группами функций х (t) £ Rn пи (£) £ Rm, связанными системой дифференциальных и конечных уравнений х = f (x,u,t), (5) g{x,u,t)= 0, (6) где / : Rn+m+l —» Rn: д : др+m+i j^jn -ракие системы сочетают в себе особенности нелинейных алгебраических или трансцендентных уравнений с особенностями обыкновенных дифференциальных уравнений. Здесь начальные или граничные условия должны быть согласованы с (6).

Неструктурированное ДАУ (1) может быть преобразовано к форме (5), (6) введением функции и (t), обозначающей производную решения х (*) = и (£).

При этом исходная система (1) принимает вид (5), (6) с конечным уравнением F(u,x,t) = 0. Такое преобразование позволяет выбирать различные варианты численных схем метода параметризации, описанные далее.

Трудности в построении теории и применении численных методов к решению начальной или граничной задачи для (1) или (5), (6) появляются, когда матрица Якоби (3) или матрица Якоби конечной подсистемы (б) dg{x,u,t)

Зи [ } вырождается на решении {.т(£)} или {х (t) ,u(i)} соответственно.

Известно, что сингулярные дифференциальные уравнения имеют ряд особенностей по сравнению с регулярными дифференциальными уравнениями [5, 47]. В частности, их решения (если они существуют) не обязательно непрерывно зависят от начальных данных, а множество решений может иметь бесконечную размерность. Примеры таких систем приведены в п. 1.2.1. Соответственно, в общем сингулярном случае, теория и численные методы решения должны разрабатываться в рамках теории некорректных задач [41, 54]. Такие методы обычно основываются на. использовании дополнительной эвристической информации о решении и на переходе от исходной проблемы к некоторой вариационной задаче.

Во многих случаях система. (1) может быть преобразована к нормальной форме Коши (2) с помощью конечного числа дифференцирований по независимой переменной и алгебраических преобразований. Минимальное число таких дифференцирований называется индексом дифференцирования [45].

В приложениях встречаются системы индекса 5 и выше [51]. Более того, известны задачи, которые нельзя привести к нормальной форме(т.е. они не имеют конечного индекса дифференцирования), например уравнение Шрёдингера с вырожденным кусочно-непрерывным потенциалом [64]. Простейший пример - уравнение tx (L) + Ъ (t) я (t) = / (t) , рассматриваемое на. отрезке, содержащем точку t = 0. Легко видеть, что оно имеет переменное вырождение и не может быть приведено к нормальной форме посредством дифференцирований по независимой переменной. Другие примеры приведены в п. 1.2.1.

Современное состояние проблемы численного решения неявных дифференциальных уравнений (1) достаточно полно отражено в монографиях Э.Хайрера и Г.Ваннера [45], P. Kunkel и V. Mehrmann [63], в докторских диссертациях М.В.Булатова [8], Г.Ю.Куликова [30], В.Ф.Чистякова [47], защищенных в 2002 году, а также в работах [11, б, 49]. Увеличившееся в последнее время количество работ [2, 7, 9], посвященных данной проблеме, говорит об ее актуальности. Кроме того, остается актуальным совершенствование методов численною решения "жестких" задач для уравнений в нормальной форме (2).

Первые результаты, полученные при изучении НОДУ, относятся к линейным уравнениям (4) с постоянными коэффициентами. В работе Н.Н. Лузина [31] 1940 г. доказан критерий совместности линейной системы с постоянными матричиымп коэффициентами Aj и приведено общее правило нахождения решения этой системы Большое влияние на развитие теории и методов решения линейных сингулярных систем (4) оказало применение Ф.Р. Гантмахером теории пучков матриц [13].

Систематическое изучение (1) и построение численных методов их решения началось в 1970-х гг. Большой вклад в исследование линейных систем (4) внесли работы Ю.Е. Боярннцева и его учеников [4, 5, 6, 7]. Ими рассматривались взаимосвязи кронекеровой структуры пучка матриц XA(t)-\-B(t) со структурой общего решения н свойствами численных методов. В работах М.В. Булатова и В.Ф. Чистякова [8, 9] предлагается понижать индекс исходной системы (4) с помощью "левого регуляризирующего оператора", определенного через полуобратные матрицы. Таким образом, для систем (4) развиты общие методы определения и понижения индекса, основанные на сложной алгебраической технике, итерационном построении иолуобрат-ных матриц и ранговых критериях, которые трудно проверять в условиях приближенных данных. Для сингулярных задач с переменным вырождением главной части A(t), несводимых к нормальной форме (2). эти методы, как правило, неприменимы

Также следует отметить работы зарубежных ученых [50, 51, 52, 53], посвященные различным аспектам теории НОДУ (1) и методам их решения. В частности, в [52] определяется понятие индекса дифференцирования, а в работе [51] представлен специальный вариант метода наименьших квадратов для неструктурированной системы высокого индекса, основанный на разрешении «массива производных» относительно всех производных решения. Этот метод требует выполнения определенных ранговых условий, проверка которых является достаточно сложной задачей.

Для сингулярных уравнений с невысоким индексом дифференцирования (до 3) существуют эффективные специализированные методы [45, 67]. Некоторые из них реализованы в известной системе тестов [68]. Отметим также метод продолжения решения по параметру [49], в котором требуется дифференцирование по независимой переменной.

Однако дифференцирование исходного уравнения в общем случае приводит к уравнению с более широким множеством решений, и такое расширение может привести к постороннему решению (см. п. 1.2.1).

Ввиду возможной сложности нахождения индекса дифференцирования и различных способов его уменьшения, очевидна ограниченность перечисленных выше подходов. Операция дифференцирования функций, определяющих реальную модель и заданных приближенно, является некорректной [3, 41]. Она может существенно увеличить ошибку задания функций. Некорректна также задача вычисления полуобратных матриц [6]. Соответственно, нормализация сингулярных НОДУ высокого индекса является плохо обусловленной вычислительной задачей, н .мера ее обусловленности растет с ростом индекса системы, а порядок аппроксимации решения нормализованной задачи снижается. Кроме того, указанные методы нельзя применять к уравнениям с неременным вырождением матрицы Якоби (3), как у приведенного выше элементарного примера. Такие уравнения не сводятся к нормальной форме, а значит, не имеют конечного индекса.

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

Для линейных уравнений (4) существует метод нормальной сплайн-коллокации, разработанный В.К. Горбуновым [15] и далее развитый его учениками в работах [22, 23, 24, 40], который может использоваться для решения задач с переменным вырождением матрицы A (t).

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

Метод параметризации был развит и применен в диссертации И.В. Лу-тошкина [32] для различных вырожденных задач оптимального управления [17, 56]. При этом были охвачены задачи классического вариационного исчисления, порождаемые НОДУ произвольной структуры. Минимизируемым функционалом здесь является интегральная квадратичная невязка части или всех уравнений решаемой системы. Производные искомых функций считаются управлениями. Такая схема позволила использовать эффективную технику сопряженных систем для вычисления первых и вторых производных функционала параметризованной вариационной задачи. Этот подход позволяет не вычислять индекс дифференцирования и решать системы любых индексов, а также системы с переменным вырождением, обычно не имеющие конечного индекса. В случае возможной неединственности решений здесь естественно доопределять вариационную задачу с использованием дополнительных условий на искомое решение [54].

В данной диссертации разрабатывается метод вариационных сплайнов, который представляет собой спецификацию метода параметризации для НОДУ. Специфика вариационной задачи, порождаемой такими уравнениями, позволяет не вводить управляющие переменные, а представлять искомое решение непосредственно в виде сплайна с подвижными узлами. Параметры сплайна определяются из условия минимизации невязки решаемой системы. Такой сплайн назван вариационным в отличие от интерполяционных и коллокационных сплайнов, параметры которых определяются решением соответствующих систем уравнений [26. 36]. Производные функционала невязки по параметрам находятся непосредственным дифференцированием без использования сопряженных систем, как в исходном варианте метода параметризации.

Диссертация состоит из 4 глав. В первой главе рассматривается современное состояние проблем, связанных с численным решением НОДУ. Во второй главе излагается метод вариационных сплайнов. Приводится общая вычислительная схема метода и формулы вычисления производных функционала невязки по параметрам. Для уравнения в нормальной форме (2) исследуются проблемы оценки точности и устойчивости, что является актуальным для "жестких" задач. Третья глава посвящена практическим и теоретическим аспектам, связанным с предлагаемым методом. Четвертая глава содержит результаты решения задач.

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

На. защиту выносятся:

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

• Оценка точности новой схемы МП и исследование ее устойчивости для уравнения в нормальной форме.

• Сравнительный анализ новой и исходной схем МП.

• Комплекс программ, реализующий новый вариант МП.

Выносимые на защиту результаты опубликованы в работах [18, 19, 20, 21, 33, 34, 35, 58, 59, 60], в том числе в журнале "Обозрение прикладной и промышленной математики" [21], входящем в список изданий, рекомендованных ВАК.

Исследование выполнено в рамках проекта, поддержанного Российским Фондом Фундаментальных Исследований, проект № 07-01-90000.

Результаты работы докладывались и обсуждались на семинарах Сред-неволжского математического общества, ряде российских и зарубежных конференций, в том числе:

1. XIII Байкальская международная школа-семинар "Методы оптимизации и их приложения" (Иркутск, Байкал, 2-8 шоля 2005 г.);

2. VI международная конференция "Математическое моделирование физических, технических, экономических, социальных систем и процессов" (Ульяновск, 19-21 октября, 2005 г.);

3. IX Всероссийский симпозиум по прикладной и промышленной математике (Кисловодск, 1-8 мая, 2008 г.);

4. The Fourth International Conference "Inverse problems: Modelling and Simulations" (Fethiye-Turkey, May 26-30, 2008 г.).

1. Алберг Дж., Нильсон Э., Уолш Дж. Теория сплайнов и её приложения. - М.: Мир, 1973.

2. Балакина. Е.А., Кузнецов Е.Б. Решение систем дифференциально- алгебраических уравнений высоких индексов. // Ж. вычисл. матем. и матем. физ. 2000. Т.40. С. 199-206.

3. Бахвалов Н.С., Жидков Н.К., Кобельков Г.М. Численные методы. -М.: Мир, 1988.

4. Бояринцев Ю.Е. Регулярные и сингулярные системы линейных обыкновенных дифференциальных уравнений. Новосибирск: Наука, 1980.

5. Бояринцев Ю.Е. Методы решения вырожденные системы обыкновенных дифференциальных уравнений. Новосибирск: Наука, 1988.

6. Бояринцев Ю.Е., Данилов В.А., Логинов А.А., Чистяков В.Ф. Численные методы решения сингулярных систем. Новосибирск: Наука, 1989.

7. Бояринцев Ю.Е. Линейные и нелинейные алгебро-дифференциальные системы. Новосибирск: Наука, 2000.

8. Булатов М.В. Методы решения дифференциально-алгебраических и вырожденных интегральных систем. Дисс. . д.ф.-м.н. Иркутск, 2002.

9. Булатов М.В., Чистяков В.Ф. Об одном численном методе решения дифференциально-алгебраических уравнений j j Ж. вычисл. матем. и матем. физ. 2002. Т. 42. № 4. С.459-4Т0.

10. Верлань А.Ф., Сизиков B.C. Интегральные уравнения: методы, алгоритмы, программы. Киев: Наукова думка, 1986.

11. Вырожденные системы обыкновенных дифференциальных уравнений / Ред. Бояринцев Ю.Е. Новосибирск: Наука, 1982.

12. Гавурин М.К. Лекции по методам вычислений. М: Наука, 1971.

13. Гантмахер Ф.Р. Теория матриц. М: Наука, 1988.

14. Горбунов В.К. Метод параметризации задач оптимального управления // Ж. вычисл. матем. и матем. физ. 1979. Т. 19. № 2. С.292-303.

15. Горбунов В.К. Метод нормальной сплайн-коллокации // Ж. вычисл. матем. и матем. физ. 1989. Т. 29. № 2. С.212-224.

16. Горбунов В.К. Регуляризация нелинейных некорректных задач с параметризованными данными // Нелинейный анализ и нелинейные дифференциальные уравнения / Под ред. В.А.Треногина и А.Ф.Филиппова. М.: Физматлит, 2003. С.418-447.

17. Горбунов В.К., Лутошкин И.В. Развитие и опыт применения метода параметризации в вырожденных задачах динамической оптимизации // Известия РАН: теория и системы управления. 2004. № 5 С. 67-84.

18. Горбунов В.К., Лутошкин И.В, Мартыненко Ю.В. Метод параметризации для сингулярных обыкновенных дифференциальных уравнений // Труды Средневолжского математического общества. Саранск: СВМО. 2006. Т.8. Вып.1. С.36-50.

19. Горбунов В.К., Мартыненко Ю.В. Метод вариационных сплайнов для неявных дифференциальных уравнений // Вестник Самарского государственного технич. ун-та. Сер. Математическая. 2007. №2(6). С.16-28.

20. Горбунов В.К., Мартыненко Ю.В. Метод вариационных сплайнов для сингулярных дифференциальных уравнений // Обозрение прикладной и промышленной математики. 2008. Т. 15. Вып. 3. С. 461-462.

21. Горбунов В.К., Петрищев В.В. Метод нормальных сплайнов в вырожденных системах дифференциальных уравнений // Учёные записки УлГУ. Сер. "Фундаментальные проблемы математики и механики" Вып. 3. Ульяновск, 1997. С. 125-132.

22. Горбунов В.К. Петрищев В.В. Развитие метода нормальной сплайн-коллокации для линейных дифференциальных уравнений. // Ж. вы-числ. матем. и матем. физ. 2003. Т. 43. № 8. С.1161-1170.

23. Горбунов В.К., Свиридов В.Ю. Метод нормальных сплайнов для численного обращения преобразования Лапласа в вещественной форме // Труды СВМО, Саранск, Т. 10. № 1. 2008. С. 46-54.

24. Дулан Э., Миллер Дж., Шилдерс У. Равномерные численные методы решения задач с пограничным слоем. М.: Мир. 1983.

25. Завьялов Ю.С., Квасов Б.И., Мирошниченко В.Л. Методы сплайн-функций. М.: Наука, 1980.

26. Крылов В.И , Бобков В.И., Монастырный П.И. Вычислительные методы. Т.2. М.: Наука, 1976.

27. Кузнецов Е.Б., Шалашилин В.И. Решения дифференциально-алгебраических уравнений с выбором наилучшего аргумента // Ж. вычисл. матем. и матсм. фпз. 1997. Т.37. № 6. С.711-722.

28. Куликов Г.Ю. Численное решение задачи Коши для систем дифференциально-алгебраических уравнений с помощью неявных методов Рунге-Кутты с нетривиальным предиктором // Ж. вычисл. матем. и матем. физ. 1998. Т. 38. № 1. С.68-84.

29. Куликов Г.Ю. Численные методы с контролем глобальной ошибки для решения дифференциальных и дифференциально-алгебраических уравнений индекса 1. Дисс. . д.ф.-м.н. Ульяновск, 2002.

30. Лузин Н.Н. К изучению матричной теории дифференциальных уравнений // Автоматика и телемеханика. 1940. №5. С.4-66.

31. Лутошкин И.В. Метод параметризации и его использование в вырожденных задачах. Дисс. . к.ф.-м.н. Ульяновск, 2000.

32. Мартыненко Ю.В. Метод вариационных сплайнов для дифференциальных уравнений в нормальной форме // Труды Средневолжского математического общества. Саранск: СВМО. 2007. Т.9. Вып.2. С.110-120.

33. Мартыненко Ю.В. Программный комплекс для решения дифференциально-алгебраических уравнений Текст] / Ю.В. Мартыненко // Инновации в науке н образовании (Телеграф отраслевого фонда алгоритмов и программ). N° 9(44). - С. 50.36