Approximation of the Multidimensional Optimal Control Problem for the Heat Equation (Applicable to Computational Fluid Dynamics (CFD))

This work is devoted to finding an estimate of the convergence rate of an algorithm for numerically solving the optimal control problem for the three-dimensional heat equation. An important aspect of the work is not only the establishment of convergence of solutions of a sequence of discrete problems to the solution of the original differential problem, but the determination of the order of convergence, which plays a very important role in applications. The paper uses the discretization method of the differential problem and the method of integral estimates. The reduction of a differential multidimensional mixed problem to a difference one is based on the approximation of the desired solution and its derivatives by difference expressions, for which the error of such an approximation is known. The idea of using integral estimates is typical for such problems, but in the multidimensional case significant technical difficulties arise. To estimate errors, we used multidimensional analogues of the integration formula by parts, Friedrichs and Poincare inequalities. The technique used in this work can be applied under some additional assumptions, and for nonlinear multidimensional mixed problems of parabolic type. To find a numerical solution, the variable direction method is used for the difference problem of a parabolic type equation. The resulting algorithm is implemented using program code written in the Python 3.7 programming language.


Introduction
The heat equation is used to find the dependence of the temperature of the medium on the spatial coordinates and time, for given coefficients of heat capacity and heat conductivity. This is a second order partial differential equation, which is a parabolic type equation. Since the need to determine temperatures in the whole space is often absent, when setting the problem, additional conditions are introduced that determine the restrictions on the solution of the problem for a given area. For example, one of these conditions is to set the temperature distribution at the boundary of the region (the Dirichlet problem).
The process of finding the temperature distribution at given times is a laborious task. Since differential problems of a continuous nature cannot be programmed due to the limited capabilities of computer technology, such problems, by discretizing them, reduce to similar difference problems. Such a transition is carried out using difference schemes.
The main task of approximation is to find such an approximate function that least, in a certain sense, deviates from a given continuous function. Due to the fact that when solving continuous problems, the differential operators are replaced by finite-difference analogues, which are written in the form of algebraic equations, problems arise for determining the convergence and approximation error.
Note that when switching from a differential operator to a finite-difference analogue, a numerical solution is obtained that differs from the original solution. In such cases, an analysis is performed that determines the approximation order. For example, in Godunov and Ryabenkii (1987) study, the one-dimensional optimal control problem of the heat conduction process and the gradient descent method are considered, on the basis of which the approximation order of the finite-difference problem was obtained [1]. The optimality criterion is based on the gradient descent method, ideas leading to the assertion of the type of maximum principle by L. S. Pontryagin [2][3][4], lead to significant complications and are not considered in this paper. Approximation of optimization problems is considered by many researchers. An important work is Serovaiskii (2013) [5], from which methods for obtaining estimates of the boundedness of the target functional are used.
In modern works, attention is paid to the convergence of functionals in optimization problems of different nature. A hyperbolic boundary value problem with a quadratic cost functional is considered in Edalatzadeh et al. (2020) study [6]. An important point is the use of a similar technique of integral estimates to obtain optimal control in an explicit form. Criteria for the existence of optimal forms in Banach spaces were established in Edalatzadeh (2016 and 2019) studies [7,8]. For a differential operator in divergent form and for an integro-differential operator in Deligiannidis et al. (2020) and Mukam and Tambue (2020) researches [9,10] using integral estimates in suitable spaces, weak convergence of the numerical method was established and the order of convergence of the functional sequence to the solution was found.
The technique developed in this paper will be transferred to parabolic problems with variable coefficients, as well as to nonlinear cases. The possibility of such a step is considered plausible due to the Guillén-González et al. (2020) and Biccari et al. (2020) works [11,12].
The aim of this work is to estimate the approximation of a finite-difference analogue for the heat equation of three spatial variables. The solution to the difference problem is constructed using the variable direction method.
Note that the original result on the convergence estimation of the sequence of the target functional in 3dimensional space is established and the constants in the O symbols are directly calculated. We can briefly formulate the sequence of actions and steps that are used in the work: • Statement of the differential problem; • Analysis of the differential problem, obtaining an estimate of the norm of the solution depending on the control function; • Building a sequence of discrete tasks; • Obtaining expressions for errors between solutions to differential and discrete problems; • Estimation of errors using the technique of Sobolev spaces and the establishment of target inequality.
Based on the discretization of the three-dimensional heat conduction problem, a numerical algorithm is developed, with the help of which a software package is created to determine the time required for uniform distribution of heat in the rod.

