Журнал вычислительной математики и математической физики, 2019, T. 59, № 7, стр. 1100-1107

1. ВВЕДЕНИЕ

Рассмотрим систему обыкновенных дифференциальных уравнений (ОДУ)

где $A(t)$, $B(t)$ есть $(n \times n)$-матрицы, $f(t)$ и $x(t)$ – заданная и искомая n-мерные вектор-функции.

В том случае, если

то систему (1.1) принято называть дифференциально-алгебраическими уравнениями (ДАУ). Всюду в дальнейшем изложении предполагается, что элементы $A(t)$, $B(t)$ и $f(t)$ обладают достаточной гладкостью, необходимой для проведения выкладок, и начальные условия (1.2) согласованы с правой частью, т.е. задача (1.1) имеет достаточно гладкое решение. Под решением задачи (1.1) будем понимать вектор-функцию, которая обращает (1.1) в тождество и удовлетворяет начальному условию (1.2).

Характеристикой сложности рассматриваемых задач является понятие индекса, которое имеет несколько определений. Приведем одно из них (см., например, [1]).

Определение 1. Если $r$ – наименьшее число аналитических дифференцирований (1.1), необходимых для того, чтобы из данной системы путем алгебраических преобразований выделить систему ОДУ, разрешенную относительно производной, то r будем называть индексом системы (1.1).

Если элементы матриц $A(t)$ и $B(t)$ – вещественно-аналитические функции, то данное определение эквивалентно следующему (см. [2], [3]). Пусть существуют неособенные для любого $t \in [0,1]$ матрицы $P(t)$ и $Q(t)$ с вещественно-аналитическими коэффициентами такие, что

