Choice of the regularization parameter for the Cauchy problem for the Laplace equation

a

= multinomial coefficient of the function of distribution of temperature T̃(w);

c

= multinomial coefficient of the function of distribution of temperature T (x, y);

E (α)

= functional, energy integral;

J_α

= regularizing functional (J_α = J + α²I);

m

= number of Chebyshev nodes on the y-axis;

n

= number of Chebyshev nodes on the x-axis;

N₁–1

= degree of the polynomial describing unknown distribution of temperature on the Γ₁ boundary;

N_max

= norm;

q

= heat flux density, W/m²K;

T

= temperature, K;

T̃

= temperature, function dependent on the Chebyshev node;

w

= Chebyshev node;

W_i

= Chebyshev polynomial of the first kind of i-th degree;

x, y

= Cartesian coordinates;

[x]_n

= integer part of the division of number x by n; and

x mod n

= remainder of the division of number x by n.

α

= regularization parameter;

δ

= error;

δ_M

= error of measurement data (Morozov principle);

γ

= multinomial coefficient, pertains to the sought temperature distribution on the boundary Γ₁;

Γ

= boundary of the domain Ω, (Γ = Γ₁ ∪ Γ₂ ∪ Γ₃ ∪ Γ₄); and

Ω

= calculation domain.

A

= analytical solution;

c

= calculated value;

m

= measured value;

ran

= randomly disturbed value; and

Γ_i

= on the boundary Γ_i (for i = 1, 2, 3 and 4).

Inverse problems are ill-posed in the Hadamard (1902) sense. It means that a slight disturbance to measurement data results in significant errors of the obtained results (Joachimiak and Ciałkowski, 2017, 2014, 2018; Nowak, 2017). Therefore, problems of such type need to be regularized. There are many methods used to regularize inverse problems. Among them, there is the Tikhonov regularization (Beck and Woodbury, 2016; Chen et al., 2019; Djerrar et al., 2017; Frąckowiak et al., 2019a; Laneev, 2018; Marin, 2010, 2016; Niu et al., 2014; Sun, 2016; Tikhonov and Arsenin, 1977; Yaparova, 2016), the Tikhonov–Philips regularization (Joachimiak et al., 2019a), the discrete Fourier transform (Frąckowiak and Ciałkowski, 2018; Wróblewska et al., 2015) and SVD algorithm (Hasanov and Mukanova, 2015). In her article, Cheruvu (2017) applied the wavelet regularization of Laplace’s equation in the arbitrarily shaped domain. The solution to the Cauchy problem for the Laplace’s equation was also sought with the use of the iterative Tikhonov-type method (Delvare and Cimetière, 2017). Han et al. (2011) in the article presented numerical tests concerning the solution to the Cauchy problem for Laplace’s equation with the use of the energy regularization method. Obtained results were compared with the Tikhonov regularization for which the regularization parameter was chosen based on the Morozov principle. In the paper of Liu and Wang (2018), the Cauchy problem for the Laplace’s equation was solved with the use of the method of fundamental solutions and the energy regularization technique to choose the source points. The Laplace’s equation was also solved with the use of iterative algorithms (Frąckowiak et al., 2015a, 2015b), of the Trefftz method (Ciałkowski and Frąckowiak, 2002; Ciałkowski and Grysa, 2010; Grysa et al., 2012; Hożejowski, 2016; Lin et al., 2018), of the method of fundamental solution (Kołodziej and Mierzwiczak, 2008; Mierzwiczak et al., 2015; Mierzwiczak and Kołodziej, 2011) and of the collocation method (Joachimiak et al., 2016). In many cases, the regularization of the inverse problem concerns the problem of choosing the regularization parameter. The regularization parameter can be chosen based on the Morozov principle (Chen et al., 2019; Han et al., 2011; Joachimiak et al., 2019a; Marin, 2016; Morozov, 1984; Sun, 2016) or using the L-curve method (Jin and Zheng, 2006; Marin and Munteanu, 2010; Marin, 2005). In the study of Marin (2011), the optimal regularization parameter was sought based on the generalized cross-validation criterion. Currently, research work focuses on finding new methods of regularization (Cheng and Feng, 2014; Zhuang and Chen, 2017) and on the modification of already known and used methods (Yang et al., 2015; Zheng and Zhang, 2018). Because of a wide application of inverse problems in engineering problems, such as the cooling of the blades in gas turbines (Frąckowiak et al., 2017; Frąckowiak et al., 2019b; Frąckowiak et al., 2011), analysis of the boiling heat transfer in minichannels (Hożejowska et al., 2009; Maciejewska and Piasecka, 2017), analysis of thermal and thermo-chemical treatment (Joachimiak et al., 2019b) or monitoring of power boilers operation (Taler et al., 2016, 2017), developing methods for regularization of these problems and investigating the process of choosing the regularization parameter are very significant.

