Jtam-A4.dvi


JOURNAL OF THEORETICAL

AND APPLIED MECHANICS

53, 1, pp. 15-26, Warsaw 2015
DOI: 10.15632/jtam-pl.53.1.15

RECURSIVE DIFFERENTIATION METHOD: APPLICATION TO THE
ANALYSIS OF BEAMS ON TWO PARAMETER FOUNDATIONS

Mohamed Taha Hassan

Cairo University, Faculty of Engineering, Giza, Egypt

e-mail: mtahah@eng.cu.edu.eg; mtaha@alfaconsult.org

Eid Hassan Doha

Cairo University, Faculty of Science, Giza, Egypt; e-mail: eiddoha@sci.cu.edu.eg

The recursivedifferentiationmethod (RDM) is introducedandemployed to obtainanalytical
solutions for static and dynamic stability parameters of beams resting on two-parameter
foundations in various different end conditions. The present analysis reflects the reliability,
efficiency and simplicity of the proposedRDM in tackling boundary value problems. In fact,
it is widely common that the critical load accompanied with the first buckling mode is the
smallest critical load, and then it is the dominant factor in the static stability analysis. In
contrast, the present analysis indicates that such a conclusion is correct only for the case
of beams without foundations or in the case of a weak foundation relative to the beam. It
is proved that critical loads accompanied with higher buckling modes may be smaller than
those accompaniedwith the lowermodes and then itmay control the stability analysis. The
same phenomenon exists for natural frequencies in the presence of an axial load. Several
illustrations are introduced to highlight the effects of both the foundation stiffness and
beam slenderness on the critical loads and natural frequencies.

Keywords: critical loads, natural frequencies, recursive differentiation method, beam on
elastic foundation

1. Introduction

Numerous analytical and numerical methods have been developed to obtain approximate so-
lutions for Boundary Value Problems (BVP). The most commonly used analytical techniques
are: Adomian Decomposition Method (ADM) proposed by Adomian (1994) and used by Taha
et al. (2012), Bahnasawi et al. (2004) and Wazwas (2001). The Variational Iteration Method
(VIM) was developed by He (2007) and used by Noor andMohyud-Din (2008). The Homotopy
Perturbation Method (HPM) used by Tan et al. and Abbasbandy (2009) and Jin (2008). The
Differential Transform (DTM) byAli (2012) and perturbation techniques byNayfeh andNayfeh
(1994) and Maccari (1999). On the other hand, numerical methods such as the Finite Element
Method (FEM) used by Mullapudi and Ayoub (2010) and Naidu and Rao (1996) and the Dif-
ferential Quadrature Method (DQM) used by Taha and Nassar (2014) and Chen (2002) offer
tractable alternative solutions for many BVPs that involve complicated mathematical formula-
tions.

Analytical methods construct the solution to BVP as a polynomial such that the coeffi-
cients of the polynomial are obtained to satisfy both the governing differential equation and the
boundary conditions. However, numerical techniques transform the differential equation into a
system of algebraic equations either on the boundary of the BVP domain or at discrete points
in the BVP domain. Indeed, the degree of success in dealing with mathematical manipulation
and accuracy determines the method efficiency.



16 M.H. Taha, E.H. Doha

In addition, techniques for finding approximate solutions for differential equations, based
on classical orthogonal polynomials, are popularly known as spectral methods. Approximating
functions in spectral methods are related to polynomial solutions of eigenvalue problems in
ordinary differential equations, known as Sturm-Liouville problems. In the last few decades,
there has been a growing interest in this subject.As amatter of fact, spectralmethods provide a
competitive alternative to other standardapproximate techniques for a large variety of problems.

Initial applications were concerned with investigations of periodic solutions to BVP using
trigonometric polynomials. Subsequently, the analysis was extended to algebraic polynomials.
The reader interested in spectral method is referred to Shen (1994, 1995), Doha and Abd-
-Elhameed (2002), Gottlieb and Orszag (1977).

