Неявные численные методы решения функционально-дифференциальных уравнений и их компьютерное моделирование тема автореферата и диссертации по математике, 01.01.02 ВАК РФ

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

V7TK .-»1Q Г,9

ск ¿а

Квон О Бок

НЕЯВНЫЕ ЧИСЛЕННЫЕ МЕТОДЫ РЕШЕНИЯ ФУНКЦИОНАЛЬНО-ДИФФЕРЕНЦИАЛЬНЫХ УРАВНЕНИЙ И PIX КОМПЬЮТЕРНОЕ МОДЕЛИРОВАНИЕ

01.01.02 - дифференциальные уравнения

Автореферат

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

Екатеринбург 2000

Работа выполнена на кафедре теоретической механики Уральского государственного университета им. А.М.Горького.

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

наук, с.н.с. A.B. Ким;

— кандидат физико-математических наук, доцент В.Г. Пименов.

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

наук, профессор А.Ф. Шориков;

— кандидат физико-математических наук, доцент А.Ю. Вдовин.

Ведущая организация — Московский государственный институт электроники и математики.

Защита диссертации состоится 2000 г. в^-^ ч.

00 м. на заседании диссертационного Совета К 063.078.03 по присуждению ученой степени кандидата физико-математических наук в Уральском государственном университете им. А.М.Горького по адресу: 620083, г. Екатеринбург, проспект Ленина 51, комн. 248.

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

Автореферат разослан

года.

Ученый секретарь диссертационного Совета кандидат физ.-мат. наук, доцент 'Г'-4 " Пименов В.Г.

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

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

Основополагающий вклад в развитие теории функционально-дифференциальных уравнений внесли В. Вольтерра, Р. Беллман, H.H. Красовский, А.Д. Мышкис, Л.Э. Эльсгольц, Н.В. Азбелев, Р.Д. Драйвер, В.Б. Колмановский, К. Кордунеану, К.Л. Кук, А.Б. Куржанский, С.Б. Норкин, В.Р. Носов, Ю.С. Осипов, Д.К. Хейл, С.Н. Шиманов и многие другие математики. Полученные при этом результаты находят значительные приложения в моделировании процессов автоматического регулирования, управления и устойчивости движений, механики, различных технологических процессов, биологии, медицины, химии, экономики и в других отраслях знаний.

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

Среди многообразия исследований в области численных методов решения ФДУ отметим следующие направления.

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

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

2. Идея непрерывных методов [4, 5] состоит в том, что численная модель задает не только значение в узлах, но и во всех промежуточных точках. Эти методы обладают большой степенью общности по отношению к различным типам ФДУ. Однако, большинство эффективных методов решения ОДУ (в том числе и применяемых в различных пакетах прикладных программ) не являются непрерывными. Особенно это касается неявных методов, среди которых непрерывные методы неизвестны.

3. В работах [6-8] был предложен новый подход к конструированию численных методов решения ФДУ. Эти методы основаны на идеях разделения конечномерной и бесконечномерной фазовых составляющих, построению по конечномерной составляющей полных аналогов известных для ОДУ дискретных алгоритмов, проведению по бесконечномерной составляющей (предыстории) интерполяции с заданными свойствами и использовании специальной техники вычисления производных функционала правой части ФДУ вдоль решения (эта техника названа г-гладким анализом) [9]. В качестве интерполяционных процедур была предложена интерполяция вырожденными сплайнами и экстраполяция продолжением.

В рамках этого подхода в данной диссертации численные конструкции распространяются на широко применяемый в ОДУ класс методов - неявные одношаговые методы типа Рунге-Кутты.

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

Цель работы. Работа посвящена разработке и изучению свойств нРК-методов для ФДУ, в том числе, получении условий порядка сходимости мсТОДОВ 15 ЗаБИСИмОСТй ОТ ПСрЛДКа ЛОКсьЛЬНСЙ аППрОКСЯМа-ции, от порядка интерполяции дискретной предыстории модели и от порядка ее экстраполяции; разработке алгоритмов с автоматическим

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

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

Научная новизна. Все существенные результаты работы являются новыми. Приведем основные из них.

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

2. Получена теорема о порядке сходимости методов.

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

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

5. Для неявных методов разработаны алгоритмы автоматического выбора шага.

6. Указанные алгоритмы протестированы на ряде жестких и нежестких систем, имеющих явное решение.

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

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

Апробация результатов работы. Результаты диссертации докладывали еь по, международном семинаре 1РАС "Негладкие и раз рывные задачи управления и оптимизации"(Челябинск, июнь 1998); на международном конгрессе по автоматике (Анкоридж, май 1998);

на международной конференции по динамическим системам и их приложениям (Атланта, май 1999); на семинарах кафедры теоретической механики Уральского государственного университета, группы функционально-дифференциальных уравнений ИММ УРО РАН, на семинарах в Сеульском национальном университете и Кёнгмун университете.

Публикации. Результаты диссертации опубликованы в 9 работах, список которых приведен в конце автореферата.

Структура и объем диссертации. Работа состоит из введения и трех глав, разбитых на 12 разделов. В работе приведено 34 рисунка. Общий объем диссертации 97 страниц, библиография содержит 128 наименований. Диссертация написана на английском языке.

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

В главе 1 изучаются НРК-методы с постоянным шагом. В разделе 1 приводится постановка задачи Коши для ФДУ вида

