Jtam.dvi JOURNAL OF THEORETICAL AND APPLIED MECHANICS 46, 2, pp. 347-366, Warsaw 2008 THE METHOD OF SOLVING POLYNOMIALS IN THE BEAM VIBRATION PROBLEM Mohammad Jehad Al-Khatib Krzysztof Grysa Artur Maciąg Kielce University of Technology, Department of Mathematics, Poland e-mail: matam@tu.kielce.pl; grysa@tu.kielce.pl The paper presents a new method of finding an approximate solution to the beam vibration problem. The problem is described with a par- tial differential equation of the fourth order. Basically, linear differential equations can be solved by means of various methods. The key idea of the presented approach is to find polynomials (solving functions) that satisfy the considered differential equation identically. In this sense, it is a variant of the Trefftz method. The advantage of the method is that the approximate solution (a linear combination of the solving functions) satisfies the equation identically. The initial and boundary conditions are then satisfied approximately. The formulas for solving functions and their derivatives are obtained. The solving Trefftz functions can be used in the whole domain or can be used as base functions in nodeless FEM. Both cases are considered. A numerical example is included. Key words: beam vibration, Trefftz method, solving functions 1. Introduction Although the method of solving functions for linear partial differential equ- ations was been developedmainly during the last 10 years, it was thoroughly addressed in the literature. The method was first described by Rosenbloom and Widder (1956) where it was applied to one-dimensional heat conduction problems. In the papers by Ciałkowski et al. (1999a,b), Futakiewicz (1999), Futakiewicz and Hożejowski (1998a,b), Futakiewicz et al. (1999), Hożewski (1999), Yano et al. (1983), the authors described heat functions in different 348 M.J. Al-Khatib et al. coordinate systems for direct and inverse heat conduction problems. Ciałkow- ski (1999) considered a highly interesting idea of using heat polynomials as a new type of finite-element base functions. In another paper Ciałkowski and Frąckowiak (2000) dealt with numero- us cases, involving other differential equations such as the Laplace, Poisson, Helmholtz ones. Solving functions for thewave equation are presented inCiał- kowski (2003), Ciałkowski and Frąckowiak (2000, 2003, 2004), Ciałkowski and Jarosławski ((2003), Maciąg (2004, 2005), Maciąg and Wauer (2005a,b). Ap- plications of wave polynomials for elasticity problems are shown in Maciąg (2007b). In Section 2, solving polynomials and their partial derivatives are consi- dered. Section 3 describes the method of solving functions. An example is discussed in Section 4. Concluding remarks are presented in Section 5. 2. Solving polynomials Let us consider the equation that describes beam vibration ∂4u ∂x4 + 1 a2 ∂2u ∂t2 =0 (2.1) where a2 = EJ/(ρS); E stands for the coefficient of elasticity, J describes the inertia moment of the beam cross section with S being the surface of the beam cross-section. Equation (2.1) describes also vibration of the tuning fork and some problems concerning rotating shafts. In dimensionless coordinates, equation (2.1) has a form ∂4u ∂x4 + ∂2u ∂t2 =0 x∈ (0,1), t > 0 (2.2) Equation (2.1) or (2.2) has to be completed with initial and boundary condi- tions depending on the type of fixation of the beam. There are two ways to obtain solving polynomials for equation (2.2). The first one is to expand the function satisfying the equation in power series and reducing some components of the expansion with the use of this equation. The second method consists in using a ”generating function”. Both ways are equivalent and lead to the same polynomials. The method of solving polynomials... 349 2.1. Expansion of the function satisfying the governing equation into Taylor series Thefirstway to get solvingpolynomials for (2.2) is to expandthe function that satisfies this equation into Taylor series (comp. Ciałkowski and Frąckowiak, 2000). Let the function u(x,t) satisfy equation (2.2) with given initial and boundary conditions. We assume u ∈ CN+1 in the neighborhood of (0,0). Then, the Taylor series for the function u is u(x,t) =u(0,0)+ du 1! + d2u 2! + . . .+ dn−1u (n−1)! +Rn (2.3) where dνu= (∂u ∂x x+ ∂u ∂t t )(ν) We denote the k-th derivative of the function u(x,t) in the point (0,0) (i.e. m-th with respect to t and (k−m)-th with respect to x, m=0,1, . . . ,k) as uxk−mtm. Hence dku k! = 1 k! k ∑ m=0 ( k m ) uxk−mtmx k−mtm Using equation (2.2) we can substitute the partial derivative of the order 4ν, ν =1,2, . . . as follows ∂4νu ∂x4ν =(−1)ν ∂2νu ∂t2ν (2.4) Formula (2.4) holds for differentials of the order greater than 3. For example du 1! = ( 1 0 ) uxx+ ( 1 1 ) utt d2u 2! = ( 2 0 ) ux2 x2 2! + ( 2 1 ) uxt xt 2! + ( 2 2 ) ut2 t2 2! d3u 3! = ( 3 0 ) ux3 x3 3! + ( 3 1 ) ux2t x2t 3! + ( 3 2 ) uxt2 xt2 3! + ( 3 3 ) ut3 t3 3! d4u 4! =− ( 4 0 ) ut2 x4 4! + ( 4 1 ) ux3t x3t 4! + ( 4 2 ) ux2t2 x2t2 4! + ( 4 3 ) uxt3 xt3 4! + ( 4 4 ) ut4 t4 4! d5u 5! =− ( 5 0 ) ut2x x5 5! − ( 5 1 ) ut3 x4t 5! + ( 5 2 ) ux3t2 x3t2 5! + ( 5 3 ) ux2t3 x2t3 5! + + ( 5 4 ) uxt4 xt4 5! + ( 5 5 ) ut5 t5 5! . . . 350 M.J. Al-Khatib et al. Let us transform formula (2.3)with the use of (2.4), and then group the coeffi- cients at consecutive derivatives. Finally, we obtain the following polynomials S0,0 =1 S1,0 = ( 1 0 ) x S1,1 = ( 1 1 ) t S2,0 = ( 2 0 ) x2 2 S2,1 = ( 2 1 ) xt 2 S2,2 = ( 2 2 ) t2 2! − ( 4 0 ) x4 4! S3,0 = ( 3 0 ) x3 3! S3,1 = ( 3 1 ) x2t 3! S3,2 = ( 3 2 ) xt2 3! − ( 5 0 ) x5 5! S3,3 = ( 3 3 ) t3 3! − ( 5 1 ) x4t 5! S4,1 = ( 4 1 ) x3t 4! S4,2 = ( 4 2 ) x2t2 4! − ( 6 0 ) x6 6! S4,3 = ( 4 3 ) xt3 4! − ( 6 1 ) x5t 6! S4,4 = ( 4 4 ) t4 4! − ( 6 2 ) x4t2 6! + ( 8 0 ) x8 8! . . . For m=0,1,2, . . .; k=0,1,2,3; m−k­ 0 the polynomials have a form Sm,m−k = [m−k 2 ] ∑ j=0 (−1)jtm−k−2jxk+4j (m−k−2j)!(k+4j)! (2.5) where [x] = floor(x). It is convenient to renumberpolynomials (2.5) according to Table 1. Table 1.New numbering of the polynomials m=0 m=1 m=2 m=3 m=4 m=5 m=6 .. . 0,0 0 1,0 1 2,0 3 3,0 6 4,1 10 5,2 14 6,3 18 . . . 1,1 2 2,1 4 3,1 7 4,2 11 5,3 15 6,4 19 . . . 2,2 5 3,2 8 4,3 12 5,4 16 6,5 20 . . . 3,3 9 4,4 13 5,5 17 6,6 21 . . . The method of solving polynomials... 351 Then we get S0 =S0,0(x,t) = 1 S1 =S1,0(x,t)=x S2 =S1,1(x,t) = t S3 =S2,0(x,t)= x2 2 S4 =S2,1(x,t) =xt S5 =S2,2(x,t)= t2 2! − x4 4! and for m­ 3 and k=0,1,2,3 Sn =S4m−k−3 =Sm,m−k = [m−k 2 ] ∑ j=0 (−1)jtm−k−2jxk+4j (m−k−2j)!(k+4j)! Let us denote k = 4[n−2 4 ]−n+5 and m = [n−2 4 ]+2. Then, for n ­ 3, we obtain Sn(x,t)= [ n−3[ n−2 4 ]−3 2 ] ∑ j=0 (−1)jtn−3[ n−2 4 ]−3−2jx(4[ n−2 4 ]−n+5)+4j (n−3[n−2 4 ]−3−2j)!((4[n−2 4 ]−n+5)+4j)! (2.6) 2.2. Partial derivatives of th solving polynomials In numerical practice, recurrent formulas are very useful. Theorem 1 presents such formulas for the derivatives of solving polynomials in the beamvibration problem. Properties shown in the theorem prove that the considered polyno- mials satisfy equation (2.2). That is why they are called the ”solving polyno- mials”. Theorem 1 For polynomials Sn, we get: a) for n­ 5 ∂Sn ∂x = { Sn−3 for n 6=4p+1, p=1,2, . . . −Sn+1 for n=4p+1, p=1,2, . . . (2.7) b) for n­ 8 ∂2Sn ∂x2 = { Sn−6 for n=4p+2 or n=4p+3, p=2,3, . . . −Sn−2 for n=4p or n=4p+1, p=2,3, . . . (2.8) 352 M.J. Al-Khatib et al. c) for n­ 11 ∂3Sn ∂x3 = { Sn−9 for n=4p+2, p=3,4, . . . −Sn−5 for n 6=4p+2, p=3,4, . . . (2.9) d) for n­ 11 ∂4Sn ∂x4 =−Sn−8 (2.10) e) for n­ 7 ∂Sn ∂t =Sn−4 (2.11) f) for n­ 11 ∂2Sn ∂t2 =Sn−8 (2.12) Properties (2.10) and (2.12) prove that the polynomials Sn satisfy equation (2.2). Proof of the Theorem 1 a) For n­ 3 Sn(x,t)= [ n−3[ n−2 4 ]−3 2 ] ∑ j=0 (−1)jtn−3[ n−2 4 ]−3−2jx(4[ n−2 4 ]−n+5)+4j (n−3[n−2 4 ]−3−2j)!((4[n−2 4 ]−n+5)+4j)! In the first summand of Sn, the variable x to the power of w=4[ n−2 4 ]−n+5 appears. It can be easily shown that for n = 4p; n = 4p+2; n = 4p+3 (p= 1,2, . . .) we obtain w ­ 1, and for n= 4p+1 we get w = 0. Then, for n=4p ∂Sn(x,t) ∂x = [ n−3[ n−2 4 ]−3 2 ] ∑ j=0 (−1)jtn−3[ n−2 4 ]−3−2jx(4[ n−2 4 ]−n+4)+4j (n−3[n−2 4 ]−3−2j)!((4[n−2 4 ]−n+4)+4j)! = = [ 4p−3[ 4p−2 4 ]−3 2 ] ∑ j=0 (−1)jt4p−3[ 4p−2 4 ]−3−2jx(4[ 4p−2 4 ]−4p+4)+4j (4p−3[4p−2 4 ]−3−2j)!((4[4p−2 4 ]−4p+4)+4j)! = = [ p 2 ] ∑ j=0 (−1)jtp−2jx4j (p−2j)!(4j)! The method of solving polynomials... 353 On the other hand, Sn−3(x,t)=S4p−3(x,t)= = [ 4p−3−3[ 4p−3−2 4 ]−3 2 ] ∑ j=0 (−1)jt4p−3−3[ 4p−3−2 4 ]−3−2jx(4[ 4p−3−2 4 ]−4p−3+4)+4j (4p−3−3[4p−3−2 4 ]−3−2j)!((4[4p−3−2 4 ]− (4p−3)+4)+4j)! = = [ p 2 ] ∑ j=0 (−1)jtp−2jx4j (p−2j)!(4j)! This proves property a) for n=4p. For n=4p+2; n=4p+3 the proof is similar. For n=4p+1we get ∂Sn(x,t) ∂x = = [ 4p+1−3[ 4p+1−2 4 ]−3 2 ] ∑ j=1 (−1)jt4p+1−3[ 4p+1−2 4 ]−3−2jx(4[ 4p+1−2 4 ]−(4p+1)+4)+4j (4p+1−3[4p+1−2 4 ]−3−2j)!((4[4p+1−2 4 ]− (4p+1)+4)+4j)! = = [ p+1 2 ] ∑ j=1 (−1)jtp+1−2jx−1+4j (p+1−2j)!(−1+4j)! = substituting j= s+1 =− [ p−1 2 ] ∑ s=0 (−1)stp−1−2sx3+4s (p−1−2s)!(3+4s)! On the other hand, for n=4p+1 −Sn+1(x,t) = =− [ 4p+2−3[ 4p 4 ]−3 2 ] ∑ j=0 (−1)jt4p+2−3[ 4p 4 ]−3−2jx(4[ 4p 4 ]−(4p+2)+5)+4j (4p+2−3[4p 4 ]−3−2j)!((4[4p 4 ]− (4p+2)+5)+4j)! = =− [ p−1 2 ] ∑ j=0 (−1)jtp−1−2jx3+4j (p−1−2j)!(3+4j)! Thiscompletes theproofof propertya).Propertiesb), c), d) followpropertya). 354 M.J. Al-Khatib et al. e)We have Sn(x,t) = [ n−3[ n−2 4 ]−3 2 ] ∑ j=0 (−1)jtn−3[ n−2 4 ]−3−2jx(4[ n−2 4 ]−n+5)+4j (n−3[n−2 4 ]−3−2j)!((4[n−2 4 ]−n+5)+4j)! The variable t appears in the derivative ∂Sn/∂t if and only if in Sn the inequality w = n − 3[n−2 4 ] − 3 − 2j ­ 1 holds. Hence, for j we get j¬ (n−3[n−2 4 ]−4)/2 and ∂Sn(x,t) ∂t = [ n−3[ n−2 4 ]−4 2 ] ∑ j=0 (−1)jtn−3[ n−2 4 ]−4−2jx(4[ n−2 4 ]−n+5)+4j (n−3[n−2 4 ]−4−2j)!((4[n−2 4 ]−n+5)+4j)! On the other hand Sn−4 = [ n−3[ n−2 4 −1]−7 2 ] ∑ j=0 (−1)jtn−3[ n−2 4 −1]−7−2jx(4[ n−2 4 −1]−n+9)+4j (n−3[n−2 4 −1]−7−2j)!((4[n−2 4 −1]−n+9)+4j)! = = [ n−3[ n−2 4 ]−4 2 ] ∑ j=0 (−1)jtn−3[ n−2 4 ]−4−2jx(4[ n−2 4 ]−n+5)+4j (n−3[n−2 4 ]−4−2j)!((4[n−2 4 ]−n+5)+4j)! This proves property e). Property f) results directly from e). 2.3. Generating function Using the method of variable separation in equation (2.2), we arrive at the function W(x,t,p) called the generating function for the solving functions W(x,t,p)= epx+ip 2t (2.13) with i= √ −1. Expanding the function W into the power series, we obtain W(x,t,p)= epx+ip 2t = =1+px+ p2 2! (2it+x2)+ p3 3! (6itx+x3)+ p4 4! (−12t2+12itx2+x4)+ . . .= =S0+pS1+p 2(iS2+S3)+p 3(−iS4+S6)+p4(−S5+iS7)+ . . . It is clear that the real and imaginary parts of coefficients at successive powers of the parameter p are the previously obtained solving polynomials. The method of solving polynomials... 355 3. The method of solving functions Themethod of solving functions discussedherebelongs to a class of theTrefftz methods. The solving polynomials may be used to obtain an approximate solution to equation (2.2) with prescribed initial and boundary conditions.We can determine a solution in thewhole domain (0,1)×(0,∆t) orwe canuse the polynomials as base functions in the Finite Element Method. In the second case, at least three approaches are applied: continuous FEM, non-continuous FEM (Ciałkowski and Frąckowiak, 2007) or nodeless FEM (Maciąg, 2007b). Unfortunately, for beamvibrationproblems, theuse of the solvingpolynomials in continuous anddiscontinuousFEMdoesnot lead to satisfactory results.The matrix obtained inFEM is then ill-conditioned. This problemdoes not appear in the nodeless FEM. Therefore, in this paper we consider two cases: how to use the solving polynomials in the whole domain (0,1)× (0,∆t) and how to use them in the nodeless FEM. In both cases we take a linear combination of the polynomials in order to findan approximate solution (in thewhole domain or in time-space elements) u≈w= N ∑ n=1 cnSn (3.1) Because the polynomials Sn satisfy equation (2.2), their linear combination satisfies this equation aswell. The coefficients cn are chosen to obtain the best fitting to the initial and boundary conditions. Additionally, in the nodeless FEM we have to minimize the difference of solutions between the elements (see example). 4. Example The goal of the example presented here is not to show all possible cases of beam vibration. The authors want only to show how to find an approximate solution to the beam vibration problem with satisfactory accuracy using the method of solving functions. In numerical practice, the ”exact (analytical) solution” very often does not lead to accurate numerical values of solution because of truncation errors. For vibrations of a beam, an analytical solution includes roots of a transcendental equation. That may cause some numerical problems. 356 M.J. Al-Khatib et al. The method is tested here for a problem for which the exact solution is known.We consider two ways of using the solving polynomials. The first one is the use of the polynomials in the whole domain. In the second approach, the solving functions stands for base functions in the nodeless Finite Element Method. 4.1. Problem formulation Consider vibrations of a beamwhich is fixed at x=0 and free at x=1: — equation of motion ∂4u ∂x4 + ∂2u ∂x2 =0 x∈ (0,1), t > 0 (4.1) — initial conditions u(x,0)=u0(x)= 1 1000 x2 (4.2) ∂u(x,0) ∂t =0 (4.3) — boundary conditions, t> 0 u(0, t) = 0 (4.4) ∂u(0, t) ∂x =0 (4.5) ∂2u(x,1) ∂x2 =0 (4.6) ∂3u(x,1) ∂x3 =0 (4.7) The exact solution to the problem reads u(x,t)= ∞ ∑ n=1 Kn(x)cos(γ 2 nt) ∫1 0 Kn(x)u0(x) dx (sinγn+sinhγn)2 (4.8) where Kn(x)= (cosγn+coshγn)[sinh(γnx)− sin(γnx)]+ −(sinγn+sinhγn)[cosh(γnx)− cos(γnx)] and γn are successive positive roots of the equation 1+cosxcoshx=0. The method of solving polynomials... 357 4.2. Approximate solution in the whole domain An approximate solution u(x,t) has a form as follows u≈w= N∑ n=1 cnSn Because the function w satisfies the equation ofmotion of the vibrating beam (with w being a linear combination of the solving polynomials), in order to find the coefficients cn weminimize the norm ‖w−u‖ for the boundary and initial conditions.We look for an approximate solution u in the time interval (0,∆t). Hence, the coefficients cn have to be appropriately chosen tominimize the functional I = 1∫ 0 [w(x,0)−u0(x)]2 dx ︸ ︷︷ ︸ cond. (4.2) + 1∫ 0 [∂w(x,0) ∂t −0 ]2 dx ︸ ︷︷ ︸ cond. (4.3) + ∆t∫ 0 [w(0, t)−0]2 dt ︸ ︷︷ ︸ cond. (4.4) + ∆t∫ 0 [∂w(0, t) ∂x −0 ]2 dt ︸ ︷︷ ︸ cond. (4.5) (4.9) + ∆t∫ 0 [∂2w(1, t) ∂x2 −0 ]2 dt ︸ ︷︷ ︸ cond. (4.6) + ∆t∫ 0 [∂3w(1, t) ∂x3 −0 ]2 dt ︸ ︷︷ ︸ cond. (4.7) The necessary condition to minimize the functional I is ∂I ∂c1 = . . .= ∂I ∂cN =0 (4.10) The linear system of equations (4.10) can be written as AC=B (4.11) 358 M.J. Al-Khatib et al. where C= [c1, . . . ,cN] ⊤ and the elements of matrices A andB are aij = 1∫ 0 Si(x,0)Sj(x,0) dx+ 1∫ 0 ∂Si(x,0) ∂t ∂Sj(x,0) ∂t dx + ∆t∫ 0 Si(0, t)Sj(0, t) dt+ ∆t∫ 0 ∂Si(0, t) ∂x ∂Sj(0, t) ∂x dt + ∆t∫ 0 ∂2Si(1, t) ∂x2 ∂2Sj(1, t) ∂x2 dt+ ∆t∫ 0 ∂3Si(1, t) ∂x3 ∂3Sj(1, t) ∂x3 dt bi = 1∫ 0 Si(x,0)u0(x) dx Note that aij = aji (thematrix A is symmetric) which simplifies the calcula- tions. From equation (4.11), we obtain the coefficients cn: C=A −1 B. In the time intervals (∆t,2∆t), (2∆t,3∆t), . . . , we proceed analogously. Here, the initial condition for the time interval ((m−1)∆t,m∆t) is the value of the function u at the end of the interval ((m−2)∆t,(m−1)∆t). All results shown in Figures 1-4 have been obtained for ∆t = 1. The number N of polynomials in formula (3.1) depends on the degree D of the polynomials. Here N =2D+1. In the whole time interval (0,∆t) a good approximation is obtained. For example, Fig.1 shows (a) the exact solution, (b) the approximation with the solving polynomials of degrees up to 30, (c) the difference between (a) and (b). It is clear that the approximation is satisfactory. The largest error has been obtained for x=1 (for x=0 the beam is fixed). Figure 2 shows the initial condition u(x,0) in the case of (a) the exact solution, (b) the approximationwith the solving polynomials of degrees up 30, (c) the difference. Figure 3 shows the solution u(x,0.5) in the case of (a) the exact solution, (b) the approximationwith the solvingpolynomials of degrees up to 30, (c) the difference. Figures 2 and 3 show that the approximation with polynomials of degrees up to 30 is satisfactory. The determination of the highest degree of the solving polynomials used in approximation (3.1) is very important. Our goal is to minimize the inaccuracy of the approximation. The solution should be ”as accurate aspossible”.Themethodof solving functions shouldbeconvergent. It means that usingmore polynomials one should obtain a better approximation. The method of solving polynomials... 359 Fig. 1. Solutions (a) exact, (b) approximated with the solving polynomials of degrees up 30, (c) difference Fig. 2. The initial condition in the case of (a) exact solution, (b) approximationwith the solving polynomials of degrees up to 30, (c) difference In general, it is true. Unfortunately, a too big number of polynomials leads to big dimensions of the matrix in (4.11). That may cause numerical problems. Figure 4 shows the exact and approximate solution for x=1with the solving polynomials of degrees up to (a) 10, (b) 25, (c) 30. 360 M.J. Al-Khatib et al. Fig. 3. Solution for t=0.5 in the case of (a) exact solution, (b) approximation with the solving polynomials of degrees up to 30, (c) difference Fig. 4. Exact and approximate solutions for x=1 with the solving polynomials of degrees up to (a) 10, (b) 25, (c) 30 Let wD(1, t) denote the approximation w(1, t) withpolynomials of degrees 0 to D. The average, relative difference between the solutions wD(1, t) and u(1, t) may be defined as follows The method of solving polynomials... 361 E(D)= √ √ √ √ √ √ √ √ ∆t∫ 0 [wD(1, t)−u(1, t)]2 dt ∆t∫ 0 [u(1, t)]2 dt Table 2 shows the difference E(D)which depends on the highest degree D of the solving polynomial. Table 2.Error in relation to the highest degree D of the solving polynomial OrderD 10 15 20 25 30 E(D) 0.0917 0.0884 0.0843 0.0741 0.0221 The difference E(D) decreases when the number of the polynomials in the approximation u increases. It means that the truncation error for the considered values of D is negligible, and that the procedure is convergent. 4.3. Nodeless Finite Element Method The secondway of using solving polynomials is the nodeless FEM.Let us divi- de the interval (0,1) into subintervals and seek the solution in the time-space subdomains Lk = ((k−1)∆x,k∆x)×∆t, k = 1, . . . ,K. Let us introduce a local co-ordinate system in each subdomain Lk andassume the approximation of the solution in the form uk ≈wk = N∑ n=1 cknSn (4.12) To find the coefficients ckn, we minimize the norm ‖w−u‖ for the boundary and initial conditions. Moreover, the norms: ‖wk−1 −wk‖, ‖∂wk−1∂x − ∂wk ∂x ‖, ‖∂ 2wk−1 ∂x2 − ∂ 2wk ∂x2 ‖ and ‖∂ 3wk−1 ∂x3 − ∂ 3wk ∂x3 ‖ are minimized simultaneously. In the time intervals (∆t,2∆t), (2∆t,3∆t), . . . , we proceed analogously. Here, the value of the function w at the endof the interval ((m−2)∆t,(m−1)∆t) stands for the initial condition for the time interval ((m− 1)∆t,m∆t). In order to compare the solution in the whole interval x ∈ (0,1) with the approximate solution obtainedwith thenodelessFEMweassume ∆t=1/2. Figure 5 shows the exact solution, approximation with the solving polynomials of degrees up to 20 and the difference for x=1, for the number of subintervals K =2 and for the number of time steps‘ TS=2. When comparingFig.4 andFig.5,wenotice that for thenodelessFEMthe approximation is better than for the solution in the whole domain.Moreover, 362 M.J. Al-Khatib et al. Fig. 5. (a) Exact and approximate solutions for x=1 (b) the difference usingonly two subintervals (i.e. K=2)allowedone to reduce thehighestdegree of the solvingpolynomials from30down to 20.Whenusing the nodeless FEM, theadvantage is such that the elements canbe relatively big.Moreover, in each element we have a solution satisfying the equation. Like in Section 4.2, wedefine the relative error of the approximate solution, E(D), as follows E(D) = √ √ √ √ √ √ √ √ TS∆t∫ 0 [wD(1, t)−u(1, t)]2 dt TS∆t∫ 0 [u(1, t)]2 dt (4.13) where TS stands for the number of time steps. Table 3 shows values of the error E(D) for TS = 2, ∆t = 1/2, which depend on the highest degree, D, and the number of subintervals, K. Table 3.Error dependence on the polynomial highest degree, D, and number of subintervals K K =2 K =3 K =4 K =5 D=10 0.0928 0.0951 0.0943 0.0932 D=15 0.06 0.0515 0.0436 0.0433 D=20 0.0073 0.011 0.0092 0.008 It is clear that the number of polynomials used in the approximation has the greatest influence on the error, E(D). Generally, the increasing number of subintervals decreases the error, but the relationship is not so distinct: probably the truncation errors increase in that case.When comparingTable 2 and 3, we see that the nodeless FEM leads to amore accurate approximation than the solution in the whole domain. It is interesting how the nodeless FEMworks for a greater number of time steps; the method should then give ”satisfactory” results. However, when we The method of solving polynomials... 363 increase the number of time steps, the error of approximation increases too because of the truncation errors. And because the solutions for each time step are not accurate, the error propagates in time. Figure 6 shows the exact solution, approximation with the solving polynomials of degrees up to 20 and the difference for x = 1 with the number of subintervals K = 2 and the number of time steps TS=10. Fig. 6. (a) Exact and approximate solutions for x=1 (b) the difference Table 4 shows E(D) for ∆t = 1/2, degree D = 20 and the number of subintervals K =2 as a function of the number of time steps TS. Table 4.Error in relation to the number of time steps TS 2 4 6 8 10 E(20) 0.0073 0.013 0.02 0.028 0.035 The error increases when the number of time steps increases too, but re- mains small. For 10 time steps, the error does not exceed 3.5%, which shows that the approximation is good. 5. Concluding remarks A new simple technique for solving the problem of beam vibration has been developed. The test example shows that the obtained approximations are very good. The solution, which is a linear combination of the solving polynomials, exactly satisfies the differential equation of motion and – approximately – the initial and boundary conditions. The solving polynomials presented here lead to good results both in the whole domain and in the nodeless FEM. In the second case, the elements can be relatively big. The main advantage of the method of solving functions is 364 M.J. Al-Khatib et al. that the approximate solutionof the consideredproblemsatisfies thegoverning equation. Acknowledgments The paper has been financed by Ministry of Science and Higher Education, Poland,Warsaw. Grant N51300332/0541. References 1. Ciałkowski M.J., 1999 Solution of inverse heat conduction problemwith use new type of finite element base functions, In: Proceedings of the Internatio- nal Symposium on Trends in Continuum Physics, Maruszewski B.T., Muschik W., Radowicz A. (Eds.), World Scientific Publishing, Singapore, New Jersey, London, Hong Kong, 64-78 2. Ciałkowski M.J., 2003, Thermal and related functions used in solving certa- in problems of mechanics. Part I – Solution of certain differential equations by means of inverse operations, Studia i materiały. Technika. Uniwersytet Zielo- nogórski, 3, 7-70 [in Polish] 3. Ciałkowski M.J., Frąckowiak A., 2000,Heat Functions and Their Appli- cation for Solving Heat Transfer and Mechanical Problems, Poznań University of Technology Publishers [in Polish] 4. Ciałkowski M.J., Frąckowiak A., 2003, Thermal and related functions used in solving certain problems of mechanics. Part II – Effective determina- tion of inverse operations applied to harmonic functions, Studia i materiały. Technika. Uniwersytet Zielonogórski, 3, 71-98 [in Polish] 5. CiałkowskiM.J.,FrąckowiakA., 2004,Applicationof symbolic operations to the determination of Trefftz functions for a heat flowwave equation,Zeszyty Naukowe Politechniki Poznańskiej, Maszyny Robocze i Transport,57 [in Polish] 6. Ciałkowski M.J., Frąckowiak A., Grysa K., 2007, Solution of a sta- tionary inverse heat conduction problems by means of Trefftz non-continuous method, International Journal of Heat Mass Transfer, 50, 2170-2181 7. Ciałkowski M.J., Futakiewicz S., Hożejowski L., 1999a, Heat polyno- mials applied to direct and inverse heat conduction problems, In: Proceedings of the International Symposium onTrends in ContinuumPhysics,Maruszewski B.T.,MuschikW., RadowiczA. (Eds.),World ScientificPublishing, Singapore, New Jersey, London, Hong Kong, 79-88 8. Ciałkowski M.J., Futakiewicz S., Hożejowski L., 1999b,Method of he- at polynomials in solving the inverse heat conduction problems, ZAMM, 79, 709-710 The method of solving polynomials... 365 9. Ciałkowski M.J., Jarosławski M., 2003, Application of symbolic calcula- tions in generating the solution of the wave equation, Zeszyty Naukowe Poli- techniki Poznańskiej, Maszyny Robocze i Transport, 56 [in Polish] 10. Futakiewicz S., 1999,Heat Functions Method for Solving Direct and Inverse Heat Conductions Problems, PhDThesis, PoznańUniversity of Technology [in Polish] 11. Futakiewicz S., Grysa K., Hożejowski L., 1999, On a problem of bo- undary temperature identification in a cylindrical layer, In: Proceedings of the International Symposium onTrends in ContinuumPhysics, Maruszewski B.T., MuschikW., Radowicz A. (Eds.),World Scientific Publishing, Singapore, New Jersey, London, Hong Kong, 119-125 12. Futakiewicz S., Hożejowski L., 1998a, Heat polynomials in solving the direct and inverse heat conduction problems in a cylindrical system of coordi- nates, In: Advanced Computational Method in Heat Transfer V, Nowak A.J., Brebbia C.A., Bialecki R., ZerroukatM. (Eds.), ComputationalMechanics Pu- blications, Southampton, Boston, 71-80 13. Futakiewicz S., Hożejowski L., 1998b, Heat polynomials method in the n-dimensional direct and inverse heat conduction problems, In:AdvancedCom- putational Method in Heat Transfer V, NowakA.J., Brebbia C.A., Bialecki R., Zerroukat M. (Eds.), Computational Mechanics Publications, Southampton, Boston, 103-112 14. Hożejowski L., 1999,Heat Polynomials and their Application for Solving Di- rect and Inverse Heat Condutions Problems, PhD Thesis, Kielce University of Technology [in Polish] 15. Maciąg A., 2004, Solution of the three-dimensional wave equation by using wave polynomials,PAMM – Proc. Math. Mech., 4, 706-707 16. MaciągA., 2005, Solution of the three-dimensionalwave polynomials,Mathe- matical Problems in Engineering, 5, 583-598 17. Maciąg A., 2007a, Two-dimensional wave polynomials as base functions for continuity and discontinuity Finite ElementsMethod, In:Współczesne techno- logie i urządzenia energetyczne, Taler J. (Eds.), Kraków, 371-381 [in Polish] 18. Maciąg A., 2007b, Wave Polynomials in elasticity problems, Engineering Transactions, 55, 2, 129-153 19. Maciąg A.,Wauer J., 2005a, Solution of the two-dimensionalwave equation byusingwavepolynomials,Journal of EngineeringMathematics,51, 4, 339-350 20. Maciąg A., Wauer J., 2005b, Wave polynomials for solving different types of two-dimensionalwave equations,Computer AssistedMechanics and Engine- ering Sciences, 12, 87-102 21. RosenbloomP.C.,Widder D.V., 1956, Expansion in terms of heat polyno- mials and associated functions,Trans. Am. Math. Soc., 92, 220-266 366 M.J. Al-Khatib et al. 22. Yano H., Fukutani S., Kieda A., 1983, A boundary residual method with heatpolynomials for solvingunsteadyheat conductionproblems,Franklin Inst., 316, 291-298 Metoda wielomianów rozwiązujących dla problemu drgań belki Streszczenie Artykułprzedstawianowąmetodęprzybliżonego rozwiązywaniaproblemówdrgań belki, które opisywane są cząstkowym równaniem różniczkowym czwartego rzędu. Generalnie takie równania mogą być rozwiązywane różnymimetodami. Główna idea prezentowanego tutaj podejścia polega na znalezieniu wielomianów spełniających w sposób ścisły dane równanie różniczkowe (funkcje rozwiązujące). Za rozwiązanie przybliżone przyjmuje się kombinację liniową tych funkcji, która również spełnia rów- nanie.Współczynniki kombinacji liniowejwyznaczane są tak, aby uzyskane rozwiąza- nie w sposób najlepszy (w sensie średniokwadratowym) było dopasowane do danych warunków początkowych i brzegowych. Jest to zatemwariantmetody Trefftza. W pracy wyznaczono efektywne wzory dla funkcji rozwiązujących i ich pochod- nych.Uzyskanewielomianymogąbyćwykorzystanedowyznaczenia rozwiązaniawca- łym obszarze lub też mogą być użyte jako funkcje bazowe w bezwęzłowej Metodzie ElementówSkończonych.Obaprzypadki rozważonowpracy. Załączono również przy- kłady numeryczne. Manuscript received December 5, 2007; accepted for print January 24, 2008