In this article, the solution to the Cauchy problem for the Laplace’s equation was investigated with the use of the Chebyshev polynomials. To regularize the solution to the inverse problem, the modification of the Tikhonov and of the Tikhonov–Philips regularizations, described in the article (Joachimiak et al., 2019a) was applied. The choice of the regularization parameter was made based on the Morozov principle, the minimum energy integral criterion and the L-curve method.

In many technical applications, it is impossible to measure temperature on the boundary of the heated component of the device or machine, such as the combustion chamber, the inner side of the body of a turbine or a boiler. Then, the distribution of temperature can be determined by finding the solution to the inverse problem. Based on the distribution of temperature on the part of the boundary $(TΓ2,TΓ3,TΓ4, fig. 1)$ and, additionally, knowing the heat flux density on the part of the boundary $(qΓ3, fig. 1)$ one can determine the distribution of temperature on the boundary, where it is impossible to measure this temperature $(TΓ1, fig. 1)$⁠. Such a posed problem is the Cauchy problem, particularly sensitive to errors in measurement and in the calculation. In the stationary thermal field, the heat equation is reduced to the Laplace’s equation (for the non-linear case, the Kirchhoff’s substitution transforms the heat equation into the Laplace’s equation).

Laplace’s equation is solved in the domain Ω = {(x, y) ∈ ℝ²: −1 ≤ x ≤ 1, −1 ≤ y ≤* 1} with the following boundary conditions (Figure 1):

It was assumed that the solution can be noted as the linear combination of the Chebyshev polynomials (Paszkowski, 1975).

To solve the Cauchy problem, the collocation method was used. It was assumed that there are n points along the x-axis and m points along the y-axis (including points on the boundary). Collocation points being inside the interval (−1, 1) are the Chebyshev nodes (Paszkowski, 1975). Nodes were renumbered, which enables the temperature function to be noted in the following equation (7):

where the coefficients a_k (k = 1, 2, …, mn) are unknown. Sought temperature distribution on the boundary Γ₁ was assumed as the linear combination of the Chebyshev polynomials (Paszkowski, 1975).

Coefficients a_k (k = 1, 2, …, mn) are expressed by the values of coefficients γ_i (i = 1, 2, …, N₁). Hence, the determination of the temperature distribution is reduced to the determination of coefficients γ_i. To do so, the functional of the following form was minimized.

where c in subscript denotes the calculated value, while m denotes the measured value. The integral in equation (9) on the boundary Γ₂ can be noted in the following equation (10):

Applying numerical integration we have:

where $Δx1=x2 − x12$⁠, $Δxi=xi+1 − xi−12$ for i = 2, 3, …, n−1 and $Δxn=xn−xn−12$⁠. Having inserted the equation (7) into the equation (11), we obtained:

Solving the direct problem, where the temperature on boundaries Γ₁, Γ₂, Γ₃ and Γ₄ was known, was reduced to solving the matrix equation

what was described in detail in the paper (Joachimiak et al., 2019a). Based on the solution to the direct problem, constants a_k [Equation (12)] are of the following equation (14):

where $Hh,k=∑j=2m−1A˜k,jnWh−1(yj)$⁠, while Ã_k,jn (k = 1, 2, …, mn; j = 2, 3, …, m − 1) are elements of the matrix A⁻¹. After substituting equation (14) into equation (12) we obtained:

We would like the integral J_Γ2 to have a value equal to zero or as close to zero as possible, hence, we equate the squared expression [Equation (15)] to zero. Hence, we have that:

We obtain n linear equations of the following equation (17):

for i = 1, 2, …, n and h = 1, 2, …, N₁. It can be reduced to the matrix equation.

Similarly, for other integrals [Equation (9)] we obtained:

After J_Γ3, J_q,Γ3 and J_Γ4 had been equated to zero, equations of the following forms were obtained:

Based on equations (16) and (22)-(24), an oversized system of linear equations was obtained as the matrix equation, which would be solved with the use of the SVD algorithm:

what can be noted in the shorter form:

Because of a great sensitivity of results to disturbances to measurement data, the Cauchy problem was regularized. The regularizing functional of the following form was assumed:

Regularization term can be noted as the sum of integrals.

where $Γ1=∪i=1m−1Γ1i$⁠. Performing numerical integration using the trapezoidal rule, we obtained:

On the boundary Γ₁ we have:

Hence,

Each of the components [Equation (32)] was equated to zero. The equation of the following form was obtained:

Equations (33)-(36) can be reduced to the following system of equations.

where α is the regularization parameter. When the equation (26) and regularization [Equation (37)] are included, the following system of equations is obtained:

The solution to the system of equations (38) was sought in the least-squares sense with the use of the SVD algorithm.

To determine unknown regularization parameter α, the Morozov principle, the minimum of energy integral criterion and L-curve method were applied. For the solution obtained with the use of the Morozov principle, the mean (Morozov_A) and the maximal (Morozov_B) errors of the heat flux density δ_M on the boundary Γ₃ were evaluated. Interval halving method was used to determine zero of the function F_M (α) defined by the following equation (39):

Unknown regularization parameter α was also sought based on the minimization of the functional (energy integral) of the following equation (40):

where $(∇T(α))2=(∂T(α)∂x)2+(∂T(α)∂y)2$⁠. The minimum of the energy integral corresponds to satisfying the Laplace’s equation (with respective boundary conditions), which is discussed in the paper (Gelfand and Fomin, 1979). Therefore, we choose the parameter α for which $min⁡αE(α)$ occurs, i.e. the derivative $E'=dEdα$ reverses the sign. Domain Ω was divided into rectangular domains with the use of equidistant nodes and next into domains being right-angled triangles. The integral value was calculated with the use of the finite element method. Value ∇T was determined based on the form of the solution equation (6). Values γ_h were obtained by solving the equation (38).

On the basis of the solution of the equation (38), the L-curve was drawn as the correlation between ‖[BM]{γ} − {BW}‖ on the x-axis and ‖[CM]{γ} − {CW}‖ on the y-axis. We sought for the regularization parameter α with which corresponded the point of the L-curve locating on the curvature of this line. To evaluate the choice of the regularization parameter α, the following norm was defined:

Calculations were made in the domain Ω for the function.

and

We assumed such values of constants a and b that singularities would be outside the calculation domain Ω and that the values of gradients would change significantly within this domain (a = 1.3, b = 1.3, b = 1.1). Values of the norm N_max [Equation (41)], not including the regularization [Equation (26)], without disturbance (δ_ran = 0) and with random disturbance to the heat flux density up to 0.01q (δ_ran = 0.01) and to 0.02q (δ_ran = 0.02) are summarized in Table I. Disturbance q is an additive function with the uniform distribution. A slight disturbance to measurement data results in a significant error of the sought temperature on the boundary Γ₁. Hence, it is necessary to regularize the inverse problem [Equation (38)] and to choose the regularization parameter α properly.

Calculations were made for the function f₁ [Equation (42)]. Heat flux density was disturbed randomly to 0.02q (δ_ran = 0.02). Regularization parameter α was chosen with the use of the Morozov principle, the minimum of energy integral criterion and the L-curve method.