¿(4) = /(*,®(*),*4(0) (!)

с начальными условиями

х(г0) = го, хи(а) = у°(в), -т<й<0. (2)

Здесь * е [¿о><о + 0] ~ время; х 6 В.1 - фазовый вектор; £((-) = {х(< + в), -т < я < 0} € Я[~т, 0) - функция-предыстория, <3[—т, 0) -пространство кусочно-непрерывных на [—т, 0) функций со значениями в В}, г > 0 - величина, характеризующая запаздывание.

Если отображение /(¿,х,у) : [¿о, ¿о + #] х Я' х (?[—т, 0) -»• Я1 непрерывно по сдвигу [9] и липшицево по второму и третьему аргументам, то решение задачи Коши существует и единственно [9].

В разделе 2 на примере простейшего метода (неявного метода типа Эйлера с кусочно-постоянной интерполяцией) приводятся основные конструкции в построении численных методов решения систем вида (1). Пусть = ¿о + пА, п = 0,1,..., N - узлы равномерного разбие-

л

ния отрезка + здесь А = —. Для простоты мы предположим, т

что — = ш целое.

Дискретной моделью назовем вектора ип, п = О,1,..., /V которые являются приближенными значениями точного решения х(^) в узлах. Предысторией дискретной модели к моменту назовем множество {"¡)п = — т < г < п}. При п = 0,1, ...,ЛГ — 1 рассмотрим

кусочно-постоянную интерполяцию предыстории модели

гл _ / * 6 &>*•+].)> п - т < г < п

и неявную пошаговую формулу

Щ — Хо,

ип+1 = ип + ип+и Щ„+, (•)) (з)

Численный метод сходится, если |[и„ — ж(£п)|| —> 0 при Д -> О для всех п = 1,..., ТУ; и имеет порядок сходимости р, если найдется постоянная С такая, что ||щ„ — г(£п)|| < СДр для всех п = 1,..., N.

С помощью оценок для локальной погрешности (невязки) и погрешности интерполяции доказывается

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

В разделе 3 изучаются более точные способы интерполяции.

В подразделе 3.1 приводится определение оператора интерполяции и его порядка. Выписан оператор кусочно-линейной интерполяции. Для произвольного натурального р проведено построение оператора интерполяции вырожденными сплайнами р-той степени - функцией, кусочно заданной многочленами р-той степени. Доказывается теорема о том, что при определённых условиях оператор интерполяции вырожденными сплайнами р-той степени имеет порядок р+ 1.

Особенность конструкций рассматриваемых в дальнейшем НРК-ме-тодов состоит в том, что в момент Ьп приходится не только производить интерполяцию дискретной предыстории модели, но и проводить экстраполяцию модели вправо за точку tn. В подразделе 3.2 приводится г^зроттопоимс г)ттог\йт1пг\я э^ТрапОЧЯЦ"И И ^О РЛпаДКЕ!. Д.псг произвольного натурального р проведено построение оператора продолжением вырожденного сплайна р-той степени. Доказывается теорема о том,

что при определённых условиях оператор экстраполяции продолжением вырожденного сплайна р-той степени имеет порядок р+ 1.

В подразделе 3.3 понятия оператора интерполяции и оператора экстраполяции объединены в единое. Для а > 0 оператором интерполяции-экстраполяции 1Е предыстории дискретной модели назовем отображение

1Е : {и4}„ и{-) £ (¡[Ъ - т, и + аД] •

Оператор 1Е имеет порядок р на точном решении х(£), если найдутся постоянные Сг, такие что

||®(0 - "(*)|| < С-1 шах ||и, - Х<|| + С2(ДР)

п—т<г<п

для всех п — 0,1,..., N — 1, и Ь € [¿„ — т, + аД].

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

Раздел 4 является основным в первой главе.

В подразделе 4.1 приводится общая схема построения неявных методов типа Рунге-Куттьг. Пусть оператор интерполяции-экстраполяции 1Е задан. НРК-методом назовем алгоритм

По =

к

ип+1 = и„. + Д £ а Ы(щГ(-)), п = 1,..., N - 1, »=1

к

= /(<п + а«Д.Ип +Д Е Ь<ЛК,(-)),«»я+в,д(-))' ^

где щп(-) = 1Е{{щ}п).

Обозначим для дальнейшего

*(*», «(О) = Е ^К(-)). п =1,.. •, N - 1.

¿=1

Как и в случае обыкновенных дифференциальных уравнений, каждый НРК-метод характеризуется набором коэффициентов, которые можно задавать в виде матрицы Бутчера. На основе её свойств проводится классификация методов, разделение их на явные и неявные.

На каждом шаге НРК-метода основная пошаговая формула (4) метода представляет собой систему нелинейных уравнений. Доказана

Теорема 4-1- Пусть Д < ттт-1 где Ь — тах Ь - кон-

кой

станта Липшица функционала х, у) по второму аргументу. Тогда существует единственное решение Я = Ь.2,..., Ик)' системы (4); функционалы /12,..., Ль и Ф липшицевы по «(•).

Основным результатом подраздела 4.2 является теорема о порядке сходимости НРК-метода.

Пусть хп = х{Ьп), п = 0,1,..., ^-значения точного решения х{£) задачи (1) - (2) узлах ¿п. Невязкой НРК-метода назовем функционалы

где

к к = Е^^п ^ = /(*п + аг-Д,з;п + Д £ д(-)).

г=1 ;=1

Будем говорить, что НРК-метод имеет порядок аппроксимации (невязки) р , если найдется константа С такая, что

Цд„(х(-))Ц < СА", п = 0,...,ЛГ-1

В ОДУ порядок сходимости НРК-метода определяется лишь порядком аппроксимации, однако в ФДУ порядок сходимости НРК-метода зависит также от порядка оператора интерполяции и экстраполяции, а именно справедлива

Теорема Пусть НРК-метод имеет порядок аппроксимации р\ > О и оператор 1Е имеет порядок р2 > 0. Тогда метод сходится и имеет порядок сходимости р = тт{р1,р2}-

В подразделе 4.3 изучаются условия, обеспечивающие порядок аппроксимации метода. Показывается, что особенность сконструированных НРК-методов для ФДУ состоит в том, что они являются полными аналогами соответствующих методов для ОДУ в смысле совпадения порядка аппроксимации методов для гладких задач ОДУ и ФДУ с одинаковыми таблицами Бутчера. В качестве примера выведено правило средней точки и приведены другие примеры, при этом используются конструкции г—гладкого анализа.

Указываются способы решения нелинейной пошаговой формулы: схема предиктор-корректор, метод Ньютона-Рафсона, особенности линейных систем.

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

