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

Изучение многих процессов, происходящих в природных и технических системах, сводится к анализу свойств их математических моделей, что приводит к необходимости исследования систем обыкновенных дифференциальных уравнений (ОДУ). Часто в приложениях встречаются системы ОДУ, у которых матрица при производной является вырожденной. Такие системы называются системами сингулярных дифференциальных уравнений или дифференциально-алгебраических уравнений (ДАУ). К решению систем ДАУ приводят многие из задач механики, кинетики химических реакций, теории управления, электрических цепей (см., например, [14, 40, 41, 42, 44]).

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

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

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

Подобные задачи возникают также и в механике. В работе [14] процесс вибросверления описывается при помощи системы ДАУ с запаздывающим аргументом.

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

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

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

1) у: К1 —> Мп, х: Ж1 —>• Мт, * Е М1 /: Еп+т+1 —> Кп, <2: Шп+т+1 —> Ет. Начальные условия имеют вид.

У (0) = Уо, ж (0) = х0.

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

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

G (y0, xq, t0) = 0.

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

Численное решение начальной задачи (1), (2) впервые, по-видимому, исследовалось в работе 1971 года (Gear C.W.) [65]. В данной работе рассматривалось решение системы разрешенных относительно производных уравнений, описывающих процессы, протекающие в электрических цепях. Для дискретизации данной системы использовались формулы дифференцирования назад (ФДН). Система нелинейных уравнений, возникающая на каждом шаге процедуры интегрирования, решалась при помощи метода Ньютона.

Первый результат по сходимости для методов ФДН был получен в 1985 году (P. Lotstedt) [74]. Доказано, что для систем индекса 1 ФДН-метод сходится с тем-же порядком, что и для ОДУ.

В работах [55, 75] получены результаты о сходимости схем, основанных на ФДН для систем индекса 2. Доказано, что /с-шаговая ФДН-схема при к < 6 сходится с порядком р — к, если погрешность стартовых значений имеет порядок 0{hv+l). В работе [67] доказательство сходимости схем, основанных на ФДН дано для случая переменной длины шага.

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

В работах R. Marz [76] рассмотрен частный случай системы (1): применялись многошаговые методы.

Исследование сходимости общих многошаговых методов проведено в работе [68].

Однако, как отмечается в монографии [54], среди многошаговых методов наилучшие результаты при решении прикладных задач показывают методы, основанные на ФДН. Причиной этому являются хорошие свойства устойчивости и точности данных методов.

Первые результаты о сходимости неявных методов Рунге-Кутты были получены в работах [62, 68] для ДАУ индеска 1. При помощи теоремы о неявной функции доказано, что если метод является жестко точным, то порядок метода будет таким-же, как и для ОДУ.

Для методов Рунге-Кутты, не являющимися жестко точными, чьи решения обычно не удовлетворяют уравнениям связи (недифференциальным где /: Ж2тг+1 —> Mn, C (t) — матрица п х п. Для решения таких задачсоотношениям), в работе [59] предложен следующий подход: решение ищется из системы нелинейных уравнений тогда как уи+х вычисляется при помощи неявного метода Рунге-Кутты. При таком подходе метод имеет такой-же порядок, как и для ОДУ. Недостатком является увеличение размерности нелинейной системы.

В [56] аналогичные исследования были проведены для ДАУ высших индексов. Были получены оценки погрешности для х и у-компоненты. В работе [69] данные результаты были улучшены (использовался другой подход — локальная и глобальная погрешности исследовались отдельно).

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

Серия работ Г. Ю. Куликова (см., например [33, 34]) посвящена рассмотрению весьма частного случая систем ДАУ индекса 1 вида.

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

Ук+Ъ хк+1, Ь+г) = 1(у, х), у (0) = уо, х^) = д (у, х), х{0) = х0. и получены оценки полной погрешности, которая складывается из погрешности дискретизации и погрешности метода Ньютона.

В работе [11] был предложен метод параметризации задач оптимального управления, который позднее (см., например, [13, 36]) был применен для решения задач вида (1), (2). Согласно данному методу приближенное решение представляется в виде сплайна с подвижными узлами, параметры которого определяются из условия минимизации невязки сингулярной части системы (1). Такой сплайн назван вариационным.

