CC BY
125
ВЕСТНИК ТОМСКОГО ГОСУДАРСТВЕННОГО УНИВЕРСИТЕТА
2007 Управление, вычислительная техника и информатика № 1
УДК 519.865
В.В. Поддубный, О.В. Романович
МОДИФИКАЦИЯ МЕТОДА ЭЙЛЕРА С УРАВНИВАНИЕМ ДЛЯ РЕШЕНИЯ ДИФФЕРЕНЦИАЛЬНЫХ УРАВНЕНИЙ С ПЕРЕМЕННЫМ ЗАПАЗДЫВАНИЕМ
Рассматривается задача численного интегрирования обыкновенных нелинейных дифференциальных уравнений, разрешенных относительно производной, с переменным, в том числе случайным запаздыванием аргумента. Получена интерполяционная модификация вычислительной схемы Эйлера с уравниванием, пригодная для решения дифференциальных уравнений с произвольными, в том числе случайно изменяющимися и как угодно малыми запаздываниями, когда становится неприменимым известный метод шагов. Показано, что сохраняется второй порядок точности этой модификации как для итерационной структуры алгоритма метода Эйлера с уравниванием, так и для ее первого приближения (метода Хьюна).
Проблема численного интегрирования обыкновенных дифференциальных уравнений с запаздывающим аргументом при постоянных или переменных, но строго положительных конечных запаздываниях решается сравнительно просто с использованием алгоритмов численного интегрирования обыкновенных дифференциальных уравнений без запаздывания на основе известного метода шагов [1]. Однако эта проблема резко усложняется в случае, если переменные запаздывания (в том числе случайные) могут принимать как угодно малые значения. В этом случае метод шагов неприменим. А поскольку любые методы численного интегрирование дифференциальных уравнений оперируют с конечным, отличным от нуля шагом интегрирования, не обязательно соизмеримым с произвольно изменяющимися запаздываниями, в процессе численного интегрирования приходится прибегать к интерполяции решения в точках, где значения решения не были вычислены.
В [4] предложен алгоритм такой интерполяции для случая произвольного, в том числе как угодно близкого к нулю, постоянного запаздывания, модифицирующий известную схему метода Эйлера с уравниванием второго порядка точности. В данной работе этот подход развивается и распространяется на случай переменных запаздываний, в том числе случайных и как угодно малых. При этом проводится оценка порядка точности модифицированного для случая переменного запаздывания итерационного метода Эйлера с уравниванием. Показывается, что модифицированный итерационный метод Эйлера с уравниванием, как и его первое приближение (метод Хьюна), сохраняют второй порядок точности.
1. Постановка задачи
Рассмотрим начальную задачу для обыкновенного (в общем случае нелинейного) дифференциального уравнения, разрешенного относительно производной, с переменным запаздыванием т(г):
= f((У()>У(-т()) *о ^ ^ T , (1)
dt
с заданными начальными условиями
У 0) = Ф О), t е [to - max т (t), t0); y (t0) = y0,
(2)
где ф (г) - начальная функция, у0 - начальное значение решения, причем ф ((0) не обязательно равно у0.
Пусть функция ф (г) непрерывна, функция / (г, у, z) непрерывна по всем аргументам и удовлетворяет условиям Липшица по второму аргументу. Тогда, как известно [1], решение начальной задачи для дифференциального уравнения (1) существует и единственно. В случае возможного скачка начального условия (ф ((0) Ф у0), а также в случае, если начальная функция ф (г) не удовлетворяет дифференциальному уравнению (1), решение у (г) при г > г0 испытывает скачки производной первого рода. Непрерывное и всюду дифференцируемое решение возможно лишь в случае, если ф (г) удовлетворяет условию
При постоянном, заранее известном т для интегрирования дифференциального уравнения с запаздывающем аргументом (1) при начальных условиях (2) можно использовать известный метод шагов [1]. Однако в случае переменного т = т (7), тем более случайно изменяющегося во времени, метод шагов не применим. В этом случае задача интегрирования нелинейного дифференциального уравнения
(1) существенно усложняется. Существует множество вычислительных схем решения начальной задачи для нелинейного дифференциального уравнения с переменным запаздыванием [2]. Однако по большей части эти методы либо носят частный характер, либо весьма сложны в реализации.
В данной работе предлагается метод численного решения начальной задачи для дифференциального уравнения (1), обобщающий метод Эйлера с уравниванием при отсутствии запаздывания [3] и при постоянном запаздывании [4] на случай переменного (в том числе случайного и как угодно малого) запаздывания т (/).
Для дальнейшего анализа необходимо конкретизировать свойства функции х(г).
Пусть х(?) случайно изменяется во времени. В простейшем случае такое случайное изменение можно описать кусочно-постоянной случайной функцией, порождаемой случайным временным рядом тк, к = 1,2,..., в котором все тк неотрицательны, статистически независимы и имеют одинаковую плотность распределения р (тк). Например, р (тк) = 1/ттах - равномерное в интервале [0, ттах > 0] распределение. Тогда при любом разбиении г0 <г1 <••• <ги = Т интервала интегрирования [ ,Т] на п примыкающих друг к другу интервалов \гк, +1), к = 1, п -1, имеем т (г) в виде случайного «меандра»:
Поскольку любое запаздывание тк с некоторой отличной от нуля вероятностью может оказаться меньше любого заданного s > 0, хорошо известный для решения дифференциального уравнения с запаздыванием метод шагов [1] оказывается, вообще говоря, неприменимым (для его применимости необходимо, чтобы min т(t) > 0 ). В связи с этим для интегрирования дифференциального уравнения
dф (t0 - 0)/dt = f (t0, ф (t0), ф (t0 -т(t0))).
(3)
k=1
со случайно изменяющимся запаздыванием приходится применять вычислительные методы, использующие интерполяцию значений решения в точках, где не были вычислены значения решения при применении той или иной вычислительной схемы интегрирования.
Методы первого порядка точности (метод ломаных Эйлера и подобные ему) не используются при интегрировании дифференциальных уравнений, в том числе дифференциальных уравнений с запаздывающим аргументом, в силу его низкой точности (медленной сходимости к точному решению с уменьшением шага интегрирования). Поэтому рассмотрим схему второго порядка точности - метод Эйлера с уравниванием [3, 5], называемый также итерационным методом Эйлера [5], методом Эйлера с пересчетом [6], симметричной схемой Эйлера [7], и его упрощенный (неитерационный) вариант - метод Хьюна [8] (он же - усовершенствованный метод Эйлера - Коши [9]). Модифицируем эту схему для решения дифференциального уравнения с произвольным, в том числе и как угодно малым запаздыванием, введя в нее возможность интерполяции.
2. Обобщение метода Эйлера с уравниванием и метода Хьюна на случай переменного запаздывания
Рассмотрим дифференциальное уравнение (1) при произвольных изменениях запаздывания во времени: т = т(г). Разобьем интервал [г0,Т] на п примыкающих друг к другу подынтервалов равной длины к = (Т -г0)/п разбиением: г0 < г1 < — < гк < — < гп = Т , где гк = г0 + кк , к = 1, п. Проинтегрируем формально дифференциальное уравнение (1) на интервале [г0, г ] с учетом начальных условий (2):
г
у () = Уо + | / ( у (г), у (-т (г))) .
Используя это выражение при г = гк и I = 1к+1, получаем рекуррентное соотношение
*к +1
У(гг+1 ) = У(гк)+ I /(у(г),у(-т(*))) . (4)
(к
Считая функцию f (г, у (г), у (г - т (г))) непрерывной функцией своих аргументов, представим приближенно при малом шаге интегрирования к значение интеграла в правой части точного соотношения (4) по формуле трапеции:
У (**+1 ) = У (**)+ -2(f (>У (**). У -т (к))) +
+/ (гк+1> У(гк+1)> У (гк+1 - т (гк+1)))). (5)
Обозначая ук = у (¡к), тк = т (гк), представим рекуррентное уравнение (5), определяющее приближенное решение уравнения (1) при начальных условиях (2), в виде
к —
Ук+1 = Ук + -(/ (к > У к > У (к -тк)) +/ (1к+1, У к+1, У (4+1 -тк+1))), к п, (6)
с начальным условием у0 в момент времени (0 и начальной функцией ф (гк - тк)
при к -тк < Iо и Ф (¡к+1 - тк+1) при гк+1 - тк+1 < го .
В случае автономного уравнения (1) функция / (г, у (г), у (г-т (г))) не зависит явно от г: f = / (у (г), у (г-т (г))). В этом случае уравнение (6) принимает более простой вид:
к —
Ук+1 = Ук + ^(f (ук>у(-тк^+/(ук+1 >у(гк+1 -тк+1))), к п. (7)
Ограничимся в дальнейшем случаем автономного уравнения (1). Это не нарушает общности рассуждений, так как любое неавтономное уравнение вида (1) можно свести к автономной системе уравнений, расширив состояние у1 (() = у (() дополнительным состоянием у2 (/) = (, так что расширенное состояние примет вид двумерного вектора состояния у (г) = (у1 (г) у2 (г)). Тогда уравнение (1) станет системой уравнений автономного вида ¿у (г)/¿1 = / (у (г), у (г -т (г)), где
/ = / (( (), У1 (-т (/))) У2 (г)) - двумерная вектор-функция.
Для автономного дифференциального уравнения в качестве начальной функция ф(/) удобно взять функцию ф(?) = у*, тождественно равную постоянному равновесному значению решения (точке покоя) у*, обращающему в нуль правую часть автономного дифференциального уравнения ЖуУ /Ж = / (у*, у*) = 0 .
Уравнение (7) обобщает известное соотношение метода Эйлера с уравниванием на случай переменного тк. Как и в схеме Эйлера с уравниванием, уравнение (6) будем решать при каждом фиксированном к относительно ук+1 методом последовательных приближений (простых итераций), используя в качестве начального приближения у^ решение по схеме метода Эйлера, использующего вычисление интеграла в правой части (3) по формуле прямоугольников:
•Ук+1 = У к + ¥ {у к, у у-\к)), (8)
Ук+10 = Ук + 2(/(Ук>У(Ь-тк)) +/(У*+[,УМ (*+1 -т*+1))], v = 0,1,2,... . (9)
2
Окончание процесса последовательных приближений происходит по достижении требуемой абсолютной точности решения е:
И;;0 - уй| <е. (10)
Упрощенным вариантом итерационной схемы (9-10) является ее первое приближение у^ (обобщение метода Хьюна [8]).
Поскольку тк может принимать любое неотрицательное значение, не обязательно кратное шагу к, моменты времени гк - тк также не обязательно будут кратными к. Для приближенного вычисления значения у {1к+1 -тк+1) через значения в точках, ближайших слева и справа к точке ^+1 - тк+1, но кратных к, проведем линейную интерполяцию решений внутри соответствующих отрезков длиной к:
у (+1- т+1) = у к +^1 --ту ) у к+1, 0 *т к+1< ь;
У (к+1 -Тк+1 ) = [1 —+1 ~%1+1 ~-т ] У т + к+1 +1 - ^ Ут+1 ,
к — тк+1 — гк+1 — г0 , %+1 — тк+1 е [ш ’ +1 ] ;
У (к+1 -тк+1 ) = У , Тк+1 < ?к - ?0 . (11)
Как видим, интерполяция производится на том (первом при 0 < т к+1< к или т-м при к — тк+1 — гк+1 - г0) интервале, куда попадает гк+1 - тк+1.
Интерполяционная схема (11) обобщает на случай переменного запаздывания тк интерполяционную схему, предложенную в [4] для постоянного запаздывания т.
3. Порядок точности обобщенного метода Эйлера с уравниванием и интерполяцией
Найдем теперь порядок точности обобщенного на случай уравнений с переменным запаздыванием метода Эйлера с уравниванием и его упрощенного варианта - метода Хьюна. Воспользуемся определениями порядка точности и порядка аппроксимации вычислительных схем решения дифференциальных уравнений, приведенными, например, в [6]. Порядком аппроксимации (соответственно порядком точности) вычислительной схемы называют минимальную степень шага к в разложении в ряд Маклорена невязки вычислительной схемы метода численного решения (соответственно невязки численного решения). Невязкой вычислительной схемы называется разность между левой и правой частями схемы метода на точном решении, деленная на величину шага к на каждом к-м шаге решения. Невязкой решения называется разность между приближенным и точным решением на каждом к-м шаге решения. Можно показать, что порядок точности совпадает с порядком аппроксимации. Поэтому порядок аппроксимации часто называют также порядком точности. Конструктивный алгоритм вычисления порядка аппроксимации вытекает из определения невязки метода.
Запишем выражение для невязки метода Эйлера с уравниванием в соответствии с определением, обозначив точное решение дифференциального уравнения для к-го шага (в к-й момент времени) через ик = и (¡к), где и (г) является точным решением уравнения (1). Очевидно, ик подчиняется точному интегральному аналогу (4) дифференциального уравнения (1). Для автономного уравнения (1) рекуррентное соотношение (4) принимает вид
*г +1 _
ик+1 = ик + | / (и (), и (-Т ()))&, к = 1, п -1. (12)
(к
Невязка метода Эйлера с уравниванием и интерполяцией при случайных запаздываниях вида (3) принимает вид
ч к(к) = К+1- ик- (к/2)[ / (ик>и (ч -тк)) +/ (ик+1-и (гк+1 1 Ш/к, (13)
где от шага к на к-м шаге зависят ик+1 через верхний предел гк+1 = гк + к точного равенства (12) и и (гк+1 -тк+1) по формулам (11), куда вместо всех у подставляет-
ся и. При этом следует иметь в виду, что от к на к-м шаге зависят только ик+1 и гк+1 = гк + к. Величины к, стоящие в знаменателях формул (11), относятся не к к-му интервалу, а к более ранним интервалам, поэтому они в расчете порядка точности участия не принимают.
Разложив в ряд Маклорена правую часть равенства (12), получим с точностью до членов 4-го порядка малости
ик+1 = ик + /кк + (к(к + ЛкСт )) + [/11 к /к + 2/12к/кст + /22с1 +
к3
+/'('к/к + /кСт )] — + О (к4 ), (14)
3!
где обозначено:
/к = /(«к>«(гк-Тк)), /1 'к = д/(«к>«(к-Тк))/5«к ,
/2к = д/ («к >«(к - Тк))/( ~хк), /1!* = д2/ (ик,и (1к - тк))/ди1 ,
/12к = д2/ («к >«(к - Тк))//икд« (к -Тк), /22к = д2/(«к >« (к - Тк))/д« (к -Тк )2 ,
Ст = (ит+1 — ит )/^т , ^к — Тк+1 ^ [т ’ ^т+1 ] . (15)
В последнем равенстве через Нт обозначена длина т-го, более раннего, чем к-й, интервала [1т, 1т+1 ].
Разлагая аналогичным образом в ряд Маклорена функцию
/ («к+1, и (1к+1 тк+1 )) во втором слагаемом выражения (13) как сложную функцию к, с учетом (12) и
(11) (с заменой у на и) и с учетом обозначений (15) получим
f (ик+1, и (к+1 -Хк +1 )) = /к + (/1к/к + /2кст) к +
+ [/пк /к + 2/12к/кст + /22ск + Л (к(к + Лкст )]) + 0( ) . (16)
Подставляя выражения (14) и (16) с учетом обозначений (15) в (13), приведем ^¥к(к) к виду
^к (А) = 12[/к2 + 2/ш:/кС„ + /^ + /' (/¡к/к + /1кст)] к1 + О(к3). (17)
Из разложения (17) видно, что порядок точности метода Эйлера с уравниванием и интерполяцией остался равным 2.
Аналогичным образом можно получить разложением в ряд Маклорена по степеням к на к-м шаге невязку метода Хьюна
^к (А) = К-и -«к -(А/2)[/(«к>«(¡к -тк)) +
+/(ик + к/ (ик,и (1к -хк)), и (1к+1 -хк+1 ))]^к . (18)
Проводя разложение второго слагаемого в (18) как сложной функции к и используя разложение (14), с учетом обозначений (15) получим
^к (А) = 12[-^к /2 - 2/2к/кс„ - /¿Ск2 +2// (/к/к + /2кСт )] А2 + О(к3) . (19)
Отсюда видно, что порядок точности метода Хьюна с интерполяцией также остается равным двум. Заметим, что в случае линейных дифференциальных уравнений, для которых функция f линейна по своим аргументам, вторые производные в выражениях (17) и (19) обращаются в нуль, и мы видим, что функция ¥k (h)
для метода Хьюна оказывается при малых h ровно в 2 раза больше, чем для метода Эйлера с уравниванием, что свидетельствует все-таки о несколько меньшей скорости сходимости к точному решению вычислений по неитерационной схеме Хьюна по сравнению с итерационной схемой Эйлера с уравниванием.
ЛИТЕРАТУРА
1. Эльсгольц Л.Э. Введение в теорию дифференциальных уравнений с отклоняющимся аргументом. М.: Наука, 1964. 128 с.
2. Солодов А.В., Солодова Е.А. Системы с переменным запаздыванием. М.: Наука, 1980. 384 с.
3. Тихонов А.Н., Васильева А.Б., Свешников А.Г. Дифференциальные уравнения. М.: Наука, 1980. 233 с.
4. Поддубный В.В. Численное решение дифференциальных уравнений с постоянным запаздыванием методом Эйлера с уравниванием и интерполяцией // Обработка данных и управление в сложных системах. Вып. 7. Томск: Изд-во Том. ун-та, 2005. С. 165 - 174.
5. Эльсгольц Л.Э. Дифференциальные уравнения и вариационное исчисление. М.: Наука, 1965. 324 с.
6. Бабенко К.И. Основы численного анализа. М.: Наука, 1986. 744 с.
7. Самарский А.А., Гулин А.В. Численные методы. М.: Наука, 1989. 432 с.
8. Ортега Дж., Пул У. Введение в численные методы решения дифференциальных уравнений. М.: Наука, 1986. 288 с.
9. Коллатц Л. Численные методы решения дифференциальных уравнений. М.: ИЛ, 1953.
Статья представлена кафедрой прикладной информатики факультета информатики Томского государственного университета, поступила в научную редакцию 29 июня 2007 г.