Practically, in static analysis of axially loaded beams (columns), the determination of axial
loads at which the beam losses its static stability is the main issue, while in analysis of forced
vibration of beams, natural frequencies of the beam are the dominant factor in the avoidance of
the resonance phenomenon which leads to unbounded response.

In the presentwork, theRecursiveDifferentiationMethod (RDM) is introduced and employ-
ed to investigate static and dynamic behaviour of beams resting on two-parameter foundations
taking into account rotational inertia of the beam. The analysis indicates that theRDM is stra-
ightforward in its mathematical formulation and very efficient in achieving accurate solutions.
Analytic expressions for the amplitude of the lateral displacement are derived, and then the
applications of boundary conditions at the beam ends yield the corresponding characteristic
equation in two parameters (Pcr,ωn). The solution to the characteristic equation yields either
the critical loads (for ω = 0) or the natural frequencies in the case (P < Pcr). The effect of
different parameters on both the critical loads and natural frequencies will be analysed.

Thepaper is organizedas follows: inSections 2and3, theRDMis introducedandemployed to
obtain analytical solutions for static and dynamic stability parameters of beams resting on two-
-parameter foundations in various different end conditions. In Section 4, some numerical results
are discussed and a comparison with other algorithms available in the literature is presented.
Some concluding remarks are given in Section 5.

2. Recursive Differentiation Method (RDM)

In this Section, we are interested in developing the RDM to solve analytically the n-th order
differential equation

y(n)(x)=
n−1
∑

i=0

Ai,0y
(i)(x)+f0(x) x0 ¬x¬x1 (2.1)

subject to the boundary conditions

gk(x,y,y
(1), . . . ,y(n−1))= bk k=1,2, . . . ,n (2.2)

where y(i)(x) is the i-th derivative, f0(x) is the source function and Ai,0, i = 0,1,2, . . . ,n− 1
and bk are known constants.

In the RDM, we seek a solution to Eq. (2.1) subject to Eq. (2.2) in the form

y(x)=
∞
∑

m=0

Tm
(x−x0)

m

m!
(2.3)

whereTm,m=0,1, . . . are unknowncoefficients to bedetermined such that differential equation
(2.1) with its boundary conditions (2.2) is to be satisfied.



Recursive differentiation method: application to the analysis of... 17

The coefficients Tm are related to the governing differential equation on the boundary as

Tm = y
(m)
∣

∣

x=x0
(2.4)

Now, if we differentiate Eq. (2.1) once, and eliminate y(n)(x) from the resulting equation, and
after some little manipulations, we can write

y(n+1)(x)=
n−1
∑

i=0

Ai,1y
(i)(x)+f1(x) (2.5)

where

A0,1 =A0,0An−1,0 Ai,1 =Ai−1,0+Ai,0An−1,0 i=1,2, . . . ,n−1

f1(x)= f
(1)
0 (x)+f0(x)An−1,0

(2.6)

Recursive differentiations of Eq. (2.4) k-times lead to

y(n+k)(x)=
n−1
∑

i=0

Ai,ky
(i)(x)+fk(x) x0 ¬x¬x1 (2.7)

where the recurrence formulae for the coefficients Ai,k and fk(x), i = 1,2, . . . ,n − 1 and
k=1,2, . . ., may be expressed as

A0,k =A0,0An−1,k−1 Ai,k =Ai−1,k−1+Ai,0An−1,k−1

fk(x)= f
(1)
k−1(x)+An−1,k−1f0(x)

(2.8)

Making use of Eq. (2.7) enables one to get a recurrence relation for the coefficients Tn+k,
k=0,1,2 . . . in the form

Tn+k =
n−1
∑

i=0

Ai,kTi+fk(x0) (2.9)

and substitution of Eqs. (2.8) and (2.9) into Eq. (2.3) yields the solution toEq. (2.1) in the form

y(x)=
n−1
∑

j=0

TjRj(x)+Rf(x) (2.10)

where the recursive functionsRj(x) and the force recursive functionRf(x) are

