Программирование, 2022, № 4, стр. 50-56

1. ВВЕДЕНИЕ

Линейные разностные уравнения возникают в различных областях математики, следовательно, поиск решений таких уравнений является одной из математических задач, имеющей многочисленные приложения в различных областях науки и техники. Например, разностные уравнения часто используются в моделях динамики с дискретным временем [1, 2], а также для приближенного решения дифференциальных уравнений [3], в комбинаторном анализе разностные уравнения в сочетании с методом производящих функций дают мощный аппарат исследования перечислительных задач (см., например, [4–6]). Алгоритм вычисления производящей функции решения задачи Коши для двумерного разностного уравнения с постоянными коэффициентами по коэффициентам разностного уравнения и начальным данным задачи Коши представлен в работе [7]. В работе [8] представлен алгоритм вычисления решения задачи Коши для двумерного разностного уравнения с постоянными коэффициентами в точке по коэффициентам разностного уравнения и начальным данным задачи Коши. В данной работе представлен алгоритм вычисления решения задачи Коши для двумерного разностного уравнения с постоянными коэффициентами по коэффициентам разностного уравнения и начальным данным задачи Коши, которые заданы в “полосе”.

2. ПОСТАНОВКА ЗАДАЧИ И ИЗВЕСТНЫЕ РЕЗУЛЬТАТЫ

Обозначим $\mathbb{Z}$ – множество целых чисел, ${{\mathbb{Z}}^{2}} = \mathbb{Z} \times \mathbb{Z}$ – двумерная целочисленная решетка и $\mathbb{Z}_{ + }^{2}$ – подмножество этой решетки, состоящее из точек с целыми неотрицательными координатами. Пусть δ₁ – оператор сдвига по переменной $x$, т.е. ${{\delta }_{1}}f(x,y) = f(x + 1,y)$, а ${{\delta }_{2}}$ – оператор сдвига по переменной $y$, т.е. ${{\delta }_{2}}f(x,y) = f(x,y + 1)$. Зададим “полосу”, $\Pi = \{ (x,y) \in {{\mathbb{Z}}^{2}},0 \leqslant x \leqslant B,y \geqslant 0\} $ в положительном октанте целочисленной решетки, число $B + 1$ будем называть шириной “полосы” $\Pi $. Рассмотрим разностный полиномиальный оператор с постоянными коэффициентами вида

где $m$ и b определяют размер схемы, P_j(δ₁) = = $\sum\limits_{i = 0}^b {{{c}_{{ij}}}} \delta _{1}^{i}$, $j = 0,\;1,\; \ldots ,\;m$.

Многочлен $P(z,w) = \sum\limits_{j = 0}^m {\sum\limits_{i = 0}^b {{{c}_{{ij}}}} } {{z}^{i}}{{w}^{j}}$ называется характеристическим. Степень $m$ многочлена $P(z,w)$ по переменной $w$ будем называть порядком разностного оператора $P({{\delta }_{1}},{{\delta }_{2}})$ и предполагать, что $b < B$.

Зафиксируем $\beta = \left( {{{x}_{\beta }},m} \right)$ такое, что ${{c}_{\beta }} \ne 0$, и рассмотрим множество Π_β = {(x, y) ∈ $\mathbb{Z}_{ + }^{2}$: $0 \leqslant x - {{x}_{\beta }} \leqslant $ B – b, y > m – 1}. Обозначим ${{L}_{\beta }} = \Pi {{\backslash }}{{\Pi }_{\beta }}$ и сформулируем следующую задачу:

найти решение разностного уравнения

удовлетворяющее условиюгде $g(x,y)$ и $\varphi (x,y)$ – заданные функции целочисленных аргументов.

Задачу (2)–(3) назовем задачей Коши для полиномиального разностного оператора (1) и приведем легко проверяемое условие ее разрешимости.

Известно (см. [9]), что задача (2)–(3) однозначно разрешима, если выполняется условие

Поставим задачу вычислить значение функции $f(x,y)$ в точке A с координатами $\left( {{{x}_{A}},{{y}_{A}}} \right)$, т.е. $f\left( {{{x}_{A}},{{y}_{A}}} \right)$.

3. ОПИСАНИЕ ВХОДНЫХ ДАННЫХ

Решение задачи Коши (2)–(3) для двумерного разностного уравнения с постоянными коэффициентами в точке $z = ({{x}_{z}},{{y}_{z}})$ представляет собой значение функции $f(x,y)$ в этой точке.

Алгоритм вычисления значения функции $f(x,y)$ в заданной точке A c координатами $({{x}_{A}},{{y}_{A}})$ имеет рекурсивный характер и сводится к вычислению значений функции $f(x,y)$ на конечном подмножестве точек $\left( {x,y} \right) \in {{L}_{\beta }} = \Pi {{\backslash }}{{\Pi }_{\beta }}$.

Начальные данные (3) задаются матрицей F размерности $(B + 1) \times ({{y}_{A}} + 1)$, содержащей конечное подмножество значений начальных данных задачи Коши.