В подразделе 4.5 представлены некоторые классы функционалов, используемых при описании правых частей ФДУ. Конкретный вид функционалов позволяет эффективно проверять условия существования и единственности решений ФДУ, исследовать свойства соответствующих операторов интерполяции и проводить процедуру распараллеливания.

В главе 2 изучаются неявные методы типа Рунге-Кутты с переменным шагом.

В разделе 5 проводится описание НРК-методы при условии, что шаг дискретизации непостоянен.

Зададим на [¿о, ¿о + неравномерную временную сетку Ьп+\ — ^ + Дп, п = 0,1,..., N — 1, где Дп > О -шаг. Обозначим Атах — тах Д„, Дтт = тш Д„. Будем предполагать, что Атах < КАт{п, К - фиксированное число.

Обозначим приближение решения з;(£г1) = хп в точке £п через ип £ К1. Расширенной предысторией дискретной модели назовем множество {«;}+ = {щ, и 6 [¿п, - 2т,*„]}.

Пусть > 0. Будем говорить, что на множестве расширенных дискретных лредысторий модели {щ}* задан оператор интерполяции и экстраполяции 1Е, если

1Е{{щ)п) = и(-) ь ч[ъп - т,ъп -Ь аДта1).

Для оператора 1Е, вводится понятие порядка, согласованности (и{и) = щ) и липшицевости. Если такой оператор задан, то аналогично

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

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

В разделе 7 исследуется порядок оператора интерполяции расширенной предыстории дискретной модели вырожденными сплайнами р—степени. При неравномерной временной сетке в отрезок [¿п — т, ¿тг], длиной равный величине запаздывания, укладывается не обязательно целое число шагов. Поэтому, чтобы не потерять точности интерполяции, нужно привлекать для построения самого левого из многочленов в конструкции сплайна узлы из предыдущего отрезка „ — 2т, ?„ — т]. Доказывается, что при такой модификации оператор интерполяции расширенной предыстории дискретной модели вырожденными сплайнами р—степени имеет порядок р+1 при условии достаточной гладкости точного решения. Оператор экстраполяции расширенной предыстории дискретной модели продолжением вырожденного сплайна р—степени также имеет порядок р + 1, откуда следует существование оператора интерполяции-экстраполяции 1Е с требуемыми в теореме о порядке сходимости предыдущего раздела свойствами.

Раздел 8 посвящен описанию процедур автоматического выбора шага. Излагается правило Рунге практической оценки погрешности применительно к рассматриваемым методам. На основе оценки локальной погрешности в зависимости от величины шага приводятся алгоритмы автоматического уменьшения или увеличения шага в зависимости от заданной погрешности. Рассматривается два варианта: оценка погреш-

игч^ттд лпипт п ттоу*л<т тла'зиитлдтт шагами оценка гтот^решног,гт'ч

двух вложенных методов (методов с почти одинаковыми матрицами Бутчера). Такие процедуры являются неотъемлемой частью современ-

пых численных методов решения ОДУ и ФДУ, применяемых в пакетах прикладных программ.

В разделе 9 рассматривается возможность применения численных методов для решения задач позиционного управления системами с запаздыванием

x(t) = f(t,x(t),xt{-),u{t),v{t)),

где u(t) € Rq вектор управляемых параметров, v(t) 6 Rr - помеха.

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

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

Жесткость решений системы ФДУ понимается в следующем (конечномерном) смысле.

Пусть x(t) - решение задачи (1). Составим матрицу Якоби функционала f вдоль решения

= df(t,x(t),xt(.))^ Ох

Обозначим через Л¿(f), г = 1,..., I - собственные числа матрицы /1(f). Задачу (1) назовем жесткой, если

1) Re A¿(f) < 0 для всех f Е [fo,fo + 0], г = 1, ..., I;

о\-------

¿j J 1HU1U

SUP {max I Re A;(f)|/ min | Re A,-(t)|}

велико.

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

В разделе 11 описаны явный и неявный алгоритмы с автоматическим выбором шага и проведено их сравнение на тестах жестких и нежестких систем.

В разделе 12 рассматривается применение построенных методов к просчету одной модели из медицины. Система

2®= S - inT(t) + rT(t) [l - Щ+^t-rMZAt-r^ _

—k\V{t)T{t)

M - k\V(t)T(t) - ß\Zi(t - n) - hZ^t - n)

k2Z,(t - n) - fi2Z2{t - r2) ^ = Nß2Z2(t - r2) - k:V(t)T(t) - А*V(t)

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

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

x{t) = A{t)x(t) + Л1(г)х(4 - г) + B(t)u.

Выбором управления и — u(t, х, xt(-)) требуется минимизировать квадратичный функционал

J = хт{\3)Nüx{d) + f (uTNi(t)u)dt.

о

Итполтип ГТ_Ol uta гчгтттда*я ттълтр гппаяпрннр "меОТ ВИД

u(t,x,xt(-)) = -N^^B^t)

и

P{t)x+ I Q(t, s)x(t + s)ds

где матрицы P(t), Q{t,s) являются решениями системы обобщенных уравнений Риккати. Таким образом, даже если исходная система имела только постоянное сосредоточенное запаздывание, при подстановке оптимального управления в систему возникает распределенное запаздывание и при реализации алгоритма необходимы численные методы решения систем ФДУ.

Цитированная литература

[1] Холл Д., Уатт Д. Современные численные методы решения обыкновенных дифференциальных уравнений. М. Мир. 1979. 312 с.

[2] Bellen A. Constrained mesh methods for functional differential equations // International Series of Numerical Mathematics, Verlag, Basel. 1985. P. 52 - 70.

[3] Baker C.T.H., Paul C.A.H. and Wille D.R. Issues in the numerical solution of evolutionary delay differential equations // Advances in Comput. Math. 1995. V. 3. P. 171 - 196.

[4] Tavernini L. One-step methods for the numerical solution of Volterra functional differential equations // SIAM J. Numer. Anal. 1971. V. 8. P. 786 - 795.

[5] Хайрер Э., Hepcemm С., Ваннер Г. Решение обыкновенных дифференциальных уравнений. Нежесткие задачи. М. Мир. 1990. 512 с.

