Application of the Homotopy Analysis Method for Determining the free Vibrations of Beam

Vibrations of the building structures can be divided into the forced ones and the free ones. The forced vibrations appear when the variable load acts on the structure. The free vibrations take place when the variable load does not act on the structure during the vibrations and the motion of the structure results from the initial conditions [14].

The vibrations registered while testing the real structures can be presented with the aid of continuous function. Such function is the polynomial of degree n inscribed into the quantized signal obtained during the research. By analyzing these polynomials one may determine the period of vibration and the damping parameters. In case of structures characterized by the constant stiffness and the damping force proportional to the vibration velocity, the period of vibration and the damping parameters are constant. If the previously mentioned conditions are not satisfied, we deal with the nonlinear vibrations [6].

Mathematical modeling of the building structures consists in deriving the differential equations including the linear and nonlinear parameters of the construction. With respect to the form of describing the vibrations of the real structures it is essential to get the solution of a mathematical model also in the form of continuous function. This function can be obtained by applying the homotopy analysis method. The homotopy analysis method has been developed in the 90s of the last century [15-19, 23]. This method gives the possibility to search for the solution of a wide range of problems described by means of the operator equations. In particular, it can be the ordinary and partial differential equations [3, 4, 9, 10, 18, 23, 28], fractional differential or integrodifferential equations [30, 31] (for some other applications of the fractional calculus see for example [11, 21, 24]) and also the integral equations [1, 2, 7, 8, 12, 13, 23]. Some theoretical results concerning the convergence of this method in case of differential equations can be found, among others, in papers [9, 16, 19, 20, 23, 25, 29], whereas in case of integral equations in papers [2, 7, 8, 23].

In the current paper we describe the application of the homotopy analysis method for determining the free vibrations of the beam in the linear and nonlinear case. The numerical examples are presented to confirm the exactness and fast convergence of the introduced method. In the linear case we know the exact solution, so we compare with it the solution obtained with the aid of homotopy analysis method. This case serves then to verify the exactness and convergence of the discussed method. In the nonlinear case the exact solution is unknown, therefore the approximate solution received by using the homotopy analysis method is compared with the approximate solution determined numerically by applying the Mathematica software [5].

We describe the concept of the homotopy analysis method by the example of its application for solving the ordinary differential equations of the second kind in the form (1) d 2 y ( t ) d t 2 + N 1 ( t , y ) + N 2 ( t , y ) = 0 , t ∈ ( a , b ) ,

with the initial conditions (2) y ( a ) = y a , y ′ ( a ) = v a ,

where N ₁ is the linear operator, N ₂ is the non-linear operator, y_a and v_a are the given numbers, while y is the unknown function. Let us denote (3) N ( y ) = d 2 y ( t ) d t 2 + N 1 ( t , y ) + N 2 ( t , y ) .

In the first step, we specify the homotopy operator ℋ as: (4) ℋ ( Φ , p ) ≡ ( 1 − p ) L ( Φ ( t ; p ) − y 0 ( t ) ) = p h N ( Φ ( t ; p ) ) ,

where p∈[0,1] is the embedding parameter, h≠ 0 denotes the convergence control parameter [17, 19, 22, 29], y _଴ means the initial approximation of the solution of problem 1 and L is the auxiliary linear operator with property L(0).

Considering the equation ℋ(Φp=0), we obtain the so-called zero-order deformation equation (5) ( 1 − p ) L ( Φ ( t ; p ) − y 0 ( t ) ) = p h N ( Φ ( t ; p ) ) .

Substituting p=0 we get, so L(Φ(t;0) – y ₀(t))=0. However, if we assume p=0, we obtain N(Φ(t; 1)) = 0. Therefore Φ(t; 1) = y(t), where y is the desired solution of equation 1. Thus, the change of parameter p from zero to one corresponds with the change from the trivial problem to the original problem (and thus the solution from y ₀ to y).

