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

Вычислительные технологии Том 15, № 3, 2010

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

Н. Г. Бандурин Волгоградский государственный архитектурно-строительный университет, Россия e-mail: bandurin_ng@mail.ru

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

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

Введение

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

^х,У (х),у'(х), (х - т )У (х - т),... ,1 (ж)] = 0. (1)

Здесь

^ — вектор-столбец размером N1

(у(х))Т = (У1(х),У2(х),.. .ум(х));

а(х)Г = (/1(х),/2(х),...,);

х

!г(х) = ! фг[х,П,у(х),у'(х), . . . ,у(п),у'(п), . . . ,у(и - Т),...] ¿П. 0

Начальная вектор-функция у(х) = ф(х) при —т < х < 0 считается известной. Ставится задача найти решение системы (1) на равномерной сетке х1 = 0, х2,..., хм = Ь отрезка [0, Ь].

1. Вывод формул интегрирования и дифференцирования

[0, т]

ной оси х задана равномерная сетка, в т узлах которой известны значения производной © ИВТ СО РАН, 2010.

к-го порядка достаточно гладкой функции у = f (х), определенной на [0,6], и значения производных низших порядков в начальном узле х = х1. Ставится задача вычисления в этих же узлах значений функции, ее производных и интегралов, выраженных через узловые значения производной порядка к и производные при х = XI = 0.

Вначале интерполируем производную у(к)(х) на отрезке [х1;хп] (п < т) многочленом степени (п — 1) [1, 2] и результат запишем в виде

у(к)(х) = хт р,

где хт = (1,х,..., хп-1), р = (р1,р2, ■ ■ ■,рп)• Полагая в (2) х = х1; х = х2,..х получим матричное соотношение

у(к) (х) = Zp, (3)

где у(к) — матрица-столбец узловых значений производной к-го порядка, % = х^-1. Так как определителем матрицы Z является определитель Вандермонда, то можно

(2)

хп

записать

^-1у(к)

(2)

У(к)(х)

-1у(к).

(4)

(5)

[х1, х1]

[х1, х2], ■ ■ ■ , [х1 ,хп], получим

у(к-1) = 1у(к-1) + BZ-1у(к). (6)

Здесь у(к-1) — вектор-столбец узловых значений производной (к — 1)-го порядка; I — единичный вектор размером п; элементы матрицы В Ь^ = (х^ — х'1)/з. (6)

у(к-1) = 1у(1к-1) + яу(к),

(7)

где в = BZ 1 — матрица интегрирования порядка п на отрезке [х1,хп].

Если т > п, то, применяя последовательно формулу (5) на отрезках [х^(п-1)+1, хлп-1)+п] (] = 0,1, 2,...), можно получить компоненты вектора производной У(к-1) во всех узлах отрезка [0, т]

Y(k-1) = 1у(к-1) + Ж(к), (8)

где I — единичный вектор, S — квадратная матрица интегрирования порядка т. Например, при Ь = 6, п = 3, т = 7 эта матрица имеет вид

5

0.0000 0.0000 0.0000 0.0000 0.0000 0.0000 0.0000

0.4167 0.6667 —0.0833 0.0000 0.0000 0.0000 0.0000

0.3333 1.3333 0.3333 0.0000 0.0000 0.0000 0.0000

0.3333 1.3333 0.7500 0.6667 —0, 0833 0.0000 0.0000

0.3333 0.3333 0.6667 1.3333 0.3333 0.0000 0.0000

0.3333 1.3333 0.6667 1.3333 0.7500 0.6667 —0.0833

0.3333 1.3333 0.6667 1.3333 0.6667 1.3333 0.3333

Выполнив интегрирование в соответствии с (8) г раз, получим формулы для вычис-лбния производных в узлах интерполяции^ выраженных через значения производных

порядка к в этих же узлах и значения производных низших порядков в начальном узле

г—1

у(к—г) = ^ 8Чу(к—г+г) + угуг < к.

(9)

i=0

Эту формулу можно использовать для решения обыкновенных ИДУ в случае, когда за неизвестные дискретной модели приняты значения производной наивысшего порядка в узлах интерполяции и значения производных низших порядков в узле х = х1.