[6] Kim А. V., Pimenov V. G. Numerical methods for time-delay systems on the basis of г-smooth analysis // Proc. of the 15th World Congr. on Scient. Computation, Modelling and Applied Mathematics. Berlin, August

1997. V. 1: Computational Mathematics. P. 193 -- 196.

[7] Ким А.В., Пименов В.Г. О применении i-гладкого анализа к разработке численных методов решения функционально-дифференциальных уравнений // Тр. Ин-та математики и механики УрО РАН.

1998. Т. 5. С. 104 - 126.

[8] Kim А. V., Pimenov V.G. Multistep numerical methods for functional differential equations //Mathematics and Computers in Simulation. 1998. V.45. P. 377 - 384.

[9] Kim A.B. Functional differential equations. Application of ¿-smooth calculus. Kluwer Academic Publishers. The Netherlands. 1998. 236 p.

[10] Андреева E.A., Колмаповский В.В., Шайхет JI.E. Управление системами с последействием. М. Наука. 1992. 336 с.

РАБОТЫ АВТОРА ПО ТЕМЕ ДИССЕРТАЦИИ

[11] Квип О.Б., Пименов Б.Г. Неявные методы типа Рунге-Кутты // Изв. УрГУ 1095 ЛГо 10. С. 69 70.

[12] Kim.A.V., Kwon О.В., Pimenov V.G. Lyapunov-Krasovski func-tionals in stability and control problems // Proceeding of the World Automation Congress (May 9-14,1998). Anchorage, Alaska, USA, p. 105-110.

[13] Kwon O.B., Kim A. V., Pimenov V.G. Numerical modeling of control time-delay systems // Proceeding of the IFAC Workshop "Nons-mooth and discontinuous problems of control and optimization"(June 1720, 1998). Chelyabinsk, Russia, p. 137-139.

[14] Song S.H, Kwon O.B., Lim C.M. System analysis and design (Editor O.B. Kwon). Daerim Publishing Co. 1998. 468 p.

[15] Kwon O.B., Kim K.B., Park H.K., Yoo G.S. Computer modeling (Editor O.B. Kwon). Seoul. Hakmun Publishing Co. 1998. 426 p. (In Korean)

[16] Kwon O.B. Numerical simulation of delay differential equations // Third International Conference on Dynamic Systems and Applications. Atlanta, USA. May 1999. P. 43.

[17] Kwon O.B., Kim A. V., Pimenov V.G. Numerical modeling of control time-delay systems //' Nonsmooth and discontinuous problems of control and optimization / Edited by V.D.Batukhtin, F.M.Kirilova and Ukhobotov. Elsevier Publishers. London. 1999. P. 155 - 160.

[18] Квон О.Б. Численные алгоритмы с переменным шагом для решения функционально-дифференциальных уравнений // УрГУ, Екатеринбург. Рук. деп. в ВИНИТИ, 24.03.2000, №792-В00, 18 с.

[19] Квоп О.Б. Тестирование явных и неявных методов типа Рунге-Кутты для систем с запаздываниями //УрГУ, Екатеринбург. Рук. деп. в ВИНИТИ, 24.03.2000, №193-В00, 32 с.

л:,/ П./ / / !У П О

V* ' I/'/' 7 / - и У- /

/

Уральский государственный университет

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

Квон О Бок

НЕЯВНЫЕ ЧИСЛЕННЫЕ МЕТОДЫ РЕШЕНИЯ ФУНКЦИОНАЛЬНО-ДИФФЕРЕНЦИАЛЬНЫХ УРАВНЕНИЙ И ИХ КОМПЬЮТЕРНОЕ МОДЕЛИРОВАНИЕ

01.01.02 - дифференциальные уравнения

Диссертация

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

Екатеринбург 2000

URAL STATE UNIVERSITY

On the rights of manuscript

KWON Oh Bok

IMPLICIT NUMERICAL METHODS FOR FUNCTIONAL DIFFERENTIAL EQUATIONS AND THEIR SOFTWARE

REALIZATION.

01.01.02 - Differential equations

THESIS

submitted for degree of doctor of philosophy

Ekaterinburg - 2000

Contents

Introduction..........................................................4

Chapter 1. Implicit one-step methods (with fixed step size)

for functional differential equations 12

1. Statement of the problem........................................12

2. Implicit Euler method with piece-wise constant interpolation 13

3. Interpolation and extrapolation of the model pre-history . . 16

3.1. Interpolating operators..................................16

3.2. Extrapolating operators..................................19

3.3. Interpolating-extrapolating operator....................21

4. Implicit Runge-Kutta-like methods..............................21

4.1. General scheme ..........................................22

4.2. Convergence order........................................25

4.3. Approximation order ....................................27

4.4. Potential for parallel computation......................31

4.5. On the structure of functionals...................32

Chapter 2. Numerical algorithms (with variable step size) for

functional differential equations 34

5. Description of IRK-methods with variable step size control . 34

6. Convergence order................................................37

7. Methods of interpolation and extrapolation of the extended

pre-history of discrete model..................................39

8. Choice of the step size............................................43

9. Numerical modeling of control time-delay system..............45

Chapter 3. Testing and comparing explicit and implicit Runge—Kutta—like methods for functional differential equations 49

10. Realization and comparing explicit and implicit numerical methods with fixed step size..................................50

10.1. Testing of explicit and implicit Euler methods .... 50

10.2. High order methods with fixed step size ..............63

11. Explicit and implicit methods with automatic step size control 66

11.1. Explicit Runge-Kutta-like method of 2(3) order with

automatic step size control............................66

11.2. Implicit trapezoidal method with interpolation and automatic step size control . . . .....................67

12. Model examples (medicine and control theory)................79

12.1. HIV-infection model......................................79

12.2. Modeling of the linear quadratic control problem with

delay....................................................81

Bibliography 86

Introduction

Investigation of the various phenomena allows one to conclude that future behaviour of many processes depends not only the present state, but also it depends on pre-history. Mathematical description of such processes can be realized on the basis of equations with different types of delays -the so called functional differential equations (FDE), equations with delays, hereditary systems.