Rj(x)=
(x−x0)

j

j!
+
∞
∑

i=0

Aj,i
(x−x0)

n+i

(n+ i)!
j=0,1,2, . . . ,n−1

Rf(x)=
∞
∑

i=0

fi(x0)
(x−x0)

n+i

(n+ i)!

(2.11)

In fact, the practical solution to Eq. (2.3) is truncated as

y(x)=
M
∑

m=0

Tm
(x−x0)

m

m!
(2.12)

whereM is the truncation index selected to achieve the pre-assigned degree of accuracy.
Now, the application of the boundary conditions yields the characteristic equation of the

systemwhichmay be solved to investigate the significance of different parameters on the system
behaviour. It is to be noted that the transformation of the solution domain to [0,1] has great
effect on enhancing the convergence of the obtained solutions.



18 M.H. Taha, E.H. Doha

3. Free vibration of beams on two-parameter foundation

The equations of motion of an infinitesimal element of an axially loaded beam resting on two-
-parameter foundations shown in Fig. 1, taking into consideration the rotational inertia of the
beam, are

∂V

∂x
+q(x,t)−k1y(x,t)+k2

∂2y

∂x2
= ρA

∂2y

∂t2

V (x,t)+p
∂y

∂x
−
∂M

∂x
= ρI
∂2θ

∂t2

(3.1)

while the slope-deflection and force-displacement relations are

θ=
∂y

∂x
M(x,t)=−EI

∂2y

∂x2
(3.2)

whereEI is the flexural stiffness of the beam, ρ is the density,A is the area of the cross section,
p is theaxially applied load,k1 andk2 are the linear andshear foundation stiffnessperunit length
of the beam, q(x,t) is the lateral excitation, E is the modulus of elasticity, I is the moment of
inertia, θ(x,t) is the rotation,V (x,t) is the shear force,M(x,t) is the bendingmoment, y(x,t) is
the lateral response of the beam, x is the coordinate along the beam and t is time.

Substitution of Eqs. (3.1)2 and (3.2) into Eq. (3.1)1 yields the equation of the beam lateral
response in the form

EI
∂4y

∂x4
+(p−k2)

∂2y

∂x2
−ρI

∂4y

∂x2∂t2
+k1y(x)+ρA

∂2y

∂t2
= q(x,t) (3.3)

Although the proposed RDM algorithm enables finding solutions for forced vibration of beams
with different spatial distributions of the excitation function f0(x), in the present work the
numerical analysis is limited to calculate the natural frequencies resulting from free vibration
analysis, i.e. q(x,t)= 0.

Assuming that the solution toEq. (3.3) is in the formy(x,t)=Y (x)exp(iωt) and introducing
the dimensionless variables ξ = x/L and w(x) = Y (x)/L, where ω is the natural frequency,
w(x) is the dimensionless amplitude of the lateral displacement andL is the beam length, then
Eq. (3.3) may be expressed as

d4w

dξ4
+
(

P −K2+
λ4

η2

)d2w

dξ2
+(K1−λ

4)w(ξ) = 0 (3.4)

where

K1 =
k1L
4

EI
K2 =

k2L
2

EI
P =
pL2

EI

λ4 =
ρAω2L4

EI
η=
L

r
r=

√

I

A

The parameters K1 and K2 are the foundation linear and shear stiffness parameters, P is the
axial load parameter, λ is the frequency parameter, η is the slenderness parameter and r is the
radius of gyration of the beam cross section.

3.1. Boundary conditions

For the beam shown in Fig. 1, the boundary conditions in the dimensionless form may be
expressed as:



Recursive differentiation method: application to the analysis of... 19

Fig. 1. (a) Beam on a two-parameter foundation, (b) element forces

— in the case of the pinned-pinned (P-P) beam

w(0)=w′′(0)= 0 at ξ=0 w(1)=w′′(1)= 0 at ξ=1 (3.5)

— in the case of the clamped-pinned (C-P) beam

