Журнал вычислительной математики и математической физики, 2021, T. 61, № 11, стр. 1825-1838

Численное решение интегроалгебраических уравнений со слабой граничной особенностью k-шаговыми методами

М. Н. Ботороева 1*, О. С. Будникова 1**, М. В. Булатов 2***, С. С. Орлов 1****

1 Иркутский государственный университет
664003 Иркутск, ул. Карла Маркса, 1, Россия

2 Институт динамики систем и теории управления им. В.М. Матросова СО РАН
664033 Иркутск, ул. Лермонтова, 134, Россия

* E-mail: masha888888@mail.ru
** E-mail: osbud@mail.ru
*** E-mail: mvbul@icc.ru
**** E-mail: orlov_sergey@inbox.ru

Поступила в редакцию 21.11.2020
После доработки 06.04.2021
Принята к публикации 07.07.2021

Полный текст (PDF)

Аннотация

В статье излагается построение $k$-шаговых методов решения систем интегральных уравнений типа Вольтерра I и II рода со слабой степенной особенностью ядер в нижнем пределе интегрирования. Матрично-векторная форма таких систем имеет вид абстрактного уравнения с вырожденной матрицей коэффициентов при внеинтегральных слагаемых, которое называют интегроалгебраическим уравнением. Предлагаемые методы основаны на экстраполяционных формулах для главной части, многошаговых методах типа Адамса и формуле интегрирования произведений для интегрального члена. Веса построенных квадратурных формул получены в явном виде. Доказана теорема о сходимости разработанных методов. Приведены численные расчеты тестовых примеров, иллюстрирующие теоретические результаты. Библ. 30. Фиг. 2. Табл. 8.

Ключевые слова: интегроалгебраические уравнения, многошаговые методы, слабая граничная особенность.

1. ВВЕДЕНИЕ

Разработка эффективных численных методов решения интегральных уравнений является одним из центральных направлений вычислительной математики. К настоящему времени созданы и обоснованы методы решения интегральных уравнений Вольтерра (ИУВ) II рода, некоторых ИУВ I рода, в том числе и с различными особенностями. Достаточно полный обзор результатов приведен в монографиях [1] и [2]. Среди статей выделим недавно опубликованные [3] и [4]. В работе [5] методы кусочно-полиномиальной коллокации применяются для численного решения ИУВ II рода с двумя слабыми степенными особенностями ядра в нижнем и верхнем пределах интегрирования, которые авторы M. Kolk и A. Pedas называют граничной и диагональной особенностями соответственно.

В представляемой статье рассматривается система взаимосвязанных линейных ИУВ I и II рода со слабой граничной особенностью ядра. Такие системы можно представить в виде абстрактного интегрального уравнения с тождественно вырожденной матрицей в главной части, которое будем называть интегроалгебраическим уравнением (ИАУ) со слабой граничной особенностью. Данные классы ИАУ ранее, по-видимому, не рассматривались. Близкими им объектами являются дифференциальноалгебраические уравнения (ДАУ) с вырожденной матрицей перед производной, имеющей нуль в начальной точке [6], [7].

Несколько большее внимание в современной математической литературе уделено ИАУ со слабой диагональной особенностью, которые также называют ИАУ типа Абеля. Прежде всего следует отметить статью [8], где получены условия разрешимости не только ИАУ типа Абеля, но и систем вырожденных слабосингулярных интегродифференциальных уравнений (ИДУ). Продолжением этих исследований стала работа [9], в которой для ИАУ типа Абеля предложены методы численного анализа, а именно, разработанные ранее многошаговые методы для ИАУ с регулярным (достаточно гладким) ядром [10] перенесены на ИАУ со слабой диагональной особенностью в интегральном слагаемом. Примерно в это же время были получены результаты для численного решения слабосингулярных ИАУ методами коллокационного типа [11]. В целом теория ИАУ в начале своего развития. По всей видимости, первым изучать ИАУ начал В.Ф. Чистяков в своей работе [12], в которой сформулированы достаточные условия существования единственного непрерывного решения ИАУ с регулярным ядром и предложен численный метод первого порядка точности. С подробным обзором результатов в области ИАУ можно ознакомиться по монографии [13] В.Ф. Чистякова и недавно опубликованным статьям [14] О.С. Будниковой, [15] V. Balakumar и K. Murugesan, [16] M.S. Farahani и M. Hadizadeh, [17] М.В. Булатова и Е.В. Чистяковой.

Родственными объектами ИАУ с особенностями ядра являются системы интегральных уравнений Вольтерра с кусочно-гладкими ядрами, которые имеют разрывы I рода. Исследованию однозначной разрешимости таких систем в классе непрерывных функций и распределений Соболева–Шварца с компактным носителем, свойств решений, а также разработке численных методов посвящена монография [18] Д.Н. Сидорова. Рассмотрен случай, когда ядра интегральных уравнений имеют конечное количество линий разрыва, заданных явно и пересекающихся в начальной точке отрезка интегрирования. Развернутая форма такой системы содержит интегральные слагаемые с функциональными пределами интегрирования, принимающими одинаковые значения в левой концевой точке отрезка, на котором задана исходная система. Более сложная ситуация имеет место, когда указанное равенство не выполняется. В этом случае в соответствующих интегральных слагаемых возникает дополнительный эффект эредитарности (запаздывания) и, как следствие, необходимость задать решение на предыстории. Качественная теория и численные методы уравнений Вольтерра I рода с функциональными пределами интегрирования детально и полно изложены в монографии [1] А.С. Апарцина. Эти классы уравнений представляют большой интерес с точки зрения математического моделирования развивающихся динамических систем. В частности, на их основе проводились численные расчеты в задачах функционирования двухсекторной экономики [19], ввода разнотипных генерирующих мощностей электроэнергетической системы [20], управление накопителями электроэнергии для возобновляемых источников [21]. ИАУ с функциональными пределами интегрирования и запаздыванием, помимо вырождения главной части, наследуют также и все трудности исследования, присущие функциональным и функционально-дифференциальным уравнениям. В статье [22] М.Н. Ботороевой разработаны аналитические и эффективные численные методы решения специальных классов таких ИАУ.

В настоящей статье предлагаются многошаговые методы решения интегроалгебраических уравнений с особенностью ядра в нижнем пределе интегрирования. Предлагаемые методы основаны на экстраполяционных формулах для главной части, на явных многошаговых методах типа Адамса и формуле интегрирования произведений для интегрального слагаемого. Выбор численных методов обусловлен тем, что они хорошо себя зарекомендовали при решении ИУВ I рода с ядром, отличным от нуля на диагонали.

2. ПОСТАНОВКА ЗАДАЧИ

Рассмотрим систему интегральных уравнений

(1)
$A(t)x(t) + \int\limits_0^t {{{s}^{{ - a}}}} K(t,s)x(s)ds = f(t),\quad t \in [0,\;1],\quad 0 < a < 1,$
где $A(t)$ и $K(t,s)$ – заданные $(n \times n)$-матрицы, $f = f(t)$ и $x = x(t)$ суть $n$-мерные известная и искомая вектор-функции такие, что $A(t) \in {{C}^{p}}([0,\;1];{{\mathbb{R}}^{{n \times n}}})$, $K(t,s) \in {{C}^{p}}(\Delta ;{{\mathbb{R}}^{{n \times n}}})$ и $f(t) \in {{C}^{p}}([0,\;1];{{\mathbb{R}}^{n}})$ (элементы соответствующих матриц обладают необходимой степенью гладкости $p \in \mathbb{N}$), и при всех $t \in [0,\;1]$ выполняется условие
(2)
$A(t) \ne {{\mathbb{O}}_{n}},\quad detA(t) = 0.$
Здесь и всюду далее $\Delta $ – компакт вида $\Delta = {\text{\{ }}(t,s) \in {{\mathbb{R}}^{2}}\,:0 \leqslant t \leqslant 1,\;0 \leqslant s \leqslant t{\text{\} }}$, символ ${{\mathbb{R}}^{{n \times n}}}$ обозначает вещественное векторное пространство квадратных матриц порядка $n$ с операторной нормой $\left\| A \right\| = sup\left\{ {{{{\left| {Ax} \right|}}_{{{{\mathbb{R}}^{n}}}}}:{{{\left| x \right|}}_{{{{\mathbb{R}}^{n}}}}} = 1} \right\}$, а ${{\mathbb{O}}_{n}}$ – нулевую матрицу порядка $n$. Систему интегральных уравнений (1) с условием (2) назовем ИАУ со слабой граничной особенностью. Под решением ИАУ (1) будем понимать любую непрерывную на отрезке $[0,\;1]$ вектор-функцию $x = x(t)$, которая обращает данное ИАУ в тождество.

