Accepted Manuscript Title: Metaheuristic algorithms for approximate solution to ordinary differential equations of longitudinal fins having various profiles Author: Ali Sadollah Younghwan Choi Joong Hoon Kim PII: DOI: Reference:
S1568-4946(15)00276-8 http://dx.doi.org/doi:10.1016/j.asoc.2015.04.049 ASOC 2939
To appear in:
Applied Soft Computing
Received date: Revised date: Accepted date:
20-10-2014 19-2-2015 24-4-2015
Please cite this article as: A. Sadollah, Y. Choi, J.H. Kim, Metaheuristic algorithms for approximate solution to ordinary differential equations of longitudinal fins having various profiles, Applied Soft Computing Journal (2015), http://dx.doi.org/10.1016/j.asoc.2015.04.049 This is a PDF file of an unedited manuscript that has been accepted for publication. As a service to our customers we are providing this early version of the manuscript. The manuscript will undergo copyediting, typesetting, and review of the resulting proof before it is published in its final form. Please note that during the production process errors may be discovered which could affect the content, and all legal disclaimers that apply to the journal pertain.
Highlights
1
•
Approximate solutions to ordinary differential equations (ODEs) in engineering.
3
•
Fourier series with the aid of metaheuristics used as suggested approximate method.
4
•
The GA, PSO and HS are utilized for optimization purposes.
5
•
Obtained approximate solutions by the proposed method confirm the results by others.
6
•
Proposed approach offers acceptable accuracy for a wide range of ODEs.
us
cr
ip t
2
Ac ce p
te
d
M
an
7
1 Page 1 of 52
Metaheuristic algorithms for approximate solution to ordinary differential
8
equations of longitudinal fins having various profiles
9
Ali Sadollah, Younghwan Choi, Joong Hoon Kim∗
10
School of Civil, Environmental, and Architectural Engineering, Korea University, 136-713, Seoul,
11
South Korea
ip t
7
Abstract
13
Differential equations play a noticeable role in engineering, physics, economics, and other
14
disciplines. Approximate approaches have been utilized when obtaining analytical (exact)
15
solutions requires substantial computational effort and often is not an attainable task. Hence,
16
the importance of approximation methods, particularly, metaheuristic algorithms is
17
understood. In this paper, a novel approach is suggested for solving engineering ordinary
18
differential equations (ODEs). With the aid of certain fundamental concepts of mathematics,
19
Fourier series expansion, and metaheuristic methods, ODEs can be represented as an
20
optimization problem. The target is to minimize the weighted residual function (error
21
function) of the ODEs. The boundary and initial values of ODEs are considered as constraints
22
for the optimization model. Generational distance and inverted generational distance metrics
23
are used for evaluation and assessment of the approximate solutions versus the exact
24
(numerical) solutions. Longitudinal fins having rectangular, trapezoidal, and concave
25
parabolic profiles are considered as studied ODEs. The optimization task is carried out using
26
three different optimizers, including the genetic algorithm, the particle swarm optimization,
27
and the harmony search. The approximate solutions obtained are compared with the
28
differential transformation method (DTM) and exact (numerical) solutions. The optimization
29
results obtained show that the suggested approach can be successfully applied for
Ac ce p
te
d
M
an
us
cr
12
∗
Corresponding author. Tel: +82-02-3290-3316. E-mail:
[email protected] (J.H. Kim)
2 Page 2 of 52
approximate solving of engineering ODEs. Providing acceptable accuracy of the proposed
31
technique is considered as its important advantage against other approximate methods and
32
may be an alternative approach for approximate solving of ODEs.
33
Keywords
34
Metaheuristics, Analytical solution, Weighted residual function, Approximate solution,
35
Fourier series, Longitudinal fins.
cr
ip t
30
1. Introduction
us
36
Mathematical formulation of most physical and engineering problems involves differential
38
equations (DEs). The DEs have applications in all areas of science and engineering. Hence, it
39
is important for engineers and scientists to know how to deal with the DEs and solve them.
40
With regard to real life problems, which are highly nonlinear, many problems in engineering
41
and science often include one or more ordinary differential equations (ODEs).
42
Analytical approaches are often inefficient in tackling ODEs. Therefore, approximate
43
analytical methods are applied to obtain estimated solutions of ODEs. A number of analytical
44
methods have been utilized to develop approximate analytical solutions for engineering
45
problems, such as the variational iteration method (VIM) [1, 2], the homotopy analysis
46
method (HAM) [3], the method of bilaterally bounded (MBB) [4], and the adomian double
47
decomposition method (ADM) [5, 6].
48
Differential transformation method (DTM), which is based on Taylor series expansion, as
49
widely used approximate analytical method, was first introduced by Zhou [7] for solving
50
linear and nonlinear initial value problems (IVPs) in electrical circuits. The DTM has been
51
widely used to obtain accurate analytical solutions for nonlinear engineering problems [8-11].
52
Recently, many studies have combined the concept of the DTM with finite difference
53
approximation for increasing the capability of their approximate solution [12-14]. Further,
Ac ce p
te
d
M
an
37
3 Page 3 of 52
approximate analytical procedures such as the VIM, homotopy perturbation method (HPM),
55
and DTM have been combined with Padé approximation technique to overcome the
56
disadvantages faced by these methods in certain cases [15].
57
In particular, many studies have applied approximation methods for solving various types of
58
integro-differential equations (linear/nonlinear) [16-20]. However, each of these numerical
59
approximation techniques has its own operational limitations that strictly narrow its
60
functioning domain. Hence, it is possible that these approximate techniques fail to overcome
61
a specific problem.
62
A few such instances are mentioned in the following. It was reported that the DTM was
63
unable to produce physically reasonable data for the Glauert-jet Problem [11]. Moreover, for
64
some specific parameter values, the HPM and VIM failed to provide accurate results for the
65
motion of a solid particle in a fluid [21, 22]. Meanwhile, most of them are based on classical
66
mathematical tools.
67
Despite there being a wide range of approximate methods for solving ODEs, there is a lack of
68
a proper approach that meets most of the engineering demands having unconventional and
69
nonlinear ODEs. It should be very interesting to solve linear and nonlinear ODEs having
70
arbitrary boundaries and/or initial values.
71
Therefore, when analytical methods are not capable of solving differential equations or other
72
types of equations in a logical given time, approximation methods are considered as the best
73
solver. Among approximation methods, metaheuristic optimization algorithms, devised by
74
observing the phenomena occurring in nature, have demonstrated their capabilities in finding
75
near-optimal solutions to numerical real-valued problems [23-26].
76
Nowadays, applications that use metaheuristic methods for finding approximate solution of
77
ODEs have increased considerably. This includes genetic algorithms (GAs) [27, 28], particle
78
swarm optimization (PSO) [4], genetic programming [29], and others [30, 31], however, their
Ac ce p
te
d
M
an
us
cr
ip t
54
4 Page 4 of 52
approaches are different with each other in terms of applied strategy and base approximate
80
function. For instance, in [4], different strategy named method of bilaterally bounded, has
81
been employed along with the PSO.
82
Recently, the concept of Fourier series expansion has been used as a base approximate
83
function for finding the approximate solution of ODEs [32]. For the sake of simplicity, in
84
their approach, the weight function was set to unit weight function [32]. However, this
85
assumption may not help us in obtaining better results for all types of ODEs.
86
Inspired by [32] and in order to improve the efficiency of the proposed approach, lately,
87
Sadollah et al. [33] successfully applied weighted residual method for approximate solving of
88
10 ODE problems including mechanical vibration problems and more test ODE problems.
89
Also, superiority of using new weight function (instead of using unit weight function) is
90
shown and two metaheuristic algorithms (i.e., the PSO [34] and water cycle algorithm [35])
91
are utilized for optimization phase.
92
In this paper, using the concept of Fourier series as the base approximate function, other
93
optimizers, including the harmony search (HS) [36], PSO [34], and GA [37] are applied for
94
approximate solving of real life ODEs longitudinal heat transfer fins having various profiles
95
(i.e., rectangular, trapezoidal, and concave parabolic profiles) and properties (i.e., 14 ODEs).
96
In addition, least square weight (LSW) function is proposed for solving the ODEs instead of
97
considering unit weight function.
98
The GA [37], known to be the most famous metaheuristic algorithms, has been applied for
99
optimization phase. The reason of choosing PSO is that to apply the same optimizer to the
100
considered engineering heat transfer problems as did others in the literature [4, 32, 33]. In
101
case of HS algorithm, in light of its simple concept and coding, we selected this optimizer for
102
having more comparisons.
Ac ce p
te
d
M
an
us
cr
ip t
79
5 Page 5 of 52
Talking about the considered ODE problems, extended or finned surfaces are extensively
104
utilized in engineering applications where there is a need to increase heat transfer between a
105
hot primary surface and an adjoining coolant. A comprehensive literature review on different
106
facets of extended surface heat transfer theory and experiment over the past several decades
107
has been carried out by Kraus et al. [38].
108
Recently, Torabi et al. [39] investigated the approximate solving of convective-radiative
109
longitudinal fins for different profiles and nonlinearities using the DTM. They compared their
110
obtained approximate results with exact and numerical solutions. The approximated
111
mathematical formulations for different profiles are given in their study [39].
112
Consequently, the nonlinear fin equations have been solved either numerically or using a
113
variety of approximate analytical methods. In this paper, we modeled convective-radiative
114
longitudinal fin problems as an optimization problem and determined their approximate
115
solution using the proposed method.
116
The remainder of this paper is organized as follows: The next section describes problem
117
formulation as a case study and its behavior in form of ODEs accompanied with its detail. In
118
Section 3, suggested approximate approach for tackling considered problems (i.e., ODEs of
119
longitudinal fins) and the applied optimization methods are explained in brief. In addition, the
120
performance criteria for quantitative assessment among other methods are given in Section 3.
121
Section 4 compares statistical optimization results using the suggested approach with the
122
other approximation method. Graphical comparisons between the exact (numerical) and
123
approximate solutions are demonstrated and the obtained approximate mathematical
124
formulations are represented in this section. Finally, conclusions are drawn in Section 5.
Ac ce p
te
d
M
an
us
cr
ip t
103
2. Longitudinal heat transfer fins: A case study
125 126
An array of rectangular fins is widely used in practice to enhance heat dissipation from a hot
127
primary surface. The fin heat transfer model must include simultaneous surface convection 6 Page 6 of 52
and radiation. Fig. 1 shows a longitudinal (straight) fin having rectangular, trapezoidal, and
129
concave parabolic profiles, with length L.
d
M
an
us
cr
ip t
128
te
130
Fig. 1. Schematic view of longitudinal fin having different profiles.
Ac ce p
131 132
Each fin draws heat from its base at temperature Tb, and transfers it by convection to the
133
surroundings at temperature Ta, and by radiation to an effective sink at temperature Ts. The
134
fin thermal conductivity k, the convective heat transfer coefficient (h) and the fin surface
135
emissivity (ε) are assumed to be functions of temperature of the forms as given follows:
136
k = ka ⎡⎣1 + α (T − Ta ) ⎤⎦ ,
137
⎡ T − Ta ⎤ h = hb ⎢ ⎥ , ⎣ Tb − Ta ⎦
(1b)
138
ε = ε s ⎡⎣1 + β (T − Ts ) ⎤⎦ ,
(1c)
(1a)
m
7 Page 7 of 52
where ka is thermal conductivity at the convection sink temperature Ta, hb is the convection
140
heat transfer coefficient corresponding to the temperature difference (Tb-Ta), and εs is the fin
141
surface emissivity at the radiation sink temperature Ts. The constants α and β are measures of
142
variation of thermal conductivity and surface emissivity with temperature, respectively.
143
Based on one-dimensional (axial) heat conduction in the fin, the general energy balance for
144
longitudinal fins per unit width may be written as [39]:
145
ka
146
⎛ x ⎞ Note that the local semi-fin thicknesses is given by t ( x ) = tb + δ ⎜⎜ ⎛⎜ ⎞⎟ − 1⎟⎟ , where the exponent L ⎝ ⎠
147
n = 0, 1, 2 gives the local semi-fin thicknesses for rectangular, trapezoidal, and concave
148
parabolic fins, respectively. The quantity tb is the semi-base thickness, while the quantity δ
149
defines the fin taper (see Fig. 1). The boundary conditions of constant base temperature and
150
adiabatic tip are given as follows:
151
T ( L) = Tb ,
152
dT dx
cr
m
(T − Ta ) − σε s ⎡⎣1 + β (T − Ts )⎤⎦ (T 4 − Ts4 ) = 0 .
(2)
us
⎡ T − Ta ⎤ d ⎡ dT ⎤ ⎡1 + α (T − Ta ) ⎤⎦ t ( x) − hb ⎢ ⎥ dx ⎢⎣ ⎣ dx ⎥⎦ ⎣ Tb − Ta ⎦
ip t
139
n
an
d
M
⎠
(3a)
te
=0.
Ac ce p
153
⎝
(3b)
x =0
Let’s assume the following dimensionless parameters: θ=
154
T Tb
Nr =
θa =
σε s L2Tb3 ka tb
Ta Tb
θs =
ψ =
Ts Tb
X =
x L
C=
δ tb
Nc =
hb L2Tbm ka tb (Tb − Ta ) m
A = α Tb
B = β Tb
tb L
.
155
Therefore, the problem formulation for all three fin profiles considered has the following
156
ODE form [39]:
157
d dX
158
with the following boundary conditions:
⎡ n ⎡ ⎤ dθ ⎢ ⎡⎣1 + A (θ − θ a ) ⎤⎦ ⎣1 + C X − 1 ⎦ dX ⎣
(
)
m +1 ⎤ 4 4 ⎥ − Nc (θ − θ a ) − Nr ⎡⎣1 + B (θ − θ s ) ⎤⎦ θ − θ s = 0 , ⎦
(
)
(4)
(5)
8 Page 8 of 52
159
θ (1) = 1 ,
(6a)
160
dθ dX
(6b)
= 0. X =0
Eq. (5) shows that thermal performance (i.e., temperature distribution, fin heat transfer rate,
162
fin efficiency, and fin effectiveness) of a longitudinal fin depends on eight parameters,
163
namely, thermal conductivity parameter A, emissivity parameter B, the exponent m associated
164
with convective heat transfer coefficient, convection-conduction parameter Nc, radiation-
165
conduction parameter Nr, the two temperature ratios, θa and θs, that characterize the
166
temperatures of convection and radiation sinks, and the fin taper ratio, C.
167
In addition, we considered a longitudinal fin of constant cross-sectional area Ac, length L, and
168
perimeter of the cross-section P. Similarly, using Eq. (1), and introducing the following
169
dimensionless parameters:
170
θ=
171
Eq. (5) transforms to the following ODE given as follows:
172
d dX
173
having different boundary conditions:
174
θ (0) = 1 ,
(9a)
175
dθ dX
(9b)
M
Ta T h L2 P σε L2T 3 P x , , θ s = s , X = , Nc = b , A = α Tb , B = β Tb , Nr = s b Tb Tb L ka Ac ka Ac
(7)
d
θa =
m
⎡θ − θa ⎤ ⎤ 4 4 ⎥ − Nc (θ − θ a ) ⎢ 1 − θ ⎥ − Nr ⎡⎣1 + B (θ − θ s ) ⎤⎦ θ − θ s = 0 ⎦ a ⎦ ⎣
(
Ac ce p
dθ ⎡ ⎢ ⎡⎣1 + A (θ − θ a ) ⎤⎦ dX ⎣
te
T Tb
an
us
cr
ip t
161
=0 .
)
,
(8)
X =1
176
In this paper, using the proposed approximation method, Eqs. (5) and (8) along with their
177
boundary conditions (i.e., Eqs. (6) and (9)) have been approximately solved and have been
178
compared with their numerical and exact solutions.
179 180 9 Page 9 of 52
181
2.1.
Exact analytical solution
In order to validate the obtained results with the exact solutions, we can consider the simple
183
case of pure convection fins (Nr = 0) of constant thermal conductivity of fin material (A = 0)
184
operating in a constant heat transfer coefficient (m = 0) environment used in Eq. (5). For
185
these fins, the exact analytical solutions with the considered boundary conditions (see Eq.
186
(6)) can be derived as follows [39]:
187
(a) For rectangular profile:
188
θ ( X ) = C1 e
189
(b) For trapezoidal profile:
190
θ ( X ) = C1 J 0 ⎜ 2 − Nc ⎜⎜
191
(c) For concave parabolic profile:
192
θ ( X ) = C1 P⎛ 1
193
where J0 are Y0 are the zero-order Bessel functions of the first and second kinds, respectively.
194
The symbols P and Q represent the Legendre polynomials of the first and second kinds,
195
respectively. The constants C1 and C2 can be obtained by applying the boundary conditions in
196
Eqs. (6a) and (6b) [39].
cr us
+ θa ,
(10b)
⎛ ⎞ CX ⎜ ⎟ +θ , ⎜ C ( C − 1) ⎟ a ⎝ ⎠
(10c)
⎛ ⎛ 1 + C ( X − 1) ⎞ ⎞ ⎛ 1 + C ( X − 1) ⎞ ⎞ ⎟⎟ ⎟ + C2 Y0 ⎜ 2 − Nc ⎜⎜ ⎟⎟ ⎟ + θ a , 2 ⎜ C C2 ⎝ ⎠ ⎟⎠ ⎝ ⎠ ⎟⎠ ⎝
⎛ ⎞ CX ⎜ ⎟ + C2 Q ⎛1 ⎜ C ( C − 1) ⎟ ⎜⎜ 2 ⎝ ⎝ ⎠
C + 4 Nc − C ⎞ ⎟⎟ C ⎠
Ac ce p
te
C + 4 Nc − C ⎞ ⎜⎜ 2 ⎟⎟ C ⎝ ⎠
(10a)
an
⎜ ⎝
Nc X
M
⎛
+ C2 e −
d
Nc X
ip t
182
197
2.2.
Numerical solution
198
Eqs. (5) and (8) along with the considered boundary conditions given in Eqs. (6) and (9) were
199
solved numerically using a finite difference technique in conjunction with Richardson
200
extrapolation using the Maple14.0. The iterative process was continued until a convergence
201
criterion of 10-6 was satisfied [39].
202 203 10 Page 10 of 52
3. Proposed approximate approach for solving engineering ODEs
204 205
In mathematics, an ODE is an equality involving a function and its derivatives. An ODE of
206
order n is an equation having the following form:
207
F ( x, y, y′,..., y ( n ) ) = 0 ,
208
where y is a function of x, y′ = dy / dx is the first derivative with respect to x, and
209
y ( n ) = d ( n ) y / dx ( n ) is the nth derivative with respect to x. Nonhomogeneous ordinary
210
differential equations can be solved if the general solution to the homogenous version is
211
known [40].
212
3.1.
us
cr
ip t
(11)
an
Approximate base function
Fourier series is an expansion of a periodic function f(x) in terms of an infinite sum of sine
214
and cosine terms. Fourier series employs the orthogonal relationships of the sine and cosine
215
functions. The computation and study of Fourier series is known as harmonic analysis. This is
216
extremely useful to decompose an arbitrary periodic function into a set of simple terms.
217
Using some basic concepts of mathematics, accompanied with Fourier expansion, and
218
metaheuristic methods, it is straightforward to solve different types of ODEs having different
219
nature (linear/nonlinear) [32].
220
In traditional method for solving DEs, most of the methods are invented to handle the first or
221
second order ODEs, initial and boundary value problems. However, using the proposed
222
approximate method, there are no such limitations. In fact, using the suggested approximate
223
method in combination with metaheuristic optimizers, there is always an acceptable solution
224
for any type of ODEs for higher orders in their implicit forms.
225
Generally, using Fourier series, one can estimate any periodic function having finite terms. In
226
this paper, we used Fourier series as our base approximate function. There are three main
227
reasons for this selection; 1) A powerful theory backs its convergence for a wide variety of
Ac ce p
te
d
M
213
11 Page 11 of 52
continuous functions [41], 2) It contains sine and cosine terms, which are differentiable up to
229
any order, and 3) A unique form of the approximation function can be utilized for different
230
kinds of ODEs.
231
To clarify further, consider the implicit form of a nonlinear ODE given in Eq. (11), which has
232
to be solved for the interval span x0 to xn, having boundary and initial conditions as defined in
233
Eqs. (12), respectively:
234
y ( x0 ) = y0 ,
y′( x0 ) = y0' .... y ( xn ) = yn ,
235
y ( x0 ) = y0 ,
y′( x0 ) = y0' .... y (n −1) ( x0 ) = y0( n −1) .
236
In general, the suggested approximate base function, which is the partial sum of the Fourier
237
series (finite simple terms of sine and cosine functions) with center x0, as follows [41]:
cr
(12a)
us
y′( xn ) = yn' .... ,
an
(12b)
NT ⎡ ⎛ mπ ( x − x0 ) ⎞ ⎛ mπ ( x − x0 ) ⎞ ⎤ y ( x) ≈ Yappx ( x) = a 0 + ∑ ⎢ am cos ⎜ ⎟ + bm sin ⎜ ⎟⎥ . L L ⎝ ⎠ ⎝ ⎠⎦ m =1 ⎣
(13)
M
238
Accordingly, the derivatives of Eq. (13) are given as follows:
d
239
ip t
228
240
#
Ac ce p
te
NT ⎡ mπ ⎛ mπ ( x − x0 ) ⎞ mπ ⎛ mπ ( x − x0 ) ⎞ ⎤ ' ( x) = ∑ ⎢ − y′( x) ≈ Yappx am sin ⎜ bm cos ⎜ ⎟+ ⎟⎥ L L L ⎝ ⎠ L ⎝ ⎠⎦ m =1 ⎣ NT ⎡ mπ 2 ⎛ mπ ( x − x0 ) ⎞ mπ 2 ⎛ mπ ( x − x0 ) ⎞ ⎤ '' ( x ) = ∑ ⎢ −( ) am cos ⎜ ) bm sin ⎜ y′′( x) ≈ Yappx ⎟−( ⎟⎥ L L L L ⎝ ⎠ ⎝ ⎠⎦ m =1 ⎣
NT ⎡ mπ n ⎛ mπ ( x − x0 ) nπ (n) ( x) = ∑ ⎢( ) am cos ⎜ + y ( n ) ( x) ≈ Yappx L 2 L ⎝ m =1 ⎣
⎞ mπ n ⎛ mπ ( x − x0 ) nπ + ⎟ + ( L ) bm sin ⎜ L 2 ⎠ ⎝
. (14) ⎞⎤ ⎟⎥ ⎠⎦
241
In order to create a weighted residual function (R(x)), Eqs. (13) and (14) are replaced in Eq.
242
(11). In these equations (Eq. (12) to Eq. (14)), xo and xn are the lower and upper bounds of the
243
interval solution; L is the length of the interval solution, which is x0 - xn; NT is the number of
244
approximation terms of the sine and cosine functions.
245
Depending on the complexity of the ODEs, different values of NT may be used. It is worth
246
pointing out that, by increasing the number of approximation terms (NT), the accuracy of the
12 Page 12 of 52
approximate solutions increases. However, in practice, using a large number, as NT increases
248
the computational time and effort are increased, which is usually unnecessary.
249
The unknown coefficients in Eq. (13) (i.e., a0, am, and bm), which are the coefficient of
250
Fourier series, can be computed using metaheuristic optimization methods. The coefficients
251
of Fourier series are considered as design variables of our optimization model for the ODEs.
252
Weighted residual function
cr
3.2.
ip t
247
An optimization model requires an objective function to estimate the cost function (error) of
254
the ODEs. In order to obtain a numerical measure of error, the weighted residual function
255
(WRF) is considered as a cost function that needs to be minimized. In fact, the WRF assumes
256
that a solution can be approximated analytically or piecewise analytically. The general
257
mathematical form of the WRF is given as follows [42]:
258
WRF = ∫ W ( x) × R(x) dx ,
M
an
us
253
D
(15)
where W(x) and R(x) are called the weight and residual functions, respectively. R(x) is
260
obtained by substituting the approximate function Y(x) and its derivatives for y(x), y`(x),…,
261
and y(n)(x) in the implicit form of the Eq. (11) as follows:
262
R( x) = f ( x, Y ( x), Y ′( x),..., Y ( n ) ( x)) .
263
The optimum value of the cost function (i.e., the WRF given in Eq. (15)) is zero. In fact, the
264
notion in the WRF is to force the residual to zero in some average sense over the domain.
265
When the value of cost function approaches zero, better approximate solutions and more
266
accuracy can be obtained.
267
For solving the cost function given in Eq. (15), any numerical method such as the
268
Trapezoidal or Simpson methods can be used. In this paper, the trapezoidal method is utilized
269
for numerically solving the cost function. In the following section, new suggested approach
270
for weight function is introduced.
Ac ce p
te
d
259
(16)
13 Page 13 of 52
3.2.1. Least square weight function
271
If the continuous summation of all the squared residuals is minimized, the rationale behind
273
the name can be seen. In other words, a minimum of the following equation needs to be
274
computed:
275
WRM = ∫ R ( x) × R ( x)dx = ∫ R 2 ( x)dx , D
ip t
272
(17)
D
In order to achieve a minimum of this scalar function, the derivatives of the WRF with
277
respect to all the unknown parameters must be zero. That is:
∂R dx = 2 ∫ R ( x) × ∂ai D
us
278
.
an
∂WRM =0 ∂ai
cr
276
Therefore, the weight functions are seen to be:
280
Wi (x) = 2
281
However, the “2” in Eq. (19) can be dropped, because it cancels out in the equation.
282
Therefore, the weight functions for the least squares method are the derivatives of the
283
residual with respect to the unknown constants:
284
Wi (x) =
285
3.3.
M
279
(19)
Ac ce p
te
d
∂R . ∂ai
(18)
∂R . ∂ai
(20)
Initial and boundary conditions
286
The ODE problems in which values are specified at a particular value of y(x0) are called
287
initial value problems (IVPs). There is another class of problems called boundary value
288
problems (BVPs) in which conditions are given at both ends rather than merely at the initial
289
point (x0).
14 Page 14 of 52
It is crucial for the approximate solution to satisfy the conditions of the IVPs and BVPs to
291
fulfill the equation, simultaneously. The constraints for the ODE problems can be in the form
292
of homogeneous or nonhomogeneous conditions. They are given as follows:
293
For homogeneous conditions:
ip t
290
y ( x0 ) = 0 ⇒ g1 ( x0 ) = Yappx ( x0 ) − 0 .
#
cr
294
' y ' ( x0 ) = 0 ⇒ g 2 ( x0 ) = Yappx ( x0 ) − 0
295
us
(n) y ( n ) ( x0 ) = 0 ⇒ g n +1 ( x0 ) = Yappx ( x0 ) − 0
For nonhomogeneous conditions:
' y ' ( x0 ) = y0' ⇒ g 2 ( x0 ) = Yappx ( x0 ) − y0'
an
y ( x0 ) = y0 ⇒ g1 ( x0 ) = Yappx ( x0 ) − y0 296
(21)
.
#
M
(n) y ( n ) ( x0 ) = y0( n ) ⇒ g n +1 ( x0 ) = Yappx ( x0 ) − y0( n )
(22)
Eqs. (21) and (22) represent the constraints of our optimization model. In this paper, all
298
equality constraints were converted into inequality ones, |h(x)| - ε ≤ 0 using the degree of
299
violation ε = 2.2e−16 that was taken from MATLAB. As a constraint handling approach, the
300
absolute violations caused by the boundary or initial conditions (homogeneous or
301
nonhomogeneous) are calculated. Afterwards, using the penalty function approach, the
302
absolute values of imposed violations as penalty term are added to the objective function.
Ac ce p
te
d
297
303
3.4.
Performance metrics
304
To perform a fair quantitative evaluation and judgment between the approximate and exact
305
(numerical) solutions, we investigated the generational distance (GD) and inverted
306
generational distance (IGD) metrics, which are considered as assessment criteria widely used
307
in the literature [43, 44]. They were first presented by Veldhuizen and Lamont and Zitzler et
308
al. [45, 46]. 15 Page 15 of 52
The main objective of these criteria is to clarify the capability of the WRF and optimizers in
310
finding an approximate solution having the least distance when compared with the exact
311
(numerical) solution. Based on this definition, it is evident that the obtained approximate
312
solution with the minimum GD and IGD has the best convergence to the exact (numerical)
313
solution. These evaluation factors are defined in mathematical forms given as follows [43]:
ip t
309
,
us
314
cr
1/2
⎛ n 2⎞ ⎜ ∑ di ⎟ GD = ⎝ i =1 ⎠ n n
∑d
i
(23)
315
IGD =
316
where n is the number of points in the approximate solution, and d is the Euclidean distance
317
between point i in approximate solution and the nearest point (minimum Euclidean distance)
318
in the exact (numerical) solution. Lower values of GD and IGD imply that obtained
319
approximate solution must be very close to the optimal Pareto front and cannot miss any part
320
of the whole Pareto front. In this paper, for comparison purposes, the chosen value for n is set
321
to 100 for the ODEs of longitudinal fins. This means that 100 candidate points in the
322
approximate and exact (numerical) solutions are compared with each other in this paper. Fig.
323
2 shows a schematic view of these performance metrics for evaluating the accuracy between
324
the exact (numerical) and approximate solutions.
,
(24)
Ac ce p
te
d
M
n
an
i =1
16 Page 16 of 52
ip t cr
325
Fig. 2. Schematic view of the GD (IGD) criterion for evaluating the accuracy of the obtained
327
approximate solution.
an
3.5.
328
us
326
Metaheuristic algorithms
In this section, different optimizers are explained in brief. For solving the ODEs using the
330
approximate WRF method, metaheuristic algorithms can be considered as alternative
331
approaches for finding unknown coefficients for the Fourier series. In this paper,
332
metaheuristic algorithms including PSO, GA, and HS are used to find an approximate
333
solution for the considered ODEs.
te
d
M
329
3.5.1. Particle swarm optimization
Ac ce p
334 335
Particle swarm optimization (PSO) is an evolutionary computation technique for solving the
336
global optimization problems developed by Kennedy and Eberhart [34]. This computation
337
technique functions through individual improvement, population cooperation, and
338
competition. The PSO is based on the simulation of simplified social models, such as bird
339
flocking, fish schooling, and the swarm theory.
340
Researchers found that the synchrony of an animal’s behavior was through maintaining
341
optimal distances between individual members and its neighbors. The PSO exhibits common
342
evolutionary computation attributes including initialization with a population of random
343
solutions and searching for optima by updating generations [47]. 17 Page 17 of 52
Potential solutions, called particles, are then flown through the problem space by following
345
the current optimum particles. Each particle keeps track of the coordinates in the problem
346
space that are associated with the best solution (fitness) it has achieved so far. This value is
347
called pBest. Another best value, which is tracked by the global version of the PSO, is the
348
overall best value and this location obtained is the best so far by any particle in the
349
population. This location is called gBest.
350
The PSO’s concept consists of changing the velocity of each particle toward its pBest and
351
gBest locations. Acceleration is weighted by a random term with separate random numbers
352
being generated for acceleration towards pBest and gBest locations. The basic swarm
353
parameters position and velocity are updated using the following equations [34]:
354
G G G G G G Vi +1 = wVi + c1r1 ( pBesti − X i ) + c2 r2 ( gBesti − X i ) ,
355
G G G X i +1 = X i + Vi +1 ,
356
where w is the inertia weight for previous velocity (between 0 and 1), Xi is the current
357
position for particle i, Vi+1 is the updated velocity of particle i, pBesti is the best solution
358
found by particle i, gBesti is the best solution found by the swarm, r1 and r2 are uniform
359
random numbers in the range of [0, 1], c1 is the cognitive component, and c2 is the social
360
component.
361
Here, c1 and c2 are the constants that influence how each particle is directed towards good
362
positions by considering their personal best and global best information, respectively. The
363
role of w is crucial for the convergence of the PSO. It is applied to control the effect of the
364
previous velocity on the current particle velocity [47].
365
Thus, the user-defined parameter w adjusts the trade-off between the global and the local
366
searches in the PSO. A high value of w performs more exploration, while a low value
367
increases the exploitation effects. A suitable value for the w provides balance between the
(25) (26)
Ac ce p
te
d
M
an
us
cr
ip t
344
18 Page 18 of 52
368
global and local phases, and, thus results in finding better solutions. Based on the practical
369
experiments, it is more desirable to primarily set the w to a large value for having more
370
exploration at early iterations, and gradually decrease it to obtain refined solutions.
3.5.2. Harmony search algorithm
ip t
371
Since the harmony search (HS) algorithm was first developed in 2001 [36], it has been
373
applied to various research areas and obtained considerable attention in different disciplines
374
[48]. The HS conceptualizes a musical process of searching for a perfect state of harmony.
375
Musical performances seek a fantastic harmony determined by aesthetic estimation, as the
376
optimization techniques seek a best state (global optimum) determined by objective function
377
value.
378
Aesthetic estimation is performed by the set of the sounds played by musical ensemble, as
379
objective function value is evaluated by the set of the values produced by adjusted variables;
380
the better aesthetic sounds can be improved by constant practice, as the minimization
381
/maximization of the objective function can mostly be improved by repeating iteration.
382
More information about the HS can be found in the work of Geem et al. [36] and Kim et al.
383
[49]. The main steps of the original HS are summarized as below:
384
Step 1: Generate random vectors (x1, x2, …, xHMS) equal to the harmony memory size (HMS).
385
Then, store them in harmony memory (HM).
386
⎡ x11 " x1n f ( x1 ) ⎤ ⎢ ⎥ # ⎥ . HM = ⎢ # % # ⎢ x hms " x hms f ( x hms ) ⎥ n ⎣ 1 ⎦
387
The usage of HM is vital, as it is similar to the choice of the best-fit individuals in the GA.
388
This will ensure that the best harmonies will be transferred to the new candidate solution. In
389
order to use this memory more excellently, it is typically assigned as a user-defined
390
parameter, so called harmony memory considering rate (HMCR, 0 ≤ HMCR ≤ 1). If this value
Ac ce p
te
d
M
an
us
cr
372
(27)
19 Page 19 of 52
is too low, only few best harmonies are selected and it may converge too slowly. If this rate is
392
extremely high (near 1), almost all the harmonies are used in the HM, then other harmonies
393
are not explored well, resulting immature convergence. Therefore, typically, proper value of
394
HMCR is set from 0.7 to 0.98 [36].
395
Step 2: Generate a new harmony. In generating a new harmony x ′ , for each component xi′ :
ip t
391
•
With probability HMCR, pick the stored value from the HM: xi′ ← xiint ( u (0,1)× HMS )+1
397
•
With probability (1-HMCR), pick a random value within the allowed range.
us
Step 3: Perform additional task if the value in Step 2 came from the HM. •
399
With probability PAR (pitch adjusting rate; 0 ≤ PAR ≤ 1), change xi′ by a small
an
398
cr
396
amount: xi′ ← xi′ + bw × Rand ( − 1,1) for continuous variable. bw is the amount of
401
maximum change in pitch adjustment. •
With probability (1-PAR), do nothing.
d
402
M
400
The second operator is the PAR determined by a pitch band-width (bw). Indeed, this operator
404
corresponds to generate a slightly different solution in the HS [36]. This step is similar to the
405
mutation operator in the GA. A low value of PAR with a narrow bandwidth can slow down
406
the convergence of HS. On the other hand, a very high value of PAR with a wide bandwidth
407
may cause the solution to scatter around some potential optima performing as a random
408
search. Recommend value of PAR can be chosen in the range [0.1-0.5] [36].
409
Step 4: Calculate the cost function of new generated harmony and update the HM.
410
Step 5: Repeat Step 2 to Step 5 until the termination criterion (e.g., maximum number of
411
function evaluation) is satisfied.
412
3.5.3. Genetic algorithm
Ac ce p
te
403
413
Genetic algorithm (GA) is a member of a collection of methodologies known as evolutionary
414
algorithm. The GA inspired by the principals of natural selection and evolution processes 20 Page 20 of 52
[37] and it is identified as robust metaheuristic tools capable of delivering efficient solutions
416
to diverse design problems. In general, GA begins its search with a population of random
417
individuals.
418
Each member of the population possesses a chromosome (individual) which is comprised of
419
genes (design variables). In original GA, a gene may take one of two allele values, either 1 or
420
0. The first step in creating an offspring population for GA is to construct a mating pool.
421
The mating pool contains N individuals which are copied from the parent population that will
422
be utilized to create the offspring population. To create the mating pool, different selection
423
operators (e.g., tournament and roulette wheel) may be utilized.
424
To begin, two solutions from the parent population are selected to compete in a tournament in
425
order to construct the mating pool. A copy of the winner of the tournament is kept in the
426
mating pool. This process is repeated until each solution in the parent population has
427
competed twice and the N spots in the mating pool have been filled. The resulting mating
428
pool contains more copies of the more desirable solutions in the parent population and fewer
429
copies of the less desirable solutions from the parent population.
430
By completion of mating pool, the task of constructing the offspring population commences.
431
To generate two offspring solutions, one begins by selecting two individuals from the mating
432
pool at random. Once selected, the two individuals undergo crossover to create two offspring
433
solutions which are placed in the offspring population.
434
For exploration purpose, mutation operator applies on a random selected individual from the
435
mating pool [37]. There are many diverse applications of GA in the literature [27].
436
Particularly, the GA applied for solving ODEs combined with the Nelder-Mead method [28].
437
However, the approach of solving ODEs using the GA in [28] differs with the proposed
438
approach in this paper. In this paper, using the concept of weighted residual method along
Ac ce p
te
d
M
an
us
cr
ip t
415
21 Page 21 of 52
439
with Fourier series, optimization methods have been employed for approximately solving
440
ODEs. 3.6.
441
Steps of proposed approximate method
In summary, in light of all the concepts mentioned in Sections 3.1 to 3.5, the objective
443
function of the ODEs (minimization of the WRF) and its constraints are as follows:
444
Minimize
ip t
442
WRF = ∫ W ( x) × R ( x) dx ,
cr
(28)
D
subject to the BVP and/or IVP constraints:
446
⎧0 giIVP ( x0 ) = ⎨ ⎩ y0
447
⎧0 g BVP ( x0 ) = ⎨ j ⎩ y0
448
where W(x) and R(x) are the weight and residual functions, respectively. Accordingly, Eqs.
449
(29a) and (29b) are considered as the constraints for the initial and boundary conditions.
450
Steps of the approximate method are given as follows:
451
Step 1: Convert the given ODE problem in its implicit form based on Eq. (11).
452
Step 2: Prepare the constraints (initial and boundary conditions) in the proper format given in
453
Eqs. (12a) and (12b).
454
Step 3: Choose an appropriate number of terms for the Fourier series (NT), as in Eq. (13).
455
Step 4: Initiate user parameters for the applied metaheuristic algorithm (e.g., PSO and GA).
456
Step 5: Launch the selected optimizer in order to find the best coefficients for the
457
approximate function (Fourier series). Cost function and constraints are embedded in the
458
optimizer for optimization purposes.
459
Step 6: Report the obtained coefficients of the Fourier series.
i = 1, 2,..., n IVPs
j = 1, 2,..., nBVPs
,
(29a)
(29b)
Ac ce p
te
d
M
j = 1, 2,..., n BVPs
,
an
i = 1, 2,..., nIVPs
us
445
22 Page 22 of 52
Step 7: Calculate the GD and IGD performance metrics to assess the goodness of the
461
approximate solution against the exact (numerical) solution.
462
Step 8: Plot the obtained approximate and exact (numerical) solutions in the solution interval
463
from x0 to xn.
ip t
460
4. Numerical results and discussion
464
In the following sections, approximate solutions of longitudinal fins with various profiles and
466
coefficients are investigated using the DTM and proposed WRF technique, respectively. It is
467
worth mentioning that the WRF can be linked with any optimization method. Particularly, in
468
this paper, the PSO, GA, and HS are used for optimization purposes (see Section 4.2). 4.1.
469
an
us
cr
465
Approximate solutions obtained by the DTM
The ODEs given in Eq. (5) along with its boundary conditions (see Eq. (6)) have been solved
471
analytically using the DTM [39]. In the literature, many studies have been utilized and
472
described the DTM in detail [10, 11]. By considering Eq. (5) and fundamental operations
473
used in the DTM, the following are the approximate mathematical formulations obtained by
474
the DTM as follows [39]:
te
d
M
470
(a) For rectangular profile for m = 0:
Ac ce p
475
⎛
⎞
k
(1− Aθa )( k + 2)( k +1) Θ( k + 2) + A ⎜ ∑Θ(v )( k + 2 −v )( k +1−v ) Θ( k + 2 −v ) ⎟ ⎝ v =0
⎠
⎛ ⎞ +A ⎜ ∑(v +1) Θ(v +1)( k +1−v ) Θ( k +1−v ) ⎟ v = 0 ⎝ ⎠ k
476
⎛ k ⎛m ⎛v ⎛w ⎞⎞⎞⎞ −NrB ⎜ ∑ Θ( k − m ) ⎜ ∑Θ( m −v ) ⎜⎜ ∑ Θ(v −w ) ⎜ ∑Θ(w −u ) Θ(u ) ⎟ ⎟⎟ ⎟ ⎟ ⎜ v =0 ⎜ m =0 ⎝ u =0 ⎠ ⎠ ⎟⎠ ⎟⎠ ⎝ w =0 ⎝ ⎝ ⎛ k ⎛m ⎛v ⎞⎞⎞ + ( NrBθs − Nr ) ⎜ ∑ Θ( k − m ) ⎜⎜ ∑Θ( m −v ) ⎜ ∑ Θ(v −w ) Θ(w ) ⎟ ⎟⎟ ⎟ ⎜ m =0 ⎝ w =0 ⎠ ⎠ ⎟⎠ ⎝ v =0 ⎝
(
)
(
.
(30a)
)
+ NrBθs4 − Nc Θ( k ) + Ncθa + Nrθs4 − NrBθs5 δ ( k ) = 0
(b) For rectangular profile for m = 2:
477
23 Page 23 of 52
⎛
⎞
k
(1− Aθa )( k + 2)( k +1) Θ( k + 2) + A ⎜ ∑Θ(v )( k + 2 −v )( k +1−v ) Θ( k + 2 −v ) ⎟ ⎝ v =0
(
)
(c) For trapezoidal profile for m = 0:
479
⎛
k
⎝
v =0
ip t
)
. (30b)
cr
(
us
478
⎠
⎛k ⎞ +A ⎜ ∑(v +1) Θ(v +1)( k +1−v ) Θ( k +1−v ) ⎟ ⎝ v =0 ⎠ ⎛ k ⎛m ⎛ v ⎛w ⎞⎞⎞⎞ −NrB ⎜ ∑ Θ( k − m ) ⎜ ∑Θ( m −v ) ⎜⎜ ∑ Θ(v −w ) ⎜ ∑Θ(w −u ) Θ(u ) ⎟ ⎟⎟ ⎟ ⎟ ⎜ v =0 ⎜ m =0 ⎝ u =0 ⎠ ⎠ ⎟⎠ ⎟⎠ ⎝ w =0 ⎝ ⎝ ⎛ k ⎛m ⎛k ⎛ v ⎞⎞⎞ ⎛v ⎞⎞ + ( NrBθs − Nr ) ⎜ ∑ Θ( k − m ) ⎜⎜ ∑Θ( m −v ) ⎜ ∑ Θ(v −w ) Θ(w ) ⎟ ⎟⎟ ⎟ − Nc ⎜⎜ ∑Θ( k −v ) ⎜ ∑ Θ(v −w ) Θ(w ) ⎟ ⎟⎟ ⎜ m =0 ⎟ ⎝ w =0 ⎠⎠⎠ ⎝ w =0 ⎠⎠ ⎝ v =0 ⎝ v =0 ⎝ k ⎛ ⎞ +3Ncθa ⎜ ∑ Θ( k − m ) Θ( m ) ⎟ + NrBθs4 − 3Ncθa2 Θ( k ) + Ncθa3 + Nrθs4 − NrBθs5 δ ( k ) = 0 = 0 m ⎝ ⎠
⎞
(1− Aθa −C +CAθa )( k + 2)( k +1) Θ( k + 2) + ( A −CA ) ⎜ ∑Θ(v )( k + 2 −v )( k +1−v ) Θ( k + 2 −v ) ⎟ ⎠
an
⎛k ⎞ + ( A −CA ) ⎜ ∑(v +1) Θ(v +1)( k +1−v ) Θ ( k +1−v ) ⎟ + (C −CAθa )( k +1) Θ( k +1) ⎝ v =0 ⎠
M
⎛ k ⎛m ⎞⎞ +AC ⎜ ∑δ ( k −1− m ) ⎜ ∑Θ(v )( m + 2 −v )( m +1−v ) Θ( m + 2 −v ) ⎟ ⎟ ⎝ v =0 ⎠⎠ ⎝ m =0
d
480
⎛k ⎞ ⎛k ⎞ +CA ⎜ ∑Θ(v )( k +1−v ) Θ( k +1−v ) ⎟ + (C −CAθa ) ⎜ ∑( l + 2)( l + 1) Θ( l + 2) δ ( k −1− l ) ⎟ ⎝ v =0 ⎠ ⎝ l =0 ⎠ k m ⎛ ⎛ ⎞⎞ +AC ⎜ ∑δ ( k −1− m ) ⎜ ∑(v +1) Θ(v +1)( m +1−v ) Θ( m +1−v ) ⎟ ⎟ ⎝ v =0 ⎠⎠ ⎝ m =0
.
(30c)
te
⎛ k ⎛m ⎛v ⎛w ⎞⎞⎞⎞ −NrB ⎜ ∑ Θ( k − m ) ⎜ ∑Θ( m −v ) ⎜⎜ ∑ Θ(v −w ) ⎜ ∑Θ(w −u ) Θ(u ) ⎟ ⎟⎟ ⎟ ⎟ ⎜ v =0 ⎜ m =0 ⎝ u =0 ⎠ ⎠ ⎟⎠ ⎟⎠ ⎝ w =0 ⎝ ⎝
(
Ac ce p
⎛ k ⎛m ⎛ v ⎞⎞⎞ + ( NrBθs − Nr ) ⎜ ∑ Θ( k − m ) ⎜⎜ ∑Θ( m −v ) ⎜ ∑ Θ(v −w ) Θ(w ) ⎟ ⎟⎟ ⎟ ⎜ m =0 ⎝ w =0 ⎠ ⎠ ⎟⎠ ⎝ v =0 ⎝
)
(
)
+ NrBθs4 − Nc Θ( k ) + Ncθa + Nrθs4 − NrBθs5 δ ( k ) = 0
(d) For trapezoidal profile for m = 2:
481
24 Page 24 of 52
⎛
⎞ Θ (k + 2 − v)⎟ ⎠
k
(1 − Aθ a − C + CAθ a )( k + 2 )( k + 1) Θ ( k + 2 ) + ( A − CA) ⎜ ∑ Θ ( v )( k + 2 − v )( k + 1 − v ) ⎝ v =0
⎛ k ⎞ + ( A − CA ) ⎜ ∑ ( v + 1) Θ ( v + 1)( k + 1 − v ) Θ ( k + 1 − v ) ⎟ + ( C − CAθ a )( k + 1) Θ ( k + 1) ⎝ v =0 ⎠ ⎛ k ⎞ ⎛ k ⎞ +CA ⎜ ∑ Θ ( v )( k + 1 − v ) Θ ( k + 1 − v ) ⎟ + ( C − CAθ a ) ⎜ ∑ ( l + 2 )( l + 1) Θ ( l + 2 ) δ ( k − 1 − l ) ⎟ ⎝ v =0 ⎠ ⎝ l =0 ⎠
ip t
⎛ k ⎛ m ⎞⎞ + AC ⎜ ∑ δ ( k − 1 − m ) ⎜ ∑ ( v + 1) Θ ( v + 1)( m + 1 − v ) Θ ( m + 1 − v ) ⎟ ⎟ ⎝ v =0 ⎠⎠ ⎝ m=0 ⎛ k ⎛ m ⎞⎞ + AC ⎜ ∑ δ ( k − 1 − m ) ⎜ ∑ Θ ( v )( m + 2 − v )( m + 1 − v ) Θ ( m + 2 − v ) ⎟ ⎟ ⎝ v =0 ⎠⎠ ⎝ m=0
cr
⎛ k ⎛ m ⎛ v ⎛ w ⎞⎞⎞⎞ − NrB ⎜ ∑ Θ ( k − m ) ⎜ ∑ Θ ( m − v ) ⎜⎜ ∑ Θ ( v − w ) ⎜ ∑ Θ ( w − u ) Θ ( u ) ⎟ ⎟⎟ ⎟ ⎟ ⎜ v=0 ⎜ m=0 ⎝ u =0 ⎠ ⎠ ⎟⎠ ⎟⎠ ⎝ w= 0 ⎝ ⎝
(
) Θ ( k ) + ( Ncθ
(e) For concave parabolic profile for m = 0:
483
⎛
k
3 a
)
+ Nrθ s4 − NrBθ s5 δ ( k ) = 0
.
(30d)
an
482
⎛ k ⎞ +3 Ncθ a ⎜ ∑ Θ ( k − m ) Θ ( m ) ⎟ + NrBθ s4 − 3Ncθ a2 ⎝ m=0 ⎠
us
⎛ k ⎛ k ⎛ m ⎛ v ⎞⎞ ⎛ v ⎞⎞⎞ + ( NrBθ s − Nr ) ⎜ ∑ Θ ( k − m ) ⎜⎜ ∑ Θ ( m − v ) ⎜ ∑ Θ ( v − w ) Θ ( w ) ⎟ ⎟⎟ ⎟ − Nc ⎜⎜ ∑ Θ ( k − v ) ⎜ ∑ Θ ( v − w ) Θ ( w ) ⎟ ⎟⎟ ⎜ m=0 ⎟ v w v w = 0 = 0 = 0 = 0 ⎝ ⎠⎠ ⎝ ⎠⎠⎠ ⎝ ⎝ ⎝
M
(1 − Aθ a − C + CAθ a )( k + 2 )( k + 1) Θ ( k + 2 ) + ( A − CA) ⎜ ∑ Θ ( v )( k + 2 − v )( k + 1 − v )
⎞ Θ (k + 2 − v) ⎟ ⎠
te
d
⎝ v =0 ⎛ k ⎞ + ( A − CA ) ⎜ ∑ ( v + 1) Θ ( v + 1)( k + 1 − v ) Θ ( k + 1 − v ) ⎟ ⎝ v =0 ⎠ k ⎛ k ⎛ ⎞ ⎛ m ⎞⎞ + ( 2C − 2CAθ a ) ⎜ ∑ ( l + 1) Θ ( l + 1) δ ( k − 1 − l ) ⎟ + 2 AC ⎜ ∑ δ ( k − 1 − m ) ⎜ ∑ Θ ( v )( m + 1 − v ) Θ ( m + 1 − v ) ⎟ ⎟ ⎝ l =0 ⎠ ⎝ v =0 ⎠⎠ ⎝ m=0 k ⎛ ⎞ + ( C − CAθ a ) ⎜ ∑ ( l + 2 )( l + 1) Θ ( l + 2 ) δ ( k − 2 − l ) ⎟ ⎝ l =0 ⎠
Ac ce p
⎛ k ⎛ m ⎞⎞ + AC ⎜ ∑ δ ( k − 2 − m ) ⎜ ∑ ( v + 1) Θ ( v + 1)( m + 1 − v ) Θ ( m + 1 − v ) ⎟ ⎟ ⎝ v=0 ⎠⎠ ⎝ m=0 ⎛ k ⎛ m ⎞⎞ + AC ⎜ ∑ δ ( k − 2 − m ) ⎜ ∑ Θ ( v )( m + 2 − v )( m + 1 − v ) Θ ( m + 2 − v ) ⎟ ⎟ ⎝ v=0 ⎠⎠ ⎝ m=0
484
⎛ k ⎛ m ⎛ v ⎛ w ⎞⎞⎞⎞ − NrB ⎜ ∑ Θ ( k − m ) ⎜ ∑ Θ ( m − v ) ⎜⎜ ∑ Θ ( v − w ) ⎜ ∑ Θ ( w − u ) Θ ( u ) ⎟ ⎟⎟ ⎟ ⎟ ⎜ v=0 ⎜ m=0 ⎝ u =0 ⎠ ⎠ ⎟⎠ ⎟⎠ ⎝ w=0 ⎝ ⎝ ⎛ k ⎛ m ⎛ v ⎞⎞⎞ + ( NrBθ s − Nr ) ⎜ ∑ Θ ( k − m ) ⎜⎜ ∑ Θ ( m − v ) ⎜ ∑ Θ ( v − w ) Θ ( w ) ⎟ ⎟⎟ ⎟ ⎜ m=0 = 0 = 0 w v ⎝ ⎠ ⎠ ⎟⎠ ⎝ ⎝
(
)
(
)
+ NrBθ s4 − Nc Θ ( k ) + Ncθ a + Nrθ s4 − NrBθ s5 δ ( k ) = 0
.
(30e)
(f) For concave parabolic profile for m = 2:
485 486
25 Page 25 of 52
⎛
k
(1 − Aθ a − C + CAθ a )( k + 2 )( k + 1) Θ ( k + 2 ) + ( A − CA) ⎜ ∑ Θ ( v )( k + 2 − v )( k + 1 − v ) ⎝ v =0
⎞ Θ (k + 2 − v)⎟ ⎠
⎛ k ⎞ + ( A − CA ) ⎜ ∑ ( v + 1) Θ ( v + 1)( k + 1 − v ) Θ ( k + 1 − v ) ⎟ ⎝ v =0 ⎠ k ⎛ k ⎛ ⎞ ⎛ m ⎞⎞ + ( 2C − 2CAθ a ) ⎜ ∑ ( l + 1) Θ ( l + 1) δ ( k − 1 − l ) ⎟ + 2 AC ⎜ ∑ δ ( k − 1 − m ) ⎜ ∑ Θ ( v )( m + 1 − v ) Θ ( m + 1 − v ) ⎟ ⎟ ⎝ l =0 ⎠ ⎝ v =0 ⎠⎠ ⎝ m=0
ip t
⎛ k ⎞ + ( C − CAθ a ) ⎜ ∑ ( l + 2 )( l + 1) Θ ( l + 2 ) δ ( k − 2 − l ) ⎟ ⎠ ⎝ l =0 ⎛ k ⎛ m ⎞⎞ + AC ⎜ ∑ δ ( k − 2 − m ) ⎜ ∑ ( v + 1) Θ ( v + 1)( m + 1 − v ) Θ ( m + 1 − v ) ⎟ ⎟ ⎝ v =0 ⎠⎠ ⎝ m=0
⎛ k ⎞ +3 Ncθ a ⎜ ∑ Θ ( k − m ) Θ ( m ) ⎟ + NrBθ s4 − 3 Ncθ a2 ⎝ m=0 ⎠
(
Considering θ(0) = a and
dθ dX
) Θ ( k ) + ( Ncθ
)
an
488
⎛ k ⎛ m ⎛ k ⎛ v ⎞⎞⎞ ⎛ v ⎞⎞ + ( NrBθ s − Nr ) ⎜ ∑ Θ ( k − m ) ⎜⎜ ∑ Θ ( m − v ) ⎜ ∑ Θ ( v − w ) Θ ( w ) ⎟ ⎟⎟ ⎟ − Nc ⎜⎜ ∑ Θ ( k − v ) ⎜ ∑ Θ ( v − w ) Θ ( w ) ⎟ ⎟⎟ ⎜ m=0 ⎟ v w v w = 0 = 0 = 0 = 0 ⎝ ⎠⎠⎠ ⎝ ⎠⎠ ⎝ ⎝ ⎝ 3 a
+ Nrθ s4 − NrBθ s5 δ ( k ) = 0
.
(30f)
= 0 for all three fin profiles, and applying the differential
M
487
us
⎛ k ⎛ m ⎛ v ⎛ w ⎞⎞⎞⎞ − NrB ⎜ ∑ Θ ( k − m ) ⎜ ∑ Θ ( m − v ) ⎜⎜ ∑ Θ ( v − w ) ⎜ ∑ Θ ( w − u ) Θ ( u ) ⎟ ⎟⎟ ⎟ ⎟ ⎜ v =0 ⎜ m=0 ⎝ u =0 ⎠ ⎠ ⎟⎠ ⎟⎠ ⎝ w=0 ⎝ ⎝
cr
⎛ k ⎛ m ⎞⎞ + AC ⎜ ∑ δ ( k − 2 − m ) ⎜ ∑ Θ ( v )( m + 2 − v )( m + 1 − v ) Θ ( m + 2 − v ) ⎟ ⎟ ⎝ v =0 ⎠⎠ ⎝ m=0
X =0
transformation, boundary conditions can be stated as follows:
490
Θ(0) = a ,
491
Θ(1) = 0 .
492
Accordingly, from a process of inverse differential transformation, Θ( k + 2) is calculated in
493
Eq. (30). Finally, the following solution by the differential inverse transform of Θ( k ) can be
494
obtained as follows:
495
θ ( X ) = ∑ X k Θ( k ) .
te
d
489
(31a)
Ac ce p
(31b)
∞
(32)
k =0
496
Using Eq. (31a), the value of a can be determined from boundary condition given in Eq. (6a)
497
using the fsolve command in Maple 14 [39].
498
4.2.
Approximate solutions obtained by the proposed approach
499
This section does not aim to investigate the ODEs of longitudinal fins from the analytical
500
point of view. In fact, an applicable approach, with the aid of optimization methods, is 26 Page 26 of 52
applied in this paper for approximate solving of ODEs of convective-radiative longitudinal
502
fins having different profiles (n) and heat transfer coefficients (m).
503
In this paper, three well-known optimization methods, named as the PSO, GA, and HS have
504
been selected to show that the proposed approximate method provides high accuracy
505
approximate solutions. Also, interested readers may apply other optimizers in the literature.
506
Choice of metaheuristic algorithm depends on the user and, therefore, choice of optimization
507
method is optional.
508
The task of optimization related to the optimization part was carried out using 30 independent
509
runs. Regarding the side constraints for unknown coefficients of Fourier series, lower and
510
upper bounds of the design variables are considered between -1 and 1, respectively. Table 1
511
represents the sensitivity analyses of applied optimization methods for the longitudinal fins.
512
Table 1. Sensitivity analyses for the user parameters of considered optimizers (n=0 and m=0).
513
“SD” stands for standard deviation.
d
M
an
us
cr
ip t
501
1 0.0084 0.0465 0.1573 0.0574
Ac ce p
Best WRF Average WRF Worse WRF SD
te
Measures
w = 0.75 c1 and c2 1.2 1.5 0.0084 0.0084 0.0414 0.0090 0.1810 0.0104 0.0590 0.0007
Measures
0.05 0.3068 Best WRF Average WRF 0.7594 Worse WRF 1.2986 0.3399 SD
PSO
2 0.3 0.0084 0.0085 0.0093 0.0791 0.0142 0.2822 0.0018 0.1012 HS
HMCR = 0.98 PAR 0.1 0.3 0.0440 0.0162 0.1170 0.0620 0.2213 0.1700 0.0603 0.0427
c1 and c2 =1.5 w 0.5 0.75 0.0084 0.0084 0.0155 0.0086 0.0422 0.0090 0.0127 0.0001
0.9 0.0087 0.0117 0.0226 0.0043
PAR = 0.3 HMCR 0.97 0.98 0.0266 0.0161 0.0657 0.0624 0.1164 0.1049 0.0328 0.0288
0.5 0.95 0.99 0.0406 0.0175 0.0226 0.1051 0.0631 0.0625 0.1841 0.1306 0.1328 0.0436 0.0368 0.0349 GA Crossover Rate (CR) = 0.8 Mutation Rate (MR) = 0.05 Measures MR CR 0.01 0.05 0.1 0.3 0.5 0.7 0.8 0.9 0.4943 0.3954 0.4134 0.4534 0.4123 0.2466 0.3912 0.2433 Best WRF Average WRF 1.8223 1.0643 1.7454 1.3354 1.8443 2.6972 1.6201 3.7203 27 Page 27 of 52
Worse WRF SD
4.6233 2.4545 5.5043 4.3534 4.9723 7.1447 3.9545 7.8076 1.2145 0.5934 1.6052 1.2023 1.4652 2.0805 1.2121 2.3685
514
Different values have been used for each user parameter for all considered optimizer as can
516
be seen in Table 1. By observing Table 1, the highlighted bold values have provided better
517
statistical results against other selected values. Therefore, the initial parameters for the PSO,
518
based on the results given in Table 1, were as follows; population size of 300, the inertia
519
weight (w) for velocity was 0.75, and the cognitive and social components (i.e., c1 and c2)
520
were 1.5.
521
Accordingly, the user parameters of HS for all considered ODEs were: harmony memory size
522
of 30, the chosen values for the HMCR, PAR, and bw were set to 0.98, 0.3, and 0.01
523
respectively. Regarding the GA, the user-defined parameters were set as follows: population
524
size of 300, and cross over and mutation rates of 0.8 and 0.05, respectively. For having fair
525
comparison, maximum number of function evaluation (NFEs) of 300,000 was assumed as
526
stopping condition for all reported optimizers.
527
The MATLAB programming software was used for coding and implementation purposes.
528
The NT (number of approximation terms of the sine and cosine in Fourier series) is set to 3
529
(i.e., 7 design variables) for all considered ODEs for the longitudinal fins (i.e., 14 ODEs).
530
Hence, approximate mathematical formulations obtained using the PSO for different profiles
531
having different m are given as follows:
Ac ce p
te
d
M
an
us
cr
ip t
515
(a) For rectangular profile (n = 0) for m = 0:
532
θ appx = (0.917338) − (0.093953) cos (π X ) − (0.054890) sin (π X ) − 533
(0.005338) cos ( 2π X ) + (0.026293) sin ( 2π X ) +
.
(33a)
(0.005953) cos ( 3π X ) + (7.67 E − 04) sin ( 3π X ) (b) For rectangular profile (n = 0) for m = 2:
534
28 Page 28 of 52
θ appx = (0.904644) − (0.113495) cos (π X ) − (0.074055) sin (π X ) − 535
(0.010344) cos ( 2π X ) + (0.034491) sin ( 2π X ) +
.
(33b)
(0.007795) cos ( 3π X ) + (0.001691) sin ( 3π X )
θ appx = (0.908804) − (0.104220) cos (π X ) − (0.057281) sin (π X ) − 537
(0.007104) cos ( 2π X ) + (0.027046) sin ( 2π X ) +
.
us
(d) For trapezoidal profile (n = 1) for m = 2:
538
(33c)
cr
(0.005920) cos ( 3π X ) + (0.001063) sin ( 3π X )
ip t
(c) For trapezoidal profile (n = 1) for m = 0:
536
θ appx = (0.896238) − (0.122556) cos (π X ) − (0.076025) sin (π X ) − 539
(0.010682) cos ( 2π X ) + (0.035729) sin ( 2π X ) +
(33d)
.
(33e)
.
(33f)
an
(0.007988) cos ( 3π X ) + (0.001522) sin ( 3π X )
.
(e) For concave parabolic profile (n = 2) for m = 0:
M
540
θ appx = (0.901403) − (0.110344) cos (π X ) − (0.057967) sin (π X ) − 541
(0.005595) cos ( 2π X ) + (0.028214) sin ( 2π X ) +
d
(0.006159) cos ( 3π X ) + (5.12 E − 04) sin ( 3π X ) (f) For concave parabolic profile (n = 2) for m = 2:
te
542
543
Ac ce p
θ appx = (0.892889) − (0.127323) cos (π X ) − (0.078611) sin (π X ) − (0.012050) cos ( 2π X ) + (0.036506) sin ( 2π X ) + (0.008084) cos ( 3π X ) + (0.001865) sin ( 3π X ) 544
Tables 2 to 4 show statistical optimization results obtained by different metaheuristic
545
algorithms for the WRF (cost function), GD, and IGD performance metrics for the series of
546
ODEs given in Eqs. (5) and (6). By observing Table 2 for the cost function (i.e., WRF), all
547
algorithms have tried to reduce (minimize) the existing error between approximate and exact
548
solutions. However, in this competition, the PSO provided better statistical results over the
549
HS and GA offering minimum values for the WRF and evaluation metrics.
29 Page 29 of 52
550
Table 2. Comparison of obtained statistical optimization results for the WRF using the
551
applied optimizers for the longitudinal fins (Eqs. (5) and (6)).
552
m=0 1.53e-02 5.43e-03 1.11e-01 3.50e-02 7.22e-03 1.22e-01 6.14e-02 2.02e-02 5.89e-01 1.27e-02 4.58e-03 2.16e-01 1.13e-05 7.60e-04 5.75e-04 4.57e-05
m=2 2.16e-02 1.57e-02 6.11e-01 5.60e-02 1.74e-02 2.3811 0.12149 3.17e-02 3.9522 3.45e-02 5.00e-03 1.0844 6.21e-04 1.57e-06 6.73e-05 4.51e-05
m=0 9.08e-03 5.20e-03 1.76e-01 4.22e-02 6.93e-03 1.6407 9.03e-02 1.33e-02 2.8048 2.80e-02 2.50e-03 1.0219 1.12e-03 1.08e-05 6.66e-04 1.12e-04
m=2 2.45e-02 1.77e-02 1.2514 5.42e-02 2.34e-02 2.3590 8.41e-02 4.86e-02 4.0512 1.80e-02 1.12e-02 9.31e-01 5.39e-06 9.97e-05 2.19e-05 4.58e-05
ip t
m=2 2.38e-02 1.60e-02 5.24e-01 5.57e-02 1.60e-02 2.3711 8.38e-02 1.62e-02 5.3152 2.15e-02 7.23e-05 1.8409 1.84e-05 1.19e-22 2.83e-03 4.44e-05
n=2
cr
m=0 2.20e-02 8.47e-03 1.36e-01 3.52e-02 1.09e-02 1.0912 5.24e-02 2.56e-02 2.7120 8.84e-03 5.22e-03 7.73e-01 5.11e-07 9.49e-05 1.62e-03 4.54e-05
n=1
us
HS Best WRF PSO GA HS Average WRF PSO GA HS Worst WRF PSO GA HS SD PSO GA HS T-test PSO (p-value=0.01) GA Friedman Test
n=0
an
Methods
M
Measures
Table 3. Comparison of obtained statistical optimization results for the GD metric using the
554
reported optimizers for the longitudinal fins (Eqs. (5) and (6)). Measures
Methods
n=0 m=0 m=2 4.57e-04 2.27e-04 3.95e-04 5.66e-04 4.96e-04 6.36e-04 5.89e-04 7.09e-04 4.85e-04 5.76e-04 3.00e-03 8.00e-03 7.41e-04 1.01e-03 7.91e-04 5.98e-04 8.70e-03 1.83e-02 1.03e-04 2.18e-04 1.30e-04 9.65e-06 2.81e-03 6.95e-03 2.30e-08 2.90e-06 8.92e-07 1.67e-17 8.31e-03 5.41e-03 3.31e-03 2.25e-04
Ac ce p
HS PSO GA HS Average GD PSO GA HS PSO Worst GD GA HS SD PSO GA HS T-test PSO (p-value=0.01) GA Friedman Test
te
d
553
Best GD
n=1 m=0 m=2 1.66e-04 1.53e-04 8.02e-05 2.40-e04 1.20e-03 1.54e-03 3.35e-04 3.00e-04 1.62e-04 3.15e-04 7.20e-03 8.06e-03 5.03e-04 5.89e-04 3.95e-04 5.67e-04 2.24e-02 1.45e-02 1.12e-04 1.38e-04 1.10e-04 9.40e-05 8.28e-03 4.56e-03 5.60e-06 7.61e-05 1.21e-03 2.17e-06 7.44e-06 3.42e-04 1.12e-04 5.00e-04
n=2 m=0 m=2 1.90e-04 1.28e-04 5.37e-05 7.53e-05 7.24e-04 3.76e-03 3.28e-04 3.35e-04 1.65e-04 2.96e-04 5.67e-03 8.43e-03 5.15e-04 7.25e-04 4.14e-04 7.93e-04 1.11e-02 1.43e-02 1.00e-04 1.97e-04 1.25e-04 1.84e-04 4.02e-03 3.52e-03 2.77e-06 4.55e-04 2.41e-03 6.65e-04 1.63e-03 3.47e-05 2.25e-04 5.53e-04
555 30 Page 30 of 52
556
Table 4. Comparison of obtained statistical optimization results for the IGD metric using the
557
considered optimizers for the longitudinal fins (Eqs. (5) and (6)).
558
m=0 1.34e-03 6.87e-04 9.08e-03 2.82e-03 1.35e-03 5.76e-02 4.60e-03 3.29e-03 1.74e-01 1.03e-03 9.61e-04 6.44e-02 1.20e-05 1.60e-03 3.21e-03 4.54e-05
m=2 1.21e-03 2.15e-03 1.42e-02 2.49e-03 2.83e-03 6.72e-02 4.67e-03 5.24e-03 1.19e-01 1.15e-03 8.92e-04 3.61e-02 7.44e-05 3.44e-06 2.35E-04 5.00e-04
m=0 1.71e-03 4.59e-04 5.78e-03 2.64e-03 1.40e-03 4.76e-02 4.21e-03 3.73e-03 9.21e-02 8.04e-04 1.09e-03 3.31e-02 2.54e-06 2.91e-03 1.41e-03 2.25e-04
m=2 1.09e-02 6.08e-04 3.04e-02 2.78e-03 2.56e-03 6.63e-02 6.24e-03 6.52e-03 1.12e-01 1.76e-03 1.50e-03 2.83e-02 7.46e-04 4.50e-04 4.14e-05 5.53e-04
ip t
m=2 1.73e-03 4.92e-03 4.32e-03 6.16e-03 5.02e-03 6.20e-02 9.17e-03 5.22e-03 1.45e-01 2.02e-03 9.47e-05 5.45e-02 4.97e-06 2.80e-06 5.81e-03 3.33e-03
n=2
cr
m=0 4.17e-03 3.30e-03 3.59e-03 5.24e-03 4.16e-03 2.44e-02 6.39e-03 6.48e-03 6.85e-02 7.96e-04 1.07e-03 2.22e-02 6.38e-09 6.50e-07 7.10e-03 3.42e-03
n=1
us
HS Best IGD PSO GA HS Average IGD PSO GA HS Worst IGD PSO GA HS SD PSO GA HS T-test PSO (p-value=0.01) GA Friedman Test
n=0
an
Methods
M
Measures
Talking about the GD and IGD metrics (Tables 3 and 4), all applied optimizers performed
560
well for all profiles (n=0, 1, and 2) and heat transfer coefficients (m). In general speaking, all
561
three optimizers performed well in finding unknown coefficients of Fourier series and
562
confirmed the obtained optimization results.
563
In addition, two statistical analyses (i.e., T-test and Friedman test have) have been included in
564
this paper with the aim of proving that differences in terms of WRF, GD, and IGD are
565
statistically significant with a p-value of 0.01. T-test performs a statistical analysis of the null
566
hypothesis that data in a sample set are a random sample from a normal distribution with
567
mean 0 and unknown variance, against the alternative that the mean is not 0. In simple terms,
568
the T-test compares the actual difference between two means in relation to the variation in the
569
data.
570
Moreover, p-values for Friedman test are included in the last row of Tables 2 to 4. The null
571
hypothesis of the Friedman test is that the calculated errors obtained from the various
Ac ce p
te
d
559
31 Page 31 of 52
optimizers are not different. Such a null hypothesis is rejected if p-value is smaller than the
573
level of significance (less than α=0.01).
574
If the p-value is less than α, we can reject the null at the α level of significance, and we can
575
say there is significance evidence against the null hypothesis. In contrast, if the p-value is
576
greater than α, then we do not have significant evidence against the null hypothesis, and
577
therefore, the null hypothesis cannot be rejected.
578
Indeed, the smaller the p-value, the greater the evidence against the null hypothesis. From
579
Tables 2 to 4, obtained solutions indicate rejections of the null hypothesis at the 1%
580
significance level for the both reported statistical analyses.
581
Fig. 3 demonstrates the reduction rate for the error function, WRF (cost function), using the
582
PSO for the two selected ODEs for the rectangular profile (n = 0) with m = 0 and concave
583
parabolic profile (n = 2) with m = 2 for only 100 iterations. Looking at Fig. 3, the PSO has
584
tried to reduce the cost function (i.e., error) at each iteration. Finally, the PSO has converged
585
to its optimized approximate solution after 1000 iterations.
Ac ce p
te
d
M
an
us
cr
ip t
572
32 Page 32 of 52
ip t cr us an M
586
Fig. 3. Error history (cost reduction) of two selected ODEs using the PSO for
588
different profiles and m: (a) n = 0 and m = 0, (b) n = 2 and m = 2.
te
d
587
In addition, Tables 5 to 7 represent obtained optimization results for the WRF and GD, and
590
IGD evaluators for special cases of a longitudinal fin (i.e., rectangular fin, Eqs. (8) and (9))
591
for different m. In order to have fair comparisons with the numerical results offered in the
592
literature [50], same values for m are considered for the proposed approach.
Ac ce p
589
593
Table 5. Comparison of obtained statistical optimization results for the WRF (Eqs. (8) and
594
(9)) having different m.
Measures Best WRF
Average WRF
Worst WRF
Methods HS PSO GA HS PSO GA HS PSO GA
m=0 8.61e-02 8.10e-02 2.57e-01 1.18e-01 8.11e-02 2.1604 1.43e-01 8.14e-02 4.3624
m=0.25 1.00e-01 7.80e-02 2.98e-01 1.25e-01 7.82e-02 2.2007 1.59e-01 7.91e-02 4.3417
n=0 (Rectangular Longitudinal Fin) m=0.5 m=1 m=1.25 m=2 7.77e-02 7.33e-02 8.10e-02 7.81e-02 7.56e-02 7.19e-02 7.04e-02 6.70e-02 2.26e-01 1.36e-01 1.80e-01 1.3739 1.12e-01 1.05e-01 8.0254 1.08e-01 7.59e-02 7.49e-02 7.12e-02 6.76e-02 2.3423 1.1164 2.0039 2.4050 1.57e-01 1.50e-01 39.9085 1.44e-01 7.79e-02 8.40e-02 7.45e-02 7.02e-02 7.3030 2.522 5.7414 3.8327
m=3 8.34e-02 6.41e-02 5.61e-01 1.08e-01 6.43e-02 2.5076 1.68e-01 6.51e-02 3.9252
m=4 7.14e-02 6.22e-02 3.86e-01 9.79e-02 6.28e-02 2.7518 1.27e-01 6.62e-02 5.3644
33 Page 33 of 52
HS PSO GA HS T-test PSO (p-value=0.01) GA Friedman Test SD
1.91e-02 1.32e-04 1.2706 1.08e-08 1.28e-26 4.46e-04 4.52e-05
1.94e-02 3.98e-04 1.3284 7.51e-09 3.70e-22 5.36e-04 4.64e-05
2.69e-02 7.02e-04 1.9998 3.33e-07 7.95e-20 4.91e-03 4.24e-05
2.55e-02 4.19e-03 8.15e-01 3.57e-07 8.51e-13 1.91e-03 3.04e-04
16.6951 1.41e-03 1.7430 5.28e-06 7.49e-17 5.41e-03 4.24e-05
2.15e-02 1.08e-03 7.96e-01 6.91e-08 1.06e-17 5.27e-06 4.47e-05
2.40e-02 3.85e-04 1.2524 1.80e-07 1.60e-21 1.36e-04 4.63e-05
1.77e-02 1.22e-03 1.3715 2.88e-08 6.28e-17 1.34e-04 4.62e-05
595
Table 6. Comparison of obtained statistical optimization results for the GD metric using the
597
applied optimizers considering various m (Eqs. (8) and (9)). n=0 (Rectangular Longitudinal Fin) m=0.5 m=1 m=1.25 m=2 4.39e-04 3.74e-04 3.63e-04 2.82e-04 5.84e-04 3.52e-04 5.30e-04 4.15e-04 3.66e-04 3.92e-04 4.70e-04 1.49e-03 6.69e-04 5.01e-04 3.94e-02 4.76e-04 6.15e-04 6.20e-04 6.00e-04 5.63e-04 6.11e-03 2.47e-03 4.76e-03 5.84e-03 8.92e-04 7.57e-04 1.98e-01 6.78e-04 6.36e-04 8.17e-04 7.49e-04 6.65e-04 1.95e-02 7.13e-03 1.60e-02 9.20e-03 1.41e-04 1.17e-04 8.17e-02 1.29e-04 1.38e-05 1.17e-04 5.83e-05 6.12e-05 5.47e-03 2.20e-03 4.82e-03 2.42e-03 1.13e-07 2.75e-07 2.88e-05 9.88e-07 2.32e-16 4.21e-08 1.20e-10 3.25e-10 6.32e-03 6.21e-03 4.82e-03 3.16e-05 5.52e-03 4.51e-02 6.08e-02 5.00e-04
m=3 2.15e-04 5.13e-04 6.60e-04 5.06e-04 5.74e-04 5.88e-03 7.54e-04 6.81e-04 1.11e-02 1.78e-04 5.28e-05 3.37e-03 8.50e-06 7.33e-11 3.74e-04 2.22e-03
m=4 1.85e-04 5.39e-04 9.33e-04 5.09e-04 5.60e-04 6.88e-03 8.11e-04 6.14e-04 1.54e-02 1.73e-04 2.52e-05 4.02e-03 6.57e-06 1.20e-13 4.26e-04 5.00e-04
d
598
m=0.25 4.55e-04 6.11e-04 1.06e-03 6.30e-04 6.28e-04 5.00e-03 8.71e-04 6.51e-04 1.17e-02 1.56e-04 9.64e-06 3.66e-03 4.49e-07 7.64e-18 1.91e-03 5.00e-04
us
HS Best GD PSO GA HS Average GD PSO GA HS Worst GD PSO GA HS SD PSO GA HS T-test PSO (p-value=0.01) GA Friedman Test
m=0 5.15e-04 6.22e-04 5.75e-04 7.21e-04 6.48e-04 5.35e-03 9.44e-04 6.99e-04 1.35e-02 1.41e-04 2.46e-05 4.06e-03 5.79e-08 2.67e-14 2.42e-03 3.31e-03
an
Methods
M
Measures
cr
ip t
596
Table 7. Comparison of obtained statistical optimization results for the IGD metric using the
600
reported optimizers considering various m (Eqs. (8) and (9)).
Ac ce p
te
599
Measures
Methods
HS Best IGD PSO GA HS Average IGD PSO GA HS Worst IGD PSO GA HS SD PSO GA HS T-test PSO (p-value=0.01) GA Friedman Test
m=0 4.39e-03 5.61e-03 4.64e-03 6.41e-03 5.84e-03 3.83e-02 8.49e-03 6.30e-03 1.01e-01 1.35e-03 2.20e-04 2.96e-02 1.09e-07 2.50e-14 2.72e-03 3.33e-03
m=0.25 3.99e-03 5.52e-03 8.24e-03 5.58e-03 5.66e-03 3.46e-02 7.83e-03 5.86e-03 8.59e-02 1.39e-03 8.52e-05 2.64e-02 4.78e-07 6.35e-18 2.51e-03 5.00e-04
n=0 (Rectangular Longitudinal Fin) m=0.5 m=1 m=1.25 m=2 3.86e-03 3.16e-03 3.14e-03 2.21e-03 5.28e-03 3.18e-03 4.79e-03 3.77e-03 2.89e-03 3.24e-03 3.81e-03 1.30e-02 6.00e-03 4.40e-03 3.67e-01 4.07e-03 5.55e-03 5.60e-03 5.41e-03 5.08e-03 4.56e-02 1.75e-02 3.54e-02 4.40e-02 7.95e-03 6.72e-03 1.84e-01 5.94e-03 5.74e-03 7.32e-03 6.65e-03 6.00e-03 1.46e-01 5.02e-02 1.19e-01 6.73e-02 1.30e-03 1.09e-03 7.63e-01 1.21e-03 1.20e-04 1.05e-03 4.94e-04 5.45e-04 4.07e-02 1.53e-02 3.61e-02 1.87e-02 1.47e-07 4.60e-07 6.26e-07 2.19e-06 1.69e-16 3.92e-08 6.91e-11 2.94e-10 6.33e-03 5.62e-03 3.22e-03 3.96e-05 5.52e-03 1.36e-02 6.08e-02 3.71e-04
m=3 1.70e-03 4.64e-03 3.96e-03 4.49e-03 5.18e-03 4.42e-02 6.91e-03 6.13e-03 8.24e-02 1.69e-03 4.64e-04 2.55e-02 1.51e-05 5.82e-11 3.90e-04 7.41e-03
m=4 1.74e-03 4.87e-03 7.77e-03 4.44e-03 5.06e-03 5.36e-02 7.41e-03 5.54e-03 1.16e-01 1.54e-03 2.26e-04 3.03e-02 7.91e-06 1.15e-13 3.40e-04 3.71e-04
601
34 Page 34 of 52
By observing Table 5, in terms of the WRF’s values, the PSO has outperformed the HS and
603
GA, while outperforming the HS over the PSO and GA for the best obtained GD and IGD for
604
almost all cases (see Tables 6 and 7). Looking at Tables 5 to 7, obtained results prove
605
rejections of the null hypothesis at the 1% significance level for the both considered statistical
606
analyses.
607
In fact, the WRF method equipped with optimization methods has tried to force the residual
608
to be zero in some average sense over the domain. When the value of cost function reaches
609
zero or close to zero, better approximate solutions and more accuracy can be achieved. As
610
can be seen in Tables 3 and 5, the obtained values of WRF are close enough to zero for all
611
profiles having different m.
4.3.
Comparison of approximate solutions
M
612
an
us
cr
ip t
602
This section aims to compare the obtained approximate solutions by the DTM [39] and
614
numerical solutions by the Maple 14.0 [50] with WRF approach using the reported optimizers
615
(i.e., WRF with PSO, WRF via HS, and WRF using GA). In terms of complexity of
616
approximate functions as seen in Eqs. (33), one can be grasped that the obtained approximate
617
functions using the WRF equipped with PSO (WRF-PSO) are in a simpler mathematical
618
forms than those suggested by the DTM (see Eqs. (30)) [39].
619
Table 8 shows the comparison between the DTM results and the WRF method for different
620
profiles (n) and m in terms of the GD and IGD metrics and values of BCs. It is worth
621
mentioning that the BC1 and BC2, shown in Table 8, stand for the first (Eq. (6a)) and second
622
(Eq. (6b)) boundary conditions.
623
However, in order to show values of these boundary conditions in a better way (accuracy for
624
satisfying boundary conditions), we considered the differences between the exact (numerical)
625
solutions and estimated boundary conditions (e.g., BC1=1 - BC1).
Ac ce p
te
d
613
35 Page 35 of 52
626
Table 8. Comparison of best obtained approximate solution using the DTM and WRF
627
approach for Eqs. (5) and (6).
WRF using HS
WRF using GA
ip t
M
628
n=2 m=0 m=2 1.63e-06 1.47e-06 1.61e-05 1.41e-05 5.20e-03 1.77e-02 5.37e-05 7.53e-05 4.59e-04 6.08e-04 4.17e-14 0 -8.17e-10 0 9.08e-03 2.45e-02 1.90e-04 1.28e-04 1.71e-03 1.09e-02 3.80e-06 3.29e-06 5.06e-06 -1.32e-05 1.76e-01 1.2514 7.24e-04 3.76e-03 5.78e-03 3.04e-02 2.28e-03 9.47e-02 -2.03e-05 -6.89e-10
cr
WRF using PSO
GD IGD Error GD IGD BC1 BC2 Error GD IGD BC1 BC2 Error GD IGD BC1 BC2
n=1 m=0 m=2 1.24e-07 4.62e-05 1.23e-06 4.45e-04 5.43e-03 1.57e-02 8.02e-05 2.40-e04 6.87e-04 2.15e-03 6.26e-10 0 -1.04e-08 0 1.53e-02 2.16e-02 1.66e-04 1.53e-04 1.34e-03 1.21e-03 1.19e-06 6.85e-07 -1.30e-06 1.32e-06 1.11e-01 6.11e-01 1.20e-03 1.54e-03 9.08e-03 1.42e-02 0 3.63e-02 6.93e-18 9.50e-12
us
DTM [39]
n=0 m=0 m=2 1.85e-11 8.93e-11 1.50e-10 8.58e-10 8.47e-03 1.60e-02 3.95e-04 5.66e-04 3.30e-03 4.92e-03 0 0 0 0 2.20e-02 2.38e-02 4.57e-04 2.27e-04 4.17e-03 1.73e-03 3.59e-06 1.52e-06 2.76e-05 4.33e-06 1.36e-01 5.24e-01 4.96e-04 6.36e-04 3.59e-03 4.32e-03 2.61e-03 2.75e-02 1.16e-10 -6.90e-05
an
Methods
Indeed, Table 8 summarizes the results and reports for estimated solutions in more detail.
630
From Table 8, in terms of GD, IGD, and BC values, the DTM has surpassed the WRF-PSO,
631
WRF-HS, and WRF-GA. However, the offered approximate solutions by the WRF equipped
632
with optimization methods have acceptable accuracies and values of BCs are satisfactory.
633
For graphical comparison and presentation, Fig. 4 depicts the comparisons of approximate
634
solutions using the WRF-PSO and the DTM [39] against exact and numerical solutions. As
635
can be seen from Fig. 4, the approximate solutions obtained by the WRF-PSO are very close
636
to the solutions attained by the DTM.
Ac ce p
te
d
629
36 Page 36 of 52
ip t cr us an M d te
637
Fig. 4. Comparison of approximate solutions using the WRF equipped with the PSO (dashed
639
colored line), DTM results [39] (solid line), and the exact analytical and numerical solutions
640
(solid box) for n = 0 (rectangular fin), n = 1 (trapezoidal fin), and n = 2 (concave parabolic
641
fin): a) m = 0, b) m = 2.
Ac ce p
638
642
Except of longitudinal fins having rectangular profile (n = 0) for both m (i.e., m = 0 and 2),
643
estimated solutions by the WRF-PSO for other profiles (i.e., trapezoidal and concave
644
parabolic profiles) are competitive with those obtained by the DTM (see Fig. 4). Similarly, in
645
order to compare with the numerical solutions for Eqs. (8) and (9), Table 9 shows the values
646
of BCs, GD, and IGD metrics for the all considered approaches.
37 Page 37 of 52
647
Table 9. Comparison of best obtained approximate solution using the considered optimization
648
methods with numerical solutions [50] for Eqs. (8) and (9) considering various m.
WRF using GA
n=0 (Rectangular Longitudinal Fin) m=0.5 m=1 m=1.25 m=2 7.56e-02 7.19e-02 7.04e-02 6.70e-02 5.84e-04 3.52e-04 5.30e-04 4.15e-04 5.28e-03 3.18e-03 4.79e-03 3.77e-03 0 0 2.22e-16 3.33e-16 0 0 1.93e-10 1.45e-15 7.77e-02 7.33e-02 8.10e-02 7.81e-02 4.39e-04 3.74e-04 3.63e-04 2.82e-04 3.86e-03 3.16e-03 3.14e-03 2.21e-03 8.58e-05 1.00e-05 1.73e-06 8.48e-06 -1.07e-05 -1.74e-05 -5.93e-06 -6.60e-05 2.26e-01 1.36e-01 1.80e-01 1.3739 3.66e-04 3.92e-04 4.70e-04 1.49e-03 2.89e-03 3.24e-03 3.81e-03 1.30e-02 9.97e-04 4.84e-04 1.68e-05 0.1100 9.07e-07 2.65e-06 0 8.67e-18
m=3 6.41e-02 5.13e-04 4.64e-03 0 0 8.34e-02 2.15e-04 1.70e-03 2.46e-06 1.49e-06 5.61e-01 6.60e-04 3.96e-03 2.88e-02 -1.73e-17
649
m=4 6.22e-02 5.39e-04 4.87e-03 4.68e-14 -4.50e-11 7.14e-02 1.85e-04 1.74e-03 3.70e-06 -1.95e-05 3.86e-01 9.33e-04 7.77e-03 1.37e-02 -1.82e-06
ip t
WRF using HS
m=0.25 7.80e-02 6.11e-04 5.52e-03 4.65e-10 -2.08e-17 1.00e-01 4.55e-04 3.99e-03 2.46e-06 -1.06e-05 2.98e-01 1.06e-03 8.24e-03 5.35e-03 2.08e-17
cr
WRF using PSO
Error GD IGD BC1 BC2 Error GD IGD BC1 BC2 Error GD IGD BC1 BC2
m=0 8.10e-02 6.22e-04 5.61e-03 3.46e-08 2.08e-17 8.61e-02 5.15e-04 4.39e-03 1.02e-05 3.84e-07 2.57e-01 5.75e-04 4.64e-03 5.23e-03 -3.73e-12
us
Methods
Fig. 5 demonstrates comparisons among numerical and approximate solutions by the WRF
651
via PSO for different heat transfer coefficients (m). From Fig. 5, the WRF with PSO has
652
obtained approximate solutions with acceptable accuracy close to the numerical solutions in
653
the literature [50].
Ac ce p
te
d
M
an
650
654 655
Fig. 5. Comparison of approximate solutions using the WRF equipped with the PSO (dashed
656
lines) and the numerical solutions (solid lines) [50] for different m for Eqs. (8) and (9).
38 Page 38 of 52
It is worth mentioning that sometimes due to the complexity of an ODE, there is not any
658
analytical way (analytical solution) to solve it, or at least it is a consuming task accompanied
659
with possessing high skills in mathematics. However, using the WRF approach, it is useful to
660
suggest a handy mathematical formulation.
661
However, the DTM has its own operational restrictions that severely narrow its functioning
662
domain. Therefore, it is possible that the DTM fails to overcome a specific problem. For
663
instance, it was reported that the DTM was unable to produce physically reasonable data for
664
the Glauert-jet Problem [11].
665
The DTM is suitable for approximate solving of IVPs, while the WRF approach using
666
metaheuristic algorithms can tackle both IVPs and BVPs. Therefore, there is no restriction for
667
the proposed technique for handling ODEs having different boundary conditions. Meanwhile,
668
the DTM is based on classical mathematical tools and in contrast, for the WRF equipped with
669
metaheuristic algorithms, there are no limitations in terms of orders and natures of a problem.
670
They can handle a wide variety of nonlinear ODEs despite of their orders and natures. The
671
ODEs problems of longitudinal fins were examples of highly nonlinear ODEs which,
672
recently, are actively investigating in the literature [39, 50]. Other ODE problems in the
673
literature can be solved using the proposed WRF approach with the aid of optimization
674
methods.
675
In general, firstly, approximate analytical methods give accurate predictions only when the
676
nonlinearities are weak. Second, the approximate analytical methods often involve a complex
677
mathematical analysis leading to analytic expressions with a large number of terms (see Eq.
678
(30) for the DTM). Finally, methods such as the HPM, HAM, ADM, and VIM, if routinely
679
implemented, can sometimes lead to erroneous results as shown and reported in the literature
680
[51].
Ac ce p
te
d
M
an
us
cr
ip t
657
39 Page 39 of 52
All in all, the WRF approach equipped with LSW and reported optimization techniques has
682
suggested acceptable approximate solution for solving ODEs of longitudinal fins having
683
different profiles and properties.
684
It is worth pointing out that when there is an exact and analytical solution for a problem in a
685
logical time, there is no need to solve the problem using approximate methods. However, real
686
life problems, as we considered in this paper (longitudinal fins having various profiles and
687
features, Eqs. (5), (6), (8), and (9)), are highly nonlinear [39, 50] and analytical approaches
688
are often inefficient for tackling ODEs.
us
cr
ip t
681
5. Conclusions
an
689
This paper introduced a novel approximate approach for solving ordinary differential
691
equations (ODEs) in engineering, particularly longitudinal fins. Finding the exact solution of
692
an ODE often demands heavy computational effort; moreover, it is not an achievable task
693
sometimes. In this situation, the need for approximate solutions arises. Metaheuristic
694
algorithms are well-known approximate methods for solving and optimizing problems. Using
695
the Fourier series as a base approximate function, the ODEs can be modeled as an
696
optimization problem.
697
The aim is to minimize the error of the weighted residual function (WRF) equipped with new
698
proposed weight function, called least square weight (LSW) function. Boundary conditions
699
and initial values are imposed as constraints. In order to compare the results obtained using
700
the WRF (error function) equipped with metaheuristic methods with the exact (numerical)
701
solutions, the generational distance (GD) and inverted GD (IGD) metrics were used.
702
Longitudinal fins having different profiles and properties were examined using the proposed
703
method. For optimizing purposes, three metaheuristic methods including the particle swarm
704
optimization (PSO), the genetic algorithm (GA), and the harmony search (HS) were applied
705
and compared in terms of statistical optimization results and statistical analyses.
Ac ce p
te
d
M
690
40 Page 40 of 52
The best approximate solutions obtained by the WRF approach were compared with the
707
DTM results. The proposed approach offers approximate solutions having acceptable
708
accuracies compared with numerical solutions. Furthermore, obtained approximate solutions
709
using the WRF approach with finite series of cosines and sins shows its performance in terms
710
of acceptable accuracy (error).
711
Due to the nonlinearity nature and behavior of ODE for the most real life problems, the
712
proposed WRF having LSW function accompanied with optimization methods may consider
713
as an alternative approximate approach for efficiently tackling ODEs particularly longitudinal
714
fins.
715
As further study, other applications of the proposed method may be carried out for solving
716
highly nonlinear ODEs in engineering, physics, economics, fluids, and so forth. Different
717
weight functions may be applied and comparative study with other metaheuristic-based
718
approximate methods given in introduction section may be considered as other contributions
719
among these types of approximation methods.
720
Acknowledgment
721
This work was supported by the National Research Foundation of Korean (NRF) grant
722
funded by the Korean government (MSIP) (NRF-2013R1A2A1A01013886).
723
References
724
[1] S.B. Coşkun, M.T. Atay, Analysis of convective straight and radial fins with temperature-
725
dependent thermal conductivity using variational iteration method with comparison with
726
respect to finite element analysis, Math. Probl. Eng. ID 42072 (2007) 1-15.
727
[2] S.B. Coşkun, M.T. Atay, Fin efficiency analysis of convective straight fins with
728
temperature dependent thermal conductivity using variational iteration method, Appl. Therm.
729
Eng. 28 (2008) 2345-2352.
Ac ce p
te
d
M
an
us
cr
ip t
706
41 Page 41 of 52
[3] G. Domairry, M. Fazeli, Homotopy analysis method to determine the fin efficiency of
731
convective straight fins with temperature-dependent thermal conductivity, Commun.
732
Nonlinear Sci. 14 (2009) 489-499.
733
[4] Z.Y. Lee, Method of bilaterally bounded to solution Blasius equation using particle
734
swarm optimization, Appl. Math. Comput. 179 (2006) 779-786.
735
[5] C.H. Chiu, C.K. Chen, Application of the decomposition method to thermal stresses in
736
isotropic circular fins with temperature-dependent thermal conductivity, Acta Mechanica.
737
157 (2002) 147-158.
738
[6] C. Arslanturk, Correlation equations for optimum design of annular fins with temperature
739
dependent thermal conductivity, Heat Mass. Transfer 45(4) (2009) 519-525.
740
[7] J.K Zhou, Differential transform and its applications for electrical circuits, Wuhan,
741
Huarjung University Press, 1986
742
[8] M.M. Rashidi, N. Laraqi, S.M. Sadri, 2010. A novel analytical solution of mixed
743
convection about an inclined flat plate embedded in a porous medium using the DTM-Padé,
744
Int. J. Therm. Sci. 49 (2010) 2405-2412.
745
[9] B. Kundu, D. Barman, 2010, Analytical study on design analysis of annular fins under
746
dehumidifying conditions with a polynomial relationship between humidity ratio and
747
saturation temperature, Int. J. Heat Fluid Fl. 31 (2010) 722-733.
748
[10] H. Yaghoobi, M. Torabi, The application of differential transformation method to
749
nonlinear equations arising in heat transfer, Int. Commun. Heat Mass 38 (2011) 815-820.
750
[11] M. Torabi, H. Yaghoobi, A. Fereidoon, Application of differential transformation
751
method and padé approximant for the Glauert-jet problem, Recent Patents on Mechanical
752
Engineering 5 (2012) 150-155.
753
[12] H.P. Chu, C.L. Chen, Hybrid differential transform and finite difference method to solve
754
the nonlinear heat conduction problem, Commun. Nonlinear Sci. 13 (2008) 1605-1614.
Ac ce p
te
d
M
an
us
cr
ip t
730
42 Page 42 of 52
[13] H.P. Chu, C.Y. Lo, Application of the hybrid differential transform-finite difference
756
method to nonlinear transient heat conduction problems, Numer. Heat Transfer: Part A 53
757
(2008) 295-307.
758
[14] H.S. Peng, C.L. Chen, Hybrid differential transformation and finite difference method to
759
annular fin with temperature-dependent thermal conductivity, Int. J. Heat Mass Tran. 54
760
(2011) 2427-2433.
761
[15] M. Torabi, H. Yaghoobi, K. Boubaker, Accurate solution for motion of a spherical solid
762
particle in plane Couette Newtonian fluid mechanical flow using HPM–Padé approximant
763
and the Boubaker polynomials expansion scheme BPES, Int. J. Heat Mass Tran. 58(1-2)
764
(2013) 224-228.
765
[16] S. Yalcinbas, M. Sezer, The approximate solution of high-order linear Volterra–
766
Fredholm integro-differential equations in terms of Taylor polynomials, Appl. Math. Comput.
767
112 (2000) 291-308.
768
[17] P. Darania, E. Abadian, Development of the Taylor expansion approach for nonlinear
769
integro-differential equations, IJCMS 1(14) (2006) 651-664.
770
[18] P. Darania, A. Ebadian, A method for the numerical solution of the integro-differential
771
equations, Appl. Math. Comput. 188 (2007) 657-668.
772
[19] P. Darania, K. Ivaz, Numerical solution of nonlinear Volterra–Fredholm integro-
773
differential equations, Comput. Math. Appl. 56 (2008) 2197-2209.
774
[20] P. Roul, P. Meyer, Numerical solutions of systems of nonlinear integro-differential
775
equations by Homotopy-perturbation method, Appl. Math. Model. 35 (2011) 4234-4242.
776
[21] M. Torabi, H. Yaghoobi, Novel solution for acceleration motion of a vertically falling
777
spherical particle by HPM-Padé approximant, Adv. Powder Technol. 22(5) (2011) 674-677.
778
[22] H. Yaghoobi, M. Torabi, Novel solution for acceleration motion of a vertically falling
779
non-spherical particle by VIM–Padé approximant, Powder Technol. 215-216 (2012) 206-209.
Ac ce p
te
d
M
an
us
cr
ip t
755
43 Page 43 of 52
[23] I.H. Osman, G. Laporte, Metaheuristics: a bibliography, Ann. Oper. Res. 63 (1996) 513-
781
623.
782
[24] F. Glover, G.A. Kochenberger, G., Handbook of metaheuristics, Kluwer Academic
783
Publishers, 2003.
784
[25] X.S. Yang, Nature-Inspired metaheuristic algorithms, second ed., Luniver Press, 2010.
785
[26] X.S. Yang, Engineering optimization: an introduction with metaheuristic applications,
786
John Wiley & Sons, 2010.
787
[27] G.D. Mateescu, On the application of genetic algorithms to differential equations,
788
Romanian Journal of Economic Forecasting 2 (2006) 5-9.
789
[28] E.N. Mastorakis, Unstable ordinary differential equations: solution via genetic
790
algorithms and the method of Nelder-Mead, In: Proceedings of the sixth WSEAS Int. Conf.
791
on Systems Theory & Scientific Computation, 119, Elounda, Greece, (2007), pp. 297-354.
792
[29] H. Cao, L. Kang, Y. Chen, Evolutionary modeling of systems of ordinary differential
793
equations with genetic programming, Genet. Program. Evol. M. 1 (2000) 309-337.
794
[30] C. Reich, Simulation of imprecise ordinary differential equations using evolutionary
795
algorithms, In: Proceedings of the 2000 ACM Symposium on Applied Computing 1, (2000)
796
pp. 428-432.
797
[31] C.L. Karr, E. Wilson, A self-tuning evolutionary algorithm applied to an inverse partial
798
differential equation, Appl. Intell. 19 (2003) 147-155.
799
[32] M. Babaei, A general approach to approximate solutions of nonlinear differential
800
equations using particle swarm optimization, Appl. Soft Comput. 13 (2013) 3354-3365.
801
[33] A. Sadollah, H. Eskandar, D.G. Yoo, J.H Kim, Approximate solving of nonlinear
802
ordinary differential equations using least square weight function and metaheuristic
803
algorithms, Eng. Appl. Artif. Intel. 40 (2015) 117-132.
Ac ce p
te
d
M
an
us
cr
ip t
780
44 Page 44 of 52
[34] J. Kennedy, R. Eberhart, Particle swarm optimization. Proc. of the IEEE international
805
conference on neural networks, Perth, Australia, (1995), pp 1942-1948.
806
[35] H. Eskandar, A. Sadollah, A. Bahreininejad, M. Hamdi, Water cycle algorithm - a novel
807
metaheuristic optimization method for solving constrained engineering optimization
808
problems, Comput. Struct. 110-111 (2012) 151-166.
809
[36] G.W. Geem, J.H. Kim, G.V. Loganathan, A new heuristic optimization algorithm:
810
harmony search, Simulation 76(2) (2001) 60-68.
811
[37] D. Goldberg, Genetic algorithms in search, optimization and machine learning, Addison-
812
Wesley, Reading, MA, 1989.
813
[38] A.D. Kraus, A. Aziz, J.R. Welty, Extended surface heat transfer, New York: John Wiley,
814
2002.
815
[39] M. Torabi, A. Aziz, K. Zhang, A comparative study of longitudinal fins of rectangular,
816
trapezoidal and concave parabolic profiles with multiple nonlinearities, Energy 51 (2013)
817
243-256.
818
[40] W.E. Boyce, R.C. Diprima, Elementary differential equations and boundary value
819
problems, Sixth ed., John Wiley & Sons, New York, 1997.
820
[41] E. Kreyszig, Advanced Engineering Mathematics, ninth ed., John Wiley & Sons, New
821
York, 2009.
822
[42] K.J. Bathe, Finite Element Procedures, second ed., Prentice Hall, New Jersey, 1996.
823
[43] C.A.C. Coello, G.T. Pulido, Multiobjective structural optimization using a micro
824
genetic algorithm, Struct. Multidiscip. O. 30(5) (2005) 388-403.
825
[44] W. Gong, Z. Cai, L. Zhu, An efficient multiobjective differential evolution algorithm for
826
engineering design, Struct. Multidiscip. O. 38(2) (2009) 137-157.
827
[45] D.A.V. Veldhuizen, G.B. Lamont, Multiobjective evolutionary algorithm research: A
828
history and analysis. Technical Report TR-98-03, Department of Electrical and Computer
Ac ce p
te
d
M
an
us
cr
ip t
804
45 Page 45 of 52
Engineering, Graduate School of Engineering, Air Force Institute of Technology, Wright-
830
Patterson AFB, Ohio, 1998.
831
[46] E. Zitzler, L. Thiele, M. Laumanns, C.M. Fonseca, V.G. da Fonseca, Performance
832
assessment of multiobjective optimizers: an analysis and review, IEEE Trans. Evol. Comput.
833
7 (2) (2003) 117-132.
834
[47] Y. Dong, J. Tang, B. Xu, D. Wang, An application of swarm optimization to nonlinear
835
programming, Comput. Math. Appl. 49 (2005) 1655-1668.
836
[48] D.G. Yoo, J.H. Kim, Z.W. Geem, Overview of harmony search algorithm and its
837
applications in civil engineering, Evol. Intel. 7 (2014) 3-16.
838
[49] J.H. Kim, Z. Geem, E. Kim, Parameter estimation of the nonlinear Muskingum model
839
using harmony search, J. Am. Water Resour. Assoc. 37(5) (2001) 1131-1138.
840
[50] A. Aziz, M. Torabi, Convective-Radiative fins with simultaneous variation of thermal
841
conductivity, heat transfer coefficient, and surface emissivity with temperature, Heat Tran.
842
Asian Res. 41(2) (2012) 99-113.
843
[51] A. Fernandez, On some approximate methods for nonlinear models, Appl. Math.
844
Comput. 215 (2009) 168-174.
Ac ce p
te
d
M
an
us
cr
ip t
829
46 Page 46 of 52
Ac
ce
pt
ed
M
an
us
cr
i
*Graphical abstract (for review)
Page 47 of 52
Ac
ce
pt
ed
M
an
us
cr
i
Figure 1
Page 48 of 52
Ac
ce
pt
ed
M
an
us
cr
i
Figure 2
Page 49 of 52
Ac ce p
te
d
M
an
us
cr
ip t
Figure 3
Page 50 of 52
Ac ce p
te
d
M
an
us
cr
ip t
Figure 4
Page 51 of 52
Ac
ce
pt
ed
M
an
us
cr
i
Figure 5
Page 52 of 52