Коэффициенты двумерного разностного уравнения задаются матрицей $C$, имеющей прямоугольный вид (размерность: $(m + 1) \times (b + 1)$).

Проиллюстрируем процедуру задания функции начальных данных на примере.

Для разностного уравнения

и $\beta = (1,2)$ матрица коэффициентов C имеет вид а матрица начальных данных F размерности, например, 4 × 5 будет иметь вид

Здесь элементы, обозначенные $*$, вычисляются при выполнении алгоритма. При этом, вычислить элемент $\varphi (1,\;2)$ не представляется возможным без вычисления элемента $\varphi (2,\;2)$. Таким образом, для нахождения неизвестных элементов необходимо решить систему линейных разностных уравнений вида (5) с использованием начальных данных $\varphi (i,j)$, где $i = 0,\; \ldots ,\;4$, $j = 0,\;1$ и $\varphi (0,\;2)$, а именно:

Все, что нам известно, перенесем вправо и получим константы:

Матрица полученной системы уравнений будет иметь вид:

Поскольку на всех диагоналях матрицы, параллельных главной диагонали и на самой главной диагонали, элементы матрицы одинаковые, то такая матрица является теплицевой [10].

Вернемся к входным данным. Входные данные конечны и имеют вид:

1. точка $\beta = \left( {{{x}_{\beta }},m} \right)$;

2. прямоугольная марица коэффициентов $C = ({{c}_{\alpha }})$, $\alpha = \left( {{{\alpha }_{1}},{{\alpha }_{2}}} \right)$, ${{\alpha }_{1}} = 0,\; \ldots ,\;b$, ${{\alpha }_{2}} = 0,\; \ldots ,\;m$ размера $(m + 1) \times (b + 1)$ из коэффициентов ${{c}_{\alpha }}$ полиномиального разностного оператора (1);

3. точка $A$ с координатами $({{x}_{A}},{{y}_{A}})$, определяющая координаты искомого значения функции $f(x,y)$ и размерность матрицы начальных данных $F$;

4. матрица начальных данных $F = (\varphi (x,y))$, размера $(B + 1) \times ({{y}_{A}} + 1)$, где $(x,y) \in {{L}_{\beta }}$, для всех остальных значений (x, y) значения $\varphi (x,y) = 0$.

Для технической реализации алгоритма целесообразно отразить матрицы коэффициентов C и начальных данных $F$ зеркально отностельно горизонтальной оси. Поскольку координаты элементов матрицы коэффициентов разностного оператора и матрицы начальных данных в декартовой системе координат $(X,Y)$ не совпадают с их координатами в зеркально отраженной матрице (строка$ \times $столбец), то следует перейти из декартовой системы координат $\left( {D\left( {{{d}_{1}},{{d}_{2}}} \right)} \right)$ в “матричную” $\left( {M\left( {{{m}_{1}},{{m}_{2}}} \right)} \right)$ по правилу: $D\left( {{{d}_{1}},{{d}_{2}}} \right) \to M\left( {{{m}_{1}},{{m}_{2}}} \right)$, где ${{m}_{1}} = 1 + {{d}_{2}}$, ${{m}_{2}} = 1 + {{d}_{1}}$.

Например, элемент ${{c}_{{00}}}$ матрицы коэффициентов C, имеющий в декартовой системе координат координаты $\left( {{{d}_{1}},{{d}_{2}}} \right)$ = (0, 0), в “матричной” системе координат будет иметь координаты $({{m}_{1}},{{m}_{2}})$ = = (1, 1), а, например, элемент $\varphi $(1, 0) матрицы начальных данных $F$ в “матричной” системе координат будет иметь координаты (1, 2).

Далее необходимо будет проверить задачу Коши (2)–(3) на разрешимость, т.е. проверить, выполняется ли условие (4) для коэффициентов разностного оператора (1).

4. ПРИМЕР

Алгоритм был реализован в среде MatLab2014 32bit. Вычисления производились на машине Intel(R) Core(TM) i5-3330S CPU 2.70 GHz, 32bit, ОЗУ 4.00 Гб под управлением Windows 7 Корпоративная SP1. Время счета для приведенного примера составило менее 1 секунды.

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

Зафиксируем $\beta = ({{x}_{\beta }},m)$, ${{x}_{\beta }} = 2$, $m = 1$, $b = 3$, $B = 5$.

Для полиномиального разностного оператора

задача Коши будет иметь вид где Π = $\{ (x,y) \in {{\mathbb{Z}}^{2}},0 \leqslant x \leqslant 5,y \geqslant 0\} $, Π_{(2, 1)} = = $\{ (x,y) \in \mathbb{Z}_{ + }^{2}{\text{:}}\;2 \leqslant x \leqslant 4,y \geqslant 1\} $ и ${{L}_{{\left( {2,1} \right)}}} = \Pi {{\backslash }}{{\Pi }_{{\left( {2,1} \right)}}}$.

Зададим матрицу коэффициентов полиномиального разностного оператора (6)