Taking the Maclaurin series of function Φ (with respect to parameter p) we get (6) ( Φ ( t ; p ) = ∑ m = 0 ∞ y m ( t ) p m ,

where y ₀(t) = Φ(t; 0) and (7) y m ( t ) = 1 m ! ∂ m Φ ( t ; p ) ∂ p m | p = 0 , m = 1 , 2 , 3 , …

If the above series is convergent for p = 1, then we obtain the required solution (8) y ( t ) = ∑ m = 0 ∞ y m ( t ) .

If we are not able to determine the sum of series in 8, then as an approximate solution of considered equation we can accept the partial sum of this series (9) y ^ n ( t ) = ∑ m = 0 n y m ( t ) .

In order to determine the form of function y_m we differentiate m times, with respect to parameter p, the left and the right hand side of relation 5. Next, we divide the result by m! and we substitute p = 0 which gives the so-called m th order deformation equation (m > 0): (10) L ( y m ( t ) − χ m y m − 1 ( t ) ) = h R ¯ m ( y ¯ m − 1 , t ) ,

where y ¯ m − 1 = { y 0 ( t ) , y 1 ( t ) , … , y m − 1 ( t ) } and (11) χ m = { 1 m > 1 0 m ≤ 1,

and (12) R ¯ m ( y ¯ m − 1 , t ) = 1 ( m − 1 ) ! ( ∂ m − 1 ∂ p m − 1 N ( ∑ i = 0 ∞ y i ( t ) p i ) ) | p = 0 .

In view of the above definition of operator N (equation 3) we have for m= 1: R ¯ 1 ( y ¯ 0 , t ) = N ( y 0 ( t ) ) = d 2 y 0 ( t ) d t 2 + N 1 ( t , y 0 ) + N 2 ( t , y 0 ) ,

whereas for m ≥ we get R ¯ m ( y ¯ m − 1 , t ) = d 2 y m − 1 ( t ) d t 2 + N 1 ( t , y m − 1 ) + 1 ( m − 1 ) ! ( ∂ m − 1 ∂ p m − 1 N 2 ( t , ∑ i = 0 ∞ y i ( t ) p i ) ) | p = 0 .

Thus, the equation 10 for m= 1 takes the form (13) L ( y 1 ( t ) ) = h ( d 2 y 0 ( t ) d t 2 + N 1 ( t , y 0 ) + N 2 ( t , y 0 ) ) ,

whereas for m ≥ 2 it takes the form L ( y m ( t ) ) = L ( y m − 1 ( t ) ) + h ( d 2 y m − 1 ( t ) d t 2 + N 1 ( t , y m − 1 ) + 1 ( m − 1 ) ! ( ∂ m − 1 ∂ p m − 1 N 2 ( t , ∑ i = 0 ∞ y i ( t ) p i ) ) | p = 0 ) .

After selecting operator L the above two equations enable to determine recursively the successive functions y_m and, in result, to obtain the solution 8.

As the operator L one can choose the operator of the highest order derivative occurring in the equation. It is the most natural and very often used procedure. This operator can also be selected in other way. In particular, Vajravelu and Van Gorder [27, 29] proposed three methods of selecting the auxiliary linear operator: method of linear partition matching, method of highest order differential matching, method of complete differential matching (see also [23]).

In the considered case we receive the simple form of equations by taking (15) L ( y ) = d 2 y d t 2 .

Then we get for m=1: (16) d 2 y 1 ( t ) d t 2 = h ( d 2 y 0 ( t ) d t 2 + N 1 ( t , y 0 ) + N 2 ( t , y 0 ) ) ,

and for m ≥ 2: (17) d 2 y m ( t ) d t 2 = d 2 y m − 1 ( t ) d t 2 + h ( d 2 y m − 1 ( t ) d t 2 + N 1 ( t , y m − 1 ) + 1 ( m − 1 ) ! ( ∂ m − 1 ∂ p m − 1 N 2 ( t , ∑ i = 0 ∞ y i ( t ) p i ) ) | p = 0 ) .

The last component in the above equation is nonlinear. We will be able to determine this nonlinear element by knowing the form of nonlinear operator N ₂.

As the initial approximation the best is to select a function satisfying the given initial conditions 2, for example in the form (18) y 0 ( t ) = y a + v a ( t − a ) .

But if one makes another choice of the linear operator (other than in the form of second order derivative), then one can connect the initial approximation with the selected linear operator and the given initial conditions. On can then determine y ₀ as the solution of linear equation L(y ₀) = 0 with the same initial conditions as in the original equation.

In order to ensure the uniqueness of solution of equations 16 and 17 we need to define the corresponding initial conditions. If we put to each equation the homogeneous initial conditions (19) y m ( a ) = 0 , y m ′ ( a ) = 0

for m ≥ 1, then the approximate solution ŷ_n ,defined by relation 9, will satisfy the initial conditions 2.

In this way the solution of considered problem is reduced to the solution of a sequence of differential equations 16 and 17 with the homogeneous initial conditions 19. In these equations the sought functions appear only once within the symbol of second derivative. Thus the solution of these equations consist in the double integration of their right hand sides.

The proper selection of the convergence control parameter h affects the area of convergence of the created series and the rate of this convergence as well [19, 20, 26]. For determining this coefficient we can use the so-called “optimization method” [19], in which we define the following squared residual of governing equation (20) E n ( h ) = ∫ a b ( N [ y ^ n ( t ) ] ) 2 d t .

And next, the optimum value of the convergence control parameter is obtained by determining the minimum of this squared residual.

Vibrations of the simply supported reinforced concrete beam of span l= 6 m and cross-section (b/h) 0.30 m/0.60 m can be described with the aid of a linear differential equation. Uniformly distributed mass m′ of the beam, under the assumption that the density of the beam material is ρ = 2500 kg/m³, is equal to (21) m ′ = b h ρ = 450 kg/m .

Such beam can be considered as the set of segments of length dx and the kinetic energy of such elementary segment, treated as the element of mass m′ dx (Figure 1a), is equal to (22) m ′ d x y ˙ 2 ( x ) 2 ,

Figure 1.

where ẏ(x) is the velocity of point of coordinate x (in the beam). In case of the simply supported beam wrapped by the mass m′ the deflection line y(x) is symmetric with respect to the center of gravity and for 0 ≤ x ≤ 1 2 l is expressed by relation (23) y ( x ) = 16 5 y m ( x l − 2 x 3 l 3 + x 4 l 4 ) ,

where y_m is the largest deflection of the beam (Figure 1a).

Velocity of every point of coordinate x is equal to (24) y ˙ ( x ) = 16 5 y ˙ m ( x l − 2 x 3 l 3 + x 4 l 4 ) ,

where y ˙ m = d y m d t and the kinetic energy of the whole beam is equal to (25) E k = 2 ∫ 0 l / 2 m ′ y ˙ 2 ( x ) 2 d x .

By substituting 24 into 25 and by calculating the integrals we get (26) E k = 3968 7875 m ′ l y ˙ m 2 2 = 0.504 m ′ l y ˙ m 2 2 .

The kinetic energy described by equation 26 corresponds to the kinetic energy of the system with one degree of freedom and the substitute mass equal to (27) m = 0.504 m ′ l = 1360.8 kg

focused in the deflection, that is in the point in which the deflection of a beam is the largest (Figure 1b). Stiffness of the substitute system is equal to (28) k = 48 E I l 3 ,

where E= 27 GPa is the modulus of elasticity of the beam material and I is the moment of inertia of the cross-section which is equal to b h 3 12 = 5.4 · 10 − 3 m 4 Thus we have k= 32.4 MN/m.

Equation describing the first free vibration of the reinforced concrete beam of span l= 6 m and cross-section (b/h) 0.30 m/0.60 m is then of the form (29) m y ¨ ( t ) + k y ( t ) = 0,

where m= 1360.8 kg denotes the substitute mass of the beam and k= 32.4 MN/m means the stiffness of the beam in the middle of its span.

Equation 29, after dividing by m, takes the form (30) y ¨ ( t ) + ω 0 2 y ( t ) = 0,

where ω 0 = k m = 140.8573 rad/s is the natural frequency of non-damping vibrations. The above equation must be completed with the initial conditions (31) y ( 0 ) = y a , y ˙ ( 0 ) = v a .

Example 1.

We show now how the described above homotopy analysis method may be used for solving the equation 30 with conditions 31. For the value ω ₀ = 140.9 rad/s and the initial conditions y_a =0.01 m and v_a =0 m/s, the exact solution is given by function y e ( t ) = 0.01 cos ( 140.9 t ) .

If we take as the initial approximation the function satisfying the initial conditions y 0 ( t ) = y a + v a t = 0.01 ,

then in the first steps of the method we get y 1 ( t ) = 99.2641 h t 2 , y 2 ( t ) = 99.2641 ( h t 2 + h 2 t 2 ) + 164223 ⋮ h 2 t 4 ,

Figure 3 presents the squared residual for n= 9. Numerically determined optimum value of the convergence control parameter is equal to −0.950692782. In Table 1 there are collected the errors in reconstructing the exact solution for the successive approximate solutions ŷ_n , n∈ {1,2, …, 15}. The presented errors are calculated in the supremum norm Λ n = ‖ y e − y ^ n ‖ = sup | y e ( t ) − y ^ n ( t ) | , t ∈ ( a , b )

Figure 2.

Figure 3.

Table 1.

Values of errors in the reconstruction of the exact solution (Δ _n = ||y_e – ŷ_n ||)

n	1	2	3	4	5
Δ n	0.187	0.390	0.287	0.101	1.767·10-2
n	6	7	8	9	10
Δ n	1.187·10-3	4.884·10-5	7.253·10-6	3.744·10-7	2.572·10-8
n	11	12	13	14	15
Δ n	3.130·10-9	8.340·10-11	1.165·1-11	1.101·10-12	3.706·10-14

where y_e means the exact solution and ŷ_n denotes the approximate solution. At first, for n= 2 the errors increase but afterwards they quickly decrease. So, for n= 5 the errors are about 10^–2, for n= 10 they are about 10^–8 and for n= 15 the errors are at the level of 10^–14. Next the error decreases a little bit slower, so for n= 20 the error is equal to 3.747 · 10^–16. In Figure 4 the absolute error |y_e (t) – ŷ_n (t)| for n= 6 and n= 10 is plotted. As revealed by the above results, starting from some moment the errors quickly decrease when the number of components in the sum 9 increases.

Figure 4.

Example 2.

In the second example we execute the calculations for the changed initial conditions by taking y_a = 0.01 m and v_a =0.2 m/s. Value of parameter ω ₀ is not changed (ω ₀=140.9 rad/s). In this case we also know the exact solution given by the function y e ( t ) = 0.01 cos ( 140.9 t ) + 0.00141945 sin ( 140.9 t ) .

As the initial approximation we take the function fulfilling the initial conditions y 0 ( t ) = y a + v a t = 0.01 + 0.2 t .

Then we obtain successively y 1 ( t ) = 99.2641 h t 2 + 661.76 h t 3 , y 2 ( t ) = 99.2641 ( h t 2 + h 2 t 2 ) + 661.76 ( h t 3 + h 2 t 2 ) + 164223 h 2 t 4 + 656890 h 2 t 5 , ⋮

Figure 5 displays the squared residual for n = 9.

Figure 5.

Numerically computed optimum value of the convergence control parameter is equal to –0.9580711. In Table 2 there are compiled the errors in reconstructing the exact solution for the successive approximate solutions ŷ_n , n∈ {1,2, …, 15}. Similarly as in the previous example, at first the errors increase but afterwards they quickly decrease. In Figure 6 there is plotted the absolute error of reconstructing the exact solution for n= 10 and n= 15.

Table 2.

Values of errors in the reconstruction of the exact solution (Δ _n = ||y_e – ŷ_n ||)

n	1	2	3	4	5
Δ n	0.236	0.456	0.333	0.122	2.388·10-2
n	6	7	8	9	10
Δ n	2.283·10-3	4.017·10-5	8.665·10-6	8.046·10-8	3.744·10-8
n	11	12	13	14	15
Δ n	7.420·10-10	1.238·10-10	7.527·10-12	1.437·10-13	2.802·10-14

Figure 6.

The initial condition (32) y ( 0 ) = 0 , y ˙ ( 0 ) = v a

of the free vibrations can be achieved by dropping on the beam from height h the body of mass m ₁ (the beater – Figure 7a). When the impact is of plastic character then the structure and the beater would perform the free vibrations (Figure 7b). Assuming m ₁ = 100 kg and h= 0.5 m velocity of the beater is equal to (33) v 1 = 2 g h = 3.13 m/s,

Figure 7.

where g is the acceleration of gravity. Basing on the momentum conservation principle the velocity v_a of the system after the strike is equal to (34) v a = m 1 v 1 m 1 + m = 0.181 m/s .

After the impact of the plastic character the motion of the oscillator is executed according to the equation (35) ( m 1 + m ) y ¨ ( t ) + k y ( t ) = m 1 g ,

the solution of which, after taking into account the initial conditions (36) y ( 0 ) = 0 , y ˙ ( 0 ) = v a ,

is given by the function (37) y ( t ) = v a ω 01 sin ( ω 01 t ) + m 1 g k ( 1 − cos ( ω 01 t ) ) ,

where ω 01 = k m 1 + m = 136.73 rad/s .

Example 3.

Now we deal with the solution of equation 35 with conditions 36, for the following parameters m ₁ = 100 kg, m = 1633 kg, k= 32.4 MN/m, g=9.81 m/s ², and the initial conditions y_a =0 m and v_a =0.181 m/s. This time by using our method we get in turn y 1 ( t ) = − 0.283035 h t 2 + 563.159 h t 3 , y 2 ( t ) = 0.283035 ( h t 2 − h 2 t 2 ) + 563.159 ( h t 3 + h 2 t 3 ) − 440.967 h 2 t 4 + 526438 h 2 t 5 .

Numerically determined optimum value of the convergence control parameter is equal to –0.77670315. Table 3 compiles the errors in reconstructing the exact solution for the successive approximate solutions ŷ_n , n ∈ {1,2, …, 15}. In Figure 8 the absolute error of reconstructing the exact solution for n= 10 and n= 15 is presented. Similarly as before the errors quickly decrease.

Table 3.

Values of errors in the reconstruction of the exact solution (Δ _n = ||y_e – ŷ_n ||)

n	1	2	3	4	5
Δ n	3.367·10-2	2.084·10-2	1.741·10-3	1.374·10-3	4.188·10-4
n	6	7	8	9	10
Δ n	5.135·10-5	1.785·10-5	9.227·10-5	2.483·10-6	4.230·10-7
n	11	12	13	14	15
Δ n	6.678·10-8	2.612·10-8	1.181·10-8	3.571·10-9	8.398·10-10

Figure 8.

Nonlinearity in the mechanical dynamic system can occur in result of nonlinearity of the elastic force S, the value of which does not depend linearly on the displacement y. Characteristic of this force can be presented in the form of the well known in literature Duffing’s proposal (38) s ( y ) = k y ± k ^ y 3 .

The plus sign preceding k ^ means that the stiffness increases together with the displacement y, whereas the minus sign means that it decreases. Nonlinearity is the greater, the greater is the displacement. Equation describing such nonlinear free vibrations with no damping takes the form (39) m y ¨ ( t ) + k y ( t ) ± k ^ y 3 ( t ) = 0 ,

with the initial conditions (40) y ( 0 ) = 0 , y ˙ ( 0 ) = v a .

Referring to the previous example when we assumed the mass m= 1633 kg and stiffness k= 32.4 MN/m, the parameter k ^ can take the values about 10⁹ with the amplitude of variable y up to 0.002 m. Certainly the value of k ^ can be increased and then the nonlinearity effect would be greater, however, for the bigger amplitudes of y, the problem loses its physical sense.

The free vibrations with the damping include additionally one more term by the first derivative denoted usually as c: (41) m y ¨ ( t ) + c y ˙ + k y ( t ) ± k ^ y 3 ( t ) = 0.

In our case the physically real value of c (determined from the engineering relations, ζ= 0.05 – damping ratio) is equal to (42) c = 2 ζ k m = 23 ⋅ 10 3 kg/s .

Example 4.

Proceeding with the computational example we discuss the solution of equation 39 with initial conditions 40 for m = 1633 kg, k= 32.4 MN/m, k ^ = 10 9 and v_a =0.2 m/s, under the assumption that together with the increase of variable y the stiffness decreases (therefore we use the minus sign in the proper equation). This time we do not know the exact solution of considered equation. That is way the approximate solution obtained by applying the homotopy analysis method will be compared with the approximate solution received numerically with the aid of Mathematica software. Assuming as the initial approximation y ₀ the function satisfying the initial conditions y 0 ( t ) = v a t = 0.2 t ,

we get, by using the homotopy analysis method, the function y 1 ( t ) = 661.359 h t 3 − 244.948 h t 5 , y 2 ( t ) = 661.359 ( h t 3 + h 2 t 3 ) − 244.948 h t 5 + 655850 h 2 t 5 − 1.27285 ⋅ 10 6 h 2 t 7 + 249998 h 2 t 9 , ⋮

The optimum value of the convergence control parameter is equal to –0.7227439. The plot of the squared residual for n= 20 is presented in Figure 9. Figure 10 shows the plot of residual R s n ( t ) = | 1633 y ^ n ( t ) + 32.4 ⋅ 10 6 y ^ n ( t ) − 10 9 ( y ^ n ( t ) ) 3 |

Figure 9.

Figure 10.

for n= 10 and n= 20, that is the plot of error of satisfying the equation 39 by the approximate solution ŷ_n . The above error decreases very quickly when the number of elements in the approximate solution increases (with the increase of number n). So, for n= 5 the maximum of this error takes the value 67.72, for n= 10 this maximum does not exceed the value 0.45, for n= 15 the value 0.081, and finally for n= 20 it does not exceed the value 0.0003. The initial conditions are fulfilled exactly by every approximate solution.

Figure 11 presents the comparison of approximate solutions obtained by applying the homotopy analysis method (solid line) and the numerical methods available in the Mathematica software (dots). The differences are slight even for the small number n of elements (see Table 4).

Figure 11.

Table 4.

Maximal absolute differences (Δ _n ) between the approximate solutions obtained by the homotopy analysis method and the numerical methods available in Mathematica software

n	1	3	5
Δ n	3.35·10-2	4.98·10-3	7.12·10-5
n	10	15	20
Δ n	1.40·10-6	2.99·10-8	3.71·10-9

Example 5.

In the next example we solve the equation 41 with initial conditions 40 for m = 1633 kg, k = 32.4 MN/m, k ^ = 10 9 , c= 23·10³ kg/s and v_a = 0.2 m/s, under the assumption that when the coordinate y increases then the stiffness decreases (so we include the minus sign in the proper equation). This time again we do not know the exact solution of investigated equation. Taking, similarly as before, the function satisfying the initial conditions as the initial approximation y ₀, that is y 0 ( t ) = v a t = 0.2 t ,

we get successively the functions y 1 ( t ) = 1.40857 h t 2 + 661.359 h t 3 − 244.948 h t 5 , y 2 ( t ) = 1.40857 ( h t 2 + h 2 t 2 ) + 661.359 h t 3 + 667.973 h 2 t 3 + 4657.87 h 2 t 4 − 244.948 h t 5 + 655850 h 2 t 5 − 4025.32 h 2 t 6 − 1.27285 ⋅ 10 6 h 2 t 7 + 249998 h 2 t 9 , ⋮

Optimum value of the convergence control parameter is equal to –0.686325. The plot of the squared residual for n= 20 is shown in Figure 12. In Figure 13 there is displayed the plot of residual R ^ s n ( t ) = | 1633 y ^ n ( t ) + 23 ⋅ 10 3 y ^ n ( t ) + 32.4 ⋅ 10 6 y ^ n ( t ) − 10 9 ( y ^ n ( t ) ) 3 |

Figure 12.

Figure 13.

for n= 10 and n= 20, it means the plot of error of satisfying the equation 41 by the approximate solution ŷ_n . Similarly as in the previous example, the above error decreases very quickly when the number of elements in the approximate solution increases. For n= 5 the maximum of this error takes the value 66.37 for n= 10 this maximum does not exceed the value 0.78, for n= 15 the value 0.12 and, in turn, for n= 20 it does not exceed the value 0.00052. The initial conditions are satisfied exactly by every approximate solution.

Figure 14 shows the comparison of approximate solutions obtained by applying the homotopy analysis method (solid line) and the numerical methods available in the Mathematica software (dots). The differences are slight even for small number n of components. For n= 5 the absolute value of the difference does not exceed the value 2.56 · 10^–4, for n= 10 it does not exceed 3.12 · 10^–6, for n= 15 the value 3.61 · 10^–8, whereas for n= 20 it is not higher than 4.2 · 10^–9.

Figure 14.

In the paper we have presented the application of the homotopy analysis method for determining the free vibrations of the reinforced concrete beam in the linear and nonlinear case. In the investigated method the series is created, the successive terms of which are derived by calculating the proper integrals obtained from the previous terms. The formed series is in general quickly convergent. For the optimal selection of the convergence control parameter computation of just several first terms of the series ensures a very good approximation of the sought solution. Performed calculations show that the proposed method is effective in solving the problems under consideration.

The presented computational examples confirm the precision and the fast convergence of investigated method. In the linear case we knew the exact solutions, so we could compare with them the solutions obtained with the aid of homotopy analysis method. This case served then to verify the exactness and convergence of the discussed method. In the nonlinear case, when the exact solution was unknown, we compared the approximate solutions received by using the homotopy analysis method with the approximate solutions determined numerically by applying the Mathematica software. Differences between the obtained solutions were slight. However, the advantage of the examined method is that we receive here the approximate solution in the form of continuous function which can be used then in a further analysis or to perform various simulations. Whereas when we apply the numerical methods, we get the discrete set of values which may be potentially used later for approximating the appropriate functions.

The executed calculations indicate that it is possible to obtain the mathematical solution of the model describing the free vibrations of the reinforced concrete beam in the form of polynomial. The obtained results can be compared with the results of experimental research also approximated by polynomials. Such comparison will give the possibility to verify the taken model of structure.