The delay (in general form - deviation of the argument in an equation) can cause interesting mathematical phenomena: instability, periodic solutions, limit cycles, appearance and disappearance of stiffness. As the consequence, at present there is the intensification of the investigation of qualitative properties of systems with delays in different directions [22, 39, 40, 54, 77, 78, 79, 100, 117, 124].

The obtained results are applied for modeling and analysis of automatic control processes, stability of motion, mechanics, technological processes, biology, medicine, chemistry, economic and so on (one can find examples [1, 5, 39, 43, 50, 54, 77, 84, 99].

Foundations and fundamental results in the theory of functional differential equations were developed by V.Volterra, R.Bellman, N.N.Krasovskii, A.D.Myshkis, S.N.Shimanov, L.E.El'sgolts, N.V.Azbelev, H.T.Banks, T.A.Burton, G.A.Kamenskii, C.Corduneanu, M.C.Delfour, R.D.Driver, J.Hale, L.Hatvani, V.B.Kolmanovskii, A.V.Kryazhimskii, A.B.Kurzhanskii, V.Laksmikantham, A.A.Martynyuk, A.Manitius, K.L.Cook, S.B.Norkin, V.R.Nosov, Yu.S.Osipov, A.L.Skubachevskii, S.B.Razumihin and many other mathematicians.

In the present work we follow the functional approach proposed and elaborated by academician N.N.Krasovskii and Sverdlovsk scientific school of mathematics and mechanics.

One of the obstacles for investigating systems with delays consists in the absence of the general analytical representations for FDE solutions. Besides, systems with delays are essentially infinite dimensional systems, and this fact causes principal difficulties for analysis of such systems.

Hence elaboration of the effective numerical methods for FDE and their software realization is the very important and urgent task. Insufficient development of the software programs for simulating and analysis of FDE

is one of the main obstacles for more wide application of delay models in applied problems.

Let us mention some directions of the theory of numerical methods for FDE.

1. Numerical schemes which are based on the step method [40]. If we consider the system with constant delays

x = f(t,x(t),x{t- r)),

then, substituting known initial function instead of the term x(t — r), we obtain on the interval [¿o, ¿o + t] ODE that can be solved by a standard method. Then we use the obtained on the interval [to, to + r] solution as the initial function to find solutions on the interval [¿o 4- r, to + 2r], and so on.

Advantage of the method consists in a greater simplicity. Unfortunately this method cannot be directly applied to systems with other types of delay (varying delay, distributed delay), besides the method is not suitable for realization of numerical procedures with automatic step size control (the basis of the most algorithms realized in contemporary toolboxes).

2. There are many papers devoted to numerical methods that use the specific form of FDE. Especially there is a big bibliography for systems for systems with lumped time-depending delays

x = f(t,x{t\,x(t-T{t)))

and for integro-differential Volterra equations

t

x = J f(t, s,x(s))ds.

to

For such systems there were elaborated the analogs of almost all known for ODE numerical methods (see, for example, reviews [33, 56, 19, 12]).

3. The idea of the continuous methods [122, 35, 51] consists in constructing approximate models not only at points tn, but also in all points of the delay interval. Such methods can be used for general systems of the form

® = /(*, *«(■)). *«(•) = {x(t + s), -T<8< 0}.

Nevertheless even for ODE continuous methods are not usually used in software packages. Especially it is necessary to note that there are no continuous methods among implicit methods.

4. Many types of FDE can be reduced to integral equations [56] and solved by standard for integral equations methods. However, in contrary to ODE, the integral equations are solved by non-positional methods, i.e. for positional methods at point t one can use information about some part of trajectory, meanwhile for solving integral equation we calculate equation on the whole interval. This fact is the obstacle of application such methods in the positional control problems for FDE.

5. Systems with delays can be approximated by ODE systems of a high dimension. The method is based on the ideas of [80, 113]. The accuracy of the method is not very high, and the is applied for low dimension control problems [14, 15].

6. Functional approach [79, 54], representation of FDE as differential equations in Banach spaces and the semi-group approach are very fruitful for theoretical investigating FDE. However the approaches do not allow one to elaborate effective numerical methods because of the difficulties of calculating derivatives of functionals.

7. In the works [70, 71, 72] a new approach to constructing numerical methods for FDE

x = Sit,x(t),£,(•))> st(.) = {ait + s),-r<s<0}. (0.1)

was elaborated. The approach is based on:

1) separation of finite dimensional and infinite dimensional components in FDE structure;

2) constructing with respect to the finite dimensional component the complete analogs known for ODE numerical algorithms;

3) interpolating (with the given properties) of the discrete model prehistory;