w(0)=w′(0)= 0 at ξ=0 w(1)=w′′(1)=0 at ξ=1 (3.6)

in the case of the clamped-clamped (C-C) beam

w(0)=w′(0)= 0 at ξ=0 w(1)=w′(1)= 0 at ξ=1 (3.7)

in the case of the clamped-free (C-F) beam

w(0)=w′(0)= 0 at ξ=0 w′′(1)= 0 w′′′(1)=−Pw′(1) at ξ=1 (3.8)

3.2. Application of the Recursive Differentiation Method

To use the RDM, the governing equation of the beam-foundation system (Eq. (3.4)) is rew-
ritten in the recursive form

w(4)(ξ)=A0,0w
(0)(ξ)+A2,0w

(2)(ξ) (3.9)

where the constants Ai,0, i=0,1,2,3 are

A0,0 =−K1+λ
4 A1,0 =0 A2,0 =−P +K2−

λ4

η2
A3,0 =0 (3.10)

MakinguseofEqs. (2.7) and(2.8), the coefficientsAi,k andT4+k for i=1,2,3 and fork=1,2, . . .
can be obtained; hence the amplitude of the lateral displacement may be expressed as

w(ξ) =
3
∑

j=0

TjRj(ξ) (3.11)

where the recursive functionsRj(ξ) are obtained as

Rj(x)=
ξj

j!
+
M
∑

i=0

Aj,i
ξn+i

(n+ i)!
j=0,1,2,3 (3.12)

Substitution of Eq. (3.11) into the boundary conditions (Eqs. (3.5)-(3.8)) leads to the correspon-
ding characteristic (frequency) equations in different cases of end conditions as follows

P −P case: R10R32−R12R30 =0

C−C case: R20R31−R30R21 =0 (3.13)



20 M.H. Taha, E.H. Doha

C−P case: R20R32−R30R22 =0

C−F case: R22(R33+PR31)−R32(R23+PR21)= 0

whereRij =R
(j)
i (1), i,j=0,1,2,3.

Using aproper iterative technique, the solution of the corresponding eigenvalue problemwith
the two parameters (Pcr,ωn) can be obtained. However, for ωn = 0, the eigenvalues represent
Pcr while the eigenvectors represent the corresponding buckling modes. On the other hand,
assigning a value for P < Pcr, then the eigenvalues represent the natural frequencies ωn while
the eigenvectors represent the mode shapes of free vibration.

3.3. Verification

To verify the analytical expressions obtained from theRDM, the critical load parameter Pcr
and the natural frequency parameter λn for a beam resting on a two-parameter foundation
calculated from the RDM (Eqs. (3.13)) for the truncation index M = 25 and those obtained
from FEM (Mullapudi and Ayoub, 2010) are presented in Tables 1 and 2. In Table 1, values
of Pcr are presented for different foundation parameters and different end conditions, while in
Table 2, values of λn are presented through the effect of the loading ratio γ. It is clearly seen
that the RDM results are in close agreement to those calculated from the FEM. The critical
load parameter Pcr, the natural frequency parameter λn and the loading ratio γ are defined as

λ4n =
ρAω2nL

4

EI
Pcr =

pcrL
2

EI
γ =

P

Pcr
(3.14)

Table 1.Critical load parameter for beams on elastic foundations (η=50)

P-P beams C-C beams
K1 K2 FEM RDM FEM RDM

Pcr

0
0 9.8696 9.8696 39.479 39.4784
π2 19.739 19.7392 49.349 49.3480

100
0 20.002 20.0020 47.007 47.0066
π2 29.871 29.8713 56.877 56.876

104
0 201.41 201.4060 233.82 233.785
π2 211.28 211.2751 243.69 243.655

Table 2. Frequency parameter for beams on elastic foundations (η=50)

P-P beams C-C beams
K1 K2 γ FEM RDM FEM RDM

λn

0
0

0 3.1415 3.1416 4.7300 4.7286
0.6 2.5097 2.4984 3.7508 3.7795