Большое число работ посвящено исследованию линейных систем ДАУ, имеющих вид.

A (t)x (t) + B (t)x (t) = f (t), х (0) = Xq, где det A{t) = 0.

В работе Ю. Е. Бояринцева [4] для систем вида (3) с постоянными матрицами, А и В предложен метод возмущения, который состоит в переходе от исходной системы с вырожденной матрицей при производной к нахождению решения системы ОДУ с матрицей невырожденной, но близкой к исходной. Данные алгоритмы предложены для случая регулярного пучка матриц Л, А + В. Дальнейшее развитие методы возмущения получили в работах [45, 68]. В работе [6] предложен алгоритм возмущения, который применим и для систем с сингулярным пучком матриц XA (t) + В (t). В данной работе система (3) сначала записывается в интегральной форме, затем системе интегральных уравнений сопоставляется возмущенная система, после чего полученное равенство дифференцируется и вновь получается задача Коши для системы ОДУ, но уже с невырожденной матрицей при производной.

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

В.К. Горбуновым в работе [12] для решения систем вида (3) предложен метод нормальной сплайн-коллокации, который заключается в переходе от системы дифференциальных уравнений к конечной системе дифференциальных соотношений в узлах коллокационной сетки и в определении элемента, удовлетворяющего этим соотношениям и начальным условиям и имеющего минимальную норму гильберто-соболевского пространства!^^.

На сегодняшний день численному решению систем вида (1), (2) посвящено очень большое количество работ. Достаточно полно современное состояние проблемы отражено в монографиях Э. Хайрера и Г. Ваннера [44], Ю. Е. Бояринцева [5], В. И. Шалашилина и Е. Б. Кузнецова [49], докторских диссертациях М. В. Булатова [8], Г. Ю. Куликова [35] и В. Ф. Чистякова [46], а также в работах [1, 13, 54].

Необходимо отметить, что решение ДАУ является более сложной задачей по сравнению с решением ОДУ. В [54, 65] отмечаются следующие трудности:

• начальные условия должны быть согласованы с недифференциальными соотношениями;

• система уравнений плохо обусловлена для мелких шагов интегрирования;

• ошибка метода чувствительна к несогласованности в начальных условиях и к резкому изменению решения;

• численное решение в большей степени зависит от точности аппроксимации, чем для ОДУ.

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

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

Перейдем к рассмотрению интегродифференциально-алгебраических уравнений с запаздывающим аргументом, которые представляют собой систему взаимосвязанных дифференциальных, нелинейных алгебраических и интегро-дифференциальных уравнений с запаздывающим аргументом (У) У: Ут) Ут) X, Хт, ХТ) ¿-г- = 0) 0(у, ут, х, хт, г,2 т, 1) = 0,.

4).

0 / г = 1, к. у: M Rn, X: R M™, z: R Rfc, t 2(n+m+fc)+l ^ jgmz q. ]g>37i+3m+4&+l > -jj^fc.

Здесь индекс г > 0 определяет запаздывание аргумента функции, т. е. 2/тСО = - т), ?r (?) = x (t — г), zT (?) = z (t — г).

