Baixe Método dos Mínimos Quadrados na Modelagem de Trocas de Calor Terra-Ar em Ductos e outras Trabalhos em PDF para Algoritmos, somente na Docsity!
A Matlab code to fit periodic data
Ruth da Silva Brum^1 Jairo Valões de Alencar Ramalho^2 Luiz Alberto Oliveira Rocha^1 Liércio André Isoldi^3 Elizaldo Domingues dos Santos^3
Abstract: This paper presents a computer method to find the best sine-based function, in the sense of least squares, to fit periodic data. Even though the least squares method is not a novelty, there is a void in the literature about its use to find trigonometric functions, particularly when it gives rise to nonlinear systems, as it is done in this article. The respective code, implemented in the Matlab programming language, is detailed and analyzed exploring experimental data from the air and soil temperatures measured along the year in earth air heat exchangers (EAHE) built in the facilities of a case study house in the Brazilian state of Rio Grande do Sul. The fitting curves have been employed by the authors in different works to define boundary conditions to study new computer models of EAHE.
Keywords: Earth air heat exchanger. Least squares. Trigonometric curve fitting.
Resumo: Este artigo apresenta um método computacional para encontrar uma função do tipo senoidal que melhor se adequa, no sentido dos mínimos quadrados, a um conjunto de dados periódicos. Embora o método dos mínimos quadrados não seja uma novidade, foi percebida uma lacuna na literatura relativa ao seu uso para a obtenção de funções trigonométricas, particularmente envolvendo sistemas não lineares, como é feito neste trabalho. O código proposto, escrito na linguagem de programação do Matlab, é detalhado e avaliado por meio da exploração de dados empíricos e numéricos das temperaturas do ar e do solo medidos em trocadores de calor solo-ar (TCSA) instalados em uma casa experimental no interior do estado brasileiro do Rio Grande do Sul. As curvas ajustadas estão sendo empregadas pelos autores em diferentes pesquisas na definição de condições de contorno para novos modelos computacionais de TCSA.
Palavras-chave: Ajuste de funções trigonométricas. Mínimos quadrados. Trocadores de calor solo-ar.
1 Introduction
To face the increasing energy consumption, it is necessary to consider not only alternative/renewable energy sources and devices, but also to review the building construction models. As a remarkable example, in the summer of 2014, there was a change in the electrical energy demand peak in Brazil. It occurred at noon, instead of late afternoon, due to a heavy use of air conditioners during a heat wave [1].
In this context, the Earth Air Heat Exchangers (EAHE) stand out as a natural solution to improve air conditioning at low energy cost. In general, EAHE are composed by a system of fans and buried ducts. The fans blow the air inside the ducts, where it exchanges heat with the soil and returns at a milder temperature [2, 3, 4].
(^1) Programa de Pós-Graduação em Engenharia Mecânica, UFRGS, Porto Alegre (RS), Brasil. {[email protected], [email protected]} (^2) Departamento de Matemática e Estatística, UFPel, Capão do Leão (RS), Brasil.
{[email protected]} (^3) Programa de Pós-Graduação em Modelagem Computacional, Furg, Rio Grande (RS), Brasil.
{liercioisoldi,[email protected]}
http://dx.doi.org/10.5335/rbca.2015.
All this happens because the soil acts as a seasonal thermal reservoir, storing heat during summer and releasing it to the atmosphere along the winter. Hence, in the winter, the soil is warmer than the external air, what allows heating the air inside buried ducts. On the contrary, in the summer, the soil is colder than the environment air, reducing the air temperature inside the ducts [2].
On the other hand, to develop new computer models for EAHE it is necessary to know experimentally the air and soil temperatures of a given place. In particular, such states as the Rio Grande do Sul (RS) in Brazil, where the four seasons are well defined, it is important to measure this data for at least one year, as it is done in [2].
After plotting the results in a graph, it is noticeable that the annual temperatures in RS are periodic following the pattern of a translated sine (or cosine) curve in the interval I = [0, 2π]. Thus, the year begins in summer when the environment air temperatures reach their peak. Between March and June, the temperatures drop slowly (which characterizes the fall) to reach a minimum between June and September (in the winter). From September to December the temperatures go rising again (in the spring), and the year ends as it started with high summer temperatures. The soil temperatures follow an analogous behavior to the air ones, although there is a lag between them which depends on the soil depth.
Given that, it is very natural to look for the “best” sine-based function which represents (models) the air or soil temperatures along the year. To find it, this work employed the experimental data from [2], where the temperatures were measured daily, hour by hour. However, choosing to adopt the term “best” in the sense of least squares, the authors found a gap in classical literature [5-10] concerning the expected trigonometric fitting.
To be more specific, there is a lack of texts exploiting the problem of finding a real continuous function f(x) = a sen(bx+c)+d (with a, b, c, and d being real constants) to fit, by the least squares method, a set of n ordered pairs U={(xi , yi), i = 1, 2, …, n}. Despite using a trigonometric function, instead of a high order polynomial, brings the drawback of solving a non-linear system of equations to find the coefficients a, b, c, and d, this research found the Newton method very effective to circumvent this matter in all the simulated cases.
In general, this paper addresses two objectives: fill a void in the literature about trigonometric fitting; illustrate its application based on a computer code which can help the study and simulation of EAHE.
2 Methodology
2.1 Least Squares Method
Considering a set of ordered pairs {(xi, yi), i = 1, 2, …, n}, relative to the temperatures values yi measured at n regular time intervals xi, one can find a representative continuous function f to fit this data, as in Figure 1, by minimizing the sum of the squared differences between yi and f(xi).
Figure 1: An example of a fitting function
u(ξ 0 ) + ua(ξ 0 )(a-a 0 ) + ub(ξ 0 )(b-b 0 ) + uc(ξ 0 )(c-c 0 ) + ud(ξ 0 )(d-d 0 ) = 0, (12)
v(ξ 0 ) + va(ξ 0 )(a-a 0 ) + vb(ξ 0 )(b-b 0 ) + vc(ξ 0 )(c-c 0 ) + vd(ξ 0 )(d-d 0 ) = 0 , (13)
w(ξ 0 ) + wa(ξ 0 )(a-a 0 ) + wb(ξ 0 )(b-b 0 ) + wc(ξ 0 )(c-c 0 ) + wd(ξ 0 )(d-d 0 ) = 0 , (14)
y(ξ 0 ) + ya(ξ 0 )(a-a 0 ) + yb(ξ 0 )(b-b 0 ) + yc(ξ 0 )(c-c 0 ) + yd(ξ 0 )(d-d 0 )=0. (15)
In matrix form, the equations (12)-(15), can be written as:
ua(ξ 0 ) ub(ξ 0 ) uc(ξ 0 )^ ud(ξ 0 ) va(ξ 0 ) vb(ξ 0 ) vc(ξ 0 )^ vd(ξ 0 ) wa(ξ 0 ) ya(ξ 0 )
wb(ξ 0 ) yb(ξ 0 )
wc(ξ 0 ) yc(ξ 0 )
wd(ξ 0 ) yd(ξ 0 ) (^) )
a b c d
ua(ξ 0 ) ub(ξ 0 ) uc(ξ 0 )^ ud(ξ 0 ) va(ξ 0 ) vb(ξ 0 ) vc(ξ 0 )^ vd(ξ 0 ) wa(ξ 0 ) ya(ξ 0 )
wb(ξ 0 ) yb(ξ 0 )
wc(ξ 0 ) yc(ξ 0 )
wd(ξ 0 ) yd(ξ 0 ) (^) )
a 0 b 0 c 0 d 0
u(ξ 0 ) v(ξ 0 ) w(ξ 0 ) y(ξ 0 ) (^) )
The matrix of partial derivatives, which is repeated in both sides of equation (16), is known as the Jacobian matrix, J, of the vector field F. Thus, the equation (16) can be written simply as:
J(ξ 0 ) ξ= J(ξ 0 ) ξ 0 – F (ξ 0 ), (17)
where ξ = (a, b, c, d), and, as before, ξ 0 = (a 0 , b 0 , c 0 , d 0 ).
Following these formulations, an approximation to a root, xr, of a vector field F can be found iteratively by choosing a suitable initial approximation ξ 0 and computing the successive approximations:
J(ξk) ξk+1= J(ξk) ξk – F (ξk), (18)
for k = 0, 1, 2, …, etc. These resolutions should continue until reaching a maximum number of iterations or a convergence is achieved, that is, the difference in a chosen norm between successive solutions ξk and ξk+1 is under a specified tolerance. For this work, it was adopted the infinity norm [6].
This presented technique is called the Newton Method [5, 6]. In fact, to reduce computational costs, the matrix-vector multiplications in equation (18) can be avoided by introducing the variable t = ξk+1 – ξk. In this way, the algorithm for the Newton method is divided in two parts:
- Find t by solving equation J(ξk) t = – F (ξk),
- Compute ξk+1 = ξk + t. Under certain conditions the Newton method can converge “very fast” or not converge at all. Particularly, the convergence is very dependent on the initial approximation, which must be chosen carefully. Following the guidelines proposed in Section 3, for all the cases simulated in this research, the Newton method converged in less than five iterations.
Before proceeding, it is important to present the 4 by 4 Jacobian matrix, J, relative to the set of equations (6)-(9), which were obtained from the least squares methods. The elements, Jij, corresponding to the i-th row and j-th column of J, are given by:
J 11 = ∑^ n i=1 sin^2 (bxi+c), (19)
J 12 = ∑ n i=1 xi cos(bxi+c)[2a sin(bxi+c)+d-yi], (20)
J 13 = ∑ n i=1 cos(bxi+c)[2a sin(bxi+c)+d-yi], (21) J 14 = ∑ n i=1 sin(bxi+c), (22)
J 21 = ∑ n i=1 xi sin(bxi+c) cos(bxi+c), (23)
J 22 = ∑^ n i=1 xi^2 {a cos^2 (bxi+c)-[a sen(bxi+c)+d-yi] sen(bxi+c)}, (24)
J 23 = ∑^ n i=1 xi {a cos^2 (bxi+c)-[a sen(bxi+c)+d-yi] sen(bxi+c)}, (25) J 24 = ∑^ n i=1 xi cos(bxi+c), (26) J 31 = ∑ n i=1 sin(bxi+c) cos(bxi+c), (27) J 32 = J 23 , (28)
J 33 = ∑ n i=1 {a cos^2 (bxi+c)-[a sen(bxi+c)+d-yi] sen(bxi+c)}, (29) J 34 = ∑^ n i=1 cos(bxi+c), (30) J 41 = J 14 , (31) J 42 = ∑ n i=1 a xi cos(bxi+c), (32) J 43 = ∑ n i=1 a cos(bxi+c), (33) J 44 = n. (34)
2.3 Commented Matlab Code
To show how this methodology was implemented computationally, this subsection presents a corresponding code written in the Matlab (R2008a) programming language [14]. For those not used to Matlab, it is worth to read the following points.
Anything written to the right of the symbol % is taken as a commentary.
A very long command statement can be continued to the next lines using ellipsis (…). The function sum (), usually written with the syntax sum(x) , returns the sum of the elements of an array x. It is useful to compute summations without the need of a loop.
The function max (), written with the syntax max ( x ), returns the largest element of an array x.
It is possible to make operations for all the elements of an array at once. For instance, if x and y are two n×1 arrays, and k is a constant, the following statements: abs( x ), kx* , x .^2, x .* y , return n× arrays with, respectively, the absolute value of all elements of x , all the elements of x multiplied by k , all the squares of x , an element by element product of x and y.
A linear system Ax=b, where A and b are, respectively, n×n and n×1, arrays, can be solved simply by the command statement x=A\b;
As usual, Matlab counts with a conditional if-else construction and also a loop while. Even though the later has a syntax while (condition) do commands end , the loop terminates not only when the specified condition holds false, but also when a break statement is executed.
After this preamble, the commented code is given below.
% .................INPUT DATA........................................... % % Initial approximation for the sine coefficients a=7.5; b=0.0185; c=0.8316; d=20.5; % Tolerance and maximum number of iterations tol = 1e-5 ; Nmax = 10; % Time intervals (x) and Temperature values (y) % Number of time intervals n=1000; % Command to specify the path to the data files folder addpath ('\files_folder'); % Reading n lines of data from a file
Besides, it is worth to say that there is a high correlation coefficient, of 0.8693, between the fitted curve and the experimental data.
Figure 2: Air temperature at the inlet of an EAHE duct
As it was mentioned in the methodology, the Newton method converged fast, that is, in less than five iterations for all cases tested. However, it should be pointed that this depends on the quality of the initial approximation (a 0 , b 0 , c 0 , d 0 ) which can be fairly estimated from a previous analysis of the experimental data.
In this regard, computing an average of the temperatures along the summer, As, and along the winter, Aw, it is possible to estimate the amplitude term by a 0 = (As-Aw)/2. Certainly, there are other valid ways to estimate the same term. For instance, taking a look at the Figure 2, one can note the peaks of temperature during summer, in average, are around 28oC. In the same way, the average minimum temperatures during winter are close to 13 oC. Therefore, it is reasonable to assume an amplitude a 0 ≈ (28-13)/2 = 7.5.
As for the term d 0 , it represents the annual average temperatures. Thus, it can be calculated taking the mean value of the experimental data. A simple way to do it in Matlab is using the function mean(), whose usual syntax is mean(A), which computes the average or mean value of a given array A. It is worth to mention that one estimate can also be deduced from Figure 2. The value of d 0 is approximately an average between the maximum summer temperatures and the minimum winter ones, in other words, d 0 ≈ (28+13)/2 = 20.5.
Finally, to estimate the terms b 0 and c 0 , it was adopted the following procedure. For the external air, it is known empirically that the peak of summer temperatures occurs, usually, at the beginning of February (about forty days after the year start). On the other hand, the minimum winter temperatures occur in July (about two hundred and ten days after the year start). Since the sine function achieves its maximum value at the point π/ and its minimum at 3π/2, it is sufficient to make a correspondence among, respectively, the days x = 40 and x = 210, and the points π/2 and 3π/2. Thus, b 0 and c 0 are found solving the linear system:
40 b 0 + c 0 = π/2, (35) 210 b 0 + c 0 = 3π/2, (36)
whose solution is b 0 ≈ 0.0185 and c 0 ≈ 0.8316.
All these values are presented in the Matlab code above. In fact, running the program, it was obtained the following function to represent the temperature at the air inlet of the EAHE duct:
T(x) = 20.4967 + 5.6634 sin( 0.0178x + 0.9787 ). (37)
Therefore, except for the amplitude term, all the initial approximations were fairly precise.
On the other hand, one should highlight that the least squares and Newton methods are robust enough to compute the coefficients even when the initial approximations are rather rough. For instance, using the same previous estimates, it was obtained the following function to represent the temperature at the air inlet of other EAHE duct:
T(x) = 21.7895 - 5.9614 sin( -0.0184x + 5.3991 ), (38)
which is very different than the one in (37).
The Figure 3 displays a graph comparing the fitting curve with the experimental data for the temperature at this second inlet. As before, it is possible to see that the sine curve also matches the behavior of the temperature along the year. Besides, for this case, the correlation between the fitted curve and the experimental data is 0.8788, which is higher than the previous result.
This methodology was also employed to obtain the temperatures in the soil for different depths. These results were obtained from numerical data available using the EAHE model validated in [15]. In the Figure 4, there is a graph comparing the temperature curves, at the depth of 4 m, for the numerical model and the fitted one. In this case, the correlation between them is 0.9976. In fact, the curves fitted from the least squares method were used by the authors to construct new EAHE models presented in the articles [3], [4] and [15].
Figure 3: Air temperature at the inlet of a second EAHE duct
Energy Research , v. 1, p. 505 – 514, 2012. ISSN 2251-9408. Available in: <http://www.grsci- techpress.org/index.php/IJARER/article/view/62>.
[4] BRUM, R. S.; et al. A new computational modeling to predict the behavior of earth-air heat exchangers. Energy and Buildings , Elsevier, v. 64, p. 395 – 402, 2013. ISSN 0378-7788. Available in: < http://www.sciencedirect.com/science/article/pii/S0378778813003186>.
[5] GERALD, C. F. Applied Numerical Analysis. Menlo Park, CA, USA: Addison-Wesley Publishing Company, 1984.
[6] BURDEN, R. L.; FAIRES, J. D. Análise numérica. São Paulo, SP, Brasil: Pioneira Thomsom Learning,
[7] FRANCO, N. B. Cálculo Numérico. São Paulo, SP, Brasil: Pearson Prentice Hall, 2006.
[8] ARENALES, S; DAREZZO, A. Cálculo Numérico, Aprendizagem com Apoio de Software. São Paulo, SP, Brasil: Pioneira Thomsom Learning, 2008.
[9] MATHEWS, J. H. Numerical Methods for Mathematics, Science and Engineering. Upper Saddle River, NJ, USA: Prentice Hall, Inc., 1992.
[10] RUGIERO, M. A. G.; LOPES, V. L. R. Cálculo numérico, aspectos teóricos e computacionais. São Paulo, SP, Brasil: Makron Books, 1996.
[11] BRUM, R. S. Modelagem Computacional de Trocadores de Calor Solo-Ar. Dissertação (Mestrado em Modelagem Computacional) – Programa de Pós-Graduação em Modelagem Computacional da Universidade Federal do Rio Grande, 2013.
[12] STEWART, J. Cálculo Volume II. São Paulo, SP, Brasil: Pioneira Thomsom Learning, 2006.
[13] BARTLE, R. G. The elements of real analysis. New York, USA: John Wiley & Sons, 1976.
[14] CHAPMAN, S. J. Programação em Matlab para engenheiros. São Paulo, SP, Brasil: Cengage Learning,
[15] BRUM, R. S.; et al. Two-Dimensional Computational Modeling Of The Soil Thermal Behavior Due To The Incidence Of Solar Radiation. EngenhariaTérmica (Thermal Engineering) , DEMEC – UFPR, v. 12, n. 2, p. 63–68, 2013. ISSN 1676-1790. Available in: < http://demec.ufpr.br/reterm/ed_ant/21/>.