Если г > к, то компонентами век тора У(к—г) в (9) будут значения интегралов, а формула примет вид

г-1

у(к—г) ''у ^ уг—к+^уМ + у г у (к)

(10)

i=г—к

В последнем равенстве учтено, что у(к г) = 0 при г > к, поскольку в соответствии с принятыми обозначениями эти величины являются значениями интегралов в узле

Для получения формулы дифференцирования матрицу интегрирования представим в виде

511 512 521 522

У

(11)

где 511 = 0 элементы матрицы-строки в12 равны нулю, а квадратная матрица з22 — не особенная. Теперь равенство (8) примет вид

у1 2

у1 2у1

+

511 в 12 521 522

у1 2

(12)

Здесь у2 , у2 — векторы-столбцы узловых значений функции и ее производных размером (т — 1) соответствен но, 12 — единичный вектор размером (т — 1). Из (12) получим

у1 2

1

— 5

—1 22

0

521

1у1 +

00

-«—1т «—1 522 2 522

или в принятых ранее обозначениях

Y' =(Е — ББ )!у[ + DY,

00

-« —:!Т « —1 522 2 522

матрица дифференцирования имеет вид

где Б

у1 2

матрица дифференцирования. При Ь = 6, п

(13) 3, т = 7

Б

0 0 0 0 0 0 0

-1.25 1 0.25 0 0 0 0

2 —4 2 0 0 0 0

-1 2 -2.25 1 0.25 0 0

2 -4 4 -4 2 0 0

-1 2 -2 2 — 2.25 1 0.25

2 4 4 4 4 4 2

Выполнив дифференцирование в соответствии с (13) г раз, получим формулу для вычисления производных г-го порядка в узлах интерполяции, выраженных через производные в начальной точке х = х1 и значения функции в этих узлах:

г—1

уМ = ^ Бг(Е - ББ)1у(Гг) + . (14)

г=0

Формула (14) применима для решения как обыкновенных ИДУ, так и уравнений с частными производными.

Для любого целого г > 0 имеют место равенства

N

= (хг - Х1)Г

3 = 1

N

= 0, (г =1,

3 = 1

которые выражают результаты интегрирования и дифференцирования функции

I (х) = 1-

Использование формул (9) и (14) [1, 2] при численном решении краевых задач для систем интегро-дифференциальных уравнений позволяет:

1) с высокой точностью аппроксимировать системы ИДУ, имеющие достаточно гладкие решения, их дискретными аналогами, так как формулы являются точными на отрезке [0, т] для всех многочленов степени не выше щ

2) решать уравнения с заданными значениями производных не только на краях, но и внутри области интегрирования, поскольку при решении краевых задач имеется возможность за начальную точку на отрезке [0, Ь] принимать любой узел хг (г =1, 2,..., ш);

3) при решении ИДУ высокого порядка с использованием интерполирования многочленами исключить необходимость выполнения алгебраических преобразований на краях области, как это рекомендуется, например, в [3], так как формулы (9) и (14) не зависят от номеров узлов, расположенных вне области интегрирования, и содержат в явном виде необходимые для постановки граничных условий производные. Используя формулу Тейлора, можно получить оценку погрешности аппроксимации функции и ее производных^ равную 0(кга+1), где к — шаг сетки.

В литературе предложено много методов численного решения систем уравнений (1)

имеют невысокий порядок точности. В работе [8] задача решается с позиции метода продолжения по наилучшему параметру, а в [9] используется способ преобразования к наилучшему аргументу.

(1)

ствии с которым весь промежуток интегрирования [0,Ь] должен быть представлен в виде ^ ^^^^^^^ ^агов длиной т, т. е. принимается Ь = Qт. Для численной реализа-

ш - 1

узловых значений функций и их производных низшего порядка через производные наи-

(9)

Ниже приводится алгоритм итерационного процесса для вычисления узловых значений функции у = I (х), являющейся решением одного нелинейного ИДУ к-то порядка