Приведем определения и утверждения, которые понадобятся в дальнейшем изложении.

Определение 2.1 (см. [12]). Пучок матриц $\lambda B(t) + C(t)$ удовлетворяет критерию ранг-степень на отрезке $[0,\;1]$ (имеет индекс один или простую структуру), если при всех $t \in [0,\;1]$ выполняется условие

$\operatorname{rank} B(t) = \deg (\det (\lambda B(t) + C(t))) = m = {\text{const}},$
где $\lambda $ – скаляр, символ $deg(.)$ означает показатель степени многочлена $(.)$, а операция $deg(0)$ не определена.

Отметим некоторые важные свойства пучков матриц.

Лемма 2.1 (см. [12]). Пусть $B(t),C(t) \in {{C}^{p}}([0,\;1];{{\mathbb{R}}^{{n \times n}}})$, и пучок матриц $\lambda B(t) + C(t)$ удовлетворяет критерию ранг-степень на отрезке $[0,\;1]$. Тогда существуют неособые матрицы $P(t),Q(t) \in {{C}^{p}}([0,\;1];{{\mathbb{R}}^{{n \times n}}})$ такие, что

$P(t)B(t)Q(t) = {\text{diag}}\left\{ {{{\mathbb{E}}_{m}},{{\mathbb{O}}_{{n - m}}}} \right\},\quad P(t)C(t)Q(t) = {\text{diag}}\left\{ {{{J}_{m}}(t),{{\mathbb{O}}_{{n - m}}}} \right\},$
где ${{J}_{m}}(t)$ есть $(m \times m)$-матрица, а ${{\mathbb{E}}_{m}}$ – единичная матрица порядка $m$.

Из этой леммы вытекает блочное представление матрицы $B(t)$

$P(t)B(t) = \left( {\begin{array}{*{20}{c}} {{{B}^{1}}(t)} \\ \mathbb{O} \end{array}} \right),$
где ${{B}^{1}}(t)$ есть $(m \times n)$-матрица такая, что $\operatorname{rank} {{B}^{1}}(t) = m$, а $\mathbb{O}$ – нулевая $((n - m) \times n)$-матрица.

Лемма 2.2 (см. [23]). Если пучок матриц $\lambda B(t) + C(t)$ удовлетворяет критерию ранг-степень на всем отрезке $[0,\;1]$, и матрицы $B(t)$ и $C(t)$ имеют блочный вид

$B(t) = \left( {\begin{array}{*{20}{c}} {{{B}^{1}}(t)} \\ \mathbb{O} \end{array}} \right),\quad C(t) = \left( {\begin{array}{*{20}{c}} {{{C}^{1}}(t)} \\ {{{C}^{2}}(t)} \end{array}} \right),$
где ${{B}^{1}}(t)$ и ${{C}^{1}}(t)$ суть $(m \times n)$-матрицы такие, что $\operatorname{rank} {{B}^{1}}(t) = m = {\text{const}}$, $\mathbb{O}$ – нулевая $((n - m) \times n)$-матрица, ${{C}^{2}}(t)$ есть $((n - m) \times n)$-матрица, то при любом $t \in [0,\;1]$ выполняется

$det\left( {\begin{array}{*{20}{c}} {{{B}^{1}}(t)} \\ {{{C}^{2}}(t)} \end{array}} \right) \ne 0.$

Лемма 2.3. Пусть $K(t,s) \in C(\Delta ;{{\mathbb{R}}^{{n \times n}}})$, $f(t) \in C([0,\;1];{{\mathbb{R}}^{n}})$, тогда система

$x(t) + \int\limits_0^t {{{s}^{{ - a}}}} K(t,s)x(s)ds = f(t),\quad t \in [0,\;1],\quad 0 < a < 1,$
линейных интегральных уравнений Вольтерра II рода со слабой граничной особенностью имеет единственное непрерывное на отрезке $[0,\;1]$ решение.

Доказательство. Повторная суперпозиция ${{\mathcal{L}}^{{n + 1}}}:C([0,\;1];{{\mathbb{R}}^{n}}) \to C([0,\;1];{{\mathbb{R}}^{n}})$

${{\mathcal{L}}^{{n + 1}}}x = {{( - 1)}^{{n + 1}}}\int\limits_0^t {{{s}^{{ - a}}}} {{K}_{{n + 1}}}(t,s)x(s)ds + \sum\limits_{k = 1}^n {{{{( - 1)}}^{k}}} \int\limits_0^t {{{s}^{{ - a}}}} {{K}_{k}}(t,s)f(s)ds + f(t),\quad n \in \mathbb{N},$
линейного оператора $\mathcal{L}:C([0,\;1];{{\mathbb{R}}^{n}}) \to C([0,\;1];{{\mathbb{R}}^{n}})$
$\mathcal{L}x = - \int\limits_0^t {{{s}^{{ - a}}}} K(t,s)x(s)ds + f(t),$
при любых заданных $K(t,s) \in C(\Delta ;{{\mathbb{R}}^{{n \times n}}})$ и $f(t) \in C([0,\;1];{{\mathbb{R}}^{n}})$ осуществляет сжимающее отображение $C([0,\;1];{{\mathbb{R}}^{n}})$ в себя, начиная с некоторого номера $n$. Это следует из верной для всех $(t,s) \in \Delta $ оценки
$\left\| {{{K}_{n}}(t,s)} \right\| \leqslant {{K}^{n}}\frac{{{{{({{t}^{{1 - a}}} - {{s}^{{1 - a}}})}}^{{n - 1}}}}}{{{{{(1 - a)}}^{{n - 1}}}(n - 1)!}},\quad K = \mathop {max}\limits_{(t,s) \in \Delta } \left\| {K(t,s)} \right\|,$
определяемых рекуррентно
${{K}_{{n + 1}}}(t,s) = \int\limits_s^t {{{\tau }^{{ - a}}}} K(t,\tau ){{K}_{n}}(\tau ,s)d\tau ,\quad {{K}_{1}}(t,s) = K(t,s),$
элементов матрично-функциональной последовательности ${{{\text{\{ }}{{K}_{n}}(t,s){\text{\} }}}_{{n \in \mathbb{N}}}}$. Тогда линейный оператор $\mathcal{L}$ имеет в $C([0;\;1];{{\mathbb{R}}^{n}})$ единственную неподвижную точку, которая и является непрерывным решением рассматриваемой системы линейных интегральных уравнений.

Справедлива следующая

Теорема 2.1. Пусть выполнено условие (2), и

1) ${\text{rank}}\left( {A(0)\,|\,f(0)} \right) = \operatorname{rank} A(0) = m < n$;

2) $A(t) \in {{C}^{1}}([0,\;1];{{\mathbb{R}}^{{n \times n}}})$; $K(t,s),{{\partial }_{t}}K(t,s) \in C(\Delta ;{{\mathbb{R}}^{{n \times n}}})$;

3) $f(t) \in {{C}^{1}}([0,\;1];{{\mathbb{R}}^{n}})$;

4) пучок матриц $\lambda A(t) + K(t,t)$ удовлетворяет критерию ранг-степень на $[0,\;1]$.

Тогда ИАУ (1) имеет единственное непрерывное на отрезке $[0,\;1]$ решение.