π2
0 3.7306 3.7360 4.9925 4.9910
0.6 2.9842 2.9711 3.9618 3.9926

100
0

0 3.7483 3.7484 4.9504 4.9489
0.6 2.9940 2.9810 3.9323 3.9612

π2
0 4.1437 4.1437 5.1823 5.1808
0.6 3.3095 3.2954 4.1193 4.1503



Recursive differentiation method: application to the analysis of... 21

4. Numerical results

4.1. Values of the foundation stiffness parameters K1 and K2

Severalmodelshavebeenproposed to simulate the foundation reactions (Vlasov andLeontev,
1960; Zhaohua andCook, 1983). The foundation parameters (K1,K2) depend on both the beam
properties (EI,η) and the foundation linear and shear stiffness factors (k1,k2). In the present
work, k1, k2 are calculated using the expressions proposed by Zhaohua and Cook (1983) for a
rectangular beam resting on a two-parameter foundation as

k1 =
E0b

2(1−ν20)

δ

χ
k2 =

E0b

4(1−ν)

χ

δ
(4.1)

where

χ=
3

√

2EI(1−ν20)

bE0(1−ν2)
E0 =

Es
1−ν2s

ν0 =
νs
1−νs

E and ν are the elastic modulus and Poisson’s ratio of the beam, Es and νs are the respective
quantities of the foundation and δ is a parameter describing the beam-foundation loading con-
figuration (it is a common practice to assume δ=1). Typical values of the elastic modulus and
Poisson’s ratio for different types of foundations are given in Table 3.

Table 3.Typical values of modulus of elasticity and Poisson’s ratio for soils

Type of No foundation Weak foundation Medium foundation Stiff foundation
foundation (NF) (WF) (MF) (SF)

Es [N/m
2] 0 1E+07 5E+07 1E+08

Poison ratio νs 0 0.40 0.35 0.30

Though theRDM solution expressions are obtained in dimensionless forms tomake analysis
of different specific configurations possible, in the present investigation however, the properties
of the beam are (concrete beam): b = 0.2m, h = 0.5m, E = 2.1E +10Pa, Poisson’s ratio
ν =0.15. A simpleMATLAB code has been assembled to obtain the presented results.

4.2. Critical loads

It is widely known that the dominant (fundamental) critical load or natural frequency of a
beam corresponds to the first mode and is used as the limiting value for the avoidance of static
or dynamic instability. In fact, this is correct only in the cases of beams without foundations
or where the supporting foundation has a weak stiffness relative to the beam. In other words,
in the case of a beam resting on a stiff foundation, the critical loads and natural frequencies
accompanied higher modes may be smaller than those corresponding to the lower modes. This
result is crucial for stability analysis to avoid unbounded deformation in a static or dynamic
case. The variations of the critical loads and natural frequency parameters with the foundation
stiffness and slenderness ratio of the beam for different bucklingmodes and free vibrationmode
shapes are obtained and represented in Fig. 2 for the P-P case as an example. It is observed
that for a beam on an elastic foundation, as the slenderness ratio increases; i.e. the foundation
stiffness increases relative to the beam stiffness, the smallest critical load or natural frequency
parameter may correspond to one of highermodes. Thus, to avoid static or dynamic instability,
the critical loads and natural frequencies for different modes should be calculated first, and the
smaller one is then used as the limiting value for stability analysis.
Variations of the critical loads against the slenderness ratio (η) for different foundation

stiffnesses (Es,µs) and different end conditions are presented in Figs. 3 and 4.



22 M.H. Taha, E.H. Doha

Fig. 2. Critical loads and natural frequencies for differentmodes (P-P beams); (a) P
cr
(E
s
=1E+08),

(b) λ
n
(E
s
=5E+07)

Fig. 3. Variation of P
cr
with η andE

s
; (a) C-C beams, (b) C-P beams, (c) P-P beams, (d) C-F beams

Fig. 4. Effect of the beam end condition onP
cr
; (a) weak foundation (WF), (b) stiff foundation (SF)