Figure 2 presents the course of the energy integral and its derivative depending on the parameter α. To solve the Cauchy problem, the authors applied such value of the regularization parameter α for which the energy integral E (α) took the minimal value, which meant that the derivative of the energy integral $E'=dEdα$ reversed the sign. The value of α was 5.13 × 10⁻⁴ (Table II). The Cauchy problem was also regularized for the regularization parameter α amounting to 4.64 × 10⁻¹, which was determined based on the L-curve course (Figure 3).

To choose the regularization parameter with the Morozov principle, the values of mean and maximal error δ_M were evaluated for the heat flux on the boundary Γ₃. The respective values were obtained: 0.008 and 0.02 (Table II). Next, zero of the function F_M (α) was calculated as per the equation (39). The respective values of the regularization parameter were obtained: 1.98 × 10⁻⁶ and 4.0799 × 10⁻² (Table II).

The lowest value of the norm N_max amounting to 6.18 × 10⁻² (Table II) for the function f₁ was obtained for the case of choosing the regularization parameter with the use of the Morozov principle for the maximal error of the heat flux δ_M (Morozov_B). This criterion brought satisfying results, as did the choice of the regularization parameter made with the use of the minimum energy integral criterion (N_max = 9.796 × 10⁻²). When the L-curve method was used, the obtained results were considerably worse (N_max = 2.22 × 10⁻¹). For the Morozov principle, for the mean error δ_M of the heat flux (Morozov_A), the highest value of the norm N_max amounting to 50.42 was obtained. Distributions of temperature on the boundary Γ₁ resulting from the analytical solution (AS) and from the solution of the Cauchy problem are presented in Figure 4.

Calculations were made for the function f₂ [Equation (43)]. Heat flux density was disturbed randomly to 0.02q (δ_ran = 0.02). The best results were obtained for the choice of the regularization parameter made with the use of the minimum energy integral criterion (N_max = 4.77 × 10⁻²), and the worse results were obtained for the L-curve method (N_max = 0.361). Values of the regularization parameter and of the norm N_max, being the measure of the quality of the parameter α choice, for the function f₂ are summarized in Table III. Distributions of temperature on the boundary Γ₁ for the AS and for the solution to the Cauchy problem with regularization are presented in Figure 5. Distribution of temperature on the boundary Γ₁ obtained with the use of the minimum of energy integral criterion slightly diverges from the AS.

To examine thoroughly the criterion for the regularization parameter selection with the use of the minimum of energy integral, calculations were performer also for the following functions:

which were chosen based on publications (Liu et al., 2018; Conde Mones et al., 2017; Fu et al., 2013; Sun, 2017). Values of the regularization parameter and of the norm N_max for functions f₃ − f₆ are summarized in Table IV. For the function f₆ and the disturbance to the heat flux density δ_ran = 0.02 and δ_ran = 0.05 the minimum of energy integral was not achieved. Distributions of temperature on the boundary Γ₁ being sought are presented in Figure 6. For functions f₃ − f₅ the disturbance δ_ran = 0.05 was taken into account, and for the function f₆ it was δ_ran = 0.01.

This paper presents the solution of the Cauchy problem for Laplace’s equation. Obtained distributions of temperature on the boundary Γ₁ were analyzed in terms of the dependence on the method for choosing the regularization parameter. The best results were obtained for the choice of the regularization parameter made with the use of the minimum of energy integral criterion and the Morozov principle (δ_M is the maximal error for the heat flux on the boundary Γ₃). The advantage of the application of the minimum energy integral criterion is a unique determination of the regularization parameter α for which E (α) has minimal value. Regularization made with the use of the minimum energy integral criterion gives satisfying results. However, its disadvantage is the fact that not for all calculation examples the minimum energy integral was determined. For the Morozov principle, the obtained results (distribution of temperature on the boundary Γ₁) depend on calculation or evaluation of the value of the heat flux δ_M error, what is a disadvantage of this method. Inexact evaluation of the measurement data error can result in obtaining the distribution of temperature on the boundary Γ₁, which is subject to great uncertainty. The choice of the parameter α made with the use of the L-curve method gave the worst results. Smooth L-curve course was obtained, what was related to the problem with the unique determination of the regularization parameter α using this method.