2022

These authors contributed equally to this work.

[2,3]Alexander Lyapin These authors contributed equally to this work.

These authors contributed equally to this work.

[1]Kransnoyarsk State Medical University, Krasnoyarsk, 660022, Russia

2]Siberian Federal University, Krasnoyarsk, 660041, Russia

3]Fairmont State University, Fairmont, 26554, WV, United States

# Solving the Cauchy problem for a three-dimensional difference equation in a
parallelepiped

Marina Apanovich<EMAIL_ADDRESS><EMAIL_ADDRESS>Konstantin
Shadrin<EMAIL_ADDRESS>*  [  [

###### Abstract

The aim of this article is further development of the theory of linear
difference equations with constant coefficients. We present a new algorithm
for calculating the solution to the Cauchy problem for a three-dimensional
difference equation with constant coefficients in a parallelepiped at the
point using the coefficients of the difference equation and Cauchy data. The
implemented algorithm is the next significant achievement in a series of
articles justifying the Apanovich and Leinartas’ theorems about the
solvability and well-posedness of the Cauchy problem. We also use methods of
computer algebra since the three-dimensional case usually demands extended
calculations.

###### keywords:

difference operator, difference equation, Cauchy problem

###### pacs:

[

MSC Classification]05A15, 37H10, 39A05, 39A70

##  1 Introduction

Difference equations arise in various fields of mathematics and have numerous
applications in science and technology. For example, in combination with the
method of generating functions they appear to be a powerful method of study
enumerative problems in the combinatorial analysis (see, for example, Stanley
(1997, 1999); Bousquet-Melou (2000)). Another source of difference equations
is discretization of differential equations. For example, the discretization
of the Cauchy-Riemann equation resulted in the theory of discrete analytical
functions (see Duffin (1956, 1968)), which is used in the theory of Riemann
surfaces and combinatorial analysis (see, for example, Danilov (2008, 2009)).
Methods of discretization of a differential problem are an important part of
the theory of difference schemes and also lead to difference equations (see,
for example, Samarskii (1971); Godunovi (1964)). Also difference equations are
used in discrete time dynamical models (see Dudgeon (1983); Isermann (1981)).

Additional conditions ("initial", "boundary", "Cauchy data") are given for the
space of solutions of a multidimensional difference equation, which allows to
select a unique solution from an infinite set of solutions, and the
corresponding problem is called the Cauchy problem for a difference equation.
A one-dimensional case is simple since the inital data of the Cauchy problem
is finite (thus a gererating function for the solution is rational). The
multidimensional case is different since the initial data is given on an
infinite set and might have a difficult structure. A significant number of
works is devoted to the development of algorithms for solving multidimensional
difference equations with constant and polynomial coefficients (see, for
example, Abramov (2011, 2015, 2011)), and in Abramov (2011) difference
equations with coefficients in a form of rational functions are considered.
The connection between the generating function for solution to the Cauchy
problem for a two-dimensional difference equation with constant coefficients
and the generating function of the initial data is studied in Leinartas
(2009); Kytmanov (2017). In Apanovich (2021) an algorithm for calculating the
solution of the Cauchy problem for a two-dimensional difference equation with
constant coefficients at a point from the coefficients of the difference
equation and the initial data of the Cauchy problem was developed and
implemented.

However, the multidimensional cases have not been investigated properly. The
three-dimensional case is of great importance in problems of thermodynamics
(thermal conductivity, see Samarskii (1971), anisotropic diffusion, see
MurrayJD (2003)), mathematical biology (distribution of morphogen in insect
wings MurrayJD (2003), distribution of population density species bratus
(2010)). This paper considers the first step – solving the Cauchy problem for
a three-dimensional difference equation with constant coefficients in a
parallelepiped at a point from coefficients of the difference equation and the
initial data of the Cauchy problem. We give a notation of the Cauchy problem
in a "parallelepiped" and describe a function of the initial data on this set
in a way that allows us to develop and implement a computer algebra algorithm.

##  2 The Cauchy problem for a polynomial difference operator in a
parallelipiped

Let $f$ be a function of integer variables
$f:\mathbb{Z}_{+}^{n}\rightarrow\mathbb{C}$ and $\delta_{j}$ be a shift
operator by $j$-th variable, i.e.
$\delta_{j}f\left(t_{1},...,t_{n}\right)=f\left(t_{1},...,t_{j-1},t_{j}+1,t_{j+1},...,t_{n}\right)$,
$j=1,2,...,n$. If $\alpha=\left(\alpha_{1},...,\alpha_{n}\right)$ is a multi-
index, then $\|\alpha\|=\alpha_{1}+...+\alpha_{n}$,
$\delta^{\alpha}=\delta_{1}^{\alpha_{1}}...\delta_{n}^{\alpha_{n}}$. For two
multi-indexes $\alpha=\left(\alpha_{1},...,\alpha_{n}\right)$ and
$\beta=\left(\beta_{1},...,\beta_{n}\right)$ the inequality
$\alpha\leqslant\beta$ means that $\alpha_{j}\leqslant\beta_{j}$ for all
$j=1,2,...,n$.

We consider a difference polynomial operator

$P\left(\delta\right)=\sum_{0\leqslant\alpha\leqslant
m}c_{\alpha}\delta^{\alpha},$ (1)

where $c_{\alpha}$ are constants coefficients of the operator
$P\left(\delta\right)$. We denote
$\alpha=\left(\alpha_{1},...,\alpha_{n}\right)$ and
${}^{\prime}\alpha=\left(\alpha_{1},...,\alpha_{n-1}\right)$, and then write
the operator $P\left(\delta\right)$ as follows

$P\left(\delta\right)=\sum_{\alpha_{n}=0}^{m_{n}}\left(\sum_{{}^{\prime}0\leqslant^{\prime}\alpha\leqslant^{\prime}m}c_{{}^{\prime}\alpha,\alpha_{n}}\delta^{{}^{\prime}\alpha}\right)\delta_{n}^{\alpha_{n}}.$
(2)

We call a polynomial

$P\left({}^{\prime}z,z_{n}\right)=\sum_{\alpha_{n}=0}^{m_{n}}\left(\sum_{{}^{\prime}0\leqslant^{\prime}\alpha\leqslant^{\prime}m}c_{{}^{\prime}\alpha,\alpha_{n}}^{\prime}z^{{}^{\prime}\alpha}\right)z_{n}^{\alpha_{n}}=\sum_{\alpha_{n}=0}^{m_{n}}P_{j}\left({}^{\prime}z\right)z_{n}^{\alpha_{n}}$

as a characteristic polynomial for the difference operator
$P\left(\delta\right)$, and its degree by a variable $z_{n}$ is an order of
the difference operator $P\left(\delta\right)$.

Let $x=\left(x_{1},...,x_{n-1},0\right)$ be a point in the integer lattice
$\mathbb{Z}_{+}^{n}$ and

$\Pi_{{}^{\prime}x}=\left\\{t\in\mathbb{R}_{+}^{n}:^{\prime}0\leqslant^{\prime}t\leqslant^{\prime}x\right\\}$

be parallelepiped of a dimension $n-1$ in the hyperplane $t_{n}=0$.

We fix a point
$\beta=\left(\beta_{1},...,\beta_{n-1},m_{n}\right)\in\mathbb{Z}_{+}^{n}$,
${}^{\prime}0\leqslant^{\prime}\beta\leqslant^{\prime}m$ such that
$c_{\beta}\neq 0$, and consider a set
$\Pi_{{}^{\prime}\beta,^{\prime}m}=\left\\{{}^{\prime}t\in\Pi_{{}^{\prime}x}:^{\prime}\beta\leqslant^{\prime}t\leqslant^{\prime}x-^{\prime}m+^{\prime}\beta\right\\}$.
Let
$L=\left(\Pi_{{}^{\prime}x}\setminus\Pi_{{}^{\prime}\beta,^{\prime}m}\right)\times\left[0,x_{n}\right]$
be a set where the initial data is given. We consider a problem: find a
solution to the difference equation

$P\left(\delta\right)f\left(x\right)=g\left(x\right),x\in\Pi=\Pi_{{}^{\prime}x}\times\left[0,x_{n}\right],$
(3)

such that the condition

$f\left(x\right)=\varphi\left(x\right),x\in L,$ (4)

is valid, where $g\left(x\right)$ and $\varphi\left(x\right)$ are given
functions of integer variables.

The problem (3) – (4) is called the Cauchy problem for a polynomial difference
operator (1).

In Fedoruk (1987) the stability of a homogeneous two-layer linear difference
scheme with constant coefficients is investigated for $m_{n}=1$. In Rogozina
(2014) the solvability of the problem 3 – 4 was investigated for $n=2$. In the
theory of difference schemes, such problems are multilayer implicit difference
scheme. In Apanovich (2017) the well-posedness of problem (3) – (4) for $n=2$
is investigated and an easily verifiable sufficient condition for correctness
is proved. In Apanovich (2018) for $n=3$, an easily verified sufficient
condition for the solvability of the Cauchy problem (3) – (4) is proved.

We consider Cauchy problem (3) – (4) for $n=3$ and denote $\mathbb{Z}^{3}$ –
the integer lattice and $\mathbb{Z}^{3}_{+}$ is a subset of this lattice
consisting of points with non-negative integer coordinates. Let
$\delta_{1},\delta_{2},\delta_{3}$ be shift operators by variable $x,y,z$
accordingly, which means that $\delta_{1}f(x,y,z)=f(x+1,y,z)$,
$\delta_{2}f(x,y,z)=f(x,y+1,z)$, $\delta_{3}f(x,y,z)=f(x,y,z+1)$. We denote a
parallelepiped $\Pi=\left\\{(x,y,z)\in\mathbb{Z}_{+}^{3},\,0\leqslant
x\leqslant B_{x},\,0\leqslant y\leqslant B_{y},\,z\geqslant 0\right\\}$ in a
positive octant of the integer lattice $\mathbb{Z}^{3}_{+}$, and $B_{x}+1$ is
a width of the parallelipiped $\Pi$ and $B_{y}+1$ is a length of the
parallelipiped $\Pi$. The difference polynomial operatop (1) is

$P(\delta_{1},\delta_{2},\delta_{3})=\sum_{j=0}^{m}\sum_{i_{y}=0}^{b_{y}}\sum_{i_{x}=0}^{b_{x}}c_{i_{x}i_{y}j}\delta_{1}^{i_{x}}\delta_{2}^{i_{y}}\delta_{3}^{j}=\sum_{j=0}^{m}P_{j}(\delta_{1},\delta_{2})\delta_{3}^{j},$
(5)

where $b_{x}$, $b_{y}$ and $m$ define a size of a difference scheme and

$P_{j}(\delta_{1},\delta_{2})=\sum_{i_{x}=0}^{b_{x}}\sum_{i_{y}=0}^{b_{y}}c_{i_{x}i_{y}}\delta_{1}^{i_{x}}\delta_{2}^{i_{y}},j=0,1,...,m.$
(6)

The characteristic polynomial is

$P(s,w,v)=\sum_{j=0}^{m}\sum_{i_{x}=0}^{b_{x}}\sum_{i_{y}=0}^{b_{y}}c_{i_{x}i_{y}j}s^{i_{x}}w^{i_{y}}v^{j},$
(7)

where $m$ is an order of the difference operator
$P(\delta_{1},\delta_{2},\delta_{3})$, $b_{x}<B_{x}$, $b_{y}<B_{y}$.

We fix $\beta=\left(x_{\beta},y_{\beta},m\right)$ such that
$c_{x_{\beta}y_{\beta}m}\neq 0$ and consider a set
$\Pi_{\beta}=\\{(x,y,z)\in\mathbb{Z}_{+}^{3}:\,0\leqslant x-x_{\beta}\leqslant
B_{x}-b_{x},\,0\leqslant y-y_{\beta}\leqslant B_{y}-b_{y},\,z>m-1\\}$, then
$L=\Pi\setminus\Pi_{\beta}$ and the Cauchy problem is

to find a solution to the difference equation

$P(\delta_{1},\delta_{2},\delta_{3})f(x,y,z)=g(x,y,z),(x,y,z)\in\Pi$ (8)

under a condition that

$f(x,y,z)=\varphi(x,y,z),(x,y,z)\in L,$ (9)

where $g(x,y,z)$ and $\varphi(x,y,z)$ are given functions of integer
arguments.

Problem (8) – (9) is the Cauchy problem for polynomial difference operator
(5).

In Apanovich (2018) it was proven that problem (8) – (9) is uniquely solvable
if the condition

$\|c_{x_{\beta}y_{\beta}m}\|>\sum_{\left(\alpha_{1},\alpha_{2}\right)\neq\left(x_{\beta},y_{\beta}\right),\alpha_{3}=m}\|c_{\alpha_{1}\alpha_{2}\alpha_{3}}\|$
(10)

takes place.

Our problem is to compute a value of function $f\left(x,y,z\right)$ at point
$A$ with coordinates $\left(x_{1},y_{1},z_{1}\right)$.

##  3 Description of the Input Data

A solution to Cauchy problem (8) – (9) for three-dimensional difference
equation with constant coefficients at a point $A$ with coordinates
$(x_{1},y_{1},z_{1})$ is a value of the function $f(x,y,z)$ at a point $A$.
The algorithm of computing a value of function $f(x,y,z)$ at a point with
given coordinates $(x_{1},y_{1},z_{1})$ is recursive and reduces to computing
values of a function $f(x,y,z)$ on the finite set of points
$\left(x,y,z\right)$ in the set $L=\Pi\setminus\Pi_{\beta}$.

Initial data (6) is giver by a matrix $F$ of a dimension three, containing a
finite set of values of the initial data of the Cauchy problem. Coefficients
of three-dimensional difference equation are given by a matrix $C$ of
dimension three. For the technical implementation of the algorithm, matrices
of coefficients $C$ and initial data $F$ are specified in layers, starting
from the lowest one. Let us illustrate the procedure for specifying the
matrices $F$ and $C$.

For a difference equation

$\begin{array}[]{c}c_{000}f\left(x,y,z\right)+c_{100}f\left(x+1,y,z\right)+c_{200}f\left(x+2,y,z\right)+c_{010}f\left(x,y+1,z\right)+\\\
+c_{110}f\left(x+1,y+1,z\right)+c_{210}f\left(x+2,y+1,z\right)++c_{001}f\left(x,y,z+1\right)+\\\
+c_{101}f\left(x+1,y,z+1\right)+c_{201}f\left(x+2,y,z+1\right)+c_{011}f\left(x,y+1,z+1\right)+\\\
+c_{111}f\left(x+1,y+1,z+1\right)+c_{211}f\left(x+2,y+1,z+1\right)=0\end{array}$
(11)

the first layer of the matrix $C$ of coefficients is a matrix $C_{0}$:

$C_{0}=\begin{pmatrix}c_{000}&c_{100}&c_{200}\\\
c_{010}&c_{110}&c_{210}\end{pmatrix},$ (12)

the second layer of $C$ is a matrix $C_{1}$:

$C_{1}=\begin{pmatrix}c_{001}&c_{101}&c_{201}\\\
c_{011}&c_{111}&c_{211}\end{pmatrix}.$ (13)

A matrix $C$ of coefficients is written as follows:

$C=\left(\begin{array}[]{c|c}\begin{array}[]{c}C_{0}\end{array}&\begin{array}[]{c}C_{1}\end{array}\end{array}\right)=\left(\begin{array}[]{c|c}\begin{matrix}c_{001}&c_{101}&c_{201}\\\
c_{011}&c_{111}&c_{211}\end{matrix}&\begin{matrix}c_{001}&c_{101}&c_{201}\\\
c_{011}&c_{111}&c_{211}\end{matrix}\end{array}\right)$ (14)

For difference equation (11) we have $b_{x}=2$, $b_{y}=1$, $m=1$. For
$\beta=\left(1,1,1\right)$, $B_{x}=4$, $B_{y}=2$, the three-dimensional matrix
$F$ of the initial data is
$F=\left(\begin{array}[]{c|c}\begin{array}[]{c}F_{0}\end{array}&\begin{array}[]{c}F_{1}\end{array}\end{array}\right)$,
where

$F_{0}=\left(\begin{array}[]{ccccc}\varphi\left(0,0,0\right)&\varphi\left(1,0,0\right)&\varphi\left(2,0,0\right)&\varphi\left(3,0,0\right)&\varphi\left(4,0,0\right)\\\
\varphi\left(0,1,0\right)&\varphi\left(1,1,0\right)&\varphi\left(2,1,0\right)&\varphi\left(3,1,0\right)&\varphi\left(4,1,0\right)\\\
\varphi\left(0,2,0\right)&\varphi\left(1,2,0\right)&\varphi\left(2,2,0\right)&\varphi\left(3,2,0\right)&\varphi\left(4,2,0\right)\end{array}\right),$
(15)
$F_{1}=\left(\begin{array}[]{ccccc}\varphi\left(0,0,1\right)&\varphi\left(1,0,1\right)&\varphi\left(2,0,1\right)&\varphi\left(3,0,1\right)&\varphi\left(4,0,1\right)\\\
\varphi\left(0,1,1\right)&*&*&*&\varphi\left(4,1,1\right)\\\
\varphi\left(0,2,1\right)&*&*&*&\varphi\left(4,2,1\right)\end{array}\right).$
(16)

The entries of $F$, denoted by $*$, are calculated when the algorithm is
executed. However, it is not possible to calculate the element
$\varphi(1,2,1)$ without calculating the elements $\varphi\left(2,2,1\right)$,
$\varphi(1,1,1)$, $\varphi(2,1,1)$, $\varphi(3,1,1)$, $\varphi(3,2,1)$. Thus,
to find the unknown elements, it is necessary to solve a system of linear
difference equations of the form (8) using the initial data $\varphi(i,j,k)$,
where $i=0,...,4$, $j=0,1,2$, $k=0$, $\varphi(i,j,k)$, $i=0,...,4$, $j=0$,
$k=1$ and $\varphi\left(0,1,1\right)$, $\varphi\left(0,2,1\right)$,
$\varphi\left(4,1,1\right)$, $\varphi\left(4,2,1\right)$.

Finally, the input data is finite:

  1. 1.

a three-dimensional $(b_{x}+1)\times(b_{y}+1)\times(m+1)$-matrix
$C=(c_{\alpha_{1}\alpha_{2}\alpha_{3}}),\alpha_{1}=0,\ldots,b_{x},\alpha_{2}=0,\ldots,b_{y}$,
$\alpha_{3}=0,...,m$ with coefficients $c_{\alpha_{1}\alpha_{2}\alpha_{3}}$ of
three-dimensional difference equation;

  2. 2.

a point $\beta_{0}=(x_{\beta},y_{\beta},m)$;

  3. 3.

a point $A$ with coordinates $(x_{1},y_{1},z_{1})$, which defines the
coordinates of the desired value of the function $f(x,y,z)$ and a number of
layers in the three-dimensional matrix $F$ of initial data;

  4. 4.

a three-dimensional $(B_{x}+1)\times(B_{y}+1)\times(z_{1}+1)$-matrix of the
initial data $F=(\varphi(x,y,z))$ for $(x,y,z)\in L$ and $\varphi(x,y,z)=0$
for $(x,y,z)\notin L$.

Since the coordinates of the elements of the three-dimensional matrix $C$ of
coefficients of the difference operator and the three-dimensional matrix $F$
of initial data in the Cartesian coordinate system $(X,Y,Z)$ do not coincide
with their coordinates in the matrix (row $\times$ column $\times$ layer),
then they have to be transformed from Cartesian coordinates
$\left(D\left(d_{1},d_{2},d_{3}\right)\right)$ into the "matrix>> coordinates
$\left(M\left(m_{1},m_{2},m_{3}\right)\right)$ as follows:
$D\left(d_{1},d_{2},d_{3}\right)\rightarrow M\left(m_{1},m_{2},m_{3}\right)$,
where ${m_{1}=d_{2}}+1$, ${m_{2}=d_{1}}+1$, $m_{3}=d_{3}+1$.

Then it will be necessary to check the Cauchy problem (8) – (9) for
solvability, which means, we have to check whether coefficients of the
difference operator (5) satisfy to condition (10).

##  4 Example

We consider the polynomial difference operator

$\begin{array}[]{c}P(\delta_{1},\delta_{2},\delta_{3})=c_{000}+c_{100}\delta_{1}+c_{200}\delta_{1}^{2}+c_{010}\delta_{2}+c_{110}\delta_{1}\delta_{2}+c_{210}\delta_{1}^{2}\delta_{2}+\\\
+c_{020}\delta_{2}^{2}+c_{120}\delta_{1}\delta_{2}^{2}+c_{220}\delta_{1}^{2}\delta_{2}^{2}+c_{001}\delta_{3}+c_{101}\delta_{1}\delta_{3}+c_{201}\delta_{1}^{2}\delta_{3}+\\\
+c_{011}\delta_{2}\delta_{3}+c_{111}\delta_{1}\delta_{2}\delta_{3}+c_{211}\delta_{1}^{2}\delta_{2}\delta_{3}+c_{021}\delta_{2}^{2}\delta_{3}+c_{121}\delta_{1}\delta_{2}^{2}\delta_{3}+c_{221}\delta_{1}^{2}\delta_{2}^{2}\delta_{3}=\\\
=1+2\delta_{1}+3\delta_{1}^{2}+4\delta_{2}+5\delta_{1}\delta_{2}+6\delta_{1}^{2}\delta_{2}+1\delta_{2}^{2}+2\delta_{1}\delta_{2}^{2}+\\\
+3\delta_{1}^{2}\delta_{2}^{2}+3\delta_{3}+4\delta_{1}\delta_{3}+5\delta_{1}^{2}\delta_{3}+6\delta_{2}\delta_{3}+\\\
+80\delta_{1}\delta_{2}\delta_{3}+8\delta_{1}^{2}\delta_{2}\delta_{3}+3\delta_{2}^{2}\delta_{3}+4\delta_{1}\delta_{2}^{2}\delta_{3}+5\delta_{1}^{2}\delta_{2}^{2}\delta_{3},\end{array}$
(17)

where $b_{x}=2$, $b_{y}=2$, $m=1$. We fix
$\beta_{0}=(x_{\beta},y_{\beta},m)=(1,1,1)$, $B_{x}=4$, $B_{y}=3$. Then the
set of initial data is

$L=\Pi\setminus\Pi_{\left(1,1\right)},$

where $\Pi=\left\\{(x,y,z)\in\mathbb{Z}^{3},\,0\leqslant x\leqslant
4,\,0\leqslant y\leqslant 3,z\geqslant 0\right\\}$ and
$\Pi_{\left(1,1\right)}=\\{(x,y,z)\in\mathbb{Z}_{+}^{3}:1\leqslant x\leqslant
3,1\leqslant y\leqslant 2,z\geqslant 1.\\}$

The three-dimensional matrix $C$ of coefficients of polynomial difference
operator (17) is given by layers, starting from the lowest one:

$C=\left(\begin{array}[]{c|c}\begin{array}[]{c}\begin{matrix}c_{000}&c_{100}&c_{200}\\\
c_{010}&c_{110}&c_{210}\\\
c_{020}&c_{120}&c_{220}\end{matrix}\end{array}&\begin{array}[]{c}\begin{matrix}c_{001}&c_{101}&c_{201}\\\
c_{011}&c_{111}&c_{211}\\\
c_{021}&c_{121}&c_{221}\end{matrix}\end{array}\end{array}\right)=\left(\begin{array}[]{c|c}\begin{array}[]{c}\begin{matrix}1&2&3\\\
4&5&6\\\ 1&2&3\end{matrix}\end{array}&\begin{array}[]{c}\begin{matrix}3&4&5\\\
6&80&8\\\ 3&4&5\end{matrix}\end{array}\end{array}\right).$ (18)

The arrangement of the elements of the matrix of coefficients $C$ in the
Cartesian coordinate system is shown in Figure 1.

Figure 1: The arrangement of the elements of the matrix of coefficients $C$

The problem is to find a value of function $f\left(x,y,z\right)$ at a point
$A$ with coordinates $\left(2,1,4\right)$.

The matrix $F$ of the initial data is given by layers, starting from the
lowest one:

$F=\left(\begin{array}[]{c|c|c|c|c}F_{0}&F_{1}&F_{2}&F_{3}&F_{4}\end{array}\right)=\\\
\left(\begin{array}[]{c|c|c|c|c}\begin{array}[]{ccccc}1&2&3&4&5\\\
4&5&6&1&3\\\ 1&2&3&8&4\\\
1&2&3&4&6\end{array}&\begin{array}[]{ccccc}3&4&5&7&6\\\ 6&0&0&0&8\\\
3&0&0&0&4\\\ 6&5&3&2&3\end{array}&\begin{array}[]{ccccc}3&4&5&7&6\\\
5&0&0&0&2\\\ 3&0&0&0&4\\\
6&8&1&4&3\end{array}&\begin{array}[]{ccccc}3&4&5&7&6\\\ 6&0&0&0&1\\\
3&0&0&0&4\\\ 6&6&7&1&4\end{array}&\begin{array}[]{ccccc}3&4&5&7&6\\\
3&0&0&0&8\\\ 3&0&0&0&4\\\ 6&2&6&7&3\end{array}\end{array}\right).$ (19)

The arrangement of the elements of the matrix $F$ if the initial data in the
Cartesian coordinates is shown in Figure 2.

Figure 2: The arrangement of the elements of the matrix $F$

Now we have to check the following items:

  1. 1.

A transition from the Cartesian coordinates into "matrix>> coordinates:
$\beta_{0}=[1,1,1]\rightarrow\beta_{1}=[2,2,2]$, $A(2,1,4)\rightarrow
A_{1}\left(2,3,5\right)$.

  2. 2.

A solvability of the problem. Since the inequality

$80=c_{111}>\sum_{\left(\alpha_{1},\alpha_{2}\right)\neq\left(1,1\right),\alpha_{3}=1}\|c_{\alpha_{1}\alpha_{2}\alpha_{3}}\|=38$

takes place, the problem is solvable.

  3. 3.

To find unknown (zero) values in the matrix $F$ on the second layer it is
necessary to solve the system of difference equations

$\begin{cases}P(\delta_{1},\delta_{2},\delta_{3})f(x,y,z)=0,(x,y,z)\in\Pi\,forz=0,1,\\\
f(x,y,z)=\varphi(x,y,z),(x,y,z)\in L\,for\,z=0,1.\end{cases}$ (20)

The matrix $T$ of the system of difference equations (20) will be a block
Töplitz matrix (for details see Iohvidov (1982); Apanovich (2020)):

$T=\left(\begin{array}[]{cc}T_{0}&T_{1}\\\ T_{-1}&T_{0}\end{array}\right),$
(21)

where

$T_{-1}=\left(\begin{array}[]{ccc}4&5&0\\\ 3&4&5\\\
0&3&4\end{array}\right),T_{0}=\left(\begin{array}[]{ccc}80&8&0\\\ 6&80&8\\\
0&6&80\end{array}\right),T_{1}=\left(\begin{array}[]{ccc}4&5&0\\\ 3&4&5\\\
0&3&4\end{array}\right).$ (22)

Solving system (20) yields

$F_{1}=\left(\begin{array}[]{ccccc}3.0000&4.0000&5.0000&7.0000&6.0000\\\
6.0000&-2.1249&-1.8591&-2.6440&8.0000\\\
3.0000&-1.4028&-1.3917&-2.4811&4.0000\\\
6.0000&5.0000&3.0000&2.0000&3.0000\end{array}\right).$ (23)

  4. 4.

To find the unknown (zero) values in the matrix $F$ on the third layer, it is
necessary to solve the system of difference equations:

$\begin{cases}P(\delta_{1},\delta_{2},\delta_{3})f(x,y,z)=0,(x,y,z)\in\Pi\,forz=1,2,\\\
f(x,y,z)=\varphi(x,y,z),(x,y,z)\in L\,for\,z=1,2.\end{cases}$ (24)

The matrix of the system of equations (24) will be equal to the matrix of the
system of equations (9). Solving the system of equations (24) yields:

$F_{2}=\left(\begin{array}[]{ccccc}3.0000&4.0000&5.0000&7.0000&6.0000\\\
5.0000&-1.1915&-0.4974&-1.9906&2.0000\\\
3.0000&-1.0084&0.0246&-1.2240&4.0000\\\
6.0000&8.0000&1.0000&4.0000&3.0000\end{array}\right).$ (25)

  5. 5.

Continuing the process yields:

$F_{3}=\left(\begin{array}[]{ccccc}3.0000&4.0000&5.0000&7.0000&6.0000\\\
6.0000&-1.4035&-0.7450&-1.5169&1.0000\\\
3.0000&-1.4785&-0.3556&-1.2873&4.0000\\\
6.0000&6.0000&7.0000&1.0000&4.0000\end{array}\right),$ (26)

$F_{4}=\left(\begin{array}[]{ccccc}3.0000&4.0000&5.0000&7.0000&6.0000\\\
3.0000&-1.1586&-0.7080&-2.1118&8.0000\\\
3.0000&-1.1695&-0.4721&-1.8794&4.0000\\\
6.0000&2.0000&6.0000&7.0000&3.0000\end{array}\right).$ (27)

Thus, the value is $f\left(1,2,4\right)=-0.7080$.

##  5 Description of the algorithm

The algorithm was implemented in the MatLab2014 32bit environment. The
calculations were performed on an Intel (R) Core (TM) i5-3330S CPU 2.70 GHz,
32bit, 4.00 GB RAM, running Windows 7 Enterprise SP1. The counting time for
the given example in Section 4 was less then a second.

1:a point $\beta_{0}$, a matrix $C$ of coefficients, a point $A$, a matrix $F$
of the initial data

2:the value of the function $f\left(x,y,z\right)$ at a point $A$

3:procedure DATA($C,\beta_{0},A,F$)

4: $\beta_{1}:=$ are coordinates of $\beta_{0}$ in "matrix" coordinates

5: if
$\|C\left(\beta_{1}\left(1\right),\beta_{1}\left(2\right),\beta_{1}\left(3\right)\right)\|<\sum\left(\sum\|C\left(:,:,\beta_{1}\left(3\right)\right)\|\right)-\|C\left(\beta_{1}\left(1\right),\beta_{1}\left(2\right),\beta_{1}\left(3\right)\right)\|$
then return The error of entering the matrix $C$

6: end if

7: $A_{1}:=$ are coordinates of $A$ in "matrix" coordinates

8: $Q:=C\left(:,:,end\right)$

9: $p=size\left(Q,2\right)$

10: $NAD=\beta_{1}\left(1\right)-1$

11: $POD=size\left(Q,1\right)-\beta_{1}\left(1\right)$

12:
$\begin{array}[]{c}sizeTblock=size\left(F,1\right)*size\left(F,2\right)-size\left(F,1\right)*\left(size\left(Q,2\right)-1\right)-\\\
-\left(size\left(F,2\right)-\left(size\left(Q,2\right)-1\right)\right)*\left(size\left(Q,1\right)-1\right)\end{array}$

13: $T:=zeros\left(sizeTblock,sizeTblock\right)$

14: $a=Q\left(\beta_{1}\left(1\right),:\right)$

15: $e=\beta_{1}\left(2\right)$

16: $column_{c}=zeros\left(p,1\right)$

17: $col=1$

18: while $e\geqslant 1$ do

19: $column_{c}\left(col\right):=a\left(e\right)$

20: $col:=col+1$

21: $e:=e-1$

22: end while

23: $e:=\beta_{1}(2)$

24: $row_{c}:=zeros\left(1,p\right)$

25: $r:=1$

26: for $m$ from $e$ to $length\left(a\right)$ do

27: $row_{c}\left(r\right):=a\left(e\right)$

28: $r:=r+1$

29: $e:=e+1$

30: end for

31: $T_{0}:=toeplitz\left(column_{c},row_{c}\right)$

32: for $i$ from $0$ to $sizeTblock/(size(Q,2))-1$ do

33: $T(((i*p)+1):p*(i+1),((i*p)+1):p*(i+1)):=T_{0}$

34: end for

35: for $w$ from $1$ to $NAD$ do

36: $a:=Q\left(\beta_{1}\left(1\right)-w,:\right)$

37: $e:=\beta_{1}\left(2\right)$

38: $column_{c}:=zeros\left(p,1\right)$

39: $col:=1$

40: while $e\geqslant 1$ do

41: $column_{c}\left(col\right):=a\left(e\right)$

42: $col:=col+1$

43: $e:=e-1$

44: end while

45: $e:=\beta_{1}\left(2\right)$

46: $row_{c}=zeros\left(1,p\right)$

47: $r:=1$

48: for $m$ from $e$ to $length\left(a\right)$ do

49: $row_{c}\left(r\right):=a\left(e\right)$

50: $r:=r+1$

51: $e:=e+1$

52: end for

53: $T_{NAD}:=toeplitz\left(colomn_{c},row_{c}\right)$

54: for $i$ from $0$ to $sizeTblock/\left(size\left(Q,2\right)\right)-1-w$ do

55: $T((((i*p)+1)):p*(i+1),((i*p)+1)+p*w:p*(i+1)+p*w):=T_{NAD}$

56: end for

57: end for

58: for $w$ from $1$ to $POD$ do

59: $a:=Q\left(\beta_{1}\left(1\right)+w,:\right)$

60: $e:=\beta_{1}\left(2\right)$

61: $column_{c}=zeros\left(p,1\right)$;

62: $col=1$;

63: while $e\leqslant 1$ do

64: $column_{c}\left(col\right)=a\left(e\right)$

65: $col=col+1$

66: $e=e-1$

67: end while

68: $e=\beta_{1}\left(2\right)$

69: $row_{c}=zeros\left(1,p\right)$

70: $r=1$

71: for $m$ from $e$ to $length\left(a\right)$ do

72: $row_{c}\left(r\right)=a\left(e\right)$

73: $r=r+1$

74: $e=e+1$

75: end for

76: $T_{POD}=toeplitz\left(column_{c},row_{c}\right)$

77: for $i$ from $0$ to $sizeTblock/\left(size\left(Q,2\right)\right)-1-w$ do

78: $T(((i*p)+1)+p*w:p*(i+1)+p*w,((i*p)+1):p*(i+1)):=T_{POD}$

79: end for

80: end for

81: $b:=zeros\left(size\left(T,1\right),1\right)$;

82: for $k$ from $\beta_{1}\left(3\right)$ to $A\left(3\right)$ do

83: $t:=1$;

84: for $i$ from $0$ to
$\left(size\left(F,2\right)-size\left(C,2\right)\right)$ do

85: for $j$ from $0$ to
$\left(size\left(F,1\right)-size\left(C,1\right)\right)$ do

86:
$f_{iter}:=F\left(j+1:size\left(C,1\right)+j,\left(i+1\right):size\left(C,2\right)+i,k-size\left(C,3\right)+1:k\right)$

87: $f_{new}:=f_{iter}.*C$

88:
$b\left(t\right)=-\left(\sum\left(\sum\left(\sum\left(f_{new}\right)\right)\right)\right)$

89: $t:=t+1$

90: end for

91: $elements:=linsolve\left(T,b\right)$

92: $m:=1$

93: for $i$ from $0$ to
$\left(size\left(F,2\right)-size\left(C,2\right)\right)$ do

94: for $j$ from $0$ to
$\left(size\left(F,1\right)-size\left(C,1\right)\right)$ do

95:
$F\left(\beta_{1}\left(1\right)+j,\beta_{1}\left(2\right)+i,k\right):=elements\left(m\right)$;

96: $m:=m+1$

97: end for

98: end for

99: end for

100: end forreturn $f\left(A\right)$

101:end procedure

The complete code of the program is available at
http://github.com/ApanovichMS/CauchyParallelepiped.git.

##  6 Conclusion

Classical methods for solving differential equations (methods of Runge-Kutta
of the fourth order, Euler, Newton, etc.) have certain difficulties in their
application in the multidimensional case. In each specific task, it is
necessary to individually select a solving method. This makes the development
of a universal approach to solving multidimensional differential equations
using these methods laborious. In this study, we expanded the possibility of
using standard symmetric difference schemes to approximate two-dimensional
differential equations with constant coefficients, where the problem is
considered "on the plane>>, in which, as a rule, time and one spatial variable
are used as independent variables. The scheme proposed by us, firstly, has an
arbitrary number of points, and secondly, it can be used to approximate
differential equations in the three-dimensional case. For the heat conduction
equation, this would mean that when describing heat transfer, it becomes
possible to take into account the anisotropy of this process by adding an
additional spatial variable.

The process of approximating differential equations by a difference scheme in
the three-dimensional case gives a system of difference equations, the matrix
of which has a block Töplitz form. The process of solving such systems for a
large number of variables and initial data is associated with high
computational costs. In our work, we proposed an algorithm for solving such a
system using values of the coefficients of the difference equation, which made
it possible to automate the solution of the Cauchy problem with the initial
data in the parallelepiped. Thus, our proposed approach defines a unified
algorithm for solving differential equations when they are approximated by a
difference scheme in a parallelepiped.

## Funding

This work is supported by the Krasnoyarsk Mathematical Center and financed by
the Ministry of Science and Higher Education of the Russian Federation in the
framework of the establishment and development of regional Centers for
Mathematics Research and Education (Agreement No. 075-02-2021-1388).

## References

  * Stanley (1997) Stanley, R. Enumerative combinatorics, Volume 1. Cambridge University Press, Cambridge, 1996.
  * Stanley (1999) Stanley, R. Enumerative combinatorics, Volume 2. Cambridge University Press, Cambridge, 1999.
  * Bousquet-Melou (2000) Bousquet-Melou, M.; Petkovsek, M. Linear recurrences with constant coefficients: the multivariate case. Discrete Mathematics. 2000, 225, 51–75.
  * Duffin (1956) Duffin, R.J. Basic Properties of Discrete Analytic Functions. Duke Math. J. 1956, 23, 335–363.
  * Duffin (1968) Duffin, R.J. Potential theory on rhombic lattice. J. Combinatorial Theory. 1968, 5, 258–272.
  * Danilov (2008) Danilov, O.A. Lagrange interpolating formula for discrete analytic function. NSU Bulletin. Series: Mathematics, Mechanics, Informatics. 2008, 8(4), 33–39. (in Russian)
  * Danilov (2009) Danilov, O.A.; Mednykh, A.D. Discrete analytical functions of several variables and Taylor expansion. NSU Bulletin. Series: Mathematics, Mechanics, Informatics. 2009, 9(2), 38–46. (in Russian)
  * Samarskii (1971) Samarskii, A.A. Introduction to the theory of difference schemes; Publisher: Nauka, Russia, 1971; 552. (in Russian)
  * Godunovi (1964) Godunov, S.K.; Ryaben’kii, V.S. The theory of difference schemes, an infroduction; Nort-Holland publishing company, Amsterdam, 1964.
  * Dudgeon (1983) Dudgeon, D; Mersereau, R.  Multidimensional Digital Signal Processing; Publisher: Prentice Hall, 1983; p. 400.
  * Isermann (1981) Isermann, R. Digital control systems, Berlin etc., 1981.
  * Abramov (2011) Abramov, S.A.; Gheffar, A.; Khmelnov, D.E. Rational solutions of linear difference equations: Universal denominators and denominator bounds. Programming and Computer Software. 2011, 37, 78–86.
  * Abramov (2015) Abramov, S.A. Search of rational solutions to differential and difference systems by means of formal series. Programming and Computer Software. 2015, 41, 65–73.
  * Abramov (2011) Abramov, S.A.; Barkatou, M.A.; van Hoeij, M.; Petkovsek, M. Subanalytic Solutions of Linear Difference Equationce and Multidimensional Hypergeometric Sequences. Journal of Symbolic Computation. 2011, 46, 1205–1228.
  * Abramov (2011) Abramov, S.; Petkovsek, M.; Ryabenko, A. Hypergeometric solutions of first-order linear difference systems with rational-function coefficients // Lecture Notes in Computer Science (including subseties Lecture Notes in Artificial Intelligence and Lecture Notes in Bioinformatics). 2015. V. 9301. P. 1–14
  * Leinartas (2009) Leinartas, E.K.; Lyapin, A.P. On the Rationality of Multidimentional Recusive Series. Journal of Siberian Federal University. Mathematics & Physics. 2009, 2(2), 449–455.
  * Kytmanov (2017) Kytmanov, A.A.; Lyapin, A.P.; Sadykov, T.M. Evaluating the rational generating function for the solution of the Cauchy problem for a two-dimensional difference equation with constant coefficients. Programming and computer software. 2017, 43(2), 105–111.
  * Apanovich (2021) Apanovich, M.S.; Lyapin, A.P.; Shadrin, K.V. Solving the Cauchy Problem for a Two-Dimensional Difference Equation at a Point Using Computer Algebra Methods. Programming and Computer Software. 2021, 47(1), 1–5.
  * Fedoruk (1987) Fedoruk, M.V. Asymptotics: Integrals and Series; Publisher: Nauka, Russia, 1987; 544. (in Russian)
  * Rogozina (2014) Rogozina, M.S. Solvability of a difference Cauchy problem for multi-layer implicit difference schemes. SibGAU Bulletin. 2014, 55(3), 126–130. (in Russian)
  * Apanovich (2017) Apanovich, M.S. Correctness of a two-dimensional Cauchy problem for a polynomial difference operator with constant coefficients. Journal of Siberian Federal University. Mathematics & Physics. 2017, 10(2), 199–205.
  * Apanovich (2018) Apanovich, M.S. On the solvability of a three-dimensional Cauchy problem. ITNOU: Information technology in science, education and management. 2018, 5, 9–13. (in Russian)
  * Iohvidov (1982) Iohvidov, I.S. Hankel and Toeplitz matrices and forms: Algebraic theory; Birkhauser, Boston,1982; 231.
  * Apanovich (2020) Apanovich, M.S.; Lyapin, A.P.; Shadrin, K.V. Calculating the sequence of main minors of the Toeplitz band matrix. Applied Mathematics & Physics. 2020, 52(1), 5–10 (in Russian).
  * MurrayJD (2003) Murray, J.D. Mathematical biology. II: Spatial models and biomedical applications; Springer-Verlag, Berlin, Heidelberg, 2003; 811.
  * bratus (2010) Bratus, A.S.; Novozhilov, A.S.; Platonov, A.P. Dynamical systems and models in biology; Fizmatlit, Moscow, 2010; 400 (in Russian).