с постоянной величиной запаздывания аргумента т на отрезке [0,фт], которое имеет в своем составе один интеграл I(х) с переменным верхним пределом

/ [х,у(х),у' (х),...,у(к)(х),у(х — т),..., у(к)(х — т ),1 (х)] = 0, (15)

х

I (х) = У ф[х,п,у(п),у' (и),... ,у(к\и),у(и — т),..., у(кк) (и — т)] ¿и. 0

х=0

у(0) — Со = 0, у' (0) = С1,.. .,у(к—1)(0) = Ск—1. (16)

Числа ^ могут быть определены как производиые функции ф(х) при х = 0, но также могут быть заданы независимо от этой функции [11].

Рассмотрим алгоритм получения решения на первом шаге, т.е. на отрезке [0,т]. В соответствии с методом шагов [10, 11] на первом шаге значения функции с запаздыванием у(х — т) и ее производных определяются из выражения начальной функции ф(х). На последующих шагах эти значения выбираются из результатов решения на предыдущем шаге. Последовательно подставляя в (15) х1 = 0, х2 = т/(т — 1),..., хт = т, можно

ту Дискретизация интегралов выполняется с использованием формулы (10). Эти уравнения и начальные условия (16) можно записать в виде системы нелинейных уравнений относительно узловых значений функции и ее производных

Е(со,С1,...,Ск—1,У,У,...,У(к)) = 0, (17)

где = , у2^),..., у™^ , 5 = 0,1, 2,..., к. В систему (17) в качестве известных

величин входят узловые значения начальной функции ф(х).

Исключая узловые значения функции и ее производных низших порядков с помо-

(9)

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

(к)

Ф(со,С1,...,Ск—1,У (к)) = 0, (18)

которая решается с высокой точностью известными итерационными методами [12, 13].

(18)

величины ИехТо1. При известном век торе У(к) производные низших порядков на шаге вычисляются по формуле интегрирования (9) Второй шаг при т < х < 2т выполня-бтся при известных узловых значениях функции и ее производных, вычисленных на первом шаге. Очевидно, что при условии независимости подынтегрального выражения х

зультат суммируется со значением интеграла, полученным в результате выполнения

х

узла х1 текущего шага интегрирование должно быть выполнено по отрезку [0,х^, что многократно увеличивает время счета.

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

точность вычисления искомых функций. Следует учитывать, что с увеличением числа узлов m нарастает погрешность округления чисел. Заметим, что решение систем ЛИН6ИНЫХ ИДУ получается в результате выполнения одной итерации. Эффективность применения разработанной программы иллюстрируется приведенными ниже результатами решения трех тестовых примеров при IterTol = 10-12.

2. Тестовые примеры

Пример 1. Существенно нелинейное интегро-дифференциальное уравнение второго порядка с запаздыванием аргумента т = п/2

y''(x) + sin [y''(x)] + y3(x) — J [y(s — т) + y'(s)] ds + exp

0

— y'(s — т)y''(s — т) ds

— exp[0.25(cos 2x — 1)] + sin(y'(x — т)) + sin x cos2 x = 0.

Известно точное решение y(x) = sin x. Поэтому начальная функция на отрезке [—п/2, 0] принимается в виде tp(x) = sinx.

y(0) = 0, y (0) — 1 = 0 Интегрирование выполнялось при различных m в области 0 < x < 10. В результате решения получены следующие значения максимального отклонения приближенного решения от точного:

Дм=7 = 7 • 10-3, Дм=ii = 3 • 10-7, Дм=з1 = 5 • 10-10, Дм=151 = 4 • 10"10.

Пример 2. Система двух существенно нелинейных интегро-дифференциальных уравнений с запаздыванием аргумента т = п/2

y(VI11 )(x)+exp

y1'(s — т)y2(s — т) ds

+ y2(x — т) — exp[0.25(1 — cos 2x)] = 0,

x

x

x

y2 (x) + [y(VI11 )(x^ 3 — y1(x — т) — cos3 x = 0.

Точное решение: y1(x) = cos x, y2(x) = sin x. ^^^^^^^ условия для функции y1(x) при x=0 очевидны и здесь не приводятся.

Значения отклонения Д при b = 3п и разных значениях m для функции y1 (x)

Дт=9 = 1.5 • 10-5, Дт=11 = 6.5 • 10-7, Дт=21 = 1.7 • 10-4, Дт=з1 = 1.0 • 10-3,

для функции y2(x)

Дт=9 = 2.1 • 10-5, Дт=И = 1.4 • 10-7, Дт=21 = 4.6 • 10-5, Дт=31 = 4.0 • 10-4.

ЧиСЛбННОб р GTIJG HUG С VTIXG С ТВ GUUO НбЛШНбШНЪТХ ITUTG Рр О-ДШ(^) cjj Gp б ЛЦШЭ.ЛЪ НЫХ..

37

Пример 3. Система существенно нелинейных интегро-ди(^)(^)еренциальных уравне-иий с запаздыванием аргумента т = 1

y2'(x) + sin [y1(x — 1)] — 2exp(—x) sin x — sin[exp(x — 1) sin(x — 1) — 10x + 10] = 0.

Точное решение: y1(x) = exp(x) sinx — 10x, y2(x) = exp(—x) cos x.

x=0 очевидны и здесь не приводятся. Значения отклонения Д при b = 8 для функции y1 (x)

Дт=9 = 1.5 • IQ"5, Дт=15 = 7.5 • IQ"8, Дт=21 = 7.6 • 1Q"11, Дт=э1 = 3 • 1Q"9,

Дт=9 = 1.3 • 10"6, Дт=15 = 4.6 • 1Q"9, Дт=21 = 1.4 • 1Q"11, Дт=э1 = 3 • 1Q"10.

Из анализа результатов решения этих примеров можно заключить, что при некоторых значениях m достигается максимальная точность решения, в то время как при уменьшении m погрешность нарастает в связи с разрежением сетки, а при увеличении m

Список литературы

[1] Бандурин Н.Г. Новый численный метод порядка n для решения интегро-дифференциальных уравнений общего вида // Вычисл. технологии. 2002. Т. 7, № 2. С. 3-10.

[2] Бандурин Н.Г., Игнатьев В.А. Пакет программ для решения систем нелинейных интегро-дифференциальных уравнений (одно-, двух- и трехмерные начально-краевые задачи) // Мат. моделирование. 2007. Т. 19, № 2. С. 105-112.

[3] Корн Г., Корн Т. Справочник по математике. М.: Наука, 1968. 720 с.

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

[5] Трауб Дж. Итерационные методы решения уравнений. М.: Наука, 1985.

[6] Ортега Дж., Пул У. Введение в численные методы решения дифференциальных уравнений. М.: Наука, 1986.

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

[8] Дмитрие в С.С., Кузн ецов Е.Б. Численное решение систем интегро-дифференциально-алгебраических уравнений с запаздывающим аргументом // Ж урн. вычисл. математики и мат. физики. 2008. Т. 48, № 3. С. 430-444.

[9] Кузнецов Е.Б., Микрюков В.Н. Численное интегрирование системы интегро-диффе-ренциально-алгебраических уравнений с запаздывающим аргументом // Там же. 2007.

x

0

2 exp(x) cos x — exp[cos 2(1 — cos 2x)] — sin 2 sin 2x = Q;

для функции y2(x)

Т. 47, № 1. С. 83-95.

[10] Каменский Г.А. Общая теория уравнений с отклоняющимся аргументом // Докл. АН СССР. 1958. Т. 120, № 4. С. 697-700.

[11] Эльсгольц Л.Э., Норкин С.Б. Введение в теорию дифференциальных уравнений с отклоняющимся аргументом. М.: Наука, 1971. 296 с.

[12] Ортегд Дж., РеЙНБОЛДТ В. Итерационные методы решения нелинейных систем уравнений со многими неизвестными. М.: Мир, 1975.

[13] Денис Дж., Шнабель Р. Численные методы безусловной оптимизации и решения нелинейных уравнений. М.: Мир, 1988.

Поступила в редакцию 4 июня 2009 г.