In Fig. 3, case (a) represents the clamped-clamped beams (C-C), case (b) clamped-pinned
beams (C-P), case (c) pinned-pinned beams (P-P) and case (d) clamped-free (cantilever) beams
(C-F). The influence of the foundation increases as the slenderness ratio increases. The influence
of end conditions decreases as the slenderness ratio increases.



Recursive differentiation method: application to the analysis of... 23

The transition between different modes is clear in (P-P) and (C-C) cases for a medium
stiffness (MF) and stiff foundations (SF). Actually, the assembled MATLAB code always picks
the smallest critical load in spite of the buckling mode.

4.3. Natural frequencies

The variations of the natural frequency parameter λn with the foundation stiffness (Es,µs)
and slenderness ratio (η) are presented inFigs. 5 and 6 for γ=0,while the effect of γ are shown
in Figs. 7-10.

Fig. 5. Varaiations of λ
n
with η andE

s
; (a) C-C beams, (b) C-P beams, (c) P-P beams, (d) C-F beams

Fig. 6. Effect of the beam end condition on λ
n
; (a) weak foundation (WF), (b) stiff foundation (SF)

The influence of the slenderness ratio increases with the increase of the foundation stiffness
and vice versa. The effect of end conditions vanishes with the increase in the slenderness ratio
for a weak foundation faster than in the case of stiff foundation (SF). The transition between
different vibration modes is not obvious in the frequency charts, as in the absence of the axial
load, the natural frequency of the first vibration mode is always the smaller one. In the case of
axially loaded beams, the transition between the first and second mode is detected for slender
beams with (P-P) and (C-C) end conditions.



24 M.H. Taha, E.H. Doha

Fig. 7. Variations of λ
n
with γ for C-C beams; (a) effect ofE

s
, (b) effect of η

Fig. 8. Variations of λ
n
with γ for P-C bemas; (a) effect ofE

s
, (b) effect of η

Fig. 9. Variations of λ
n
with γ for P-P beams; (a) effect ofE

s
, (b) effect of η

Fig. 10. Variations of λ
n
with γ for C-F beams; (a) effect ofE

s
, (b) effect of η

In the case of axially loaded beams (Figs. 7-10), the effect of the loading ratio γ may be
approximated by a linear relation up to γ = 0.5. Also, the influences of both the slenderness
ratio and the foundation stiffness are more noticeable in the (C-F) case.



Recursive differentiation method: application to the analysis of... 25

Furthermore, it is found that taking the rotational inertia of the beam into consideration
decreases thenatural frequencies of shortbeams, and the effectmaybe ignoredas the slenderness
ratio η> 30.

5. Conclusions

The RDM is introduced and employed in the investigation of the static and dynamic stability
parameters of axially loadedbeams resting on two-parameter foundations.Recursive functions of
theproblemarederivedfirst, and thenafter applying the endconditions, the frequency equations
accompanied with different end conditions are obtained. However, it is found that the accuracy
of the obtained RDM expressions is greatly enhanced when the solution domain is transformed
to the domain [0,1].
The critical loads required in static stability analysis and natural frequencies required in

dynamic stability analysis are obtained and investigated.
It is observed that in the case of beams resting on elastic foundations, the critical load of the

first buckling mode is not always the smallest critical load in contrast to that common fact in
the case of beams without foundation. The critical load of a higher mode may be smaller than
the critical load of a lower buckling mode. This phenomenon is also observed for the natural
frequency, but in the presence of an axial load.
It is also concluded that both the influence of the foundation stiffness and the slenderness

ratio are more noticeable for the (C-F) case.
Although the proposedRDMsolution is applicable for forcedvibration, thenumerical results

are limited to free vibration to calculate the natural frequencies which are required in stability
analysis.

References

1. Adomian G., 1994, Solving Frontier Problems of Physics: the Decomposition Method, Vol. 60 of
the Fundamental Theories of Physics, Kluwer Academic Publishers, Boston

