Журнал вычислительной математики и математической физики, 2022, T. 62, № 1, стр. 105-112

ВВЕДЕНИЕ

Численные методы решения различных видов интегральных уравнений являются классическим разделом вычислительной математики. К настоящему времени вышло и продолжает выходить огромное число публикаций по данной тематике. Среди этих классов выделим интегральные уравнения Вольтерра. Прекрасный исторический обзор и обширная библиография приведены в статье [1] и монографиях [2]–[8] и др.

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

Для численного решения ИУВ II рода с достаточно гладкими входными данными можно применять ряд алгоритмов, в основу которых положены известные квадратурные формулы (см., например, [2], [5], [9]). Иное дело состоит с ИУВ I рода (ИУВ I). Для создания и обоснования методов их решения важную роль играет понятие их степени неустойчивости [8].

К настоящему времени разработаны методы решения ИУВ I только для частных случаев (у которых степень неустойчивости не превосходит двух).

В первом разделе статьи показано, что понятие степени неустойчивости для ИУВ I [8] и понятие индекса для дифференциально-алгебраических уравнений (ДАУ) [10], [11], [14] близки. Поэтому, на наш взгляд, некоторые подходы, развитые для численного решения ДАУ, можно применять и для численного решения ИУВ I.

Этот факт инициировал проведение исследований данной статьи. Здесь предложены новые одношаговые алгоритмы численного решения ИУВ I со степенью неустойчивости, равной единице.

1. ПОСТАНОВКА ЗАДАЧИ И ОПРЕДЕЛЕНИЯ

Рассмотрим интегральное уравнение Вольтерра I рода вида

Характеристикой сложности уравнения (1) является понятие степени неустойчивости [8].

Определение 1 (см. [8, c. 14]). Пусть существуют взаимно однозначно ограниченный из ${{C}_{{[0,1]}}}$ в ${{C}_{{[0,1]}}}$ оператор $\mathcal{B}$ и $\mathcal{D}$ – линейный дифференциальный оператор порядка $\xi $ такие, что $\mathcal{B}\mathcal{D}$ является обратным оператором для исходного уравнения $Vx = f$, т.е. $\mathcal{B}\mathcal{D} = {{V}^{{ - 1}}}$. Число $\xi $ называется степенью неустойчивости решения (1).

В частности, если ядро и правая часть (1) достаточно гладкие функции, то после $\xi $-кратного дифференцирования (1) мы получим уравнение вида

где функция ${{\left. {K_{{{{t}^{{\xi - 1}}}}}^{{(\xi - 1)}}(t,\tau )} \right|}_{{\tau = t}}} \ne 0$ $\forall t \in [0,1]$, т.е. уравнение (2) является уравнением II рода.

Если в (1) $K(t,\tau ) = {{(t - \tau )}^{{ - \alpha }}}k(t,\tau )$, где $\alpha \in (0,1)$ и функция $K(t,\tau ):k(t,t \ne 0)$ $\forall t \in [0,1]$, то степень неустойчивости равна $\xi = 1 - \alpha $. Такие уравнения принято называть уравнениями со слабой особенностью.

Создание алгоритмов численного решения уравнения (1) имеет специфические особенности. Так, например, ряд многошаговых методов будет неустойчивым для случая $\xi = 1$ и устойчивым для $\xi = 2$. К настоящему времени развиты многошаговые методы для решения уравнения (1) с $0 \leqslant \xi \leqslant 1$ и для $\xi = 2$.

Отметим, что ИУВ I с конечным, целым показателем степени неустойчивости $\xi $ по своей сути близки к ДАУ индекса $\xi + 1$. Понятие индекса для линейных ДАУ вида

где $A(t)$, $B(t)$ суть $(n \times n)$-матрицы, $y(t)$ – искомая, $\varphi (t)$ – заданная $n$-мерные вектор-функции в предположении, что элементы входных данных достаточно гладкие и начальное условие задано корректно, введено различными авторами по-разному [11]. В монографии [12] показано, что в случае вещественно-алгебраических коэффициентов матриц $A(t)$ и $B(t)$ в (3) эти определения эквивалентны.