Расположение элементов матрицы $C$ в декартовой системе координат представлено на рис. 1.

Поставим задачу: найти значения функции $f(x,y)$ в точке $A$ с координатами (4, 4), т.е. $f$(4, 4).

Зададим матрицу начальных данных:

Распределение температуры в образце в начальный момент времени, т.е. последняя строка матрицы $F$, имеет форму, представленную на рис. 2.

Расположение элементов матрицы $F$ в декартовой системе координат представлено на рис. 3.

1. Переход из декартовой системы координат в “матричную” осуществляется путем отражения матриц начальных данных и коэффициентов зеркально относительно горизонтальной оси, поэтому точка $\beta $ с координатами (2, 1) переходит в точку с координатами (2, 3) в “матричной” системе координат.

Пусть $p$ – количество столбцов матрицы коэффициентов, а $i$ – количество столбцов матрицы начальных данных. Тогда для рассматриваемого примера $p = 4$ и $i = 6$.

“Матричная” система координат для матрицы коэффициентов разностного оператора (6) имеет вид, представленный в табл. 1:

“Матричная” система координат для матрицы начальных данных $F$ имеет вид, представленный в табл. 2:

В табл. 2 знак $ \times $ – начальные данные. Точка $A$ имеет координаты [5, 5 ] в “матричной” системе координат.

2. Проверка на разрешимость задачи Коши (7), т.е. выполняется ли условие (4) для коэффициентов полиномиального разностного оператора (6). Так как

то условие (4) выполняется и задача Коши разрешима.

3. Зеркальное отражение относительно горизонтальной оси матрицы коэффициентов C и матрицы начальных данных $F$:

4. Последняя строка в матрице ${{C}_{{work}}}$:

5. Получение общей Теплицевой матрицы [10, 11], размер которой зависит от искомого значения.

Для построения теплицевой матрицы необходимы два вектора $colum{{n}_{c}}$ и $ro{{w}_{c}}$. Пусть матрица $neiz{{v}_{1}}$ – несимметричная Теплицева матрица размера $\left( {i - p + 1} \right) \times \left( {i - p + 1} \right)$ (в рассматриваемом примере 3 × 3), построенная на основе правила, что элементы вектора $colum{{n}_{c}}$ – ее первый столбец, а элементы вектора $ro{{w}_{c}}$ – ее первая строка, причем первые элементы векторов $colum{{n}_{c}}$ и $ro{{w}_{c}}$ совпадают. В случае, если количество элементов в векторах $colum{{n}_{c}}$ и $ro{{w}_{c}}$ меньше $\left( {i - p + 1} \right)$, то вектор меньшего размера дополняется нулями. Элементы векторов $colum{{n}_{c}}$ и $ro{{w}_{c}}$ берутся из вектора $a$. В вектор $colum{{n}_{c}}$ попадают элементы, начиная с ${{c}_{\beta }} = {{c}_{{21}}}$ в порядке уменьшения номера элемента, т.е. $colum{{n}_{c}} = {{\left( {\begin{array}{*{20}{c}} {{{c}_{{21}}}}&{{{c}_{{11}}}}&{{{c}_{{01}}}} \end{array}} \right)}^{{\text{т}}}} = {{\left( {\begin{array}{*{20}{c}} {10}&3&1 \end{array}} \right)}^{{\text{т}}}}$, а в вектор $ro{{w}_{c}}$ попадают элементы, начиная ${{c}_{\beta }} = {{c}_{{21}}}$ в порядке увеличения номера элемента, т.е. row_c = = $\left( {\begin{array}{*{20}{c}} {{{c}_{{21}}}}&{{{c}_{{31}}}} \end{array}} \right)$ = (10 1). Поскольку длина вектора $ro{{w}_{c}}$ меньше $\left( {i - p + 1} \right)$, то дополним его нулем: $ro{{w}_{c}}$ = = (10 1 0). Таким образом, матрица $neiz{{v}_{1}}$ имеет вид

6. Решение системы линейных уравнений с целью нахождения значений $f\left( {x,y} \right)$, которые могут потребоваться для вычисления искомого значения в точке $A\left( {{{x}_{A}},{{y}_{A}}} \right)$.

Количество неизвестных для каждого значения $m$ (т.е. номера строки матрицы $F$ в декартовой системе координат), показывающего номер “слоя”, на каждом шаге алгоритма одинаково. Поэтому размер матрицы Теплица остается постоянным. На каждом шаге алгоритма происходит запись найденных значений функции $f\left( {x,y} \right)$ в матрицу начальных данных $F$.

На рис. 4 представлено итоговое распределение температуры в образце в течение четырех временных шагов.

7. Переход к декартовой системе координат $f(4,4)$ = F[5, 5].

Output: значение функции $f\left( {x,y} \right)$ в искомой точке $A\left( {4,\;4} \right)$ равно 0.19.

Входными данными для работы алгоритма в этом случае будут:

1. $\beta = (2,\;1)$;

2. матрица коэффициентов

3. $A = (4,\;4);$

4. матрица начальных данных