2. Ali J., 2012, One dimensionless differential transform method for some higher order boundary
value problems in finite domain, International Journal of ContemporaryMathematical Sciences, 7,
6, 263-272

3. Bahnasawi A.A., El-TawilM.A., Abdel-NabyA., 2004, SolvingRiccati differential equation
using Adomian’s decompositionmethod,Applied Mathematics and Computation, 157, 503-514

4. Chen C.N., 2002, DQEMVibration analysis of non-prismatic shear deformable beams resting on
elastic foundations, Journal of Sound and Vibration, 255, 5, 989-999

5. DohaE.H.,Abd-ElhameedW.M., 2002,Efficient spectral-Galerkinalgorithmfordirect solution
of second-order equationsusingultrasphericalpolynomials,SIAMJournal on ScientificComputing,
24, 548-571

6. Doha E.H., Abd-Elhameed W.M., Bassuony M.A., 2013, New algorithms for solving high-
even-order differential equations using third and fourth Chebyshev-Galrkin methods, Journal of
Computational Physics, 236, 563-579

7. Gottlieb D., Orszag A., 1977, Numerical Analysis of Spectral Methods: Theory and Applica-
tions, SIAM, Philadelphia, Pennsylvania

8. HeJ.-H., 2007,Variational iterationmethod: Some recent results andnew interpretations,Journal
of Computational and Applied Mathematics, 207, 3-17

9. Jin L., 2008,Homotopyperturbationmethod for solvingpartial differential equationswith variable
coefficients, International Journal of Contemporary Mathematical Sciences, 3, 28, 1395-1407



26 M.H. Taha, E.H. Doha

10. Maccari A., 1999, The asymptotic perturbationmethod for nonlinear continuous systems,Non-
linear Dynamics, 19, 1-18

11. Mullapudi R., Ayoub A., 2010, Nonlinear finite element modeling of beams on two-parameter
foundations,Computers and Geotechnics, 37, 334-342

12. Naidu N.R., Rao G.V., 1996,Vibrations of initially stressed uniform beams on a two-parameter
elastic foundation,Computers and Structures, 57, 5, 941-943

13. Noor M.A., Mohyud-Din S.T., 2008,Modified variational iterationmethod for heat andwave-
like equations,Acta Applicandae Mathematicae, 104, 3, 257-269

14. Nayfeh A.H., Nayfeh S.A., 1994,On nonlinearmodes of continuous systems, Journal of Vibra-
tion and Acoustics, 116, 129-136

15. Shen J., 1994, Efficient spectral-Galerkin methods. I: Direct solvers of second and forth-order
equations using Legendre polynomials, SIAM Journal on Scientific Computing, 15, 1489-1505

16. Shen J., 1995, Efficient spectral-Galerkin methods. II: Direct solvers of second and forth-order
equations using Chepyshev polynomials, SIAM Journal on Scientific Computing, 16, 74-87

17. Taha M.H., NassarM., 2014,Analysis of axially loaded tapered beamswith general end restra-
ints on two parameter foundation, Journal of Theoretical And Applied Mechanics, 52, 1, 215-225

18. Taha M.H., Omar A., Nassar M., 2012, Dynamics of Timoshenko beam on nonlinear soil,
International Journal of Civil Engineering, 3, 2, 93-103

19. Tan Y., Abbasbandy S., 2008, Homotopy analysis method for quadratic Riccati differential
equation,Communications in Nonlinear Science and Numerical Simulation, 13, 3, 539-546

20. Vlasov V.Z., Leontev U.N., 1960, Beams, Plates And Shells On Elastic Foundations, Gos.
Izdat. Fiz.-Math. Lit., Moskva

21. Wazwas A.M., 2001, The numerical solution of fifth-order boundary value problems by decom-
position method, Journal of Computational and Applied Mathematics, 136, 259-270

22. Zhaohua F., Cook R.D., 1983, Beam elements on two parameter elastic foundation, Journal of
Engineering, ASCE, 109, 6, 1390-1402

Manuscript received March 14, 2014; accepted for print June 17, 2014