The Problem Statement
The following is a third-order differential heat equation that describes the process of heating a body in space: For which the following boundary conditions are given; 3 | = ( , , , ), 0 < < ; Where ( , , , )is the solution of the boundary value problem, ( , , , ) − is a control function that shows the temperature at the point ( , , ), at the moment of time , 2 − is the thermal conductivity coefficient, ( , , ) − is the temperature of the rod at the initial moment of time at each point, ( , , , ) − is a given function from 2 [(0, ) × (0, ) × (0, )]. Questions of representation of solutions, existence and uniqueness are stated in Vladimirov (1981) and Shubin (2003) works [13,14].
We denote that the control belongs to the following set: Such a problem is called the Dirichlet problem or the first boundary value problem. We find a numerical solution to this problem using numerical methods, namely, the finite difference method. By expanding the function in a Taylor series, the first and second partial derivatives are expressed, and the boundary conditions are used to determine the value of the nodes on the boundary region.
The task is to find a function ( , , , ; ), such that on the whole region 2 [(0, ) × (0, ) × (0, )] by the time we get the distribution function heat close to the given function ( , , ). The criterion for this difference problem has the form: And the boundary conditions are rewritten as follows: In this case, it is necessary to go to the finite-difference analogue of the function ( , , , ; ) and evaluate the approximation order.

Equation of a Parabolic Type
In this paper, we consider the process of temperature distribution over a three-dimensional rod with a length, height, and width equal to , , , respectively, for the time interval , which is described by the heat equation. An inhomogeneous equation is considered: Which has coefficient 2 = 1 and boundary conditions (5).
( ) = 1 cos ( ) , = 1,2 … To obtain a complete orthonormal system, we define . To do this, take the scalar product from (9), equate it to 1 and find the integral. We get that 1 = √ 2 and: Similarly, we find a generalized solution for ( ) and ( ).
We apply the variational constant method. We solve the corresponding homogeneous equation and find a generalized solution in which is an arbitrary constant on .
We put this in Equation (9) and we obtain that for the unknown function ( ) the equality ′ ( ) − 2 = ( ) must be satisfied. We get that.
Whence the solution of the Cauchy problem is given by the formula: We obtain a formula for calculating the expansion coefficients ( , , , ) in eigenfunctions. Given the orthogonality of the Sturm-Liouville problem, we obtain: Hence, on the basis of (7), (9), (11) and (12), we obtain: , and ( ) is equal to (14).

Discretization of the Problem
The difference minimization problem has the following form. It is necessary to minimize the objective function Following works [15] and [16] we perform discretization and obtain difference problems.

Theoretical Information
In the course of performing mathematical analysis, a number of theorems, equations, and inequalities were used that play a fundamental role or are often used in mathematical calculations and simplifications.

Analysis of the Differential Problem
We begin the analysis of the differential problem by deriving two estimates for sufficiently smooth classical solutions to problem (5), (6), which will be emphasized in future work. Further actions are based on functional inequalities, which are sufficiently developed in Vasilev (2002) study [17].
We multiply equation (1.6) by ( , , , ; ) and integrate the resulting equality over the rectangle In view of conditions (2), we transform the first term from the left-hand side: To estimate the second term, we introduce each term of the Laplace operator under the differential sign, after which we apply the boundary conditions (5). As a result, we have: We use the Cauchy-Bunyakovsky formula (21) for the right-hand side of equality (27), after which we pass to the maximum in time for the classical solution of problem (5), (6). We have: We replace the terms in Equation (27) in accordance with formulas (28), (29) and (30), we have: Let us estimate this inequality. To do this, we remove each term from the right-hand side in turn. Based on this, we evaluate the first term.
Therefore, if we take the integral of the square of the function with respect to the maximum on the interval [0, ], square and extract the square root, and then use the estimate for the first term, we obtain the following inequality: From (31) we make an estimate for the second term, taking into account the estimate (32), we have: Based on inequality (31) and estimates (32) and (33), we obtain the first estimate for a sufficiently smooth solution to problem (31) and (32): Multiply equation (32) by and integrate over the domain : We estimate the first scalar product from the right-hand side. To do this, we introduce each term of the Laplace operator under the differential sign. As a result, we get: Taking into account the boundary conditions (5), only the last integral does not vanish. If we introduce the derivative of the function with respect to each variable under the differential sign and use the main theorem of mathematical analysis, we get: As a result, we obtain the following equality: Based on formula (36) and the elementary inequality of paragraph 1.4 for the product from formula (35), we have; 1 2 Or: Hence we have 2 inequalities: We use the fact that takes any value on the interval [0, ], we get: In addition, if we integrate equation (6) over the domain taking into account (37), we have: Adding inequalities (37) and (38) we obtain the second estimate for a sufficiently smooth solution: We use the Friedrichs inequality and inequalities (37), and also taking the maximum in time for differentials with respect to spatial variables, we estimate the square of the solution with respect to the control function: ∀( , , , ) ∈ , max = max{ , , }.
From here we get the energy estimate:

Analysis of a Discrete Task
Using analogues with estimates (34) and (39), we derive the corresponding estimates for the discrete problem. We multiply Equation (19) by ℎ ℎ ℎ = ℎ and sum over , , from 1 to ℎ − 1 = ℎ ̅̅̅ , from 1 to ℎ − 1 = ℎ ̅ , from 1 to ℎ − 1 = ℎ ̅̅̅ respectively: It is easy to verify that; Similarly, the formula is applicable to the spatial variables and .
We substitute formulas (43) and (44) in the formula (41): Inequality (45) is summed over from 1 to some , where on the interval 1 ≤ ≤ . We use the boundary condition 0 = 0, = 0. . ℎ ̅̅̅̅̅̅̅ , = 0. . ℎ ̅̅̅̅̅̅̅ , = 0. . ℎ ̅̅̅̅̅̅̅ . Then, if we expand the right-hand side of inequality (45) according to the Cauchy-Bunyakovsky formula (21) and make the maximum transition in time for a discrete solution, we obtain a difference analogue of the inequality (30): If we add inequalities (47) and (48) We use the summation formula in parts (20) in accordance with formula (44) and the boundary conditions (19) to estimate the second term from the left-hand side. We get: The left side of inequality (53) is summed over from 1 to some , where is in the interval 1 ≤ ≤ . Given

Evaluation of the Difference of Differential Decision and Discrete Analogue
We introduce the Hilbert space 2ℎ ℎ ℎ = 2ℎ , which is the difference analogue of the space 2 ( ). The elements of this space will be the grid functions ℎ ℎ ℎ = ℎ = By ℎ ℎ we denote the piecewise constant continuation of the grid function ℎ according to the rule; The domain of the function ℎ ℎ is denoted by ℎ = {( , , , ): ℎ ≤ ≤ , ℎ ≤ ≤ , ℎ ≤ ≤ , 0 ≤ ≤ }. We note that: Based on (60), (61), we rewrite the difference equation (16): ℎ ̅ ℎ − ℎ ( ̅ ℎ + ̅ ℎ + ̅ ℎ ) = ℎ ℎ , ( , , , ) ∈ ℎ (62) Subtract (62) from equation (6), п multiply the resulting equality by − ℎ ℎ and integrate over the domain ℎ : We estimate the first term from the left side of the equality (63). We replace the integration over the entire domain with summation in accordance with the formula For the first sum, we go through the cycle in the time variable, opening the squares of the difference and using the boundary condition 0 = 0, and for the second, we calculate the time integral for the third term, substitute the limit values in the first and third elements of the term bracket. Then the final inequality takes the following form: We transform the second term from the left-hand side of (64). We note that; ∑ ∫ (( We transform the first term from the right-hand side of (65). To do this, we take out the differential from the first bracket, we apply integration by parts. For the part of the expression in which the limit values are substituted, we will go through the cycle in the variable . Then, having completed the mathematical operations, we arrive at the following inequality: The third term from the right-hand side of (65), and using the formula for summing by parts (17), can be represented as follows: Performing similar mathematical operations (65)-(67), we can obtain estimates for the variables and , replacing the variable with another spatial variable, taking into account the limits of summation and integration.
We substitute the obtained inequality (64) and equality (65) taking into account (66) and (67) into (61). We get: Before estimating | |, it is necessary to introduce some more auxiliary inequalities.
If we take the function ̅ ( +1 , , , ) and its analogue ̅ , is discrete, expanding them with respect to the variable , and summing from to , we get the following 2 equalities: Where 1 ≤ ≤ ≤ ℎ ̅̅̅ ; for = y definition, we consider the sum in any of the equalities to be 0.
We perform the same operations to determine the lower limit, but change the arguments for the criteria of the differential and discrete problems. We get: * − ℎ * ≥ ( * ) − ℎ ( ℎ * ) ≥ − √(ℎ + ℎ + ℎ + Based on the formula (113), the first term is neglected due to its lesser influence on the right side in comparison with others. For the remaining elements of the right-hand side, the conclusions remain similar to the conclusions for the formula (114).

Conclusion
The problem of determining the approximation order of the optimal control problem for the spatial process of heat conduction is considered in the paper. Using the methods of integral inequalities and the method of difference approximation, a difference problem is obtained, an algorithm for finding its solution is described, and an estimate is obtained for the deviation of the value of the difference functional from the continuous functional. The established inequality Equation (118) gives an idea of the time complexity of the process of calculating an approximate solution when the accuracy of calculations is given in advance. The time steps and spatial variables are not independent, additional restrictions must be imposed to ensure stability. The methodology for obtaining an approximation estimate can be used for implicit approximations, hybrid schemes. The methods used in this article can be successfully applied for similar parabolic problems with bounded coefficients, as well as for problems of large dimensions.

Conflicts of Interest
The authors declare no conflict of interest.