4) application of the special technique of calculating derivatives of functionals (along FDE solutions) (these methods are called i-smooth calculus

The applied interpolational procedure is based on the degenerate splines and extrapolation by continuation.

The aim of the present dissertation is to construct for FDE on the basis of the approach [TO, 71, 72] the analogs of implicit one-step Runge- Kutta methods (which are widely used for ODE, see, for example, [51]).

Implicit numerical Runge-Kutta methods have some advantages compare to other methods - they can be used for solving stiff problem [56, 52]. Besides these methods have high accuracy and potential for parallelization, compare to explicit methods.

One of the problems in practical realization of implicit methods consists in necessity of solving at every step a special nonlinear equation with respect to the current state of the model un. Approach of separation finite and infinite dimensional components in FDE structure and, respectively, in the model structure, allows us to solve this problem in the same manner as in ODE case.

Let us give the overview of the results, presented in the thesis.

In Chapter 1 we investigate implicit Runge-Kutta-like (IRK) methods with constant step size.

In Section 1 the statement of the initial value problem for FDE (0.1) is given. We also describe some conditions which guarantee existence and uniqueness of solutions (the mapping f(t,x,y) should be shift-continuous and should satisfy the Lipschitz condition with respect to the second and the third arguments.

In Section 2 we describe (on the example of the simplest implicit Euler method) basic constructions of the numerical approach for solving (0.1). The notions of the discrete model, the pre-history of the discrete model, the interpolation of the discrete model are introduced. We describe the formula of the Euler method and introduce the notions of the residual (local error), convergence of the method, the residual order and the convergence order of the method.

In the section we investigate the residual order and properties of the

piece-wise constant interpolation. Also we prove the following result: if the solution x(t) of system (0.1) is sufficiently smooth, then the implicit Euler method with piece-wise constant interpolation has the first convergence order.

In Section 3 we investigate more accurate methods of interpolating.

In Subsection 3.1 the notion of interpolating operator and its order are defined. As an example the formula of the operator of piece-wise linear interpolating is derived. For arbitrary natural p we describe the procedure of constructing interpolating operator by degenerate splines of the degree p. It is proved that under some assumptions the interpolating operator by degenerate splines of order p has the order p 4-1.

The specific features of the described IRK-methods is the following: at moment tn it is necessary to realize not only an interpolation of the model pre-history, but also the extrapolation to the right-side of point tn. In subsection 3.2 we introduce the notion of the extrapolating operator and its order. For arbitrary natural p we describe the interpolating operator by continuation of the degenerate spline of order p. We prove that under some assumptions the extrapolating operator by continuation of the degenerate spline of order p has the order p + 1.

In Subsection 3.3 we introduce the unified notion of the interpolating-extrapolating operator IE. On the example of the IE operator, constructed on the basis of degenerate splines and continued extrapolation, we investigate properties of such operators: the order, consistence, the Lipschitz condition.

The Section 4 is the main section in the first chapter.

In Subsection 4.1 we describe general scheme of constructing implicit Runge-Kutta-like methods for FDE. It is supposed that an interpolating-extrapolating operator IE is given. Similar to ODE, each IRK-method for FDE is characterized by a set of coefficients which can be presented in the form of the Butcher tableau. We introduce the classification of explicit and implicit methods for FDE basing on the structure of the Butcher tableau.

For each step of IRK-method, the main advancing formula is the system of nonlinear equations with respect to the next model state. We prove that under some conditions (small size of A) this system has a unique solution.

The principle result of the subsection is the theorem on the convergence order of IRK-methods: IRK-method has the convergence order equals to the minimal value among the approximation order (residual) and the order of interpolating-extrapolating operator IE.

In the next Subsection 4.2 we derive conditions which guarantee the corresponding approximation order of the method.

It is important to note that elaborated IRK-methods for FDE are complete analogs of the corresponding ODE numerical methods, i.e. approximation orders of ODE and FDE numerical methods are the same for similar Butcher tableaus. As an example, we consider the middle point rule, and also describe some other models. In this subsection some elements of ¿-smooth calculus are used.

Besides we describe some methods of solving nonlinear advancing formula: predict-correct scheme, Newton-Raphson method, specific features of linear systems.

In Subsection 4.4 we discuss the possibility of parallelization of computing procedures of IRK-methods for FDE of high dimension. Beside well known for ODE methods of parallelization (special structure of the Butcher tableaus, factorization of the Jacobi matrix for solving the system of nonlinear equations, dimensional parallelization) FDE have some specific potential for organization of parallel computation. These are: parallelization of interpolating and extrapolating procedures of discrete model prehistory; parallelization of computing integrals (for FDE with distributed delays); parallelization of calculating right-hand part of system (0.1) of a complicated structure.

In Subsection 4.5 we describe the structure of functionals which are used in right-hand sides of FDE. Concrete forms of right-hand sides of FDE allow us effectively to verify existence and uniqueness of solutions, to elaborate constructive interpolating procedures and parallelization procedures.

Chapter 2 is devoted to implicit Runge-Kutta-like methods with variable step size.

In Section 5 we describe general scheme of IRK-methods with variable step size. We use non-uniform partition of time interval, ratio of the maximal step to the minimal step should be less than some fixed number. Then

we introduce the discrete model, pre-history of the model and extended model pre-history (i.e. continued, if it is necessary, over the left point of the interval [tn — r, tn}. For the extended pre-history one can introduce an interpolating-extrapolating operator IE, and the corresponding notions of its order, consistence and the Lipschitz property. If such operator is given, then, similar to the case of the constant step size, we introduce the advancing formula of implicit Runge-Kutta-like method for FDE. It is proved that for sufficiently small step size there exists a solution of the system defined next state of the model.

In Section 6 we investigate convergence of IRK-methods with variable step size. The notions of residual and its order, convergence order of the method are introduced.

We prove (similar to the case of the constant step size) that IRK-method with variable step size has the convergence order equals to the minimal value among the approximation order (residual) and the order of interpolating-extrapolating operator IE. In this result we essentially use the constant that limits ratio of the maximal step to the minimal step.

In Section 7 we investigate the order of the interpolating operator (of the extended discrete model) constructed by degenerate splines of degree p. Note, if we use non-uniform grid then the interval [in — r, tn] can contains not natural number of steps of a numerical method. So, in order to preserve the accuracy of the interpolation it is necessary to use (for constructing the left polynomial in the spline) some points of the grid of the previous interval [tn — 2r, tn — r\. If we use such modification, then under sufficient smoothness of the exact solution the interpolating operator of the extended pre-history of the discrete model by degenerate splines of degree p has the order p + 1. Because the extrapolating operator by continuation of the degenerate spline of the degree p has the order p + 1, hence there exists interpolating-extrapolating operator IE with required in the theorem on convergence order (see previous section) properties.

In Section 8 we describe the procedure of automatic step size control. We base on the Runge rule of practical error estimation modified according to the specific features of the proposed approach. On the basis of the local error (depending on the step size) we derive algorithms of automatic decreasing and increasing of steps depending on a given accuracy. Two variants are considered: error estimate of one method with two different

step sizes, and error estimate of two embedded methods (methods with similar Butcher matrices). Such procedures are the bases of modern numerical methods applied in software programs for ODE and FDE.

In Section 9 we investigate the possibility of application ' of numerical methods for solving positional control problems for systems with delays. We gave the problem statement of the positional (conflict) control in the framework of differential game theory, and consider the reduction of this problem to the discrete scheme. The notion of the optimal strategy (positional control) having the accuracy order with respect to the step discretization. It is shown that under some assumptions such strategy (in combination with a suitable numerical method) guarantee necessary (a priori given) accuracy, moreover in the framework of this scheme one can organize the procedure of automatic step size control.

Chapter 3 is devoted to comparison on some models of computing properties of explicit and implicit one-step methods of solving FDE. Main goal of the chapter is to to show that implicit numerical methods are more effective than explicit methods for solving stiff FDE.

In Section 10 we describe two simplest methods: explicit Euler method and implicit Euler method with a constant step size and piece-wise constant interpolation. The results of realization of these methods for stiff and non-stiff problems are presented. We describe the technique of realization of the implicit Euler method for linear systems and for nonlinear systems (using Newton-Raphson method). In this section we also compare explicit and implicit methods of high orders with a constant step size.

In Section 11 explicit and implicit algorithms with automatic step size control are described. The results of their comparison for stiff and non-stiff problems are presented.

In Section 12 we applied the elaborated methods to simulation of a medical model - model of the human immunodeficiency (HIV) virus. In this section we also investigate and simulate an optimal control problem (linear quadratic control problem).

CHAPTER 1

Implicit one-step methods (with fixed step size) for functional differential equations

1. Statement of the problem

We consider the system with bounded delay

x(t) = f(t,x(t),xt(-)) (1.1)

with the initial conditions

x{t0) = x0, (1.2)

^W = !/°W,-r<s<0. (1.3)

Here f(t,x,y(-)) : [t0,t0 + 6} x R* x Q[-t, 0) -> Rz; [t0,t0 + 8] is the time interval (6 > 0); r > 0 characterizes the time-delay interval; R* is the /-dimensional Euclidean space with the scalar product (•, •) and the norm || • || = (•, -)1/2 ; Q[—r, 0) is the space of /-dimensional functions y(-) continuous everywhere on [—r, 0) except, may be, a finite set of points of discontinuity of the first kind (in which y(-) is continuous from the right) with the norm ||?/(-)||q — SUP {||z/(s)ll ' —t < s < 0} x g Kl is the phase vector of the system; xt(-) = {x(t -f s), — r < s < 0} € Q[—r, 0) is the pre-history, i.e. the function which characterizes delays.

We consider system (1.1) in the phase space H = Hl x Q[—r, 0), so the phase state of the system is the pair {x(t); £<(•)} G H [66].

We will make the following assumptions: the mapping f{t,xyy(-)) is shift-continuous [66] and satisfies the Lipschitz condition

<L\\xW-xW\\ + M\\yW(.)-yW(-)\\Q

for all t 6 [*o,*o + 0] and k {z(2), 2/(2)( )> € H\ where L and

M are constants. Under the above assumptions the initial value problem (1.1) - (1.3) has unique solution on the interval + 0} [66]. It is necessary to note, for continuous initial functions the conditions of existence

the solutions of problem (1.1) - (1.3) are less restricted [55].

2. Implicit Euler method with piece-wise constant interpolation

The aim of this Section is to demonstrate the basic idea of the approach on the example of the simple implicit Euler-like numerical scheme. Let us

consider the uniform partition tn = ¿o + ^A, n = 0,1,..., N, of the interval

g

[t0, to 4- 9}] here A = —. For the sake of simplicity we suppose that the

T

ratio — = m is a positive integer.

Our aim is to obtain on the interval [to, 6] an approximation to the solution x(t) of the initial value problem (1.1) - (1.3). In this section we describe implicit Euler-like method of evaluating the values un G RK n = 0,1,..., iV, which approximate x(t) at the points ¿o,- • -Jn', that is un ~ x(tn), n = 0,1,.. .,N. The sequence {un} which approximates the solution x(t) will be called the discrete model (numerical model, approximate model) of system (1.1). To find at time moment tn the next approximation un+1 using implicit Euler scheme it is necessary to calculate the right-part f{t,x,y(-)) of system (1.1) on the pre-history

{ui,n — m < i < n} (2.1)

of the discrete model. At time-moment tn the pre-history (2.1) of the discrete model is the finite set of vectors un_m,.. .,un, meanwhile the functional / in the right part of system (1.1) is defined on functions of H. Hence, to calculate the value of the functional / on the pre-history of the discrete model (2.1) it is necessary to make an interpolation (and extrapolation beyond tn) of the approximate solution un. One can use simple piece-wise constant interpolation

_ / t e fo' n-m<i<n ,

{)~\y\k-t), t£[t0-T,t0), {2-2)

and the corresponding discrete model

uq = Xq , (2.3)

un+\ = un + A/(in+i, un+i, utn+1(-)) (2.4)

will be called Implicit Euler method with piece-wise constant interpolation of the discrete pre-history (of the model).

Remark. Note, formula (2.4) is the equation with respect to un+\. In Section 1.4 we will show existence and uniqueness of the solutions of such equations for a sufficiently small step.

Definition 2.1. Numerical method

1) converges, if ||un — a;(in)j| —> 0 as A —> 0 for all n = 1,..., iV;

2) has the convergence order p, if there exists a constant C such that IK - x{tn)|| < CAP for all n = 1,..., N.

Show that implicit Euler's method (2.2) - (2.4) converges and has the convergence order p = 1.

Lemma 2.1. Let the solution x(t) of the initial value problem (1.1) -(1.3) is the continuous differentiate function. Then for the discrete model (2.3) - (2.4) with piece-wise constant interpolation (2.2) we have

ll®tn(-)-«J-)llg< m^x \\xi - Ui\\ + CiA, n = 0,...,N. (2.5)

