CC BY
5
УДК 519.62/.642 Вестник СПбГУ. Прикладная математика. Информатика... 2019. Т. 15. Вып. 4 МБС 65Ь03, 65Ь06
Комбинированный функционально-непрерывный метод для дифференциальных уравнений с запаздыванием
А. С. Еремин
Санкт-Петербургский государственный университет, Российская Федерация, 199034, Санкт-Петербург, Университетская наб., 7—9
Для цитирования: Еремин А. С. Комбинированный функционально-непрерывный метод для дифференциальных уравнений с запаздыванием // Вестник Санкт-Петербургского университета. Прикладная математика. Информатика. Процессы управления. 2019. Т. 15. Вып. 4. С. 425-441. https://doi.org/10.21638/11702/spbu10.2019.402
В работе предлагается комбинированный метод решения дифференциальных уравнений с дискретным запаздывающим аргументом, представляющий собой вложенную пару: непрерывный шестиэтапный и функционально-непрерывный (этапно-непрерывный) семиэтапный методы Рунге—Кутты четвертого порядка, которые совместно позволяют организовать эффективное решение дифференциальных уравнений с дискретными запаздываниями (ДУЗА). Комбинированный метод остается полностью явным для любых значений запаздываний. Выбор между шести- и семиэтапным методами происходит автоматически в зависимости от попадания запаздывания в текущий шаг. Тем самым на каждом шаге выполняется лишь необходимое число оценок правой части. В этих методах также реализовано повторное использование последнего этапа на новом шаге, что дополнительно снижает вычислительные затраты на каждом шаге. Представлены условия порядка и доказательство их разрешимости. Приводятся результаты тестовых расчетов, демонстрирующие эффективность построенного метода.
Ключевые слова: запаздывание, ДУЗА, непрерывный метод, функционально-непрерывный метод, этапно-непрерывный метод.
1. Введение. Рассмотрим дифференциальное уравнение с дискретным запаздывающим аргументом (ДУЗА)
и(*) = /{ь,и(г),и(г — т)), * > г0, (1)
с предысторией
и(г) = ф), * < *о. (2)
Величина т отклонения аргумента функции и(Ь — т) называется запаздыванием. В общем случае функция правой части / может зависеть от значения решения в нескольких запаздывающих значениях времени (*—т\, *—Т2 и т. д.). Однако рассматриваемый в рамках настоящей работы численный метод формулируется для одного запаздывания, так что для упрощения записи ограничимся формой (1). Формулировка метода легко распространяется на случай нескольких запаздываний.
Математические модели, использующие ДУЗА, встречаются во многих отраслях: биологии, популяционной динамике, медицине, иммунологии и эпидемиологии, экономике, нелинейной оптике, химии, механике, в задачах теории управления и во многих других (см., например, обобщающие монографии [1, 2] и литературу в них).
Если запаздывание т постоянно, то численное решение задачи (1), (2) может быть выполнено произвольным известным методом решения обыкновенных дифференциальных уравнений (ОДУ) с постоянным шагом, целое число раз укладывающимся
© Санкт-Петербургский государственный университет, 2019
в т. Особый интерес для решения численными методами представляют уравнения с запаздыванием общего вида: зависящим от времени и/или самого решения в текущий момент времени, т. е. т = т(t,u(t)). Для полноты описания приведем условия существования и единственности решения в этом случае [3].
Теорема 1. Рассмотрим уравнение (1) с запаздыванием общего вида т = т(t,u(t)) и предысторией (2). Пусть U и V — окрестности у(t0) и у(t0 — т(t0, y(t0))) соответственно. Если функция f (t,u,v) непрерывна по первому аргументу и лип-шецева по второму и третьему в [to,to + h] х U х V для некоторого h > 0, и если, кроме того, предыстория y(t) липшецева при t ^ t0, а запаздывание т(t, u) ^ 0 непрерывно по t и липшецево по второму аргументу в [t0,t0 + h] х U, то задача (1), (2) имеет единственное решение на интервале [to, t0 + ó) для некоторого ó > 0, и это решение непрерывно зависит от начальных данных.
В такой ситуации требуется применять численные методы, дающие приближение к решению не только в конечных точках совершаемых шагов, но и в любой точке промежутка решения. Эти методы называются непрерывными. Так называемому стандартному подходу к решению ДУЗА на примере непрерывных методов Рунге—Кутты посвящен п. 2.
В тех случаях, когда запаздывание становится слишком маленьким или в некоторых точках исчезает, запаздывающий аргумент попадает в текущий шаг (ситуация, называемая перекрытием) и реализация любого классического метода решения ОДУ становится полностью неявной. Само собой, применение изначально неявных методов, например полностью неявного метода Радо, лежащего в основе программы RADAR [4] — одного из наиболее популярных кодов для численного решения ДУЗА, практически не усложняется в случае перекрытия. В то же время основное применение неявные методы находят в случае жестких задач, но большое число возникающих при моделировании систем уравнений с запаздываниями можно отнести к нежестким. Их можно с успехом решать и явными методами, скорость которых несравнимо выше. Само собой, потеря явности даже на нескольких шагах существенно снижает скорость работы методов, поэтому появились методы, сохраняющие явную реализацию и для малых значений запаздывания. В п. 3 приводится их описание, а в п. 4 формулируется теорема о минимальном необходимом количестве вычислений правой части в зависимости от возникновения перекрытия на этапах.
Поскольку метод Рунге—Кутты четвертого порядка в случае перекрытия требует семи этапов, а без него — шести, эффективная реализация метода решения задачи (1), (2) должна в зависимости от наличия перекрытия выбирать, какой метод использовать. В п. 5 предложена вложенная пара методов, один из которых имеет шесть этапов и применяется в отсутствие перекрытия, а второй, семиэтапный, когда оно есть. В начале шага всегда предполагается, что перекрытия нет, и расчет начинается шестиэтапным методом. Вложенность проявляется в том, что при обнаружении перекрытия на любом из этапов все выполненные ранее вычисления используются семиэтапным методом, а общее число вычислений правой части на шаге не превышает семи.
В п. 6 приводится решение тестовых задач, подтверждающее заявленный порядок сходимости предложенного метода, а также он сравнивается с полностью семи-этапным методом, предложенным ранее в [5], по критерию «затраты/точность».
2. Стандартный подход на основе непрерывных методов. В простых случаях решение ДУЗА можно проводить с помощью непрерывных численных методов для ОДУ. Их отличие от классических дискретных методов, которые дают приближе-
ние к решению на некоторой сетке, заключается в том, что они дополнительно содержат способ приближения решения в произвольной точке отрезка решения. В частности, в многошаговых методах можно строить интерполяционный полином через точки сетки. В одношаговых методах интерполянт обычно делают локальным, т. е. опирающимся только на те же данные (оценки правой части уравнения) внутри шага, по которым находится приближение в его конце. Например, неявные коллокационные методы (см. [3, 6]) содержат такой интерполянт по определению. Другим вариантом является построение локальных эрмитовых сплайнов по значениям решения в точках сетки и их производных, что для явных методов Рунге—Кутты доступно всегда, если в задаче в явной форме выделена производная.
Рассмотрим более общий вариант непрерывного метода Рунге—Кутты (НРК) применительно к (1). Предположим, что выполнено n шагов численным методом и решение приближено на отрезке \t0,tn] с помощью интерполянта r/(t) (n(t) = V>(t), t ^ to). Явный s-этапный НРК дает приближение к решению уравнения (1) в любой точке (n + 1)-го шага через доопределение интерполянта n(t) по формуле
(tn + вН) « n(tn + вН) = n(tn) + ^bi(e)Ki, в e \0,1],
(3)
i=i
где
i-1
Ki = f(Ti,Yi,n{Ti - T (Ti ,Yi)j), Ti = tn + Cih, Yi = + Kj. (4)
j=i
Коэффициенты метода, как правило, представляются в виде так называемой таблицы Бутчера
(5)
0 0
С2 «21
сз «31 аз2
cs as 1 as 2 a>s,s-1
bi(0) ъ2(в) ■ • т ьа{в)
В ней отмечено, что Ъ^{в) — это функциональные коэффициенты, позволяющие построить непрерывное приближение.
При описании непрерывного приближения необходимо рассмотреть вопрос о его точности. В отличие от дискретных методов, в которых нас интересует погрешность в точках сетки, для непрерывных методов следует оценивать погрешность в каждой точке интервала вычислений.
После выполнения п шагов методом (3), (4) на следующем шаге решаем задачу
U(t) = f(t,u(t),u(t - T)) , t > tn, u(t) = n(t), t < tn,
(6)
где непрерывное приближение к решению п(^) уже определено вплоть до Ьп. Пусть ¿п+г^) — точное решение задачи (6). Рассмотрим отрезок от ~Ъп до ¿п+1 — ^п + к, h> 0.
Дадим несколько определений, которые потребуются в дальнейшем.
Определение 1. Локальной погрешностью (погрешностью на шаге) метода (3), (4) при решении задачи (6) назовем величину
e(t) = zn+i(t) - n(t), t e [tn,tn+i]-
Определение 2. Метод (3), (4) имеет дискретный локальный порядок p, если p ^ 1 — наибольшее целое число такое, что для любой правой части задачи (6) f e Cp и достаточно малого шага h найдется константа C > 0 такая, что
\\e(tn+i)\\ < Chp+1.
Определение 3. Метод (3), (4) имеет равномерный локальный порядок q, если q ^ 1 — наибольшее целое число такое, что для любой правой части задачи (6) f e Cq и достаточно малого шага h найдется константа C' > 0 такая, что
max \\e(t)\\ < C'hq+1.
[tn,tn+i]
Отметим, что, в силу определения, p ^ q, т. е. равномерный порядок q гарантирует и дискретный порядок q. Обратное может быть неверно. Нередко в практически применяемых методах непрерывное приближение имеет порядок, меньший, чем приближение в точках сетки. Например, использование упомянутой выше эрмитовой интерполяции, если метод Рунге—Кутты имеет порядок не ниже третьего, дает равномерный порядок три. В пакете MATLAB используется как раз такой вариант построения непрерывного приближения для классического метода Рунге—Кутты дискретного порядка четыре [7].
Определение 4. Глобальной погрешностью метода (3), (4) при решении задачи (1) в точке t ^ to шагов назовем величину
E(t) = u(t) - n(t).
Пусть на отрезке решения задачи [to,tf] сделано N шагов длин hn, n = 1,..., N.
Определение 5. Метод (3), (4) имеет дискретный глобальный порядок p (дискретный порядок сходимости), если p ^ 1 — наибольшее целое число такое, что для любой правой части f задачи (1), для которой на каждом шаге f e Cp, и достаточно малого максимального шага h = maxn hn, n = 1,..., N, найдется константа C > 0 такая, что
max_||i?(tn)|| < Chp.
n=l,N
Определение 6. Метод (3), (4) имеет равномерный глобальный порядок q (равномерный порядок сходимости), если q ^ 1 — наибольшее целое число такое, что для любой правой части f задачи (1), для которой на каждом шаге f e Cq, и достаточно малого максимального шага h = maxn hn, n = 1,..., N, найдется константа C' > 0 такая, что
max \\E(t)\\ < C'hq.
te[to,tf ]
В работе [8] была изучена общая процедура построения явных непрерывных методов Рунге—Кутты и установлено, что для построения явного метода равномерного
порядка три требуются как минимум четыре этапа, для четвертого — шесть, а для пятого — восемь (если речь не идет о системах со структурными особенностями, в которых возможно снижение числа этапов, см., например, [9-12]). Рассмотрим метод четвертого порядка с шестью этапами:
0 0
1 1
6 6
11 44 363
37 1369 1369
11 3388 8349 8140
17 4913 4913 4913
13 36764 767 32708 210392
15 408375 1125 136125 408375
1 1697 0 50653 299693 3375
18876 116160 1626240 11648
bi(0) 0 Ьз(в) Ь4(0) Ь5 (в)
bi(0) = _ 104217 и — - 37466 2 1806901 з 866577 4 618189 824252 '
Ьз(в) = 861101 2 230560 2178079 3 12308679 4 -в Н--в , 380424 5072320 '
Ъ±(в) = 63869 2 293440 6244423 3 7816583 Н--в--в , 5325936 10144640 '
Ьъ(в) = 1522125 --6 762944 2 982125 о 624375 4 + в в , 190736 217984 '
Ьа(в) = 165,2 _ 461 з 296^ 131 131 131
Следует обратить внимание на то, что в представленном методе приближение Уб на (п + 1)-м шаге оказывается равным r/(tn + h). Действительно, £>¿(1) = a^i, г = 1,6. Это возможно из-за того, что четвертый дискретный порядок может быть построен за четыре этапа. Таким образом, для получения приближения в конечной точке шага используются пять первых этапов, а шестой нужен только для непрерывного четвертого приближения порядка. Очевидно, что в таком случае Yi на новом шаге совпадает с Уб на предыдущем, и метод по своей реализации оказывается пятиэтапным. Такое повторное использование (англ. reuse) иногда называют FSAL (от англ. First Same as Last — первый, как последний) [6]. Мы будем строить метод на основе этой идеи.
Для ОДУ известно, что если на каждом шаге метод демонстрирует дискретный локальный порядок p, то он сходится с дискретным порядком p [6]. При этом равномерный локальный порядок q позволяет получить равномерный порядок сходимости max(p,q +1) [3]. Для ДУЗА, однако, глобальный порядок сходимости, как равномерный, так и дискретный, не превосходит max(p, q +1), потому необходимо строить непрерывные методы с равномерным порядком не ниже p - 1.
Это утверждение верно в том случае, если решаемая на каждом шаге начальная задача является достаточно гладкой, т. е. если ее решение на каждом шаге имеет непрерывную производную порядка p. Для уравнений с запаздывающим аргументом характерно наличие точек разрыва производных (см., например, [3]), которые должны быть включены в сетку по времени так, чтобы они не являлись внутренними точками ни на одном из шагов. Задача обнаружения и достаточно точного прибли-
жения точек разрыва описана, например, в [13, 14] и выходит за рамки настоящей работы, в которой мы рассматриваем локальное поведение конструируемого метода.
Отметим также, что предпочтительнее использовать методы с полным локальным равномерным порядком, т. е. при q = p. В этом случае контроль длины шага (как через оценку погрешности, так и через оценку невязки) оказывается более удобным и точным. Кроме того, такие методы без потери порядка можно применять к нейтральным уравнениям (таким, где ù{t) зависит и от ù{t — а) для некоторого запаздывания а). В связи с этим далее строится метод равномерного порядка четыре.
3. Функционально-непрерывные методы. Упомянутый стандартный подход на основе непрерывных методов применим в том случае, когда длина шага h < то = min t{t,u{t)). Однако в наиболее общем случае переменное запаздывание т может быть очень маленьким, что влечет за собой слишком малую длину шага для сохранения порядка метода, или вообще становится нулевым (так называемый случай исчезающего запаздывания). В последнем случае при практической реализации даже явные методы становятся полностью неявными, так как из-за попадания запаздывания внутрь совершаемого шага величины Y в (4) зависят от интерполянта n{t), который строится по ним самим и всем последующим. Ситуацию, в которой запаздывающее решение u{t — т) необходимо вычислить внутри текущего шага, будем называть перекрытием (англ. overlapping).
Еще в начале 1970-х годов появилась работа, в которой был предложен способ сохранить полностью явную реализацию одношаговых методов для уравнений с произвольным запаздыванием [15]. В ней рассмотрено расширение НРК. В [5, 16] методы такого типа получили название функционально-непрерывных, так как они сконструированы специально для решения функциональных дифференциальных уравнений запаздывающего типа. В них была построена теория и выведены условия равномерного и дискретного порядков для функционально-непрерывных методов Рунге—Кутты (ФНРК). Следует отметить, что в некоторых других статьях методы с подобными свойствами называются этапно-непрерывными (англ. continuous-stage methods), что точнее описывает их структуру (см. например, [17]).
Замечание 1. Вообще говоря, ФНРК разработаны для случая функционально-дифференциальных уравнений (ФДУ), которые являются более общими, чем уравнения с дискретным запаздыванием (см., например, [18]). Численные методы для ФДУ, в том числе обобщение методов Рунге—Кутты, разрабатываются уже длительное время, в частности они широко представлены в [19]. Заметим, что метод, предложенный в данной статье, содержит в себе часть, которая является полноценным методом ФНРК для функционально-дифференциальных уравнений. Однако в случае ФДУ предполагается, что перекрытие происходит всегда, и получить какой-то выигрыш в числе этапов для общего случая ФДУ по сравнению с семиэтапным методом с повторным использованием, описанным в [20], не получится. Именно сужение области применения ФНРК на случай уравнений с дискретным запаздыванием дает возможность скомбинировать его с НРК и получить эффективный инструмент решения ДУЗА. Для полноты изложения в п. 4 приведем условия порядка НРК и ФНРК, которым удовлетворяют построенные нами методы.
Поскольку при совершении {n + 1)-го шага интерполянт n{t) еще не определен для t ^ tn, явная реализация метода типа Рунге—Кутты требует для каждого этапа построения своего интерполянта, использующего только уже явно доступную к моменту его применения информацию.
Как и в случае НРК (3), (4), предположим, что n шагов численным методом уже
выполнено и на отрезке [to,tn] построен интерполянт r/(t). Как и раньше, n(t) = v(t), t ^ to. На (n + 1)-м шаге интерполянт n(t) доопределяется следующим образом:
t(tn + Oh) « n(tn + Oh) = n(tn) + hY^bi(0)Ki, в e [0,1],
i=i
где
i-1
Ki = f(Ti,Yi,r/i{ Ti - T (Ti,Yt)fj, Ti = tn + cih, Yi = n(tn) +a^ K3,
j=i
а этапные интерполянты вычисляются как
s
ni(tn + Oh) = n(tn) + hJ2aij (O)Kj, O e [0,Ci].
i=i
(8)
(9)
(10)
Помимо функциональных весовых коэффициентов Ъ^(в) в ФНРК присутствуют функциональные этапные коэффициенты а^(в). При этом естественным является требование совпадения У^ со значением интерполянта гц в отсутствие запаздывания, т. е. при в = е^. Равенство щО = У^ позволяет представить ФНРК с помощью таблицы Бутчера (см. (5)), лишь заменив постоянные коэффициенты а^ на их непрерывные аналоги для метода (8)—(10). Приведем коэффициенты ФНРК четвертого равномерного порядка из [5]:
0
1
1
2
0
в
2в I-
I*2 2 V 2
>2 + -О3 3 0 2в2 -з* V + 2 3*
>2 + -О3 3 0 2в2 ~2в + 2 3*
)2 + -в3 3 0 2в2 -V 3 2 + и 3
(11)
bi(0)
bl (0) = в - Ц.в2 + зо3 - У, ьб(0) = д-в2 - Що3 + 27
Ь6(в) = -%2-6в3-Щ-в\
Ьь(в) Ъв(9) Ъ7(в) в4
8
8
Условия порядка ФНРК сложнее, чем для НРК (они приведены в п. 4), и построение метода четвертого непрерывного порядка требует, как минимум, семи этапов [5, 15]. При этом этапные интерполянты, очевидно, не могут иметь тот же порядок, что и результирующий интерполянт, построенный по полному числу этапов. Тем не менее, как и в случае НРК, для сохранения дискретного глобального порядка р достаточно обеспечить локальный равномерный порядок p — 1 и локальный дискретный — р. В [15], в частности, представлен метод дискретного порядка четыре с интерполянтом порядка три, требующий всего шести этапов.
1
3
3
2
3
3
1
0
Это позволяет реализовать идею повторного использования и построить метод четвертого равномерного порядка с сохранением последнего этапа с прошлого шага и лишь шестью новыми вычислениями правой части на новом шаге [20].
Далее приводятся условия четвертого порядка ФНРК и доказывается утверждение о минимальном необходимом для него числе этапов.
4. Условия порядка ФНРК. Покажем, что условия четвертого равномерного порядка ФНРК, выведенные в [5], в случае отсутствия перекрытия на четвертом этапе можно разрешить с шестью этапами вместо семи для общего случая.
Согласно [5], явный ФНРК (8)—(10) обладает четвертым локальным равномерным порядком, если его коэффициенты удовлетворяют следующим равенствам:
¿2ъг(
1=1 в
-а, а € [0, 1 ац(а) — а, а € [0, с
э=1
в а2
ЕМ«)С» = —, а € [0, 1
¿=1
¿1
Е Ъ,(а) ( ]Г ау(/%-£)=(>, а € [0, 1
¿=1 сг = Ст
\3=1
¿1
^¿(а)с3 —
в3
€ [0,1
€ [0, 1
Е ь<(<*) Е - т = °> а 6 [°>1
г — ..., s,
в € [0,стЪ т —1,...,в*
в € [0,0^], т — 1,...,в*
(12)
(13)
(14)
(15)
(16)
(17)
(18)
\3 = 1
¿-1
(3-1
Е о, ае[0,1],/?е[0,<4], 7е[0,с?], (19)
¿=1
Сг = Ст
3=1
\к=1
для т,1 — 1,..., в*. Здесь о*, г — 1,..., в*, — все различные значения коэффициентов o¿. Условия четвертого дискретного порядка такие же, кроме (17), которое достаточно удовлетворить при а — 1.
Условия (12), (14), (15) и (17) дают четвертый равномерный порядок для случая чистой квадратурной формулы. Выполнение условия (13) дает первый равномерный порядок для каждого этапа. Условия (16) и (18) означают, что равномерный порядок приближения на этапах, используемых в финальной сумме (8) (т. е. на тех, для которых Ъ¿(а) ф 0), равен двум и трем соответственно. Последнее равенство означает, что приближения для которых должен быть обеспечен третий этапный порядок, не должны быть основаны на приближениях К, не обладающих вторым порядком (т. е. a¿j(в) ф 0, если для соответствующих г и ] не выполняется условие в скобках).
3
а
3
4
а
4
С =С
и
С =С
В [5] доказывается, что ФНРК четвертого равномерного порядка может быть построен как минимум за семь этапов, в то время как четвертый дискретный порядок достижим за шесть. Однако в том случае, если мы применяем метод к дифференциальному уравнению с дискретным запаздыванием (1), можно сформулировать следующее утверждение.
Теорема 2. При применении ФНРК (8) —(10) к ДУЗА (1) в том случае, если нет перекрытия на четвертом этапе, можно построить метод четвертого равномерного порядка с шестью этапами.
Доказательство. Наиболее общим случаем является ситуация, когда все о^, г — 1,..., 6, различны. Действительно, условия (16), (18) и (19) распадаются в этом случае на независимые равенства. Если выполняются они, очевидно, будут выполняться и равенства, включающие в себя их суммы, которые возникают при совпадении некоторых значений о¿ друг с другом. Поэтому будем считать, что все о¿ различны, в* — 6, а набор о* — о¿, г — 1,..., 6.
Поскольку рассматриваем явный метод, коэффициент 01 — 0. Таким образом, коэффициенты a¿l, г — 1,..., 6, не входят ни в одно уравнение, кроме группы условий (13), которая может быть решена за счет них при любых значениях других коэффициентов a¿j, г — 1,..., 6, ] — 2,..., г — 1. В частности,
а21(в)— в, в € [0,02].
Очевидно, что условие (16) при т — 2 влечет
Ъ2(а) ф 0, а € [0,1],
так как в скобках пропадает сумма. Совместное рассмотрение условий (16) и (18) при т — 3 дает к тому же
Ъ3(а) ф 0, а € [0,1],
так как одновременное выполнение условий аз2(в) — в2/202 и аз2(в) — в3/2о2 при произвольном в невозможно.
Условия (12), (14), (15) и (17) разрешимы, если 01 — 0, 04, 05 и 06 различны и Ъ1(а), Ъ4(а), Ъ5^(а) и Ъб(а) не нулевые.
Рассмотрим оставшиеся условия групп (16), (18) и (19) при т — 4, 5, 6. В силу отсутствия перекрытия на четвертом этапе уравнения с т — 4 должны быть удовлетворены только при значении в — 04, чтобы нужный порядок имело приближение 14. Приближение П4 при этом будет находиться по уже совершенным шагам (или предыстории) и потому будет иметь необходимый порядок точности.
Заметим, что условия (19) при т — 4, 5, 6 и I — 2 могут быть выполнены только при условии
а42 (04)—0, а52 (в) ф аб2(в) ф 0.
С учетом вышесказанного условия (16) примут вид
02
в4з(с4)сз =
2
а5з(/3)сз + а54(/3)с4 = у, /3 € [0, с5], (16')
2
абз(/?)сз + а64(/3)с4 + а65((3)с5 = —, (3 е [0,с6],
условия (18) —
е3
< \ 2 е4 а43(с4)с3 = —,
а5з(в)е! + а54(в)е24 =
3
в € [0,е5],
в3
абз(/?)с3 + а64(/3)с4 + а65(/3)с5 = —, /3 € [0, с6]
а условия (19)
а4з(е4) аз2(7)е2 -
а5з(в) аз2(7)с2 -
авз(в) аз2(7)е2 -
а54(в) а4з(е4)ез -
аб4(в) а4з(е4)ез -
= 0, 0, 0,
= 0, = 0,
1 € [0,ез], в € [0,е5], 7 € [0,ез], в € [0,ее], 7 € [0,ез], в € [0,С5], в € [0,ее],
(18')
(19')
а&ъ{Р) ^а53(7)с3 + а54(7)с4 - = 0, /?е[0,с6], 7б[0,с5]. Одновременно все условия (16'), (18') и (19') выполняются, если выполнены усло-
а32(7)с2 = у, е2
а4з(с4)с3 =
е3
2 е4
а43(с4)с3 = а5з(в)ез + а54(в)е4 = а5з(в)е3 + а54(в)е4 =
7 € [0,ез],
2 ' р
3 '
в € [0,С5], в € [0,е5],
в2
(20)
вбз(/3)сз +а64(/3)с4 + а65(/3)с5 = у, /3 £ [0,с6],
в3
абз(/?)с| +а64(/3)с4 + а,65(/3)с1 = у, Р € [0,с6].
Очевидно, что линейные по а^ условия (20) могут быть разрешены, если ез = С4. Таким образом, теорема доказана.
Замечание 2. Доказательство теоремы дает конструктивный алгоритм построения метода ФНРК для задач без перекрытия на четвертом этапе. Такие методы можно построить и с некоторыми совпадающими абсциссами, главное, чтобы выполнялось ез = С4 ив наборе е\, е4, е5, еб не было совпадений.
2
7
7
2
2
7
2
2
е
4
2
2
е
4
2
2
2
5. Комбинированный метод. В силу того, что трудоемкость применения методов типа Рунге—Кутты определяется числом этапов, которое определяет количество вычислений функции ] для явных методов и размерность систем алгебраических уравнений для неявных, естественно желание минимизировать число этапов метода при заданном порядке. Как было отмечено выше, явный НРК четвертого равномерного порядка при применении FSAL требует пяти новых этапов на шаге, а ФНРК — шести. При этом следует отметить, что в задачах с дискретным запаздыванием достаточно часто перекрытие случается лишь на некотором числе шагов, когда запаздывание становится достаточно малым (например, в уравнении пантографа лишь на нескольких первых шагах возникает перекрытие; при исчезающем запаздывании перекрытие появляется вблизи тех точек, где запаздывание равно нулю), большая же часть шагов не требует этапной интерполяции. В то же время потеря локального порядка на шагах с перекрытием, которая неизбежна при использовании непрерывного метода без этапной интерполяции, может привести к недопустимой потере глобальной точности.
Таким образом, лучшим решением будет применение двух методов: НРК с пятью этапами для шагов без перекрытия и ФНРК с шестью этапами для шагов с перекрытием. Очевидная проблема заключается в том, что обычно определить, появится перекрытие на новом шаге или нет, заранее невозможно. А в силу того, что на пограничных шагах перекрытие возникает лишь для части этапов (например, легко представить это в случае монотонного запаздывания), часть совершенных НРК вычислений может быть потеряна при обнаружении запаздывания на поздних этапах.
Однако, если обеспечить совпадение некоторых этапов НРК и ФНРК, то при обнаружении перекрытия на более поздних этапах станет возможным использовать уже выполненные вычисления.
Приведем коэффициенты метода Рунге—Кутты четвeртого равномерного порядка с шестью этапами и FSAL, в котором все этапы, кроме четвертого, снабжены ин-терполянтом и позволяют обрабатывать случай перекрытия (обозначим его cFCRK-I, от комбинированного ФНРК):
0 2 5 0 в
16 51 в--в2 4
8 2 0 6
17 17 17
19 20 ,_ 85 289^ 32 128 0 15302 32 - 867,3 128 17 „2 289 3 --в2 Н-- 8 64
1 а61(в) 0 0 аб4(в) Ябб(в) 0
(в) 0 0 &4 (в) Ьь (в) Ьб(в)
(ал а 483 85 3 (ал 5491 2 1445 3
аъ1\0) = 0--О Ч--0 , ацл(О) = -0--0 ,
п ' 304 114 ' 1 ' 2608 978 '
п ' 304 456 152 '
(ал 1600 2 6800 з
3097
9291
3097
9291
3097
23472
11736
11736
и коэффициенты полного ФНРК равномерного порядка четыре с семью этапами и FSAL ^СБК-П):
0 2
5
16
51
_8_
17
_8_
17
19
20 1
в--в2 4
в- -в2 4
,_ 85 289,3 32 128
32 128
0.71(6)
62
3
0 —в2 - — в3 32 128
0
32 128
0
3
0
17 Л2 289 Л 3
■—в Н--в3
8 64
0
17 „2 289 Л 3
--в2 Н--О3
8 64
075 (6) 076(6) 0
0
Ь1 (6) 0 0 0 65(6) Ьб (6) 67(6)
483 „2 , 85 лз _ ,пх 5491 2
0.71(0) =0--Г +-в, а7Ъ(в) = -Г--в , атв(в) =--Г +
304 114 ' 2608 978 ' 3097 9291
304 456 152 ' к ' 3097 9291 3097 '
Ыв) = 9-^о2 - 6-^в3 + Ъ7(в) = - ^в3 +
к ' 23472 11736 11736 ' п ' 9 9 9
Метод cFCRK-I построен согласно алгоритму из доказательства теоремы 2. Подробности построения cFCRK-II, который удовлетворяет полным условиям (12)—(19), можно найти в [5, 20].
В случае, если на четвертом этапе появляется перекрытие, следует применять cFCRK-II, если же нет — cFCRK-I. При этом в силу того, что первые три этапа этих методов совпадают, переход с cFCRK-I на cFCRK-II происходит без потери уже полученных результатов.
Фактически, метод cFCRK-II является тем же cFCRK-I с дополнительным четвертым этапом, который используется для построения на пятом этапе приближения необходимого порядка, в то время как в отсутствие перекрытия такое приближение получается уже на четвертом этапе. Таким образом, оба метода имеют очень близкие значения локальной погрешности и могут применяться совместно.
Дополнительно отметим, что интерполянт третьего порядка, использующийся на предпоследнем этапе обоих методов, может быть продлен до конца шага и применяться для вычисления приближения вложенным методом третьего порядка, что позволяет получить оценку локальной погрешности, на основе которой можно управлять длиной шага подобно тому, как это делается в методах Рунге—Кутты для ОДУ.
6. Численное тестирование. Для численного тестирования порядка сходимости метода и демонстрации необходимости использования этапных интерполянтов, а также для сравнения вычислительных затрат предлагаемого комбинированного метода с семиэтапным методом из [5] рассмотрим две задачи.
Задача 1. Она представляет собой начальную задачу (запаздывающий аргумент никогда не становится отрицательным), в которой в нулевой момент
времени нет запаздывания, и перекрытие возникает на нескольких первых шагах (пример 1.2.6 из [21]):
( t \(1+2i)2 ù(t) = и ( --— ) , t > 0,
V ; V(l+2tf) ' " '
u(0) = 1.
Аналитическим решением является функция
u(t) = e, t > о.
Будем решать задачу с постоянным шагом на отрезке 0 ^ t ^ 3. Число шагов обозначим через N. Покажем сравнение для трех методов: НРК (7) (обозначим его OZ), семиэтапного ФНРК (11) (MTV) и метода cFCRK.
В табл. 1, 2 через max \\E(t)|| обозначена максимальная глобальная погрешность численного решения на промежутке t G [0, 3]; p — полученный порядок сходимости, т. е. двоичный логарифм отношения погрешностей на вдвое различающихся длинах шагов (в пределе эта величина должна равняться теоретическому порядку сходимости метода); Nf — число обращений к функции правой части f, определяющее итоговую трудоемкость метода.
Таблица 1. Результаты решения Задачи 1 методом OZ в зависимости от использования экстраполяции
N Без экстраполяции С экстраполяцией
max ||E(i)II P max ||E(i)II Р
8 1.452374508 3.813702286 • 10~4
16 1.394858759 0.06 1.317968198 • 10~Б 4.85
32 1.349344295 0.05 2.481964536 • Ю-7 5.73
64 0.654110042 1.04 4.468434511 • Ю-9 5.80
128 0.290127444 1.17 8.613127989 • 10~10 2.38
256 0.104284152 1.48 4.323652547 • 10~п 4.32
512 0.036355681 1.52 2.692956969 • 10~12 4.00
1024 0.012659970 1.52 1.847411113 • 10~13 3.86
2048 0.004452930 1.51 2.842170943 • Ю"14 2.70
Таблица 2. Результаты решения Задачи 1 методами MTV и cFCRK
N MTV cFCRK
max ||E(i)II P Nf max ||E(i)Il P Nf
8 1.125584924 îo- -4 57 4.652127631 io- -3 42
16 1.565238520 10~ -b 2.84 113 6.052372897 10- -b 6.26 82
32 2.412819331 10- -ö 2.69 225 4.762033306 10- -6 3.67 162
64 3.398839610 10- -Y 2.82 449 5.764573281 10- -Y 3.04 323
128 2.005178956 10- -Я 4.08 897 2.203978511 10- -8 4.71 643
256 9.357421504 10- lu 4.42 1793 9.029577086 10- 1U 4.61 1284
512 3.989697461 10- 11 4.55 3585 3.499778245 10- 11 4.69 2566
1024 1.421085472 10- 12 4.81 7169 1.140421091 10- 12 4.94 5128
2048 1.953992523 КГ 14 6.18 14337 1.776356839 10~ 14 6.00 10250
Замечание 3. Непрерывный метод OZ можно применить к решению задач с перекрытием, если приближать решение внутри совершаемого шага посредством экстраполяции решения с предыдущего шага (или истории на первом шаге). Для рассмотренных задач метод демонстрирует сходимость с порядком, близким к четырем
(в пределе, пока вычислительная погрешность не оказывает влияния, как в последних строках табл. 1). Однако следует понимать, что если началом шага будет являться точка разрыва некоторой производной, то экстраполяция с предыдущего шага будет выполняться недостаточно точно, и порядок все же будет потерян. Это особенно характерно для начальной точки, в которой история почти всегда негладко соединяется с решением. В случае Задачи 1, в которой история не задана вовсе, мы взяли в качестве истории точное решение и использовали его значения для вычисления запаздывающего решения на первом шаге — очевидно, практически неприменимый метод, тем не менее на тестовом примере позволяющий продемонстрировать недостаточность приближения решения внутри текущего шага начальным значением на нем (случай «Без экстраполяции» в табл. 1).
Видно, что метод cFCRK требует существенно меньше вычислений правой части, чем MTV, при практически сравнимой величине глобальной погрешности. Метод OZ без экстраполяции демонстрирует порядок сходимости полтора, а не четыре. Порядок полтора возникает из-за того, что на шагах с перекрытием он падает до единицы, но число таких шагов среди всех мало, а именно, оно пропорционально 1 /л/N. Это наблюдение было описано в [22] и более детально объяснено Э. Хамфрисом на конференции «Recent trends in delay differential equations: models, theory and numerics», проходившей в Кортоне в Италии в 2012 г. *) (к сожалению, полноценная публикация до сих пор не вышла).
Задача 2. Запаздывание является периодически исчезающей функцией на всем промежутке решения:
u(t) = -u(a(t))u(t)ea(t), t > 0, u(t) = e-t, t < 0,
где a(t) = t — тщ cos(1007rf)2. Аналитическое решение u(t) = t ^ 0, является
гладким продолжением предыстории.
Выберем промежуток решения 0 ^ t ^ 0.5 и проведем те же расчеты, что и для Задачи 1. Полученные данные (табл. 3, 4) качественно совпадают с результатами для Задачи 1.
Таблица 3. Результаты решения Задачи 2 методом OZ в зависимости от использования экстраполяции
N Без экстраполяции С экстраполяцией
max ||E(i)II Р max ||E(i)|| Р
1 0.08268457982 1.017116087 • 10~4
2 0.03251547193 1.34 1.535256994 • 10~Б 2.73
4 0.01480079960 1.13 3.909773693 • 10~7 5.30
8 0.00657970202 1.16 1.808724093 • 10~8 4.43
16 0.00266893155 1.30 8.185840894 • КГ10 4.47
32 0.00097382734 1.45 4.219324889 • КГ11 4.28
64 0.00030453601 1.67 2.441713498 • КГ12 4.11
128 0.00010319090 1.56 1.482147738 • КГ13 4.04
256 0.00003579997 1.53 9.103828802 • КГ16 4.03
*) Humphries A. R. Singly diagonally implicit Runge—Kutta methods for state-dependent DDEs with overlapping // Recent trends in delay differential equations: models, theory and numerics. June 4—8, 2012. Cortona, Italy, 2012. Unpublished.
Таблица 4- Результаты решения Задачи 2 методами MTV и cFCRK
N MTV cFCRK
max ||E(i)II p Nr max ||E(i)II P Nr
1 2.662228310 10- -3 8 8.446918382 io- -4 7
2 1.229983383 10- -4 4.43 15 3.224687468 10- -b 4.71 13
4 6.158101201 10- -ö 4.32 29 1.446756357 10- -ö 4.48 25
8 2.974281966 10- -7 4.37 57 5.825843386 10- -8 4.63 49
16 1.183688003 10- -8 4.65 113 2.143614064 10- -У 4.76 97
32 4.980709178 10- 1U 4.57 225 9.249112587 10- 11 4.53 183
64 2.014843847 10- 11 4.62 449 3.962274953 10- 12 4.54 347
128 1.050270981 10- 12 4.26 897 1.965094754 10- 13 4.33 677
256 5.828670879 10- 14 4.17 1793 1.065814104 io- 14 4.20 1331
7. Заключение. Предложен метод, сочетающий в себе непрерывный метод Рунге—Кутты с пятью новыми этапами и функционально-непрерывный метод с шестью новыми этапами на шаге. Используемое число этапов — минимальное для случаев без и с перекрытием соответственно. Кроме того, за счет вложенности двух предложенных схем все расчеты, выполненные для непрерывного метода, используются и функционально-непрерывным, переход на который происходит в случае обнаружения перекрытия на любом из этапов. Это позволяет методу на шагах без перекрытия быть таким же эффективным, как и существующие непрерывные методы, а на шагах с перекрытием — полностью обеспечивать заявленный четвертый порядок точности.
Отметим, что непрерывное приближение предпоследнего этапа (у каждого метода) может быть использовано для оценки локальной погрешности и управления длиной шага, но данный вопрос требует дополнительной проработки и его рассмотрение выходит за рамки настоящей статьи.
Литература
1. Erneux Th. Applied delay differential equations. New York: Springer Science+Business Media, LLC, 2009. 204 p.
2. Smith H. An introduction to delay differential equations with applications to the life sciences. New York: Springer Science+Business Media, LLC, 2011. 172 p.
3. Bellen A., Zennaro M. Numerical methods for delay differential eqautions. Oxford: Oxford University Press, 2013. 396 p.
4. Guglielmi N., Hairer E. Implementing Radau IIA methods for stiff delay differential equations // Computing. 2001. Vol. 67. N 1. P. 1-12.
5. Maset S., Torelli L., Vermiglio R. Runge—Kutta methods for retarded functional differential equations // Mathematical Models and Methods in Applied Sciences. 2005. Vol. 15. N 8. P. 1203-1251.
6. Hairer E., N0rsett S. P., Wanner G. Solving ordinary differential equations I: Nonstiff problems, 2 ed., 3rd corr. print. Berlin; Heidelberg: Springer-Verlag, 2008. 528 p.
7. Shampine L. F. Solving ODEs and DDEs with residual control // Applied Numerical Mathematics. 2005. Vol. 52. P. 113-127.
8. Owren B., Zennaro M. Derivation of efficient, continuous, explicit Runge—Kutta methods // SIAM Journal on Scientific and Statistical Computing. 1992. Vol. 13. P. 1488-1501.
9. Еремин А. С., Олемской И. В. Вложенный метод интегрирования систем структурно разделенных обыкновенных дифференциальных уравнений // Журн. вычисл. математики и матем. физики. 2010. Т. 50. №3. C. 434-448.
10. Олемской И. В., Коврижных Н. А. Семейство шестиэтапных методов шестого порядка // Вестн. С.-Петерб. ун-та. Прикладная математика. Информатика. Процессы управления. 2018. Т. 14. Вып. 3. С. 215-229.
11. Eremin A. S., Kovrizhnykh N. A. Continuous extensions for structural Runge—Kutta methods // Computational Science and Its Applications — ICCSA 2017 / eds. O. Gervasi et al. Lecture Notes in Computer Science. 2017. Vol. 10405. P. 363-378.
12. Eremin A. S., Olemskoy I. V. Functional continuous Runge—Kutta methods for special sys-
tems // Proceedings of the International Conference on Numerical Analysis and Applied Mathematics 2015 (ICNAAM-2015). AIP Conference Proceedings. 2016. Vol. 1738. Art. no. 100003.
13. Feldstein A., Neves K. W. High order methods for state-dependent delay differential equations with nonsmooth solutions // SIAM Journal of Numerical Analysis. 1984. Vol. 21. N 5. P. 844-863.
14. Eremin A. S., Humphries A. R. Efficient accurate non-iterative breaking point detection and computation for state-dependent delay differential equations // Proceedings of the International Conference on Numerical Analysis and Applied Mathematics 2014 (ICNAAM-2014). AIP Conference Proceedings. 2015. Vol. 1648. Art. no. 150006.
15. Tavernini L. One-step methods for the numerical solution of Volterra functional differential equations // SIAM Journal of Numerical Analysis. 1971. Vol. 8. N 4. P. 786-795.
16. Bellen A., Guglielmi N., Maset S., Zennaro M. Recent trends in the numerical solution of retarded functional differential equations // Acta Numerica. 2009. P. 1-110.
17. Miyatake Yu. An energy-preserving exponentially-fitted continuous stage Runge—Kutta method for Hamiltonian systems // BIT Numerical Mathematics. 2014. Vol. 54. P. 777-799.
18. Hale J. K. Theory of functional differential equations. New York: Springer-Verlag, 1977. 366 p.
19. Kim A.V. i-Smooth analysis: Theory and applications. Beverly, MA: Scrivener Publishing LLC; Hoboken, NJ: John Wiley & Sons Inc., 2015. 296 p.
20. Eremin A. S. Functional continuous Runge—Kutta methods with reuse // Applied Numerical Mathematics. 2019. Vol. 146. P. 165-181.
21. Paul C.A.H. A test set of functional differential equations. Tech. Report N 243. Manchester: Manchester Centre for Computational Mathematics, University of Manchester, Feb. 1994. 41 p.
22. Magpantay F. M. G. On the stability and numerical stability of a model state dependent delay differential equation: PhD thesis. Montreal, Quebec, Canada: McGill University, 2011. 189 p.
Статья поступила в редакцию 30 апреля 2019 г.
^атья принята к печати 7 ноября 2019 г.
Контактная информация:
Еремин Алексей Сергеевич — канд. физ.-мат. наук; a.eremin@spbu.ru
Combined functional continuous method for delay differential equations
A. S. Eremin
St. Petersburg State University, 7—9, Universitetskaya nab., St. Petersburg, 199034, Russian Federation
For citation: Eremin A. S. Combined functional continuous method for delay differential equations. Vestnik of Saint Petersburg University. Applied Mathematics. Computer Science. Control Processes, 2019, vol. 15, iss. 4, pp. 425-441. https://doi.org/10.21638/11702/spbul0.2019.402 (In Russian)
In the paper a combined numerical method for discrete delay differential equations is presented. The method is an embedded pair of type two explicit Runge—Kutta methods of order four: a continuous method with six stages and a stage-continuous method with seven stages. Their combination provides an effective solution of discrete delay differential equations. The combined method remains explicit for any values of the delay: for small values the stage-continuous scheme is used while for large delays a faster continuous scheme is applied. The scheme to use is chosen automatically based on whether the delay falls into the current step and a switch to the stage-continuous scheme can be made at any stage when required. The embedding of the methods lets to minimize the required number of the right-hand side function computations. The order conditions and the proof of their resolvability with the stated number of stages are presented. Tests, confirming the effectiveness of the proposed methods, are made.
Keywords: delay differential equations, continuous methods, functional continuous method, stage-continuous method.
References
1. Erneux Th. Applied delay differential equations. New York, Springer Science+Business Media, LLC Press, 2009, 204 p.
2. Smith H. An introduction to delay differential equations with applications to the life sciences. New York, Springer Science+Business Media, LLC Press, 2011, 172 p.
3. Bellen A., Zennaro M. Numerical methods for delay differential eqautions. Oxford, Oxford University Press, 2013, 396 p.
4. Guglielmi N., Hairer E. Implementing Radau IIA methods for stiff delay differential equations. Computing, 2001, vol. 67, no. 1, pp. 1—12.
5. Maset S., Torelli L., Vermiglio R. Runge—Kutta methods for retarded functional differential equations. Mathematical Models and Methods in Applied Sciences, 2005, vol. 15, no. 8, pp. 1203—1251.
6. Hairer E., N0rsett S. P., Wanner G. Solving ordinary differential equations I: Nonstiff problems, 2 ed., 3rd corr. print. Berlin, Heidelberg, Springer-Verlag Press, 2008, 528 p.
7. Shampine L. F. Solving ODEs and DDEs with residual control. Applied Numerical Mathematics, 2005, vol. 52, pp. 113-127.
8. Owren B., Zennaro M. Derivation of efficient, continuous, explicit Runge—Kutta methods. SIAM Journal on Scientific and Statistical Computing, 1992, vol. 13, pp. 1488-1501.
9. Eremin A. S., Olemskoy I. V. An embedded method for integrating systems of structurally separated ordinary differential equations. Computational Mathematics and Mathematical Physics, 2010, vol. 50, no. 3, pp. 414-427.
10. Olemskoy I. V., Kovrizhnykh N. A. Semeistvo shestietapnyh metodov shestogo poriadka [A family of sixth-order methods with six stages]. Vestnik of Saint Petersburg University. Applied Mathematics. Computer Science. Control Processes, 2018, vol. 14, no. 3, pp. 215-229. (In Russian)
11. Eremin A. S., Kovrizhnykh N. A. Continuous extensions for structural Runge—Kutta methods. Computational Science and its Applications — ICCSA 2017. Eds. O. Gervasi et al. Lecture Notes in Computer Science, 2017, vol. 10405, pp. 363-378.
12. Eremin A. S., Olemskoy I. V. Functional continuous Runge—Kutta methods for special systems. Proceedings of the International Conference on Numerical Analysis and Applied, Mathematics 2015 (ICNAAM-2015). AIP Conference Proceedings, 2016, vol. 1738, art. no. 100003.
13. Feldstein A., Neves K. W. High order methods for state-dependent delay differential equations with nonsmooth solutions. SIAM Journal of Numerical Analysis, 1984, vol. 21, no. 5, pp. 844-863.
14. Eremin A. S., Humphries A. R. Efficient accurate non-iterative breaking point detection and computation for state-dependent delay differential equations. Proceedings of the International Conference on Numerical Analysis and Applied Mathematics 2014 (ICNAAM-2014). AIP Conference Proceedings, 2015, vol. 1648, art. no. 150006.
15. Tavernini L. One-step methods for the numerical solution of Volterra functional differential equations. SIAM Journal of Numerical Analysis, 1971, vol. 8, no. 4, pp. 786-795.
16. Bellen A., Guglielmi N., Maset S., Zennaro M. Recent trends in the numerical solution of retarded functional differential equations. Acta Numerica, 2009, pp. 1-110.
17. Miyatake Yu. An energy-preserving exponentially-fitted continuous stage Runge—Kutta method for Hamiltonian systems. BIT Numerical Mathematics, 2014, vol. 54, pp. 777-799.
18. Hale J. K. Theory of functional differential equations. New York, Springer-Verlag Press, 1977, 366 p.
19. Kim A. V. i-Smooth analysis: Theory and applications. Beverly, MA, Scrivener Publishing LLC, Hoboken, NJ, John Wiley & Sons Inc. Press, 2015, 296 p.
20. Eremin A. S. Functional continuous Runge—Kutta methods with reuse. Applied Numerical Mathematics, 2019, vol. 146, pp. 165-181.
21. Paul C. A. H. A test set of functional differential equations. Tech. Report no. 243. Manchester: Manchester Centre for Computational Mathematics, University of Manchester Press, Feb. 1994, 41 p.
22. Magpantay F. M. G. On the stability and numerical stability of a model state dependent delay differential equation. PhD thesis. Montreal, Quebec, Canada, McGill University Press, 2011, 189 p.
Received: April 30, 2019.
Accepted: November 07, 2019.
Author's information:
Alexey S. Eremin — PhD in Physics and Mathematics; a.eremin@spbu.ru