где ${{E}_{k}}$ – единичная матрица размером $k$, $N$ – верхнетреугольная матрица с нулевыми квадратными блоками на диагонали и индексом нильпотентности $r \geqslant 1$, т.е. ${{N}^{r}}(t) \equiv 0,\;{{N}^{{r - 1}}}(t)\not { \equiv }0$ для $t \in [0,1]$. Умножая (1.1) на $P(t)$ и производя замену $x(t) = Q(t)y(t)$, получаем где ${{(g_{1}^{{\text{т }}}(t),g_{2}^{{\text{т }}}(t))}^{{\text{т }}}} = P(t)f(t)$, $(y_{1}^{{\text{т }}}(t),y_{2}^{{\text{т }}}(t)) = {{Q}^{{ - 1}}}(t)x(t)$. В этом случае ${{y}_{1}}(t)$ является общим решением ОДУ вида $y_{1}^{'} + J(t){{y}_{1}} = {{g}_{1}}(t)$, а ${{y}_{2}}$ находят по формуле [4]: где действие оператора $T$ на вектор-функцию ${{g}_{2}}(t)$ определено по правилу

Другие определения индекса приведены в [1]–[5].

Далее приведем примеры, иллюстрирующие сложности численного решения ДАУ.

Пример 1 (см. [3]):

где $\alpha $ – числовой параметр. При любых значениях $\alpha $ данный пример имеет единственное решение индекса два $u(t) = g(t) - \alpha tv,$ $v(t) = f(t) - g{\text{'}}(t)$, не зависящее от начальных условий.

Матрицы $P(t)$ и $Q(t)$, приводящие данный пример к виду (1.4), равны $P(t) = E$, $Q(t) = \left( {\begin{array}{*{20}{c}} { - \alpha t}&1 \\ 1&0 \end{array}} \right)$ соответственно.

Если мы применим неявный метод Эйлера, то анализ, который мы опускаем для краткости, показывает, что метод неустойчив для $\alpha < - \tfrac{1}{2}$. Если $\alpha = - 1$, данная схема неприменима, потому что $\det ({{A}_{{i + 1}}} + h{{B}_{{i + 1}}}) \equiv 0$.

Пример 2 (см. [5]):

$u(0) = \tfrac{1}{{\mu - 1}}$, $v(0) = 1$, $\mu < 0$ – скалярный параметр. Точное решение данного примера: При $\mu \ll 0$ данный пример становится жестким. Применяя для данного примера неявную схему Эйлера, которая весьма успешно справляется с решением жестких ОДУ, получаем где $z = \mu h$. Таким образом, данный алгоритм для устойчивости требует ограничения на шаг интегрирования $h \in \left( {0; - \tfrac{2}{\mu }} \right)$.

Другие примеры линейных ДАУ, численное решение которых требует, в отличие от жестких ОДУ, существенного ограничения на шаг интегрирования, можно найти в [5], [6].

В статье [7], используя результаты [8], приведен пример жесткой линейной системы ОДУ, для численного решения которой ряд неявных методов РК (неявный метод Эйлера, Лобатто IIIC и др.) также требовал ограничения на шаг интегрирования.

Таким образом, разработка алгоритмов численного решения жестких ДАУ является достаточно актуальной задачей. Именно это обстоятельство инициировало продолжение исследований, начатых в [9]–[11]. В данной статье мы предлагаем переписать (1.1) в виде $(A(t)x(t)){\kern 1pt} {\text{'}} + (B(t) - A{\text{'}}(t))x(t)$ = f(t) и уже для такой записи применять блочные многошаговые методы.

2. БЛОЧНЫЕ МЕТОДЫ

Перепишем исходную систему (1.1) в виде

где $C(t) = B(t) - A{\kern 1pt} '(t)$.

Приведем построение интерполяционного варианта таких схем для системы (2.1) с условием (1.3). Пусть $L({{y}_{0}},{{y}_{1}},...,{{y}_{{k + s - 1}}},t)$ – интерполяционный полином степени $k + s - 1$, проходящий через точки $({{y}_{0}},{{t}_{0}})$, $({{y}_{1}},{{t}_{1}})$, … , $({{y}_{{k + s - 1}}},{{t}_{{k + s - 1}}})$, $t \in [{{t}_{0}},{{t}_{{k + s - 1}}}]$, ${{t}_{j}} = jh$, $j = 0,1,...,k + s - 1$.

Тогда справедливо представление

Пусть заранее заданы (вычислены) стартовые значения

Основываясь на формулах (2.2), выпишем блочные разностные схемы для решения задачи (2.1) с условием (2.3):

Приведем частные случаи схем вида (2.4). При $s = 1$ получим методы, основанные на формуле дифференцирования назад [11]. Известно, что при $k > 6$ эти методы неустойчивы (см. [1], [3]). При $s = 2$ и $k = 2$ будем иметь

Обоснование данных схем проведем для частного случая рассматриваемой задачи, а именно, когда матрицы $A(t)$ и $C(t)$ имеют блочный вид

где элементы матриц $\mathcal{A}(t)$, $\mathcal{C}(t)$ и вектор-функции $f(t)$ предполагаются достаточно гладкими. Обозначим ${{\epsilon }_{j}} = x({{t}_{j}}) - {{x}_{j}}$, тогда уравнение ошибки схемы (2.4) имеет вид где представляет собой погрешность численного дифференцирования. В силу предположения о достаточной гладкости входных данных и формул численного дифференцирования [12] справедлива оценка

Положим $k$ кратным $s$. В противном случае искусственно увеличиваем $k$ до минимального числа, кратного $s$, и полагая $\rho _{m}^{l} = 0$, $m = k + 1,...,q$, $l = 1,2,...,s$, вводим обозначения

где $p = 0,1,...,q,$ С учетом данных обозначений перепишем (2.7) в виде $i = q,q + s,...,N - s$. В силу того, что матрицы $\mathcal{M}_{i}^{p}$ состоят из квадратных матриц с двумя нулевыми блоками на диагонали (см. формулу (2.6)), то где ${{\bar {S}}_{i}}$ – блочно-диагональная матрица $\operatorname{diag} (\bar {S}_{i}^{1},\bar {S}_{i}^{2},...,\bar {S}_{i}^{s})$, и Учитывая этот факт, формулы (2.6), (2.8), (2.11) и (2.12), получаем Обозначим ${{\zeta }_{i}} = {{(\vartheta _{i}^{{\text{т }}},\vartheta _{{i - 1}}^{{\text{т }}},...,\vartheta _{{i - q + 1}}^{{\text{т }}})}^{{\text{т }}}}$ и перепишем (2.14) в виде одношагового соотношения ${{\zeta }_{i}} = {{\mathcal{R}}_{i}}{{\zeta }_{{i - 1}}} + {{\omega }_{i}},$ где ${{\omega }_{i}} = {{({{\eta }_{i}},0,0,...,0)}^{{\text{т }}}}$. Здесь 0 есть $(sn \times sn)$ – нулевые матрицы.

Хотя $\left\| {{{\mathcal{R}}_{i}}} \right\|$ неограничена, опуская выкладки, можно показать, что

Таким образом, для погрешности получим

Отсюда, с учетом формул (2.13), (2.15), оценок (2.16), формул погрешности численного дифференцирования (2.9), вытекает $\left\| \zeta \right\| = O({{h}^{{k + s - 1}}})$.

Таким образом, доказана

Лемма 1. Пусть матрицы $A(t)$ и $C(t)$ в задаче (2.1) имеют вид (2.6), элементы матриц $\mathcal{A}(t)$, $\mathcal{C}(t)$ и $f(t)$ принадлежат классу $C_{{[0,1]}}^{{k + s}}$ и начальные условия заданы в виде $\left\| {{{x}_{j}} - x({{t}_{j}})} \right\| = O({{h}^{{k + s}}})$, $j = 0,1,...,k$. Тогда разностные схемы (12) сходятся к точному решению исходной задачи со скоростью $O({{h}^{{k + s - 1}}})$.

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

где $B(t)$ есть $(n \times n)$-матрица и предполагается, что элементы $B(t)$ и $f(t)$ и обладают той гладкостью, которая необходима для дальнейших рассуждений.

Если в разностных схемах (2.4) зафиксировать $s$ – число стадий, то возникает вопрос: при каких ограничениях на $k$ данные алгоритмы будут устойчивы для задачи (2.17)?

Приведем стандартную схему исследования на предмет устойчивости методов (2.4) для задачи (2.17). Для этого достаточно рассмотреть очень простой пример:

для которого выпишем матрицу перехода от $i$-го к $i + s + 1$-му шагу. Введем обозначения Положим $\det {{\mathcal{U}}_{0}}$ не равным нулю и представим данные методы как одношаговые. Тогда матрица перехода от шага к шагу будет иметь вид

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

Определение 2. Блочный метод (2.4) называется устойчивым для задачи (2.17), если все собственные значения матрицы $\mathcal{K}$ лежат внутри единичного круга и на границе круга нет кратных собственных чисел.

Лемма 2. Пусть элементы матрицы $B(t)$ и вектор-функции $f(t)$ принадлежат классу $C_{{[0,1]}}^{{k + s}}$ и стартовые значения заданы в виде $\left\| {x({{t}_{j}}) - x(j)} \right\| = O({{h}^{{k + s}}})$, $j = 0,1,\;...,\;k - 1$. Если методы (2.4) устойчивы, то данные схемы сходятся к точному решению (2.17) и справедлива оценка $\left\| {x({{t}_{j}}) - x(j)} \right\| = O({{h}^{{k + s - 1}}})$, $j = k,k + 1,\;...,\;N$, $N = 1{\text{/}}h$.

Мы проводим доказательство леммы 2, используя стандартную технику для многошаговых схем, записанных как одношаговые методы, и теорему Лакса. Для краткости изложения доказательство не приводится.

Например, для схемы (2.5) матрица перехода от $({{x}_{{i - 1}}},{{x}_{i}})$ к $({{x}_{{i + 1}}},{{x}_{{i + 2}}})$ имеет вид

Собственные числа данной матрицы $\left| {{{\lambda }_{1}}} \right| = 1$, $\left| {{{\lambda }_{1}}} \right| = \tfrac{1}{{23}}$, т.е. метод (2.5) будет устойчивым для задачи $x{\kern 1pt} {\text{'}} + B(t)x = f(t)$ и сходится к точному решению с третьим порядком, если элементы матриц $B(t)$ и вектор-функции $f(t)$ принадлежат классу $C_{{[0,1]}}^{4}$ и $\left\| {{{x}_{1}} - x(h)} \right\| = O({{h}^{4}})$.

Относительно сходимости данных схем к точному решению (1.1), (1.2) справедлива

Теорема 1. Пусть для задачи (1.1), (1.2) выполнены условия.

1) ДАУ имеет индекс не выше двух;

2) элементы матриц $A(t)$, $B(t)$ и $f(t)$ принадлежат классу $C_{{[0,1]}}^{{k + s}}$;

3) матрица $P(t) = P$ – постоянная в формуле (1.4);

4) блочные схемы (2.4) устойчивы и собственные числа матрицы $\mathcal{K}$ лежат в единичном круге и на границе круга нет кратных собственных чисел;