n-m<i<n

Proof. Let x(t) is the piece-wise constant function defined by the rule

x{t) = {x{U), t G [th U+1), i = 0, • ■ •, N} , = (A/(•)}• Then

IK+i(0 - ^„+1(-)IIq = sup Mt) - «(i)|| <

tn~T<t<tn+1

< sup ||ic(£) — x(t) || + sup POO — w(i)|| <

tn—T<t<t„+1 tn — T<t<tn+1

< max ||X{ — «¿11 + CiA ,

n-m<i<n

where C\ = max ||i(i)||. The theorem is proved. t0<t<t0+8

The function

= Xn+1 ~ f(tn+1, ¡Cn+1, Xtn+1(-))

will be called the residual of implicit Euler's method.

Lemma 2.2. If the solution x(t) of the initial value problem (1.1) -(1.3) is twice continuous differentiable function then the residual of implicit Euler's method (2.2) - (2.4) has the first order with respect to A. Proof. Because

A2

x(tn) = z(£n+i) - A/(in+i, xn+i, xtn+1 (•)) + x(c)— , then for any n = 0,..., N — 1, we have j|^(tn)H < C2A, where

C2 = \ max ||i(i)ll-2 t0<t<t0+eu w"

Theorem 2.1. Let the solution x(t) of the initial value problem (1.1) - (1.3) is twice continuous differentiable function. Then implicit Euler's method (2.2) - (2.4) converges with the convergence order p = 1.

Proof. Denote en = ||a;n — un||. Using the definition of the residual and formula (2.4) we can estimate:

en+i — H^n+i - un+i|| = — \\xn +A f{tn+1 ) %Tl+1} 1