Доказательство. При $t = 0$ ИАУ (1) имеет вид системы линейных алгебраических уравнений $A(0)x(0) = f(0)$, которая совместна в силу первого условия теоремы. Далее в условиях 2, 3 и 4 теоремы с учетом лемм 2.1 и 2.2 исходное ИАУ допускает редукцию к системе интегральных уравнений Вольтерра II рода

$\left( {\begin{array}{*{20}{c}} {{{x}_{1}}(t)} \\ {{{x}_{2}}(t)} \end{array}} \right) + \int\limits_0^t {{{s}^{{ - a}}}} \mathop {\left( {\begin{array}{*{20}{c}} {{{A}^{1}}(t)} \\ {{{K}^{2}}(t,t)} \end{array}} \right)}\nolimits^{ - 1} \left( {\begin{array}{*{20}{c}} {{{K}^{1}}(t,s)} \\ {{{t}^{a}}{{\partial }_{t}}{{K}^{2}}(t,s)} \end{array}} \right)\left( {\begin{array}{*{20}{c}} {{{x}_{1}}(s)} \\ {{{x}_{2}}(s)} \end{array}} \right)ds = \mathop {\left( {\begin{array}{*{20}{c}} {{{A}^{1}}(t)} \\ {{{K}^{2}}(t,t)} \end{array}} \right)}\nolimits^{ - 1} \left( {\begin{array}{*{20}{c}} {{{f}_{1}}(t)} \\ {{{t}^{a}}f_{2}^{'}(t)} \end{array}} \right),$
которая имеет единственное непрерывное на отрезке $[0,\;1]$ решение, согласно лемме 3.

Отметим, что условие 4 теоремы 2.1 является достаточным. Если оно не выполняется, то, как показывают следующие ниже примеры, исследование однозначной разрешимости ИАУ (1) сталкивается с большими трудностями.

Пример 2.1. Рассмотрим систему интегральных уравнений (1) с компонентами

$A(t) = \left( {\begin{array}{*{20}{c}} {\psi (t)}&0&0 \\ 0&1&0 \\ 0&0&0 \end{array}} \right),\quad K(t,s) = \left( {\begin{array}{*{20}{c}} { - 1}&0&0 \\ 0&{ - 1}&0 \\ 0&0&1 \end{array}} \right),\quad f(t) = \left( {\begin{array}{*{20}{c}} { - {{{\left| {\psi (0)} \right|}}^{{\tfrac{1}{{1 - a}}}}}} \\ 1 \\ t \end{array}} \right),$
где $\psi :[0,\;1] \to \mathbb{R}$ – заданная функция. Здесь $A(t) \ne {{\mathbb{O}}_{3}}$ и $detA(t) = 0$ при всех $t \in [0,\;1]$. Заметим, что $\operatorname{rank} A(t) = 1$ в точках $t \in [0,\;1]$, где $\psi (t) = 0$, и $\operatorname{rank} A(t) = 2$ в точках $t \in [0,\;1]$ таких, что $\psi (t) \ne 0$. Таким образом, в данном примере критерий ранг-степень не выполняется в точках $t \in [0,\;1]$, где $\psi (t) = 0$. Это, в частности, может приводить к неединственности решения. Например, если $\psi (t) = t$, то рассматриваемое ИАУ имеет семейство непрерывных решений
$x(t) = \mathop {\left( {C\frac{{{{e}^{{ - \tfrac{{{{t}^{{ - a}}}}}{a}}}}}}{t},{{e}^{{ - \tfrac{{{{t}^{{1 - a}}}}}{{1 - a}}}}},{{t}^{a}}} \right)}\nolimits^{\text{т}} .$
При $\psi (t) = {{t}^{{1 - a}}} - {{\alpha }^{{1 - a}}}$, где $0 < \alpha < 1$, критерий ранг-степень нарушается в точке $t = \alpha $. В этом случае ИАУ имеет множество непрерывных решений вида
$x(t) = \mathop {\left( {\left\{ \begin{gathered} {{({{\alpha }^{{1 - a}}} - {{t}^{{1 - a}}})}^{{\tfrac{a}{{1 - a}}}}},\quad 0 \leqslant t \leqslant \alpha ; \hfill \\ C{{({{t}^{{1 - a}}} - {{\alpha }^{{1 - a}}})}^{{\tfrac{a}{{1 - a}}}}},\quad \alpha < t \leqslant 1; \hfill \\ \end{gathered} \right.,\quad {{e}^{{ - \tfrac{{{{t}^{{1 - a}}}}}{{1 - a}}}}},{{t}^{a}}} \right)}\nolimits^{\text{т}} .$
На отрезке $[0,\alpha ]$ решение оказывается единственным, а в точке $t = \alpha $ имеет место ветвление решения. В приведенных выше формулах $C \in \mathbb{R}$ – произвольная постоянная.

Пример 2.2. Рассмотрим систему интегральных уравнений (1) с компонентами

$A(t) = \left( {\begin{array}{*{20}{c}} t&0&0 \\ 0&1&0 \\ 0&0&0 \end{array}} \right),\quad K(t,s) = \left( {\begin{array}{*{20}{c}} 0&0&{ - 1} \\ 0&{ - 1}&0 \\ 1&0&0 \end{array}} \right),\quad f(t) = \left( {\begin{array}{*{20}{c}} 0 \\ 1 \\ t \end{array}} \right).$
Здесь $A(t) \ne {{\mathbb{O}}_{3}}$ и $detA(t) \equiv 0$, причем $\operatorname{rank} A(t) = 1$ в точке $t = 0$, и $\operatorname{rank} A(t) = 2$ при $t \in (0,\;1]$. Матричный многочлен $det(\lambda A(t) + K(t,t)) = \lambda - 1$ имеет первую степень на отрезке $[0,\;1]$, значит, критерий ранг-степень нарушается при всех $t \in (0,\;1]$. Тем не менее рассматриваемое ИАУ имеет единственное непрерывное решение вида

$x(t) = \mathop {\left( {{{t}^{a}},{{e}^{{ - \tfrac{{{{t}^{{1 - a}}}}}{{1 - a}}}}},(a + 1){{t}^{{2a}}}} \right)}\nolimits^{\text{т}} .$

3. ЧИСЛЕННЫЙ МЕТОД

Зададим на отрезке $[0,\;1]$ равномерную сетку

${{t}_{i}} = ih,\quad i = 0,\;1,\; \ldots ,\;N - 1,\quad h = \frac{1}{N},$
и введем обозначения ${{A}_{i}} = A({{t}_{i}}),$ ${{K}_{{i,j}}} = K({{t}_{i}},{{t}_{j}}),$ ${{f}_{i}} = f({{t}_{i}}),$ ${{x}_{i}} \approx x({{t}_{i}}).$

Эффективными способами борьбы со слабыми особенностями являются выделение весовой функции и применение квадратурных формул, соответствующих данной весовой функции [24], [25]. Тогда возникает вопрос о выборе подходящих квадратурных формул. В данной статье предлагаются алгоритмы, основанные на явных методах типа Адамса, так как они зарекомендовали себя при численном решении интегральных уравнений Вольтерра I рода с ядром на диагонали, не равным нулю [26]. Явные методы Адамса, описание которых можно найти, например в [10], [26], [27], с весовой функцией ${{p}_{a}}(t,s) = {{s}^{{ - a}}}$, $a \in (0,\;1)$, примут следующий вид. Для заданной функции $g = g(t)$ справедлива цепочка равенств

(3)
$\begin{gathered} \int\limits_0^{{{t}_{{i + 1}}}} {{{s}^{{ - a}}}} g(s)ds = \int\limits_0^{{{t}_{{k + 1}}}} {{{s}^{{ - a}}}} g(s)ds + \sum\limits_{j = k + 1}^i {\int\limits_{{{t}_{j}}}^{{{t}_{{j + 1}}}} {{{s}^{{ - a}}}} } g(s)ds \approx \\ \approx \;\int\limits_0^{{{t}_{{k + 1}}}} {{{s}^{{ - a}}}} L_{{k + 1}}^{k}({{g}_{0}},{{g}_{1}},\; \ldots ,\;{{g}_{k}},s)ds + \sum\limits_{j = k + 1}^i {\int\limits_{{{t}_{j}}}^{{{t}_{{j + 1}}}} {{{s}^{{ - a}}}} } L_{{k + 1}}^{j}({{g}_{{j - k}}},{{g}_{{j - k + 1}}},\; \ldots ,\;{{g}_{j}},s)ds = \\ = \;{{h}^{{1 - a}}}\sum\limits_{l = 0}^k {{{d}_{l}}} {{g}_{l}} + \sum\limits_{j = k + 1}^i {{{h}^{{1 - a}}}} \sum\limits_{l = 0}^k {{{b}_{l}}} (j){{g}_{{j - l}}} = {{h}^{{1 - a}}}\sum\limits_{l = 0}^i {{{\omega }_{{i + 1,l}}}{{g}_{l}}} . \\ \end{gathered} $
Здесь $L_{{k + 1}}^{j}({{g}_{{j - k}}},{{g}_{{j - k + 1}}},\; \ldots ,\;{{g}_{j}},s)$ – интерполяционный полином степени $k$, проходящий через точки $({{g}_{{j - k}}},{{s}_{{j - k}}})$, $({{g}_{{j - k + 1}}},{{s}_{{j - k + 1}}}),\; \ldots ,\;({{g}_{j}},{{s}_{j}})$, где $j = k + 1,k + 2,\; \ldots ,\;i$, и $g({{s}_{j}}) = {{g}_{j}}$.

Прямое применение явной квадратурной формулы (3) приведет к решению вырожденной СЛАУ в силу условия (2). Поэтому будем находить значение выражения ${{A}_{{i + 1}}}{{x}_{{i + 1}}}$ по аналогии с [10] следующим образом:

${{A}_{{i + 1}}}{{x}_{{i + 1}}} \approx {{\left. {{{A}_{{i + 1}}}L_{{k + 1}}^{i}\left( {{{x}_{{i - k}}},{{x}_{{i - k + 1}}},\; \ldots ,\;{{x}_{i}},t} \right)} \right|}_{{t = {{t}_{{i + 1}}}}}} = {{A}_{{i + 1}}}\sum\limits_{j = 0}^k {{{\alpha }_{j}}} {{x}_{{i - j}}},$
где $L_{{k + 1}}^{i}({{x}_{{i - k}}},{{x}_{{i - k + 1}}},\; \ldots ,\;{{x}_{i}},t)$ – интерполяционный полином степени $k,$ проходящий через точки $({{x}_{{i - k}}},{{t}_{{i - k}}})$, $({{x}_{{i - k + 1}}},{{t}_{{i - k + 1}}}),\; \ldots ,\;({{x}_{i}},{{t}_{i}})$. Коэффициенты ${{\alpha }_{j}}$ для различных $k = 1,$ $2,\; \ldots ,\;5$ приведены в табл. 1. С учетом отмеченного выше, предлагаемые многошаговые методы имеют вид
(4)
${{A}_{{i + 1}}}\sum\limits_{j = 0}^k {{{\alpha }_{j}}} {{x}_{{i - j}}} + {{h}^{{1 - a}}}\sum\limits_{l = 0}^i {{{\omega }_{{i + 1,l}}}} {{K}_{{i + 1,l}}}{{x}_{l}} = {{f}_{{i + 1}}},\quad i = k,\;k + 1,\; \ldots ,\;N - 1.$
Предполагается, что начальные значения ${{x}_{0}},{{x}_{1}},\; \ldots ,\;{{x}_{{k - 1}}}$ заранее вычислены с достаточной точностью. Отметим, что точность решения, получаемого по формуле (4), зависит от погрешности аппроксимации и нулей характеристического полинома
$\sum\limits_{l = 0}^k {{{b}_{l}}} (i){{\lambda }^{{k - l}}} = 0,$
связанного с коэффициентами ${{b}_{l}}(i)$ в формуле (3). Корни данного многочлена по модулю меньше единицы при $k < 6$. В этом можно убедиться непосредственной проверкой.

Таблица 1.  

Значения коэффициентов ${{\alpha }_{j}}$

$k$ ${{\alpha }_{0}}$ ${{\alpha }_{1}}$ ${{\alpha }_{2}}$ ${{\alpha }_{3}}$ ${{\alpha }_{4}}$ ${{\alpha }_{5}}$
1 2 –1
2 3 –3 1
3 4 –6 4 –1
4 5 –10 10 –5 1
5 6 –15 20 –15 6 –1

Приведем многошаговые алгоритмы при различных значениях $k$. При $k = 0$ формула (4) принимает вид

${{A}_{{i + 1}}}{{x}_{i}} + {{h}^{{1 - a}}}\sum\limits_{l = 0}^i {{{\omega }_{{i + 1,l}}}} {{K}_{{i + 1,l}}}{{x}_{l}} = {{f}_{{i + 1}}},$
где
${{\omega }_{{i + 1,l}}} = \frac{{{{{(l + 1)}}^{{1 - a}}} - {{l}^{{1 - a}}}}}{{1 - a}}.$
При $k = 1$ получим метод
${{A}_{{i + 1}}}(2{{x}_{i}} - {{x}_{{i - 1}}}) + {{h}^{{1 - a}}}\sum\limits_{l = 0}^i \,{{\omega }_{{i + 1,l}}}{{K}_{{i + 1,l}}}{{x}_{l}} = {{f}_{{i + 1}}},$
где веса имеют вид
${{\omega }_{{i + 1,l}}} = \left( {\begin{array}{*{20}{c}} {{{d}_{0}}}&{{{d}_{1}}}&{}&{}&{} \\ {{{d}_{0}}}&{{{d}_{1}} + {{b}_{1}}(2)}&{{{b}_{0}}(2)}&{}&{} \\ {{{d}_{0}}}&{{{d}_{1}} + {{b}_{1}}(2)}&{{{b}_{0}}(2) + {{b}_{1}}(3)}&{{{b}_{0}}(3)}&{} \\ {{{d}_{0}}}&{{{d}_{1}} + {{b}_{1}}(2)}&{{{b}_{0}}(2) + {{b}_{1}}(3)}&{{{b}_{0}}(3) + {{b}_{1}}(4)}&{{{b}_{0}}(4)} \\ \ldots & \ldots & \ldots & \ldots & \ldots \end{array}} \right),$
здесь
${{d}_{0}} = \frac{{{{2}^{{1 - a}}}}}{{1 - a}} - \frac{{{{2}^{{2 - a}}}}}{{2 - a}},$
${{d}_{1}} = \frac{{{{2}^{{2 - a}}}}}{{(2 - a)}},$
${{b}_{1}}(l) = - \frac{{{{{(l + 1)}}^{{1 - a}}}}}{{1 - a}} + \frac{{{{{(l + 1)}}^{{2 - a}}} - {{l}^{{2 - a}}}}}{{(1 - a)(2 - a)}},$
${{b}_{0}}(l) = \frac{{2{{{(l + 1)}}^{{1 - a}}} - {{l}^{{1 - a}}}}}{{1 - a}} - \frac{{{{{(l + 1)}}^{{2 - a}}} - {{l}^{{2 - a}}}}}{{(1 - a)(2 - a)}}.$
Приведем еще один алгоритм при $k = 2$
${{A}_{{i + 1}}}(3{{x}_{i}} - 3{{x}_{{i - 1}}} + {{x}_{{i - 2}}}) + {{h}^{{1 - a}}}\sum\limits_{l = 0}^i {{{\omega }_{{i + 1,l}}}} {{K}_{{i + 1,l}}}{{x}_{l}} = {{f}_{{i + 1}}},$
где
$\begin{gathered} {{\omega }_{{i + 1,l}}} = \\ = \left( {\begin{array}{*{20}{c}} {{{d}_{0}}}&{{{d}_{1}}}&{{{d}_{2}}}&{}&{}&{}&{} \\ {{{d}_{0}}}&{{{d}_{1}} + {{b}_{2}}(3)}&{{{d}_{2}} + {{b}_{1}}(3)}&{{{b}_{0}}(3)}&{}&{}&{} \\ {{{d}_{0}}}&{{{d}_{1}} + {{b}_{2}}(3)}&{{{d}_{2}} + {{b}_{1}}(3) + {{b}_{2}}(4)}&{{{b}_{0}}(3) + {{b}_{1}}(4)}&{{{b}_{0}}(4)}&{}&{} \\ {{{d}_{0}}}&{{{d}_{1}} + {{b}_{2}}(3)}&{{{d}_{2}} + {{b}_{1}}(3) + {{b}_{2}}(4)}&{{{b}_{0}}(3) + {{b}_{1}}(4) + {{b}_{2}}(5)}&{{{b}_{0}}(4) + {{b}_{1}}(5)}&{{{b}_{0}}(5)}&{} \\ {{{d}_{0}}}&{{{d}_{1}} + {{b}_{2}}(3)}&{{{d}_{2}} + {{b}_{1}}(3) + {{b}_{2}}(4)}&{{{b}_{0}}(3) + {{b}_{1}}(4) + {{b}_{2}}(5)}&{{{b}_{0}}(4) + {{b}_{1}}(5) + {{b}_{2}}(6)}&{{{b}_{0}}(5) + {{b}_{1}}(6)}&{{{b}_{0}}(6)} \\ \ldots & \ldots & \ldots & \ldots & \ldots & \ldots & \ldots \end{array}} \right). \\ \end{gathered} $
Коэффициенты ${{d}_{j}},{{b}_{j}},j = 0,\;1,\;2$, вычислены по формулам

${{d}_{0}} = \frac{{{{3}^{{1 - a}}}}}{{1 - a}} - \frac{{{{3}^{{3 - a}}}}}{{2(2 - a)(3 - a)}},$
${{d}_{1}} = - \frac{{{{3}^{{2 - a}}}}}{{2 - a}} + \frac{{{{3}^{{3 - a}}}}}{{(2 - a)(3 - a)}},$
${{d}_{2}} = \frac{{{{3}^{{2 - a}}}}}{{2 - a}} - \frac{{{{3}^{{3 - a}}}}}{{2(2 - a)(3 - a)}},$
${{b}_{2}}(l) = \frac{{{{{(l + 1)}}^{{1 - a}}}}}{{1 - a}} - \frac{{3{{{(l + 1)}}^{{2 - a}}} - {{l}^{{2 - a}}}}}{{2(1 - a)(2 - a)}} + \frac{{{{{(l + 1)}}^{{3 - a}}} - {{l}^{{3 - a}}}}}{{(1 - a)(2 - a)(3 - a)}},$
${{b}_{1}}(l) = - \frac{{3{{{(l + 1)}}^{{1 - a}}}}}{{1 - a}} + \frac{{4{{{(l + 1)}}^{{2 - a}}} - 2{{l}^{{2 - a}}}}}{{(1 - a)(2 - a)}} - 2\frac{{{{{(l + 1)}}^{{3 - a}}} - {{l}^{{3 - a}}}}}{{(1 - a)(2 - a)(3 - a)}},$
${{b}_{0}}(l) = \frac{{3{{{(l + 1)}}^{{1 - a}}} - {{l}^{{1 - a}}}}}{{1 - a}} - \frac{{5{{{(l + 1)}}^{{2 - a}}} - 3{{l}^{{2 - a}}}}}{{2(1 - a)(2 - a)}} + \frac{{{{{(l + 1)}}^{{3 - a}}} - {{l}^{{3 - a}}}}}{{(1 - a)(2 - a)(3 - a)}}.$

Алгоритмы для $k = 3,\;4,\;5$ не приведены из-за громозкости формул весовых коэффициентов. Отметим, что из формул приближенного вычисления следует рекуррентное соотношение

$\begin{array}{*{20}{c}} {{{\omega }_{{i + 1,l}}} = {{\omega }_{{i,l}}},\quad l = 0,\;1,\; \ldots ,\;i - k - 1,} \\ {{{\omega }_{{i + 1,i - k}}} = {{\omega }_{{i,i - k}}} + {{b}_{k}}(i),\; \ldots ,\;{{\omega }_{{i + 1,i - 1}}} = {{\omega }_{{i,i - 1}}} + {{b}_{1}}(i),\quad {{\omega }_{{i + 1,i}}} = {{b}_{0}}(i).} \end{array}$
Приведем результат о сходимости предложенных методов (4).

Теорема 3.1. Пусть выполнено условие (2), и

1) $A(t) \in {{C}^{{k + 1}}}([0,\;1];{{\mathbb{R}}^{{n \times n}}})$; $K(t,s) \in {{C}^{{k + 2}}}([0,\;1];{{\mathbb{R}}^{{n \times n}}})$; $x(t),f(t) \in {{C}^{{k + 1}}}([0,\;1];{{\mathbb{R}}^{n}})$;

2) пучок матриц $\lambda A(t) + K(t,t)$ удовлетворяет критерию ранг-степень на $[0,\;1]$;

3) $\operatorname{rank} A(0) = {\text{rank}}(A(0)\,|\,f(0))$;

4) для начальных значений справедливо

${{\left| {{{x}_{j}} - x({{t}_{j}})} \right|}_{{{{\mathbb{R}}^{n}}}}} \leqslant R{{h}^{{k + 1}}},\quad 0 < R < + \infty ,\quad j = 0,\;1,\; \ldots ,\;k - 1,\quad k \leqslant 5.$
Тогда справедлива оценка

${{\left| {{{x}_{i}} - x({{t}_{i}})} \right|}_{{{{\mathbb{R}}^{n}}}}} = O({{h}^{{k + 1}}}),\quad i = k,\;k + 1,\; \ldots ,\;N - 1.$

Доказательство данного утверждения не приводится, ввиду его чрезвычайной технической громоздкости. Оно основано на лемме 2.1 с использованием приема перехода от многошаговых методов к одношаговым и оценке остаточного члена погрешности интерполяционного полинома по методике работы [10].

4. ЧИСЛЕННЫЕ ЭКСПЕРИМЕНТЫ

Численные расчеты проводились на нескольких тестовых ИАУ со слабой граничной особенностью. Введем обозначение ${\text{err}}$ для погрешности вычислений по евклидовой норме ${{\left| {\, \cdot \,} \right|}_{{{{\mathbb{R}}^{n}}}}}$ в ${{\mathbb{R}}^{n}}$, и $p$ для порядка метода, определяемого по формуле

$p = lo{{g}_{2}}\frac{{\mathop {{\text{err}}}\nolimits_h }}{{\mathop {{\text{err}}}\nolimits_{\tfrac{h}{2}} }}.$

Пример 4.1:

$\left( {\begin{array}{*{20}{c}} 1&t \\ t&{{{t}^{2}}} \end{array}} \right)x(t) + \int\limits_0^t {{{s}^{{ - a}}}} \left( {\begin{array}{*{20}{c}} {{{e}^{{t - s}}}}&0 \\ { - 2s{{e}^{{ - s}}}}&{{{e}^{{t + s}}}} \end{array}} \right)x(s)ds = \left( {\begin{array}{*{20}{c}} {{{e}^{t}} + t{{e}^{{ - t}}} + \frac{{{{t}^{{1 - a}}}}}{{1 - a}}{{e}^{t}}} \\ {t{{e}^{t}} + {{t}^{2}}{{e}^{{ - t}}} - \frac{{2{{t}^{{2 - a}}}}}{{2 - a}} + \frac{{{{t}^{{1 - a}}}}}{{1 - a}}{{e}^{t}}} \end{array}} \right),\quad a \in (0,\;1).$

Точное решение: $x(t) = {{({{e}^{t}},{{e}^{{ - t}}})}^{{\text{т}}}}$. Результаты расчетов для $k = 1,\;2,\;3$ приведены в табл. 2–4.

Таблица 2.  

Результаты численных расчетов примера 4.1 одношаговым методом

$h$ $a = 0.1$ $a = 0.5$ $a = 0.999$
${\text{err}}$ $p$ ${\text{err}}$ $p$ ${\text{err}}$ $p$
$0.1$ $0.0094$ 2 $0.0061$ 2 0.0031 2
$0.05$ $0.0023$ 2 0.0013 2 0.00080 2
$0.025$ $0.00055$ 2 0.00027 2 0.00021 2
$0.0125$ $0.00014$ 2 $0.61 \times {{10}^{{ - 4}}}$ 2 $0.53 \times {{10}^{{ - 4}}}$ 2
$0.00625$ $0.34 \times {{10}^{{ - 4}}}$ 2 $0.14 \times {{10}^{{ - 4}}}$ 2 $0.14 \times {{10}^{{ - 4}}}$ 2
Таблица 3.  

Результаты численных расчетов примера 4.1 двухшаговым методом

$h$ $a = 0.1$ $a = 0.5$ $a = 0.999$
${\text{err}}$ $p$ ${\text{err}}$ $p$ ${\text{err}}$ $p$
$0.1$ $0.00080$ 3 $0.00060$ 3 0.00048 3
$0.05$ $0.95 \times {{10}^{{ - 4}}}$ 3 $0.67 \times {{10}^{{ - 4}}}$ 3 $0.32 \times {{10}^{{ - 4}}}$ 3
$0.025$ $0.11 \times {{10}^{{ - 4}}}$ 3 $0.70 \times {{10}^{{ - 5}}}$ 3 $0.39 \times {{10}^{{ - 5}}}$ 3
$0.0125$ $0.14 \times {{10}^{{ - 5}}}$ 3 $0.74 \times {{10}^{{ - 6}}}$ 3 $0.52 \times {{10}^{{ - 6}}}$ 3
$0.00625$ $0.17 \times {{10}^{{ - 6}}}$ 3 $0.81 \times {{10}^{{ - 7}}}$ 3 $0.67 \times {{10}^{{ - 7}}}$ 3
Таблица 4.  

Результаты численных расчетов примера 4.1 трехшаговым методом

$h$ $a = 0.1$ $a = 0.5$ $a = 0.999$
${\text{err}}$ $p$ ${\text{err}}$ $p$ ${\text{err}}$ $p$
$0.1$ $0.00011$ 4 $0.92 \times {{10}^{{ - 4}}}$ 4 $0.68 \times {{10}^{{ - 4}}}$ 4
$0.05$ $0.63 \times {{10}^{{ - 5}}}$ 4 $0.46 \times {{10}^{{ - 5}}}$ 4 $0.21 \times {{10}^{{ - 5}}}$ 4
$0.025$ $0.36 \times {{10}^{{ - 6}}}$ 4 $0.22 \times {{10}^{{ - 6}}}$ 4 $0.12 \times {{10}^{{ - 6}}}$ 4
$0.0125$ $0.22 \times {{10}^{{ - 7}}}$ 4 $0.12 \times {{10}^{{ - 7}}}$ 4 $0.80 \times {{10}^{{ - 8}}}$ 4
$0.00625$ $0.13 \times {{10}^{{ - 8}}}$ 4 $0.64 \times {{10}^{{ - 9}}}$ 4 $0.51 \times {{10}^{{ - 9}}}$ 4

Наиболее близкими объектами к рассматриваемым ИАУ (1) со слабой граничной особенностью являются сингулярные дифференциально-алгебраические уравнения (ДАУ). Для них в работах [6] и [7] предложены коллокационные методы приближенного решения и указаны различные приложения. В отличие от систем обыкновенных дифференциальных уравнений, разрешенных относительно производных, начальные условия для ДАУ должны быть согласованы с правой частью. Одним из подходов к такому согласованию является переход к ИАУ – интегральной форме исходного ДАУ. Приведем численные расчеты для интегрального аналога тестового примера, рассмотренного в статье [7].

Пример 4.2:

$\left( {\begin{array}{*{20}{c}} 1&0 \\ 0&0 \end{array}} \right)x(t) + \int\limits_0^t {{{s}^{{ - a}}}} \left( {\begin{array}{*{20}{c}} 1&0 \\ 0&{coss} \end{array}} \right)x(s)ds = \left( {\begin{array}{*{20}{c}} {tsint + \int\limits_0^t {{{s}^{{1 - a}}}} sinsds} \\ { - \int\limits_0^t {{{s}^{{ - a}}}} {{e}^{{2s}}}ds} \end{array}} \right),\quad a \in (0,\;1).$

Точное решение: $x(t) = {{\left( {tsint, - \tfrac{{{{e}^{{2t}}}}}{{cost}}} \right)}^{{\text{т}}}}$. Результаты расчетов для $k = 1,\;2,\;3$ приведены в табл. 5–7.

Таблица 5.  

Результаты численных расчетов примера 4.2 одношаговым методом

$h$ $a = 0.3$ $a = 0.6$ $a = 0.9$
${\text{err}}$ $p$ ${\text{err}}$ $p$ ${\text{err}}$ $p$
$0.1$ $0.204223$ 2 $0.203065$ 2 0.201907 2
$0.05$ $0.053732$ 2 0.053575 2 0.053419 2
$0.025$ $0.013819$ 2 $0.013799$ 2 0.013778 2
$0.0125$ $0.003507$ 2 $0.003504$ 2 $0.003502$ 2
$0.00625$ $0.000883$ 2 $0.000883$ 2 $0.000883$ 2
Таблица 6.  

Результаты численных расчетов примера 4.2 двухшаговым методом

$h$ $a = 0.3$ $a = 0.6$ $a = 0.9$
${\text{err}}$ $p$ ${\text{err}}$ $p$ ${\text{err}}$ $p$
$0.1$ $0.033079$ 3 $0.032876$ 3 0.032675 3
$0.05$ $0.004598$ 3 $0.004583$ 3 $0.004569$ 3
$0.025$ $0.000607$ 3 $0.000606$ 3 $0.000605$ 3
$0.0125$ $0.78 \times {{10}^{{ - 4}}}$ 3 $0.78 \times {{10}^{{ - 4}}}$ 3 $0.78 \times {{10}^{{ - 4}}}$ 3
$0.00625$ $0.99 \times {{10}^{{ - 5}}}$ 3 $0.99 \times {{10}^{{ - 5}}}$ 3 $0.99 \times {{10}^{{ - 5}}}$ 3
Таблица 7.  

Результаты численных расчетов примера 4.2 трехшаговым методом

$h$ $a = 0.3$ $a = 0.6$ $a = 0.9$
${\text{err}}$ $p$ ${\text{err}}$ $p$ ${\text{err}}$ $p$
$0.1$ $0.005652$ 4 $0.005522$ 4 $0.005476$ 4
$0.05$ $0.407 \times {{10}^{{ - 3}}}$ 4 $0.406 \times {{10}^{{ - 3}}}$ 4 $0.405 \times {{10}^{{ - 3}}}$ 4
$0.025$ $0.28 \times {{10}^{{ - 4}}}$ 4 $0.27 \times {{10}^{{ - 4}}}$ 4 $0.27 \times {{10}^{{ - 4}}}$ 4
$0.0125$ $0.18 \times {{10}^{{ - 5}}}$ 4 $0.18 \times {{10}^{{ - 5}}}$ 4 $0.18 \times {{10}^{{ - 5}}}$ 4
$0.00625$ $0.11 \times {{10}^{{ - 6}}}$ 4 $0.11 \times {{10}^{{ - 6}}}$ 4 $0.11 \times {{10}^{{ - 6}}}$ 4

Результаты численных экспериментов подтверждают полученные выше теоретические результаты: при уменьшении шага в 2 раза погрешность уменьшается в ${{2}^{{k + 1}}}$ раза.

Для исследования устойчивости разработанных алгоритмов к возмущению входных данных были проведены численные расчеты тестового примера 4.1 с пилообразным возмущением правой части, т.е. вместо точного значения $f({{t}_{i}})$ вектор-функции $f$ задано ее возмущенное значение

$\tilde {f}({{t}_{i}}) = f({{t}_{i}}) + {{( - 1)}^{i}}\delta .$
Результаты численных экспериментов иллюстрируют устойчивость разработанных алгоритмов к возмущениям правой части, измеряемым в норме ${{\left| {\, \cdot \,} \right|}_{C}}$ векторного пространства $C([0,\;1];{{\mathbb{R}}^{n}})$:
$A(t)\tilde {x}(t) + \int\limits_0^t {{{s}^{{ - a}}}} K(t,s)\tilde {x}(s)ds = \tilde {f}(t),$
где ${{\left| {\tilde {f}(t) - f(t)} \right|}_{C}} = {{\left| {\delta (t)} \right|}_{C}} \leqslant \delta $, и коррелируют с результатами, полученными в [14], [28]: для $k = 1$ при уменьшении возмущения в 1000 раз, шаг уменьшается в 10 раз, а погрешность решения в 100 раз, т.е. ${{\left| {{{{\tilde {x}}}_{i}} - x({{t}_{i}})} \right|}_{{{{\mathbb{R}}^{n}}}}} = O({{\delta }^{{\tfrac{2}{3}}}})$. В табл. 8, следуя работе [1], значение шага $h$, при котором погрешность ${{\left| {{{{\tilde {x}}}_{i}} - x({{t}_{i}})} \right|}_{{{{\mathbb{R}}^{n}}}}}$ принимает минимальные значения, называется квазиоптимальным шагом и обозначено ${{h}_{{k.o.}}}$.

Таблица 8.  

Результаты численных расчетов одношаговым методом примера 4.1 с пилообразным возмущением правой части

$a = 0.1$ $a = 0.5$ $a = 0.9$
$\delta $ ${{h}_{{k.o.}}}$ $er{{r}_{{k.o.}}}$ $\delta $ ${{h}_{{k.o.}}}$ $er{{r}_{{k.o.}}}$ $\delta $ ${{h}_{{k.o.}}}$ $er{{r}_{{k.o.}}}$
${{10}^{{ - 2}}}$ $0.2154$ $0.040507$ ${{10}^{{ - 2}}}$ $0.2154$ $0.029492$ ${{10}^{{ - 2}}}$ $0.2154$ $0.020908$
${{10}^{{ - 3}}}$ $0.1$ $0.009690$ ${{10}^{{ - 3}}}$ $0.1$ $0.007398$ ${{10}^{{ - 3}}}$ $0.1$ $0.003327$
${{10}^{{ - 4}}}$ $0.0464$ $0.002272$ ${{10}^{{ - 4}}}$ $0.0464$ $0.001485$ ${{10}^{{ - 4}}}$ $0.0464$ $0.000920$
${{10}^{{ - 5}}}$ $0.0215$ $0.000476$ ${{10}^{{ - 5}}}$ $0.0215$ $0.000374$ ${{10}^{{ - 5}}}$ $0.0215$ $0.000197$
${{10}^{{ - 6}}}$ $0.01$ $0.000105$ ${{10}^{{ - 6}}}$ $0.01$ $0.000083$ ${{10}^{{ - 6}}}$ $0.01$ $0.000043$
${{10}^{{ - 7}}}$ $0.0046$ $0.000022$ ${{10}^{{ - 7}}}$ $0.0046$ $0.000015$ ${{10}^{{ - 7}}}$ $0.0046$ $0.000009$

Далее проиллюстрируем область применения полученных результатов на примере решения содержательной задачи определения долгосрочной стратегии ввода генерирующих электромощностей на основе интегральных уравнений Вольтерра с учетом естественного износа оборудования. Следует отметить, что впервые интегральные модели для описания долгосрочной стратегии ввода мощностей электроэнергитической системы (ЭЭС) России были успешно применены А.С. Апарциным и А.М. Тришечкиным [29] в 1986 г. Позднее последовала серия работ по детальному исследованию таких моделей (см., например, статьи [19], [30] и приведенную там библиографию).

Построение долгострочного прогноза ввода мощностей ЭЭС России осуществляется на основе различных типов математических моделей, различающихся предположениями о показателях эффективности и механизмах старения элементов системы вплоть до исключения нерентабельного оборудования. В статье [30] приводится неклассическое уравнение Вольтерра I рода, описывающее упомянутую задачу электроэнергетики. В данной модели выделяются четыре возрастные группы оборудования: оборудование, функционирующее в системе менее 30 лет, считается новым и работает с полной мощностью; оборудование возраста от 31 до 50 лет работает с коэффициентом эффективности $0.97$; оборудование возраста от 51 до 60 лет эксплуатируется на 90 процентов; оборудование, проработавшее более 60 лет, считается устаревшим и выводится из системы. Коэффициенты эффективности каждой возрастной группы были определены методом экспертных оценок. За нуль (год зарождения системы) принят 1950 г.

Пример 4.3. Сформулируем задачу определения долгосрочной стратегии ввода генерирующих электромощностей $x(t)$ в ЭЭС России для достижения роста располагаемой мощности $p(t)$ на 2% на основе интегральных уравнений Вольтерра I рода

$\int\limits_0^t \mu (t,s)\beta (t,s)x(s)ds = p(t),\quad t \in [0,\;100],$
где $\mu (t,s)$ – коэффициент интенсивности использования в момент времени $t$ единицы мощности, введенной ранее в момент времени $s$, а $\beta (t,s)$ – относительная скорость создания новых мощностей в момент времени $t$ в расчете на единицу мощности, введенной в момент времени $s$. Предполагается, что электростанции используются на 100%, т.е. $\mu (t,s) = 1$ при всех $(t,s) \in [0,\;100] \times [0,\;t]$, а безразмерная величина $\beta (t,s)$ изменяется по релаксационному закону
$\beta (t,s) = \mathop {\left( {\frac{t}{s}} \right)}\nolimits^a ,\quad a \in (0,\;1),$
а именно, она бесконечна только в пределе $\beta (t,s) \to + \infty $ при $s \to + 0$ и любом $t \in (0,\;100]$, конечна и убывает с ростом $s$ при каждом фиксированном $t$, выходя на стационарный режим при $s = t$, т.е. $\beta (t,t) = 1$ при всех $t \in [0,\;100]$. Предлагаемый постулат модели соответствует свойствам реального объекта, в частности, правильно описывает как мгновенное поведение ЭЭС (скорость введения генерирующих мощностей конечна), так и долговременный процесс увеличения ее мощности (учитываются затраты энергии на введение новых мощностей). Таким образом, рассматриваемая математическая модель имеет вид уравнения
(5)
$\int\limits_0^t {{{s}^{{ - a}}}} {{t}^{a}}x(s)ds = p(t),\quad t \in [0,\;100]\,,$
которое является частным случаем интегроалгебраического уравнения (1) со слабой граничной особенностью при $n = 1$, $A(t) = 0$, $K(t,s) = {{t}^{a}}$, значит, для проведения расчетов можно применить разработанные выше многошаговые алгоритмы. В условиях $p(t) \in C_{{[0,100]}}^{1}$ и
$\mathop {lim}\limits_{t \to + 0} \frac{{p(t)}}{{{{t}^{a}}}} = 0,$
точным непрерывным на $[0,\;100]$ решением этого уравнения является функция

$x(t) = p{\text{'}}(t) - a\frac{{p(t)}}{t}.$

Для максимального согласования численных результатов с результатами, представленными в [30], при расчетах было выбрано значение $0.0257$ параметра $a$, определяющего коэффициент потери энергии в качестве затрат на введение в ЭЭС новых мощностей. На фиг. 1 отражено сравнение результатов, полученных вводов генерирующих мощностей в ЭЭС России для достижения ежегодного роста располагаемой мощности на 2% с прогнозами, полученными с помощью моделей на основе интегральных уравнений Вольтерра I рода с переменными пределами интегрирования простейшими квадратурными формулами [30] и многошаговыми методами [20]. Видно, что, согласно расчетам на основе модели (5), потребность в ежегодном вводе генерирующих мощностей значительно меньше. Это вполне ожидаемый результат, поскольку (5) подразумевает использование всего ранее введенного оборудования, работающего с меньшей эффективностью, но, тем не менее, вырабатывающего значительные мощности, которые необходимо восполнять в условиях отказа от старого оборудования. Данная ситуация, как раз, описывается интегральными моделями с переменными верхним и нижним пределами интегрирования и соответственно предполагает ежегодно большего ввода генерирующих мощностей. С другой стороны, при кажущейся наиболее выгодной ситуации с необходимостью ввода меньшего количества генерирующих мощностей возникает необходимость учета финансовых затрат на введение нового и обслуживание старого малоэффективного оборудования. На фиг. 2 представлено сравнение расчетных данных с реальными данными вводов генерирующих мощностей до 2014 г., которое наглядно отображает адекватность модельного уравнения (5) и предложенных алгоритмов численного решения для поставленной содержательной задачи.

Фиг. 1.

Стратегия вводов генерирующих мощностей для достижения ежегодного роста располагаемой мощности на 2%, МВт.

Фиг. 2.

Сравнение численных вводов генерирующих мощностей ЭЭС России с реальными данными, МВт.

Авторы статьи выражают огромную благодарность Е.В. Марковой и И.В. Сидлер, сотрудникам Института систем энергетики имени Л.А. Мелентьева СО РАН, за предоставленные значения вводов генерирующих мощностей до 2010 г., что позволило реализовать численные эксперименты для содержательной задачи (5), провести сравнение результатов численных расчетов с прогнозами ввода генерирующих мощностей с 2010 г. по 2050 г. (см. фиг. 1) и с реальными данными вводов генерирующих мощностей до 2014 г. (см. фиг. 2), а также осуществить анализ эффективности использования рассмотренных в статье уравнений.

5. ЗАКЛЮЧЕНИЕ

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

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

  1. Апарцин А.С. Неклассические уравнения Вольтерра I рода: теория и численные методы. Новосибирск: Наука. Сибирская издательская фирма РАН, 1999. 193 с.

  2. Brunner H. Collocation Methods for Volterra Integral and Related Functional Differential Equations. Cambridge: Cambridge Unversity Press, 2004. 597 p.

  3. Biazar J., Ali Asadi M. RBFs for integral equations with a weakly singular kernel // American Journal of Applied Mathematics. 2015. V. 3. 6. P. 250–255.

  4. Mokhtary P. Discrete collocation method for Volterra type weakly singular integral equations with logarithmic kernels // Iranian Journal of Numerical Analysis and Optimization. 2018. V. 8. 2. P. 95–117.

  5. Kolk M., Pedas A. Numerical solution of Volterra integral equations with weakly singular kernels which may have a boundary singularity // Math. Model. and Analysis. 2009. V. 14. P. 79–89.

  6. Korh O., März R., Praetorius D., Weinmüller E. Collocation methods for index 1 DAEs with a singularity of the first kind // Math. of Comput. 2010. V. 79. 269. P. 281–304.

  7. Auzinger W., Lehner H., Praetorius D., Weinmüller E. An efficient asymptotically correct error estimator for collocation solutions to singular index-1 DAEs // Bit Numerical Math. 2011. V. 51. P. 43–65.

  8. Bulatov M.V., Lima P.M., Weinmüller E.B. Existence and uniqueness of solutions to weakly singular integral-algebraic and integro-differential equations // Central European Journal of Math. 2014. V. 12. 2. P. 308–321.

  9. Булатов М.В., Будникова О.С. Численное решение интегро-алгебраических уравнений со слабой особенностью в ядре $k$-шаговыми методами // Известия Иркутского гос. ун-та. Серия Математика. 2015. Т. 13. С. 3–15.

  10. Будникова О.С., Булатов М.В. Численное решение интегро-алгебраических уравнений многошаговыми методами // Ж. вычисл. матем. и матем. физ. 2012. Т. 52. № 5. С. 829–839.

  11. Pishbin S., Ghoreishi F., Hadizadeh M. The semi-explicit Volterra integral algebraic equations with weakly singular kernel: the numerical treatments // J. of Comput. and Appl. Math. 2013. V. 245. № 1. P. 121–132.

  12. Чистяков В.Ф. О сингулярных системах обыкновенных дифференциальных уравнений и их интегральных аналогах. Функции Ляпунова и их применения. Новосибирск: Наука, 1987. С. 231–239.

  13. Чистяков В.Ф. Алгебро-дифференциальные операторы с конечномерным ядром. Новосибирск: Наука. Сибирская издательская фирма РАН, 1996. 279 с.

  14. Будникова О.С. Численное решение интегро-алгебраических уравнений многошаговыми методами. Дис. … канд. физ.-матем. наук. Иркутск, 2015. 125 с.

  15. Balakumar V., Murugesan K. Numerical solution of Volterra integral-algebraic equations using block pulse functions // Appl. Math. and Comput. 2015. V. 263. P. 165–170.

  16. Farahani M.S., Hadizadeh M. Direct regularization for system of integral-algebraic equations of index-1 // Inverse Problems in Science and Engng. 2018. V. 26. 5. P. 728–743.

  17. Farahani M.S., Hadizadeh M., Bulatov M.V., Chistyakova E.V. Adaptive iterative regularization schemes for two-dimensional integral-algebraic systems // Math. Meth. in the Appl. Science. 2019. V. 42. 18. P. 6635–6647.

  18. Sidorov D. Integral dynamical models. Singularities, signals and control. World Scientific Series on Nonlinear Science. Series A: Monographs and Treatises, Vol. 87. Hackensack, New Jersey: World Scientific Publishing, 2015. 243 p.

  19. Маркова Е.В., Сидлер И.В., Труфанов В.В. О моделях развивающихся систем типа Глушкова и их приложениях в электроэнергетике // Автоматика и телемехан. 2011. № 7. С. 20–28.

  20. Ботороева М.Н. Моделирование развивающихся систем на основе интегральных уравнений Вольтерра. Дис. … канд. физ.-матем. наук. Иркутск, 2019. 128 с.

  21. Sidorov D., Muftahov I., Karamov D., Tomin N., Panasetsky D., Dreglea A., Liu F., Foley A. A Dynamic analysis of energy storage with renewable and diesel generation using Volterra equations // IEEE Transactions on Industrial Informatics. 2020. V. 16. 5. P. 3451–3459.

  22. Ботороева М.Н., Булатов М.В. Приложения и методы численного решения одного класса интегро-алгебраических уравнений с переменными пределами интегрирования // Известия Иркутского гос. ун-та. Серия Математика. 2017. Т. 20. С. 3–16.

  23. Булатов М.В. О преобразовании алгебро-дифференциальных систем уравнений // Ж. вычисл. матем. и матем. физ. 1994. Т. 34. № 3. С. 360–372.

  24. Бахвалов Н.С. Численные методы. М.: Наука, 1975. 632 с.

  25. Weiss R., Anderssen R.S. A product integration methods for a class of singular first kind Volterra equations // Numerische Mathematik. 1972. V. 18. 5. P. 442–456.

  26. Тен Мен Ян. Приближенное решение линейных интегральных уравнений Вольтерра I рода. Дис. … канд. физ.-матем. наук. Иркутск, 1985. 215 с.

  27. Linz P. Analytical and Numerical Methods for Volterra Equations. Philadelphia: SIAM, 1985. 227 p.

  28. Булатов М.В., Будникова О.С. Об устойчивых алгоритмах численного решения интегро-алгебраических уравнений // Вестник Южно-Уральского гос. университета. Серия Матем. моделирование и программирование. 2013. Т. 6. № 4. С. 5–14.

  29. Апарцин А.С., Тришечкин А.М. Применение моделей В.М. Глушкова для моделирования долгосрочных стратегий развития ЕЭЭС // Тезисы докл. Всесоюзной конференции “Курс-4”. Рига: ЛГУ, 1986. С. 17–19.

  30. Апарцин А.C., Сидлер И.В. Применение неклассических уравнений Вольтерра I рода для моделирования развивающихся систем // Автоматика и телемехан. 2013. № 6. С. 3–16.

Дополнительные материалы отсутствуют.