ЖУРНАЛ ВЫЧИСЛИТЕЛЬНОЙ МАТЕМАТИКИ И МАТЕМАТИЧЕСКОЙ ФИЗИКИ, 2013, том 53, № 9, с. 1448-1459
УДК 519.624.2
ИССЛЕДОВАНИЕ МНОГОШАГОВЫХ МЕТОДОВ ДЛЯ РЕШЕНИЯ ИНТЕГРОАЛГЕБРАИЧЕСКИХ УРАВНЕНИЙ: ПОСТРОЕНИЕ ОБЛАСТЕЙ УСТОЙЧИВОСТИ^
© 2013 г. М. В. Булатов*, О. С. Будникова**
(*664033 Иркутск, ул. Лермонтова, 134, Ин-т динамики систем и теории управления СО РАН; **664011 Иркутск, ул. Нижняя-Набережная, 6, Восточно-Сибирская гос. академия образования)
e-mail: mvbul@icc.ru, osbud@mail.ru Поступила в редакцию 14.02.2013 г.
В работе рассмотрены системы интегральных уравнений Вольтерра с тождественно вырожденной матрицей перед главной частью (интегроалгебраические уравнения). Для численного решения выделенного класса задач предложены многошаговые методы и проведено их исследование. Библ. 10. Фиг. 8. Табл. 6.
Ключевые слова: интегроалгебраические уравнения, многошаговые методы, область устойчивости.
БО1: 10.7868/80044466913070077
1. ВВЕДЕНИЕ
Многие математические модели, описывающие различные природные процессы, включают в себя интегральные уравнения Вольтерра I и II родов. Если эти уравнения взаимосвязанные, то, объединяя их, мы получаем систему интегральных уравнений Вольтерра с вырожденной матрицей перед главной частью, которые принято называть интегроалгебраическими уравнениями (ИАУ) (см. [1]). Качественная теория и численные методы решения ИАУ находятся в самом начале своего развития. Если работ по качественной теории интегральных уравнений Вольтерра II рода и некоторых интегральных уравнений Вольтерра I рода вышло и продолжает выходить достаточно большое количество (см., например, [1]—[4] и приведенную там библиографию), то по ИАУ этого сказать нельзя. Вышло всего несколько публикаций на эту тему (см. [5]—[8], там же приведен и обзор).
Если ИАУ описывают разномасштабные по времени процессы, то такие задачи будем называть, по аналогии с обыкновенными дифференциальными уравнениями (ОДУ), жесткими ИАУ.
В [5] приведены многошаговые методы высокого порядка точности для решения выделенного класса ИАУ. Цель данной работы — исследование этих методов для ИАУ, содержащих жесткие и быстроосциллирующие компоненты, а также построение областей устойчивости предложенных алгоримов. В [9] выписаны полиномы для многошаговых методов, примененных к жестким ОДУ первого порядка, записанных в интегральной форме.
2. ПОСТАНОВКА ЗАДАЧИ И ПРЕДЛАГАЕМЫЕ МЕТОДЫ
Рассмотрим систему интегральных уравнений
г
А(г)х(г) + г, ¿)х(*)й* = /(г), 0 < ^ < г < 1, (1)
0
где Л(1) и К(?, ж) суть (п х п)-матрицы,/(?) и х(0 п-мерные известная и искомая вектор-функции. Предполагается, что элементы Л(0, К(?, ж), /(?) обладают необходимой степенью гладкости. Под
1) Работа выполнена при финансовой поддержке РФФИ (коды проектов 11-01-00639-а, 13-01-93002-вьет-а).
решением исходной задачи (1) будем понимать любую непрерывную вектор-функцию х(0, обращающую (1) в тождество.
В работе рассмотрен случай, когда матрицаЛ(1) тождественно вырожденная, т.е.
ёеЫ (^ = 0. (2)
Системы (1) с условием (2) принято называть интегроалгебраическими уравнениями (ИАУ). Исследование таких задач на предмет существования и единственности решения — достаточно сложная проблема. Первая статья вышла в 1987 году (см. [10]), и с той поры опубликовано всего несколько работ по этой теме, в основном для полуявных ИАУ, когда задача (1) имеет вид
0V ^ '.......■ ^ ^
V 0 0
У (t) V Z (t) J
+
Ku (t, s) K12(t, s) V K21 (t, s) K22( t, s )J
У ( s) V z ( s )
ds =
Ф( t)
V V(t)
0 < s < t < 1, (3)
где Em — единичная матрица размерности m, K11(t, s), K12(t, s), K21(t, s), K22(t, s) — матрицы размерности (m x m), (m x (n — m)), ((n — m) x m), ((n — m) x (n — m)) соответственно, 0 — нулевые матрицы подходящих размерностей, y(t) и ф(0 суть m-мерные, z(t) и y(t) суть (n — m)-мерные вектор-функции.
В [6] для задачи (3) с условием detK22(t, t) Ф 0 Vt е [0, 1] предложены методы типа коллокаци-онных. Достаточно полная библиография по численному решению ИАУ приведена в [5], [7], [8]).
Приведем некоторые необходимые для дальнейшего изложения определения и утверждения.
Определение 2.1 (см. [10]). Пучок матриц XA(t) + B(t) удовлетворяет критерию "ранг-степень" на отрезке [0, 1] (имеет индекс один, имеет простую структуру), если
rank A (t) = deg( det(^A (t) + B( t))) = m = const Vt е [ 0, 1 ],
где X — скаляр, символ deg(.) обозначает показатель степени многочлена (.), а операция deg(0) не определена.
В [10] приведены достаточные условия существования единственного непрерывного решения задачи (1).
Теорема 2.1 (см. [10]). Пусть для задачи (1) выполнены следующие условия:
1 1 2
1) A(t) е С^ 1 ] , f(t) е Cje, 1 ] , K(t, s) е C2A , А = {0 < s < t < 1};
2) пучок XA(t) + K(t, t) удовлетворяет критерию "ранг-степень" на всем отрезке [0, 1];
3) rank A(0) = rank(A(0)| f(0)).
Тогда исходная система имеет единственное непрерывное решение.
Замечание 2.1. 1. Если A(t) = 0, то второе условие теоремы означает, что detK(t, t) Ф 0 Vt е [0, 1], а третье условие примет вид f(0) = 0.
2. Второе условие теоремы для полуявных ИАУ вида (3) означает, что detK22(t, t) Ф 0 Vt е [0, 1].
Так как интегроалгебраические уравнения включают в себя интегральные уравнения Вольтер-ра (ИУВ) II рода
t
y(t) + JK(t, s)y(s) ds = g( t), t е [ 0, 1 ], (4)
0
и ИУВ I рода с условием detK(t, t) Ф 0 Vt е [0, 1] (см. формулу (3)), то приведем некоторые результаты по численному решению таких задач.
В [3, с. 96] для численного решеного решения ИУВ II рода (4) предложены многошаговые методы вида
i
У1 - h^ юi,K¡У] = gi, i = k, k + 1, N. (5)
j = 0
Формулы (5) основаны на формуле дифференцирования назад. Известно, что данные методы устойчивы при значениях к < 6. Сложность применения данных методов заключается в том, что веса квадратурной формулы со, j, i = к, к + 1, ..., N, выписаны в явном виде только для к = 2, при остальных значения к > 3 коэффициенты можно получить только численно, через реккурентные соотношения.
0
Опишем предлагаемые в [5] методы для численного решения задачи (1). Зададим на отрезке [0, 1] равномерную сетку = ¡к, I = 1, 2, ..., N к = 1/Ы, и обозначимЛ1 = Л(^), К,] = К(?;, ?у),^ = /(?,), X « х(?,).
Общие многошаговые методы имеют вид
к I + 1
А, +1 ^а^ +1 + к^и, +1,К +1 х = /1 +1, / = к, к + 1, ..., N (6)
1 = о г = о
где Еук= 0 аух, +1 -у — аппроксимация х(?; + 1), к 2;' + 0 и, +1, гК,\ +1,1х1 — аппроксимация интегрального слагаемого, и,- +1,1 называются весами квадратурной формулы, а начальные значения х0, х1, ..., хк-1 заранее вычислены с достаточной точностью.
Формулы (6) могут быть:
1) явными при а0 Ф 0, и, +1,, +1 = 0;
2) неявными при а0 Ф 0, и, +1,, +1 Ф 0.
В силу вырожденности матрицы Л(?) можно применять только неявные методы. Однако ряд таких методов порождают неустойчивые процессы для интегрального уравнения Вольтерра I рода с ядром на диагонали, не равным нулю (см. [2], [4]). Поэтому поступим следующим образом: интегральное слагаемое в (1) будем считать по квадратурной формуле основанной на явном методе Адамса (см., например, [2], [4]), а х(1) в точке + 1 будем находить по экстраполяционной формуле по заранее вычисленным значениям х, х1 -1, ..., х I _к.
Таким образом, предлагаемые методы имеют вид
к,
А, +1 ^ ах, -1 + к ^ и +1, К +1, ,х, = +1. (7)
1 = о г = о
Опишем вычисление х I + 1. Пусть Р(х, х I -1, ..., х I - к, ?) — интерполяционный вектор-полином, проходящий через точки (х, ?,), (х,-1, ^-1), ..., (х,-к, -к). Тогда х^ + 1 будем находить по формуле
х, + 1 = Р(х1, х, - 1, —, х, - к, 0\t = I, +1 = ^ а]Х1 -
1 = 0
Приведем таблицу коэффициентов ау для к = 1, 2, ..., 5 (см. табл. 1) и значения коэффициентов и I+ 1 1 для к = 1, 2, 3 (см., [2]):
(
и
I +1, г
4
3 3
3 2 3
3 2 2 3
3 2 2 2 3
(
и,
,+1,г
12
9 0 27
9 5 11 23
9 5 16 7 23
9 5 16 12 7 23
9 5 16 12 12 7 23
(
+1, г =
24
64 -32 64 55 5 5 55 55 -4 42 -4 55 55 -4 33 33 -4 55 55 -4 33 24 33 -4 55 55 -4 33 24 24 33 -4 55 V........................у
Относительно сходимости методов (7) справедлива
Теорема 2.2 (см. [5]). Пусть для задачи (1) выполнены следующиеусловия:
1) элементы х(0, Л(0,ДО е с(0к+^ , Щ, 5) е +2), А = {0 < 5 < I < 1};
к
2) пучок XA(t) + K(t, t) удовлетворяет критерию "ранг-степень"на всем отрезке [0, 1];
3) rank A(0) = rank(A(0)| /(0));
4) для начальных значений справедливо ||ху- — x(ty)|| < Rhk +R < да, j = 0, 1, ..., k — 1.
Тогда метод (7) при k < 5 сходится к точному решению с порядком k + 1, т.е. справедлива оценка
IIX, - x(ti)|| = O(hk +1), i = к, к + 1.....N - 1.
3. ХАРАКТЕРИСТИЧЕСКИЕ УРАВНЕНИЯ МНОГОШАГОВЫХ МЕТОДОВ
Перед тем как проводить исследование методов (7), приведем некоторые известные результаты из теории численного решения обыкновенных дифференциальных уравнений (ОДУ).
На протяжении многих десятилетий для исследования свойств численных методов решения обыкновенных дифференциальных уравнений, содержащих жесткие и быстроосциллирующие компоненты, применяются тестовые задачи
x'(t) = Xx(t), x(0) = x0, t e[0, 1 ], ReX < 0, (8)
x"(t) = nx(t), x(0) = x0, x'(0) = x0, n< 0. (9)
Классические k-шаговые методы, примененные к уравнениям (8) и (9), имеют вид
к к
X axt +1-j = Zi X bjx, +1 - j, (10)
j = 0 j = 0 кк
X cjxi + 1 -j = Z2 Z dx + 1 -j, (11)
j = 0 j = 0
где z1 = Xh и z2 = nh2.
Важнейшую роль при исследовании таких схем играют характеристические полиномы
к
X (a - Z1 bj)p -j = 0, (12)
j = 0 к
X (Cj - Z2dj)qk-j = 0. (13)
j = 0
Говорят (см., например, [11]), что разностные схемы (10)((11)) устойчивы, если корни полинома zj= 0 ajpk j = 0 (zj= 0 Cjqk j = 0) лежат на (или внутри) единичной окружности и кратность корней на единичной окружности не превосходит единицы (двух). Те значения z1 (z2), при которых корни характеристического уравнения (12) ((13)) лежат в единичном круге и на границе нет кратных корней (кратность корней не превосходит двух), принято называть областью устойчивости разностной схемы (10) ((11)).
Хорошо известно, что явные и неявные (порядка выше двух) методы Адамса для ОДУ (8) имеют ограниченные области устойчивости (см., например, [11]).
Для полиномов (12
Для дальнейшего прочтения статьи необходимо приобрести полный текст. Статьи высылаются в формате PDF на указанную при оплате почту. Время доставки составляет менее 10 минут. Стоимость одной статьи — 150 рублей.