Начальные условия для системы (4) имеют вид y (t) = y (t), № = Ш x (t) = x (t), ±(t) = x (t), t G [¿-о — T, to), [z (t) = z (t), = y (tо) = y (tо) = 2/0, z (i0) = = х0, z (t0) = z (t0) = z0.

Начальные условия должны быть согласованными, т. е. должно выполняться равенство.

G (y (to), yT (to), x (to), xT (t o), z (tQ), zT (t0), t0) = 0.

В отличие от задачи (1), (2), численному решению которой посвящено множество работ, численное решение систем вида (4), (5) практически не изучено.

В работе [71] система, названная авторами системой интегрально-дифференциально-алгебраических уравнений, состоящая из интегро-дифференциальных уравнений и интегральных уравнений, решается при помощи метода «waveform relaxation», основанного на расщеплении нели-нейностей и методе Ньютона, а также при помощи многошаговых методов.

В работе [80] для системы, состоящей из дифференциальных, алгебраических уравнений и уравнений Фредгольма, предложен метод сведения такой системы к системе ДАУ.

При отсутствии интеграла, система (4), (5) представляет собой систему ДАУ с запаздывающим аргументом. Количество работ, посвященным таким системам, также невелико. В работе [37] для линейных дифференциально-алгебраических систем управления с запаздыванием, доказана экспоненциальная оценка роста решений. В работе [32] решение систем ДАУ с запаздывающим аргументом исследовано с позиции метода продолжения решения по наилучшему аргументу и доказано, что таковым является длина дуги интегральной кривой задачи.

Некоторое количество работ посвящено рассмотрению так называемых вырожденных интегро-дифференциальных уравнений (ИДУ), имеющих вид где йеЬ А (&euro-) = 0.

В работе [7] выделен класс систем вида (6), которые имеют единственное непрерывное решение. Для таких задач предложен численный метод решения, основанный на неявном методе Эйлера и квадратурной формуле левых прямоугольников.

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

Следует отметить, что все трудности, присущие численному решению систем вида (1), (2), перечисленные выше, остаются справедливы и для систем вида (4), (5). Кроме того, добавляются трудности, связанные с наличием запаздывания и интеграла. Для преодоления этих трудностей в данной работе предлагается метод продолжения по наилучшему аргументу, которым является длина дуги интегральной кривой задачи.

Диссертация состоит из 4 глав. В первой главе приводится обзор основных методов решения систем ДАУ и вырожденных интегро-дифференциальных уравнений. Вторая глава посвящена численному реб) о ж (0) = х0, шению систем ДАУ при помощи метода продолжения по наилучшему аргументу. Излагаются общие идеи метода и алгоритм, названный непрерывным продолжением [30], основанный на дифференцировании недифференциальных соотношений. Предлагается метод, названный дискретным продолжением [27], в котором уравнения продолжения получаются без дифференцирования недифференциальных соотношений. Получена оценка погрешности решения системы нелинейных уравнений, решаемой на каждом шаге процедуры интегрирования. Эффективность метода дискретного продолжения демонстрируется на тестовых примерах. В третьей главе рассматривается численное решение систем интегродифференциально-алгебраических уравнений с запаздывающим аргументом. Доказано, что для того, чтобы преобразовать такую систему к наилучшему аргументу, необходимо и достаточно выбрать в качестве такового длину дуги интегральной кривой. Предложены методы непрерывного и дискретного продолжения для таких систем. Эффективность предложенного преобразования демонстрируется на ряде тестовых примеров. Глава 4 посвящена численным результатам решения прикладных задач. Рассмотрено численное решение уравнений модели процесса вибросверления, а также численное решение систем уравнений, моделирующих линейную и нелинейную нестационарные электрические цепи. Продемонстрированы преимущества предложенного преобразования, а также произведено сравнение эффективности методов непрерывного и дискретного продолжения.

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

• XII международный симпозиум «Динамические и технологические проблемы механики конструкций и сплошных сред» (13−17 февраля 2006 г., Ярополец).

• VIII Харитоновские чтения по проблемам физики высоких плотностей энергии (21−24 марта 2006 г., Саров).

VII Международная научная конференция «Дифференциальные уравнения и их приложения» (17−19 мая 2006 г., Саранск).

VI Международная конференция по неравновесным процессам в соплах и струях NPNJ-2006 (26 июня — 1 июля 2006 г., Санкт-Петербург).

General Linear methods and Differential Equations 2008 (GLADE'08) (Auckland, New Zealand, 14−25 July 2008).

XVI международная конференция по вычислительной механике и современным прикладным программным системам ВМСППС'2009. (2531 мая 2009 г., Алушта).

12th Seminar NUMDIFF on Numerical Solution of Differential and Differential-Algebraic Equations (Halle, Germany, 14−18 September 2009).

XX крымская осенняя математическая школа-симпозиум «Спектральные и эволюционные задачи» (18−29 сентября 2009 г., Севастополь, Украина).

Совместный семинар кафедр теоретической механики и дифференциальных уравнений Московского авиационного института (г. Москва, 2009 г.).

Заключение

.

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

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

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

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

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