—un — A/(in+i, un+i, Win+1(-))|| <

< en + A\\f(tn+1 , un+1, utn+1

Taking into account Lemma 2 and properties of the functional / we have

+ A(L\\xn+i - «B+1|| + M\\xtn+1 (•) - utn+1 (•)||q) + C2A2.

Prom this inequality and Lemma 1 it follows that

(1 - AL)en+1 < en + AM max ||€i|| 4- A2(MCt + C2). (2.6)

n-m<i<n

Let step A is sufficiently small, that 2AL < 1.

Let us prove by the mathematical induction method the following inequality

e„<(l±^±ilr(MC1 + C2) A. (2.7)

For n = 0 the estimate (2.7) is satisfied, because eo = 0. Let us suppose that (2.7) is valid for all indices less than n + 1, and prove (2.7) for n + 1.

Assume that the maximum in the right side of (2.6) is achieved at some index io < n then

en + AMeio + A2(MCi + C2) . en+i < l^AL -

(l + A(M + l))n(MCi + C2)A (1 - AL)n+1 +

AM(1 + A (M + l))io(MCi + C2) A A2(MCi + C2) + (1 - AL)io+1 + 1 - AL -

(1 + A(M + l))"(MCi + C2)A(1 + MA + A) - (1 - AL)n+1

Thus the estimate (2.7) is proved.

q

Because N = — hence for all n — 0,..., N, A

e„<(i±^±i))i(MC1 + C2)A<

<

(1 + + < + ca)A .

This inequality completes the proof of the theorem.

The described implicit Euler's method with piece-wise constant interpolation is the simplest of the implicit converging methods. To obtain more accurate methods it is necessary to use high orders interpolational procedures and more complicated discrete models. Such methods will be discussed in the next sections.

3. Interpolation and extrapolation of the model pre-history

3.1. Interpolating operators

In this section we describe methods of interpolation and extrapolation of the pre-history of the discrete model un using functions composed by polynomials of p-th degree.

Let us consider the same partition of the time interval [¿o,£o + 0) as in the previous section. Remember, thiis partition is uniform only for the sake of simplicity.

Also remember that the discrete pre-history {wi}n of the discrete model {ui}-m at time moment tn is the set of m + 1 vectors:

{wi}n = {v>i G R*, n — m <i <n} .

This set of vectors defines at time moment tn the future dynamics of the discrete model.

Definition 3.1. Interpolational operator I of the discrete model prehistory is a mapping I : {ui}n ->• u(-) G Q[t n — T, tn I.

Definition 3.2. We say that an interpolational operator I has an approximation order p at a solution x(t) if there exist constants C\, C2 such that

\\x{t) - u(t)\\ < _ max. \\Ui - x{\\ + C2AP (3.1)

i>0,fi—NT<i<n

for all n = 0,1,..., N and t G [£„ — r, ¿J.

Example 3.1. The following mapping uses the piece-wise linear interpolation and is the interpolational operator of the order p = 2:

I ■ «W =

((£ - U)ui+1 + (£¿41 - t)ui)^ , t€[ti, ti+1], ^ 2j

y°{to-t), t e [t0-T,t0). The general interpolational operators I can be constructed using splines

171

of a degree p. Without loss of generality we can suppose that — = k is

P

a natural, otherwise we can take m divisible p. Let us divide the interval [tn ~ T,tn] = [tn-m,tn] by k subintervals [tHi+l, tUi], i = 0,1,..., k - 1, of the length pA in such a way that tno = tn, tHl = tn-pr... At every subinterval [¿nj+1, tnJ we construct an interpolational polynomial Lp(t) = L*p(t) according to the values un.-p% uni-p+1, ...,uni:

= n t^zt- (3-3)

/=0 j—rii—p; j^rii-l lrii-l Ij

Then we can define the following interpolational operator I (of the discrete pre-history)

Theorem 3.1. Let the solution x(t) of the initial value problem (1.1) -(1.3) is (p + l)-times continuous differentiable on the interval [¿o — t, to + 0]. Then the interpolational operator (3.4) has an approximation error p + 1.

Proof. Fix arbitrary n E {0,..., N} and t E [tn — t, tn]. Let t E [tni+l, tni] then

|x(t) - «(t)|| = IIx(t) - Vp(t)II < IIx(t) - Vp(t)II + ||LJ,(f) - l;(«)|| where

p m t — t

= n j-*-,

1=0 j=ni+i,j^ni-llni-llj

i.e. Lp(t) is the interpolational polynomial of the degree p constructed by the values Xj at the interval [tni+1,tni\- The interpolation error can be calculated according to the following formula [115, p. 133]

x(p+ i)ft\ p

II^PWII -1i*m - ^wn = n-^nyf n c - *»i-i)ii

as £ E [tni+1,tni\. Thus, if t E [i„j+1,i„J th