Metaheuristic algorithms for approximate solution to

0 downloads 0 Views 920KB Size Report
Feb 19, 2015 - Differential equations play a noticeable role in engineering, physics, economics, and other. 13 ... With the aid of certain fundamental concepts of mathematics,. 18 ..... for any type of ODEs for higher orders in their implicit forms.
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 





Approximate solutions to ordinary differential equations (ODEs) in engineering.  





Fourier series with the aid of metaheuristics used as suggested approximate method. 





The GA, PSO and HS are utilized for optimization purposes. 





Obtained approximate solutions by the proposed method confirm the results by others. 





Proposed approach offers acceptable accuracy for a wide range of ODEs. 

us

cr

ip t



Ac ce p

te

d

M

an



1    Page 1 of 52

Metaheuristic algorithms for approximate solution to ordinary differential



equations of longitudinal fins having various profiles



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



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