5) стартовые значения заданы следующим образом $\left\| {x({{t}_{j}}) - x(j)} \right\| = O({{h}^{{k + s}}})$.

Тогда блочные схемы (2.4) сходятся, и справедлива оценка $\left\| {x({{t}_{j}}) - x(j)} \right\| = O({{h}^{{k + s - 1}}})$.

Доказательство. Умножим (2.4) на матрицу $P$ из формулы (1.4) и произведем замену переменной ${{x}_{j}} = {{Q}_{j}}{{y}_{j}}$, где ${{Q}_{j}} = Q({{t}_{j}})$. Будем иметь

В силу условий теоремы матрица $P{{A}_{{i + s - j}}}{{Q}_{{i + s - j}}}$ имеет блочный вид:

Таким образом, из данной формулы относительно ${{y}_{q}} = {{(y_{{1q}}^{{\text{т }}},y_{{2q}}^{{\text{т }}})}^{{\text{т }}}}$ будем иметь разностные схемы где ${{\Psi }_{{i + m}}}$ – вектор, состоящий из первых $d$ элементов $P{{f}_{{i + m}}}$, $m = 1,2, \ldots ,s$.

Оценка $\left\| {{{y}_{{1j}}} - {{y}_{1}}({{t}_{j}})} \right\| = O({{h}^{{k + s - 1}}})$ вытекает из второго, четвертого и пятого условия теоремы и леммы 2. Оценка $\left\| {{{y}_{{2j}}} - {{y}_{2}}({{t}_{j}})} \right\| = O({{h}^{{k + s - 1}}})$ вытекает из второго и пятого условия теоремы и леммы 1. Вспоминая что, ${{y}_{j}} = {{(y_{{1j}}^{{\text{т }}},y_{{2j}}^{{\text{т }}})}^{{\text{т }}}} = Q_{j}^{{ - 1}}{{x}_{j}}$, получаем