Приведем наиболее близкое к определению степени неустойчивости (1) понятие индекса для ДАУ (3).

Определение 2 (см. [13]). Минимальное $r$, при котором существует линейный дифференциальный оператор

где ${{P}_{j}}(t)$ суть $(n \times n)$-матрицы такой, что где $D(t)$ есть $(n \times n)$-матрица, называют индексом ДАУ (3).

Уравнения (1) и (3) в некотором смысле близки. Например, задачу восстановления $r$-й производной функции $\psi (t)$ можно задать как в виде уравнения (1):

где $\psi (t)$ задана корректно, так и в виде ДАУ (3): где $y$ – вектор-функция размерности $r + 1$. Решением задачи (4) является $x(t) = {{\psi }^{{(r)}}}(t)$, где операторами $\mathcal{B}$ и $\mathcal{D}$ из определения 1, в частности, являются $\mathcal{B} = E$ – единичный оператор, $\mathcal{D} = \mathop {\left( {\tfrac{d}{{dt}}} \right)}\nolimits^r $. Решением (5), которое не зависит от начальных данных, является вектор-функция, определенная по правилу [14]

Оператор ${{\mathcal{P}}_{r}}$ из определения 2 имеет вид

Если (1) имеет степень неустойчивости $\xi = 1$, то, действуя на это уравнение оператором $\frac{1}{{K(t,t)}}\frac{d}{{dt}} + \alpha $, где $\alpha $ – скалярный параметр, получаем ИУВ II рода.

Если у ДАУ (3) индекс равен единице, то, действуя на эту систему оператором $(E - A(t){{A}^{ + }}(t))\frac{d}{{dt}} + E$, где $E$ – единичная матрица, ${{A}^{ + }}(t)$ – псевдообратная матрица к $A(t)$, получим систему ОДУ первого порядка с невырожденной матрицей перед производной [11], [12], [14]. Отметим, что слагаемое, содержащее вторую производную, будет тождественно нулевым в силу свойства псевдообратной матрицы: $(E - A(t){{A}^{ + }}(t))A(t) \equiv 0$.

Эти факты, а также примеры (4) и (5) показывают некоторое сходство ИУВ I (1) и ДАУ (3). Поэтому можно надеяться, что методы, предложенные для численного решения (3), можно применять и для (1).

В работах [15]–[19] для численного решения ДАУ индекса не выше двух предложены, и в ряде случаев обоснованы, коллокационно-вариационные методы. Данные методы имеют целый ряд преимуществ над ранее разработанными:

а) они не требуют вычисления проекторов и производных (см., например, [10]–[14], [20] и приведенную там библиографию);

б) просты в программной реализации;

в) в отличие от методов Рунге–Кутты не требуют весьма жестких ограничений на структуру матричного пучка $\lambda A(t) + B(t)$.

Приведем описание, исследование и анализ семейства коллокационно-вариационных методов второго порядка с использованием одной точки коллокации.

2. КОЛЛОКАЦИОННО-ВАРИАЦИОННЫЙ ПОДХОД

В данном разделе мы ограничимся частным случаем коллокационно-вариационного подхода для уравнения (1), имеющего степень неустойчивости $\xi = 1$.

