ЖУРНАЛ ВЫЧИСЛИТЕЛЬНОЙ МАТЕМАТИКИ И МАТЕМАТИЧЕСКОЙ ФИЗИКИ, 2014, том 54, № 3, с. 496-502
УДК 519.642.5
НЕКОТОРЫЕ ОСОБЕННОСТИ ПОВЕДЕНИЯ ЧИСЛЕННЫХ МЕТОДОВ РЕШЕНИЯ ИНТЕГРАЛЬНЫХ УРАВНЕНИЙ ВОЛЬТЕРРА II РОДА1)
© 2014 г. М. В. Булатов*, М. Н. Мачхина**
(* 664033 Иркутск, ул. Лермонтова, 134, ИДСТУ СО РАН; ул. Лермонтова, 83, НИ Ир ГТУ;
** 664011 Иркутск, ул. Нижняя Набережная, 6, ФГБОУВПО "ВСГАО") e-mail: mvbul@icc.ru; masha888888@mail.ru Поступила в редакцию 14.03.2013 г.
Переработанный вариант 27.08.2013 г.
Рассматриваются системы интегральных уравнений Вольтерра II рода, содержащие жесткие и осциллирующие компоненты. Для численного решения таких задач предложен неявный, безытерационный метод второго порядка. На ряде конкретных примеров продемонстрирована эффективность предложенного метода. Библ. 10. Фиг. 1. Табл. 5.
Ключевые слова: интегральные уравнения Вольтерра II рода, жесткие задачи, разностные схемы, квадратурные формулы.
DOI: 10.7868/S0044466914030028
1. ВВЕДЕНИЕ
Интегральные уравнения Вольтерра (ИУВ) II рода (ИУВП) часто встречаются в различных прикладных задачах. Качественной теории и численным методам решения таких уравнений по-свещены многочисленные работы (см., например, [1]—[5] и приведенную в них библиографию). В данной работе рассмотрен численный метод решения уравнения
t
x(t) = JK(t, s, x(s))ds +f(t), 0 < t < 1, 0 < s < t, (1)
0
где K(t, s, x(s)) и f(t) — известные, достаточно гладкие «-мерные вектор-функции, x(t) — искомая «-мерная вектор-функция.
Всюду в дальнейшем будем предполагать, что уравнение (1) имеет единственное непрерывное решение и входные данные обладают той степенью гладкости, которая необходима для дальнейших рассуждений.
К настоящему времени разработана достаточно полная теория численных методов решения задачи (1): одношаговые методы типа Рунге-Кутты и их модификации, многошаговые, общие многошаговые методы, блочные методы и др. Достаточно полная библиография по таким методам приведена в монографиях [3]—[5].
Формально все эти методы можно разделить на явные и неявные. При реализации явных методов на каждом шаге интегрирования приближенное решение находят из рекуррентных соотношений. Реализация неявных методов требует значительно больших вычислительных затрат: на каждом шаге интегрирования приходится решать системы нелинейных уравнений (или линейных, если K(t, s, x(s)) = K(t, s)x(s), где K(t, s) есть (n x п)-матрица) размерности ns x ns, где s — число стадий. Как правило, такие системы решают методом Ньютона, или его модификациями. При этом возникают следующие проблемы: выбор начального приближения для запуска итерационного процесса, выбор числа итераций и критерия остановки счета. В настоящей статье предложен алгоритм, при реализации которого на каждом шаге интегрирования требуется решать линейные уравнения.
1) Работа выполнена при финансовой поддержке РФФИ (коды проектов 11-01-00639а, 13-01-93002вьет-а).
2. ТЕСТОВОЕ УРАВНЕНИЕ
На протяжении многих десятилетий для предсказания свойств численных методов решения задачи (1) используют тестовое уравнение (см. [6]):
t
x(t) = J(X + ц(t - s))x(s)ds + a1 + a2t, (2)
0
где X и ц — вещественные числа и X < 0, ц < 0.
Уравнение (2) эквивалентно обыкновенному дифференциальному уравнению (ОДУ) второго порядка
x" (t) = Xx' (t) + цх (t), x (0) = a1, x' (t) |t = 0 = X a1 + a2, (3)
при ц Ф 0 и ОДУ первого порядка
x'( t) = Xx (t), x (0) = a1, (4)
при ц = 0, a2 = 0.
Введем на отрезке [0,1] равномерную сетку tl = lh, l = 0, 1, 2, ..., N, h = 1/N, и пустьf +1 = f(tl +1), xl +1 ~ x(tl +1),
2
z = Xh, r = цк .
Любой одношаговый метод для уравнения (4) можно представить в виде xi +1 = R(z)xi, где R(z) принято называть функцией устойчивости (см. [9]) (полином или дробно рациональная функция). Если |R(z)| < 1 при z < 0, то метод называется А-устойчивым, а если при этом lim R(z) = 0,
z
то метод называется L-устойчивым.
Остановимся на одном семействе одношаговых методов, а именно, на 0-методах, которые для задачи (1) имеют вид
i
x, + ! = h [0K(tl + ь 0, x0) + (1 - 0)K(tl + ь t, + ь x, +!)] + h £ K(t, + ь tp, xp) + ft +:, x„ = x(0) = f(0) ,(5)
p = 0
где 0 — скалярный параметр: 0 < 0 < 1.
При 0 = 0, 0 = 1, 0 = 1 получим методы, основанные на квадратурных формулах правых прямоугольников, левых прямоугольников и трапеций соответственно.
Для тестового примера (2) данные методы дают рекуррентные соотношения
xi +1 - 2x, + xi -1 = z(x, +1 - x,) + rx,, 0 = 1, (6)
xi +1 - 2x, + x,-1 = z(x,-x,-1) + rx,-1, 0 = 0, (7)
z1
x, +1 - 2x, + x,-1 = 2(x, +1 - x,-1) + rx,, 0 = 2, (8)
где z = Xh, r = ц^.
В дальнейших рассуждениях нам потребуется
Определение (см. [4]). Областью устойчивости методов (6)—(8) назовем те вещественные значения z и r, при которых корни р характеристических уравнений
(1 - z )р2 - 2 (1 - 2 + 2)р + 1 = 0, (9)
р2 - 2(1 + -)р + (1 + z- r) = 0, (10)
1 - ^р2 - 2(1 + 2)р + (1 + ^ = 0, (11)
лежат в единичном круге.
Если ц = 0 в уравнении (2), то 0-методы при 0 = 0, 0 = 1, 0 = 1 совпадают с неявным методом
Эйлера, явным методом Эйлера и методом трапеций для ОДУ (4). Функция устойчивости в этом случае
1 + 0 г
* (г) =
1 - (1 - 0) г'
т.е. 0-методы будут ^-устойчивыми при 0 < 0 < 1 и Х-устойчивым при 0 = 0.
Таким образом, для численного решения задачи (1), содержащей жесткие и быстроосцилли-рующие компоненты, явные методы не эффективны, так как требуют "малого" шага интегрирования. Как уже говорилось, реализация неявных методов требует решения систем нелинейных уравнений, что влечет за собой ряд дополнительных проблем: выбор подходящего итерационного метода, начального приближения и критерия остановки итерационного процесса.
Ниже предложен безытерационный алгоритм, который объединяет "хорошие" качества метода, основанного на квадратурной формуле правых прямоугольников (Х-устойчивость) и метода, основанного на квадратурной формуле трапеции (второй порядок). Отметим, что для жестких ОДУ такие методы были предложены в [7], [8].
3. ЧИСЛЕННЫЕ АЛГОРИТМЫ
Выпишем для задачи (1) простейший метод, основанный на квадратурной формуле левых прямоугольников (0 = 1):
х, +1 = к £ К( г1+1, у х) +Л +1, I — 0,1,2.....N -1,
У = о
(12)
Хо = Л 0),
и линеаризованный неявный метод, основанный на квадратурной формуле правых прямоугольников (0 = 0):
( е - м(+1, +1,х)) х +1 = к
£ К( +1, у х) + К( +1, ^ +1, х) - /(^ +1, ^ +1, х) х +Л+1
У — 1
(13)
д К( г ^ х)
где / = —---1—-—-, I = 0, 1, 2, ..., N — 1, Е — единичная матрица порядка п. дх
Объединяя соотношения (12) и (13) в переопределенную систему линейных алгебраических уравнений (СЛАУ) относительно +1, получаем
(
Х1 +1 —
( е - ^(+1, +1, Х)) .
к £ К( +1, у ху)
у — о
Г Л ^ Л +1
( Л+1
£ К( +ь у ху) + К( +1, +1, х) - /(+ь +1, х) х
—1
(14)
+
-20 -15
10
Фигура. Область устойчивости заштрихована.
СЛАУ (14) в общем случае не имеет решения. Найдем ее "решение" следующим образом. Умножим обе части (14) слева на матрицу (А\Б) размерности п х 2п такую, чтобы матрица Л + Б(Е — М(11 +1, +1, х)) была неособенной. Если (Л\Б) = (Е\Е — М(11 +1, +1, х)), то имеем
(Е + (Е - а,)2)х, +! = ( 2Е - 0,)1 к £ К(и + ъ у х) +
+1
+
У = 1
(15)
+ к[К( +1, X) + (Е- 0,)(К( +1,+1, X,.) - 0;Х;)] , где 01 = М(Ь + 1, ^ + 1, х), I = 0, 1, 2, ..., N - 1.
Если (Л\Б) = (Е\Е), то получим линеаризованный метод трапеций.
Метод (15) и линеаризованный метод трапеций имеют второй порядок и являются безытерационными.
Для модельной задачи (2) при ц = 0 метод (15) дает рекуррентное соотношение
( \
1
Х, +1 -
V' - с +
X,
т.е. является (0, 2) паде-аппроксимацией экспоненты (см. [9]).
На фигуре показана область устойчивости метода (15) для примера (2).
4. ЧИСЛЕННЫЕ ЭКСПЕРИМЕНТЫ
Предлагаемые методы численного решения (15) были протестированы на ряде задач. Приведем результаты расчетов некоторых примеров.
г
0.1 0.05 0.025
егх 0.249 0.101 0.0440
ег2 0.0341 0.00932 0.00236
ег3 0.0228 0.00646 0.00137
Пример 1. Рассмотрим систему ИУВ11
' ( х( г) = Йо
а(г)и(г, s) + Ь (г) V(г, ^) Ь(г)й(г)(30 + 127(г - 5)) V с(г)Ь(г)( 10 + 75(г - 5)) Ь2(г)и(г, 5) + а(г) г, 5)
X ( 5 ) + /(г) ,
где
г) = 200 + 20(СО8200 + -2-0а), а(г) = ((10 + со8220о) -О,
Ь (г) = 10 + бШ ——, с(г) = 9 + соБ ——, й(г) = 11 + СОБ—, 200 200 200
и(г, 5) = - 20 -101 (г - 5), V(г, 5) = -10 - 2б(г - 5),
/(г) =
( - й (г)(г + 1) + Ь (г)Л
Точное решение:
х (г) =
V Ь(г)(г + 1) + с(г)
Ь( г) е 5г( соБ г - г) - й( г) е 10г ( соб г - г) V Ь(г)е ш(собг-9б1иг) + с(г)е 5г(собг-г) ^
В табл. 1 приведены евклидовы нормы погрешности в точке I = 1 (последняя точка отрезка интегрирования) с шагом интегрирования к.
Здесь и в дальнейших расчетах приняты следующие обозначения:
ег1 — норма погрешности в точке I = Т метода, основанного на квадратурной формуле правых прямоугольников с дальнейшей линеаризацией;
ег2 — норма погрешности в точке I = Т метода, основанного на квадратурной формуле трапеций с последующей линиаризацией;
ег3 — норма погрешности в точке I = Т метода (15).
Норма в расчетах евклидова, Т — правая граница отрезка интегрирования. Пример 2. Рассмотрим нелинейное ИУВ11 (см. [4, с. 127])
г
х( г) = ехр (-г) + |ехр (5-г)[х (5) + ехр (-х (5))] .
0
Точное решение: х(1) = 1п(1 + ехр(1)).
В табл. 2 приведенные результаты численных расчетов на отрезке интегрирования [0; 6]. Пример 3. Рассмотрим задачу (см. [4, с. 467])
г
х(г) = 1 + 1 А.( 1 - ехр(-г2))-А.|(г-5)ехр(-(г-5)2)х(5)
0
0
0.25 0.125 0.0625
егх 1.95 0.68 0.29
ег2 0.03 0.77 х 10-2 0.19 х 10-2
ег3 0.07 0.18 х 10-1 0.42 х 10-2
Таблица 3
0.1 0.05 0.025
ег^ * * 0.87
ег2 * 0.041 0.20
ег3 * 0.126 0.35
Таблица 4
Точное решение: х(1) = 1.
В табл. 3 приведены результаты для случая X = 1000, Т = 20.
Символ * означает, что норма погрешности в I = 20 больше 100.
В [4, с. 467] приведены результаты применения методов, основанных на явных квадратурных формулах Грегори (см. [4, с. 75]) и формуле дифференцирования назад. Все эти методы неустойчивы при к = 0.1 и к = 0.05.
Из приведенных расчетов видно, что метод (15) и метод, основанный на квадратурной формуле трапеций, даю
Для дальнейшего прочтения статьи необходимо приобрести полный текст. Статьи высылаются в формате PDF на указанную при оплате почту. Время доставки составляет менее 10 минут. Стоимость одной статьи — 150 рублей.