Теорема доказана.

Отметим, что третье условие теоремы 1 нельзя ослабить. В случае переменной матрицы $P(t)$ нельзя гарантировать сходимость схем (2.4) . Для иллюстрации вышесказанного приведем простой пример ДАУ индекса два с матрицами

Умножим это уравнение на матрицу

где $\alpha $ – числовой параметр, и полученную систему запишем в виде (2.1).

Имеем

Если $\alpha = 1$ и $s = 1$ (что соответствует основанной на ФДН-методе многошаговой схеме в [11]), то тогда получим систему линейных алгебраических уравнений относительно ${{x}_{{i + 1}}}$ со следующей матрицей коэффициентов:

Очевидно, что эта матрица является вырожденной для всех ${{t}_{{i + 1}}}$, т.е. схема неприменима.

Можно также показать, что мы получаем тождественно вырожденную систему линейных алгебраических уравнений относительно ${{(x_{{i + 2}}^{{\text{т }}},x_{{i + 1}}^{{\text{т }}})}^{{\text{т }}}}$ для схемы (2.5) при $\alpha = 1{\text{/}}2$.

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

В данном разделе приведем результаты численных расчетов для предложенных схем.

Для ДАУ (1.7) в примере 1 на интервале $[0,1]$ применим метод (2.5), который согласно теореме 1 сходится c третьим порядком.

Мы взяли $u(t) = exp(t) - \alpha t(exp(0.5t) - exp(t))$, $v(t) = exp(0.5t) - exp(t)$ в качестве точного решения. В этом случае правая часть (1.7) определена следующим образом $f(t) = exp(0.5t)$, $g(t) = exp(t)$. Погрешности численного решения для значения параметра $\alpha = - 0.6$ и различных значений шага интегрирования приведены в табл. 1. Достаточно очевидно, что результаты численного эксперимента подтверждают порядок сходимости метода.

Здесь, как и в следующем примере, ${\text{er}}{{{\text{r}}}_{u}} = {{\max }_{{1 \leqslant i \leqslant N}}}\left| {{{u}_{i}} - u({{t}_{i}})} \right|$, ${\text{er}}{{{\text{r}}}_{v}} = {{\max }_{{1 \leqslant i \leqslant N}}}\left| {{{v}_{i}} - v({{t}_{i}})} \right|$.

Далее мы определяем $\mu = - {{10}^{5}}$ для ДАУ (1.8) в примере 2, который представляет собой жесткое ДАУ. Напомним, что неявный метод Эйлера, примененный к этому примеру, устойчив только для $h < 2 \times {{10}^{{ - 5}}}$. Результаты численных расчетов по методу (2.5), которые подтверждают сходимость, отражены в табл. 2.

Таким образом, мы можем сделать вывод, что численные результаты согласуются с теоретическими результатами.

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

В данной статье предложено семейство блочных многостадийных многошаговых разностных схем, которые обобщают известные методы ФДН для численного решения линейных ДАУ индекса не более $2$. Сформулированы условия сходимости предложенных схем, которые, в отличие от ранее рассмотренных авторами методов ФДН, могут иметь более высокий порядок. Кроме того, блочные схемы при $k = 1$ являются одношаговыми, что открывает возможности построения алгоритмов с автоматическим выбором шага, разработанных для жестких ОДУ. Преимущества предлагаемых методов продемонстрированы на численных расчетах тестовых примеров индекса два и жесткого ДАУ индекса один. Для данных примеров ранее разработанные алгоритмы либо были неустойчивыми, либо принципиально неприменимы, либо требовали существенного ограничения на шаг интегрирования.