Это означает, что $K(t,t) \ne 0$ $\forall t \in [0,1]$, $f(0) = 0$, $K(t,\tau )$, $K_{t}^{'}(t,\tau )$, $f{\kern 1pt} '(t)$ – непрерывные функции в области определения.

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

Зададим на отрезке $[0,1]$ равномерную сетку ${{t}_{i}} = ih,$ $i = \overline {0,N} ,$ $h = 1{\text{/}}N$, положим $N$ четным, ${{x}_{0}} = x(0)$ – заданным и обозначим ${{K}_{{ij}}} = K({{t}_{i}},{{t}_{j}})$, ${{f}_{i}} = f({{t}_{i}})$, ${{x}_{i}} \approx x({{t}_{i}}),$ $j = \overline {0,i} $. Для вычисления $\int_{{{t}_{i}}}^{{{t}_{i}} + 2} {\varphi (t)dt} $ будем использовать квадратурную формулу

где коэффициенты ${{a}_{0}}$, ${{a}_{1}}$, ${{a}_{2}}$ удовлетворяют условиям порядка (являются решением системы) [21]

Формула (7) для уравнения (1) дает

Соотношения (9) дают недоопределенную систему линейных алгебраических уравнений (СЛАУ) относительно ${{({{x}_{1}},{{x}_{2}}, \ldots ,{{x}_{{i + 2}}})}^{{\text{т}}}}$ вида

где $X = {{({{x}_{1}},{{x}_{2}}, \ldots ,{{x}_{{i + 2}}})}^{{\text{т}}}}$, $\mathcal{B} = {{({{f}_{2}} - {{a}_{0}}{{K}_{{20}}}{{x}_{0}},{{f}_{4}} - {{a}_{0}}{{K}_{{42}}}{{x}_{2}}, \ldots ,{{f}_{{i + 2}}} - {{a}_{0}}{{K}_{{i + 2,i}}}{{x}_{i}})}^{{\text{т}}}}$,

Размерность СЛАУ (10) на каждом шаге интегрирования $(i{\text{/}}2 + 1) \times (i + 2)$, $i = 0,2,4, \ldots ,N - 2$. Следовательно, формулы (10) применять для численного решения исходной задачи (1) нельзя. Поэтому поступим следующим образом: на каждом шаге интегрирования дополняем уравнения (10) условием минимума квадрата нормы приближенного решения на отрезках $[0,2h]$, $[2h,4h]$, …, $[(N - 2)h,Nh]$. Под нормой мы будем понимать некоторый дискретный аналог нормы в пространстве Соболева. Выбор данной нормы приведем ниже.

Итак, предложенный подход сводится к решению задачи математического программирования (МП), где целевая функция – это квадрат нормы приближенного решения, а ограничения типа равенств – это квадратурные формулы (10). Методы будут существенным образом зависеть от выбора этой нормы приближенного решения квадратурной формулы (7) (от выбора коэффициентов в (8)).

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

Будем под квадратом нормы приближенного решения (целевой функцией) понимать норму в некотором аналоге пространства Соболева. Пусть ${{L}_{2}}(t,{{x}_{i}},{{x}_{{i + 1}}},{{x}_{{i + 2}}})$, $t \in [{{t}_{i}},{{t}_{{i + 2}}}]$, – интерполяционный полином второй степени, проходящий через точки $({{t}_{i}},{{x}_{i}})$, $({{t}_{{i + 1}}},{{x}_{{i + 1}}})$, $({{t}_{{i + 2}}},{{x}_{{i + 2}}})$. Квадрат нормы приближенного решения на отрезке $[{{t}_{i}},{{t}_{{i + 2}}}]$ будем определять по формуле

Для вычисления определенного интеграла (11) применим одну из известных квадратурных формул. Получим

где ${{\alpha }_{m}}$, ${{\beta }_{m}}$, $m = 0,1,2$, зависят от выбора квадратурной формулы.

Полагая ${{x}_{i}}$ известным (${{x}_{0}}$ задано), на каждом шаге будем иметь задачу МП:

при ограничениях типа равенств (см. формулу (9))

Задачи (12)–(14) (задачи на условный экстремум) можно решить известным методом множителей Лагранжа.

Упростим задачу (12)–(14). Так как постоянный множитель не влияет на поиск аргумента условного минимума, то вместо целевой функции (12) можно минимизировать функцию

В силу того, что множитель $h$ является малым, то в выражении (15) можно отбросить первое или первое и второе слагаемые. В этих случаях получим целевые функции

Замечание. Для единичного ядра в случае задачи (13)–(14), (17) при условии (8) решение не зависит от выбора коэффициентов.

Покажем это для первого шага. Пусть ${{K}_{{ij}}} \equiv 1,\;i,j = \overline {0,N} $. Положим в (14) $i = 0$ и решим систему

Выразим из второго уравнения ${{x}_{2}} = 2{{x}_{1}} - {{x}_{0}}$, подставим его в первое уравнение, учитывая, что ${{x}_{0}}$ известно (например, ${{x}_{0}} = f'(0)$) и ${{a}_{0}} = {{a}_{2}}$, получим и из (8)

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

В следующем разделе приведем результаты расчетов для модельных примеров при выборе целевой функции в виде (17).

3. ЧИСЛЕННЫЕ ПРИМЕРЫ

Приведем разностную схему для задачи (13)–(14), (17) в удобном для программной реализации виде. Используем введенную равномерную сетку узлов и на каждом шаге $i = \overline {1,N{\text{/}}2} $ решаем следующую систему уравнений относительно ${{x}_{{2i - 1}}}$ и ${{x}_{{2i}}}$:

Таким образом, полагая ${{x}_{0}}$ известным (например, ${{x}_{0}} = \tfrac{{f{\kern 1pt} '(0)}}{{K(0,0)}}$), на первом шаге находим ${{x}_{1}}$ и ${{x}_{2}}$. Используя найденное значение ${{x}_{2}}$, на втором шаге находим ${{x}_{3}}$ и ${{x}_{4}}$ и т.д.

Рассмотрим применение предлагаемого метода на известных из литературы примерах. В качестве коэффициентов возьмем два варианта, удовлетворяющих (8): ${{a}_{0}} = {{a}_{2}} = 1{\text{/}}3$, ${{a}_{1}} = 4{\text{/}}3$ (I) (аналог метода Симпсона) и ${{a}_{0}} = {{a}_{2}} = 1{\text{/}}2$, ${{a}_{1}} = 1$ (II).

Пример 1 (см. [2, c. 149]).

Пусть задано уравнение вида

с точным решением $x(t) = {{e}^{{\alpha t}}}$, $T = 1$. В табл. 1 и 2 приведены результаты численных расчетов для двух вариантов выбора коэффициентов ${{a}_{i}}$, где ${{\varepsilon }_{i}} = \mathop {max}\limits_{1 \leqslant j \leqslant N} \left| {x({{t}_{j}}) - {{x}_{j}}} \right|$ – значения погрешностей численного решения для $i$-го набора коэффициентов, ${{x}_{j}}$ – сеточное решение в точке ${{t}_{j}}$, $Nh = T$, ${{p}_{i}}$ – порядок точности, измеренный с использованием шагов $h$ и $h{\text{/}}2$ ($h{\text{/}}10$), для $i$-го набора коэффициентов:

Из таблиц видно, что предлагаемый метод сходится со вторым порядком.

Приведем расчеты с этими же наборами коэффициентов для другого примера.

Пример 2 (см. [2, c. 468]).

Пусть задано уравнение вида

с точным решением $x(t) = 1$, $T = 1$. В табл. 3 и 4 приведены результаты численных расчетов.

Из таблиц видно, что предложенные алгоритмы имеют скорость сходимости не ниже квадратичной.

Пример 3.

Пусть задано уравнение вида

с точным решением $x(t) = {{e}^{{ - \lambda t}}}$, $T = 1$. В табл. 5 и 6 приведены результаты численных расчетов при значениях параметра $\lambda = 1$ и $\lambda = 10$ для двух наборов коэффициентов ${{a}_{i}}$. Данное уравнение при достаточно больших положительных $\lambda $, или при положительных $\lambda $, но на большом отрезке интегрирования $[0,T]$ ($T$ достаточно большое), является жестким. Предложенные алгоритмы для такого вида уравнений требуют существенного ограничения на шаг интегрирования, что, в свою очередь, приводит к увеличению вычислительных затрат.

ЗАКЛЮЧЕНИЕ

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

Для жестких задач предполагается создать коллокационно-вариационные алгоритмы, в основу которых заложены экстраполяционные квадратурные формулы. В работах [22], [23] были исследованы экстраполяционные методы (без условия минимизации квадрата нормы приближенного решения). Эти алгоритмы прекрасно справлялись с жесткими задачами типа (20).

Теоретические исследования и тестирование многочисленных примеров показали перспективность данного подхода.